随着多旋翼无人飞行器的广泛应用，对其飞行综合性能的要求也日益提高，促使了多种新多旋翼飞行器或多旋翼-固定翼复合构型的诞生。同时，多旋翼飞行器气动效率与操稳性之间的矛盾性也逐步体现出来。为了提高多旋翼飞行器气动效率，可以采用加大旋翼直径与降低旋翼转速的设计，但工程实践表明，这会显著影响飞行器的配平能力与飞行动态品质。此外，针对多旋翼飞行器的先进导航、制导与飞控算法也逐渐成为热门领域。上述发展对多旋翼飞行器的全机理论建模与分析工作提出了客观需求。
目前，学界关于多旋翼飞行器详细飞行力学建模与特性分析的研究数量相对较少，有****在研究飞控算法过程中针对飞控系统设计需求已经建立了一些多旋翼飞行力学模型[1-8]。然而，已有模型往往较简单，目前建模工作中普遍采用的简化假设包括：① 在旋翼产生的3个力分量与3个力矩分量中，只考虑了旋翼所产生的拉力与扭矩，忽略其他分量[1, 9]；② 假设旋翼的拉力和扭矩与旋翼转速的平方成正比，且飞行器飞行状态不会影响旋翼气动力[2]；③ 忽略了入流随拉力的变化与非均匀分布情况[10]；④ 未考虑旋翼旋转自由度上的动态特性。
基于上述假设建立的模型可以在一定程度上反映多旋翼飞行器的特性，在一定的飞行条件下，采用这类模型是合理的。随着飞行速度的提高，旋翼气动力不仅是旋翼转速的函数，同时也与来流条件紧密相关[11-12]。单个旋翼除产生拉力与扭矩外，还会产生相对较大的水平、侧向气动力与俯仰、滚转力矩[11-12]。在来流角与旋翼俯仰、滚转力矩的作用下，根据动态入流理论[13]，旋翼桨盘上的诱导速度非均匀分布情况将加剧，从而对旋翼气动力造成影响。由此可见，随着用户对多旋翼飞行器飞行性能要求的不断提高，更加全面的建模分析工作显得尤为必要。
与此同时，已有的模型大多直接用于飞控系统的设计或简单的飞行模拟验证，例如姜军等[14]针对可投放多旋翼无人飞行器控制系统的设计研究。目前对多旋翼飞行器飞行性能与飞行力学特性的深入分析研究还不多见，对多旋翼飞行器操稳特性、模态特性、气动导数特点等还缺乏系统的认识。
因此，本文首先基于已有的传统直升机飞行力学模型，建立了多旋翼飞行器飞行力学模型。然后基于此模型在悬停状态对一典型四旋翼无人飞行器进行了配平与线性化处理。最后通过线化小扰动方程，分析了多旋翼飞行器建模中旋翼旋转自由度、入流模型与气动力矩的重要性，以及多旋翼飞行器部分模态与气动导数的特点。
1 飞行力学模型 本文所采用的多旋翼飞行力学模型是由荷兰代尔夫特理工大学的DelCopter通用直升机飞行力学模型发展而来的。DelCopter模型是基于叶素理论的传统解析模型，考虑了刚性桨叶的挥舞与摆振自由度，包含了动态入流模型，已经经过了飞行测试数据验证[15]。
1.1 气动模型 本文选用了DelCopter模型中的旋翼气动力建模部分，并忽略了桨叶挥舞与摆振运动。通过对桨叶叶素上的气动力进行积分，可以得到单个旋翼在与其旋转轴重合的坐标系内的3个力分量与3个力矩分量，即[Xi, Yi, Zi, Li, Mi, Ni]。其中，[Xi, Yi, Zi]分别为x-、y-和z-轴方向上的力；[Li, Mi, Ni]分别为x-、y-和z-轴方向上的力矩。由于篇幅限制，本文不再给出上述力与力矩的具体表达式。详细的推导过程与公式，读者可参考文献[15]的附录D。
在获得单个旋翼气动特性之后，可对所有旋翼气动力求和得出旋翼总气动力。假设旋翼无安装角，即旋翼旋转轴与机体z-轴平行，可得到所有旋翼总气动力[Xr, Yr, Zr]与总气动力矩[Lr, Mr, Nr]：
 | (1) |
 | (2) |
式中：xi、yi和zi为旋翼在机体坐标系中的位置；Nr为旋翼个数。
对于机身气动力，假设机身仅产生一个沿来流方向的阻力Df，且气动中心与重心重合，即不会在重心处产生力矩：
 | (3) |
式中：ρ为空气密度；F0为机身阻力常数；u、v和w分别为机体x-、y-和z-轴方向上的线速度。
1.2 入流模型 旋翼入流模型采用了传统的Pitt-Peters入流模型，桨盘上入流速度vi分布为
 | (4) |
式中：vi0、vis和vic分别为平均、横向和纵向入流分布；re为叶素径向位置；R为旋翼半径；ψ为桨叶方位角。
因入流动态响应速度较快，本文忽略了入流的动态特性，采用了准定常假设，即在任意时刻旋翼上的气动力将直接决定入流分布情况：
 | (5) |
式中：
1.3 旋转自由度 为了简化分析，专注多旋翼飞行器本身特性，本文采用了简化的电机模型。假设电机功率输出能够实时跟随输入信号，根据功率与转速可以得到电机扭矩为
 | (6) |
式中：Qi为第i个电机扭矩；Pmax为电机最大功率；Ui为第i号电机无量纲化输入信号(0≤Ui≤1)；Ωi为当前转速。
结合旋翼气动力，可以得到旋翼与电机转速变化率为
 | (7) |
式中：Ir和Im分别为旋翼和电机转动惯量。
1.4 动力学模型 将式(1) 与式(3) 确定的气动力代入传统飞行器刚体运动方程，可以得到多旋翼飞行力学模型的机体动力学部分：
 | (8) |
式中：Eb代表机体运动非线性方程组；Xb为机体的刚体运动自由度，包括3个平动速度、3个转动角速度与3个欧拉角；Ω为旋翼转速自由度，包括所有旋翼旋转角速度。
由于多旋翼飞行器的飞行操纵是通过调节各旋翼转速实现，其操纵输入量不会直接影响机体运动自由度，因此式(8) 中没有显式地包含操纵输入。
联立所有旋翼的旋转动力学方程(式(7))，可以得到多旋翼飞行力学模型的旋翼旋转动力学部分：
 | (9) |
式中：Er代表所有旋翼旋转运动非线性方程组；U为操纵输入向量。
结合式(8) 与式(9)，可以得到多旋翼飞行器非线性飞行力学模型，总计包括9+Nr个方程。
2 计算算例 为了分析对比多旋翼飞行器与传统旋翼飞行器的特性区别，本文分别选取了四旋翼无人飞行器与传统单旋翼直升机进行对比计算分析。其中，四旋翼为某公司无人飞行器产品，如图 1所示，总体参数如表 1所示。
 |
| 图 1 算例四旋翼无人飞行器 Fig. 1 Study case of quad-rotor UAV |
| 图选项 |
表 1 四旋翼无人飞行器总体参数 Table 1 Main parameters of quad-rotor UAV
| 设计参数 | 数值 |
| 起飞质量/kg | 5 |
| 布局形式 | X4布局 |
| 总功率/W | 2 000 |
| 旋翼直径/m | 0.457 2 |
| 旋翼轴距/m | 0.85 |
| 最大速度/(m·s-1) | 20 |
表选项
常规直升机选取了Bo-105直升机算例。Bo-105为传统单旋翼带尾桨直升机，起飞质量约2 200 kg，是旋翼飞行器研究中广泛采用的计算算例。DelCopter模型对Bo-105的分析计算已经经过验证，可作为对比参照。
针对上述算例，本文利用多旋翼模型与DelCopter模型在悬停状态进行了配平计算，并在平衡点利用数值差分对全机模型进行了线性化处理，得到了对应的线化小扰动模型。根据线化模型稳定性导数矩阵，可以求得飞行器模态特征值。算例典型模态对比如表 2所示。
表 2 算例悬停模态对比 Table 2 Comparison of hovering modes between two study cases
| 四旋翼无人飞行器 | Bo-105 | |||
| 模态 | 特征值 | 模态 | 特征值 | |
| 俯仰收敛 | -1.699 | 俯仰收敛 | -5.449 | |
| 滚转收敛 | -1.699 | 滚转收敛 | -17.14 | |
| Phugoid | 0.550 6+1.259 2i | Phugoid | 0.264 9+0.474 9i | |
| 横向Phugoid | 0.550 6+1.259 2i | 荷兰滚 | -0.274 8+0.559 4i | |
| 沉浮收敛 | -0.500 2 | 沉浮收敛 | -0.419 3 | |
| 螺旋 | -0.033 5 | 螺旋 | -0.284 0 | |
表选项
通过上述对比可以发现，多旋翼飞行器在模态特性上与传统直升机有较大差别，有必要进一步分析其飞行力学特性及成因。在悬停状态下，多旋翼飞行器模态主要特点如下：
1) 由于纵横向的对称性，悬停状态下四旋翼无人飞行器的纵向模态与横向模态是相同的。
2) 由于没有尾桨与垂直尾翼，多旋翼偏航气动力矩受纵横向运动影响很小，因此多旋翼不存在横航向耦合的荷兰滚模态，取而代之的是横向Phugoid模态。
3) 多旋翼模态稳定性差，振荡频率高，特别是Phugoid模态的发散速度快，使得手动飞行非常困难。
3 旋转自由度必要性 为了明确旋翼旋转自由度动态特性对于多旋翼飞行器建模分析的重要程度，本节利用线化模型，基于弱耦合系统理论进行计算分析，分别对无增稳状态与有增稳状态情况进行了讨论。
3.1 无增稳状态 对于第1节中获得的四旋翼线化模型，可以按机体运动自由度与旋翼旋转自由度划分为2个子系统，即
 | (10) |
式中：Ab和AΩ分别为机体运动自由度和旋翼旋转自由度对各自状态量变化率的影响；AbΩ为旋翼旋转自由度对机体运动自由度状态量变化率的影响；AΩb为机体运动自由度对旋翼旋转自由度状态量变化率的影响；BΩ为操纵量对旋翼旋转自由度状态量变化率的影响。
若上述系统满足弱耦合系统条件，则分析机体运动自由度时可对旋翼旋转自由度做准定常假设，即旋翼旋转动态特性对外部扰动与控制的响应速度很快，忽略其动态特性对响应较慢的刚体运动自由度影响较小。弱耦合系统判据包括[16-17]：
 | (11) |
 | (12) |
式中：λb和λΩ分别为两独立子系统(Ab，AΩ)的所有特征值；l为Ab阶数；EbΩ和EΩb分别为AbΩ和AΩb中绝对值最大的元素(即两系统间的最大耦合项)。
代入线性化计算结果，在悬停状态下，对于算例四旋翼无人飞行器：
 | (13) |
通过上述计算可知，对于无增稳悬停状态四旋翼无人飞行器，其旋翼旋转运动自由度的动态特性对机身运动自由度影响较小，可以做准定常处理，其全系统模态与近似系统模态对比如表 3所示。由对比可知，忽略旋翼旋转动态特性对系统模态的影响很小。
表 3 无增稳全系统与近似系统部分模态对比 Table 3 Comparison of some modes between entire system and approximate system without stability augmentation
| 模态 | 特征值 | |
| 全系统 | 近似系统 | |
| 俯仰收敛 | -1.725 9 | -1.698 9 |
| Phugoid | 0.553 9+1.266i | 0.550 6+1.259i |
| 沉浮收敛 | -0.513 7 | -0.500 2 |
| 螺旋 | -0.033 52 | -0.033 52 |
表选项
3.2 有增稳状态 由于多旋翼飞行器自身稳定性的不足，一般无法实现无增稳状态下的飞行，因此，本节将对增稳状态下的四旋翼特性进行分析。
对于算例四旋翼无人飞行器，在增稳控制模式下采用了简单的角度与角速率反馈。例如，对于俯仰通道：
 | (14) |
式中：Us为俯仰增稳控制输出；Kq为俯仰角速率反馈系数；q为俯仰角速率；KΘ为俯仰角反馈系数；Θ为俯仰角。
写作矩阵形式的增稳控制量为
 | (15) |
式中：Kb为刚体运动自由度对应控制系数矩阵。由于旋翼转速为开环控制，旋转自由度对应控制系数矩阵为0。
将式(15) 代入式(10)，可得增稳状态下的等效系统稳定性导数矩阵为
 | (16) |
由式(16) 可知，在增稳状态下，控制系统的反馈操纵不会影响机身刚体运动与旋翼旋转运动自身特性，即判据C1数值不变。而控制系统将改变两子系统间的耦合响应项。对于算例无人飞行器：
 | (17) |
由此可见，在有飞控增稳情况下，由于机身刚体运动与旋翼旋转自由度间的耦合作用，算例四旋翼无人飞行器不再满足弱耦合系统判据，其旋翼旋转自由度对整机动态特性有着较大影响，其全系统模态与近似系统模态对比如表 4所示，在建模仿真分析中应当考虑旋翼旋转动态特性。
表 4 有增稳全系统与近似系统部分模态对比 Table 4 Comparison of some modes between entire system and approximate system with stability augmentation
| 模态 | 特征值 | |
| 全系统 | 近似系统 | |
| 俯仰收敛 | -2.168+0.901 6i | -1.978 |
| Phugoid | -0.094 19+1.540i | -0.095 52+1.283i |
| 沉浮收敛 | -2.209+0.299 2i | -1.239 |
| 螺旋 | -0.117 7 | -0.115 4 |
表选项
为了验证上述结果，对算例四旋翼不同线化模型在前飞速度扰动(du=1 m/s)情况下的系统响应进行了模拟，其俯仰运动响应如图 2所示。从模拟结果可以发现，在无增稳状态下，旋翼旋转自由度对机身运动响应影响很小，旋转动态特性对模拟结果的影响可以忽略。然而，在有增稳情况下，若忽略旋转动态特性则会产生较大误差，影响系统响应模拟精度。同时，图 2也显示了多旋翼飞行器自身稳定性差，而通过控制系统可以显著改善其稳定性。
 |
| 图 2 线化模型扰动响应 Fig. 2 Responses of linearized mode to state perturbation |
| 图选项 |
4 俯仰模态分析 对于悬停状态的无增稳多旋翼飞行器，其俯仰与横向Phugoid模态稳定性差，模态发散速度快，与传统直升机有着较大区别，本节将针对俯仰Phugoid模态进行分析。
旋翼飞行器悬停状态下Phugoid模态分析可以简化为只考虑[u, Θ, q]3个自由度的动力学方程[16]：
 | (18) |
式中：g为重力加速度；Xu、Mu和Mq为气动导数。
飞行器Phugoid模态特性主要由Xu、Mu与Mq 3个气动导数决定。算例的相关气动导数与模态如表 5所示。
表 5 算例俯仰模态与稳定性导数 Table 5 Pitch mode and stability derivatives of study case
| 气动导数/模态 | 特征值 | |
| 四旋翼无人飞行器 | Bo-105 | |
| Xu | -0.032 71 | -0.024 29 |
| Mu | 0.327 5 | 0.130 5 |
| Mq | -0.564 9 | -4.763 |
| Phugoid模态 | 0.550 6+1.259i | 0.015 48+0.514 4i |
表选项
对于传统旋翼飞行器，其Phugoid模态特性估算公式为[16]
 | (19) |
式中：ωp为飞行器Phugoid模态频率；ζp为飞行器Phugoid模态阻尼。
虽然估算公式(19) 对多旋翼飞行器的适用性还需要证明，但其一定程度上提供了多旋翼飞行器模态特性变化的物理背景。由于Mu的增大与Mq的减小，导致多旋翼飞行器Phugoid模态的频率与负阻尼(即不稳定性)迅速增加。
对于Xu导数，其主要来源是飞行器水平飞行速度变化时带来的机身与旋翼阻力变化，该导数主要取决于飞行器形阻系数，因此多旋翼飞行器与传统旋翼飞行器在该气动导数数值上差别不大。同时，由于Mu/Mq数值较大，多旋翼飞行器的Phugoid模态受到Xu导数的影响较小。
Mq表征的是飞行器在俯仰运动过程中产生的俯仰阻尼力矩大小。对于传统直升机，俯仰运动带来的阻尼力矩主要来源于旋翼桨叶的桨毂力矩。由于科氏力的存在，传统直升机在有俯仰角速率情况下，其桨叶会产生相反方向的随动挥舞，进而产生较大的俯仰阻尼力矩。对于多旋翼飞行器，由于其采用刚性的螺旋桨，在俯仰运动中不会产生上述俯仰力矩，其Mq导数的主要来源是旋翼拉力变化。
多旋翼飞行器进行俯仰运动时，会在各旋翼处产生垂直运动速度。该速度等效于旋翼沿旋转轴方向的来流速度变化，进而改变旋翼转速与诱导速度，引起旋翼拉力变化。该拉力变化产生的俯仰力矩效果相对旋翼挥舞产生的俯仰力矩弱很多，因此，算例多旋翼飞行器的Mq远小于Bo-105直升机。
Mu导数取决于飞行器在前飞速度变化时所产生的俯仰力矩。对于传统直升机，在前飞速度改变时，由于来流条件变化，会引起桨叶吹风挥舞。在前飞速度增加时，会产生桨盘后倾，形成抬头力矩。对于算例Bo-105直升机，其桨叶挥舞对Mu导数的贡献为
 | (20) |
式中：βc1为桨叶俯仰挥舞系数。
对于多旋翼飞行器，由于不存在挥舞引起的桨毂力矩，其Mu主要由旋翼气动力矩的变化提供。对于算例四旋翼无人飞行器，Mu的主要贡献来源于前飞速度变化时引起的非均匀入流：
 | (21) |
式中：vic(i)为第i个旋翼纵向诱导速度。
在旋翼飞行器前飞过程中，由于来流影响，旋翼桨盘平面内的入流分布是非均匀的。根据Pitt-Peters入流模型，前飞速度将造成桨盘前部入流减小，后部入流增大, 如图 3所示。该入流变化等效于增加了桨叶叶素在90°~270°方位角范围内的迎角，减小了叶素在-90°~90°方位角范围内的迎角。因此，桨盘前部产生的升力较大，后部较小，形成一个抬头气动力矩。
 |
| 图 3 前飞状态桨盘入流分布示意图 Fig. 3 Sketch map of inflow distribution on rotor disk in forward flight |
| 图选项 |
对于多旋翼飞行器的刚性旋翼，该力矩将直接传递到飞行器机体，形成抬头力矩。而对于传统旋翼，该力矩将引起桨叶的挥舞运动，进而平衡气动力矩作用。同时，由于挥舞运动的滞后，纵向入流的变化将引起横向挥舞运动，对机体俯仰力矩影响微弱。因此，虽然Bo-105算例模型中也采用了Pitt-Peters入流模型，但其计算结果显示，旋翼吹风挥舞运动是Mu导数的主要来源。
5 结论 本文基于传统直升机飞行力学模型，建立了包含动态入流模型与旋翼旋转自由度的多旋翼飞行器飞行力学模型，并基于此模型对四旋翼无人飞行器算例进行了模态特性与建模考虑因素重要性分析。
1) 在无飞控增稳情况下，多旋翼飞行器模态稳定性显著弱于传统直升机，并且其悬停模态分布与传统直升机有所不同，即在纵横向上都有Phugoid模态，且不存在荷兰滚模态。
2) 多旋翼飞行器旋翼旋转自由度对于无增稳状态下的多旋翼飞行力学特性影响很小，可以做准定常处理。
3) 在增稳飞行状态下，由于飞控系统的操纵，多旋翼飞行器机体运动自由度与旋转自由度间的耦合加剧，两系统无法独立开来，在建模与分析中有必要考虑旋翼旋转动态特性。
4) 刚性桨叶的采用明显改变了多旋翼飞行器纵横向力矩产生机制与影响因素。旋翼入流分布情况(动态入流模型)对旋翼气动力矩有很大影响。同时，多旋翼飞行器旋翼不存在挥舞引起的桨毂力矩，桨盘上的气动俯仰、滚转力矩对飞行器模态特性影响较大。上述因素都需要在建模分析中予以考虑。
利用本文建立的飞行力学模型，在后续研究工作中可以进一步对多旋翼飞行器飞行力学特性与设计趋势进行更加详细的探讨。一方面，可以系统地分析多旋翼飞行器相对传统旋翼飞行器气动导数上的特点与成因；另一方面，将会针对多旋翼飞行力学特性进行设计参数敏感性分析，理解不同设计参数对多旋翼动态特性的影响，明确大型多旋翼飞行器设计趋势。
参考文献
| [1] | 王史春. 四旋翼飞行器力学模型与控制系统设计[J].中北大学学报(自然科学版), 2014, 35(2): 218–224. WANG S C. Design of mechanical model and control system of quadrotor[J].Journal of North University of China(Natural Science Edition), 2014, 35(2): 218–224.(in Chinese) |
| [2] | 杨成顺. 多旋翼飞行器建模与飞行控制技术研究[D]. 南京: 南京航空航天大学, 2013. YANG C S.Research on modeling and flight controltechnology of multi-rotor aircraft[D].Nanjing:Nanjing University of Aeronautics and Astronautics, 2013(in Chinese). |
| [3] | 韩志凤, 李荣冰, 刘建业, 等. 小型四旋翼飞行器动力学模型优化[J].控制工程, 2013(S1): 158–162. HAN Z F, LI R B, LIU J Y, et al. Model optimization of small quad-rotor in the case of center of mass offset[J].Control Engineering of China, 2013(S1): 158–162.(in Chinese) |
| [4] | 李俊, 李运堂. 四旋翼飞行器的动力学建模及PID控制[J].辽宁工程技术大学学报(自然科学版), 2012, 31(1): 114–117. LI J, LI Y T. Modeling and PID control for a quadrotor[J].Journal of Liaoning Technical University(Natural Science), 2012, 31(1): 114–117.(in Chinese) |
| [5] | 聂博文, 马宏绪, 王剑, 等. 微小型四旋翼飞行器的研究现状与关键技术[J].电光与控制, 2007, 14(6): 113–117. NIE B W, MA H X, WANG J, et al. Study on actualities and critical technologies of micro/mini quadrotor[J].Electronics Optics & Control, 2007, 14(6): 113–117.(in Chinese) |
| [6] | BOUABDALLAH S, MURRIERI P, SIEGWART R.Design and control of an indoor micro quadrotor[C]//Proceedings of IEEE International Conference on Robotics and Automation.Piscataway, NJ:IEEE Press, 2004:4393-4398. |
| [7] | MOKHTARI A, BENALLEGUE A.Dynamic feedback controller of Euler angles and wind parameters estimation for a quadrotor unmanned aerial vehicle[C]//Proceedings of IEEE International Conference on Robotics and Automation.Piscataway, NJ:IEEE Press, 2004:2359-2366. |
| [8] | MISTLER V, BENALLEGUE A, M'SIRDI N K.Exact linearization and noninteracting control of a 4 rotors helicopter via dynamic feedback[C]//Proceedings of 10th IEEE International Workshop on Robot and Human Interactive Communication.Piscataway, NJ:IEEE Press, 2001:586-593. |
| [9] | 曾小勇, 彭辉, 吴军. 四旋翼飞行器的建模与姿态控制[J].中南大学学报(自然科学版), 2013, 44(9): 3693–3700. ZENG X Y, PENG H, WU J. Modeling and attitude control for a quad-rotor aircraft[J].Journal of Central South University(Science and Technology), 2013, 44(9): 3693–3700.(in Chinese) |
| [10] | NAIDOO Y, STOPFORTH R, BRIGHT G. Quad-rotor unmanned aerial vehicle helicopter modelling & control[J].International Journal of Advanced Robotic Systems, 2011, 8(4): 139–149. |
| [11] | PROUTY R. Helicopter aerodynamics[M].Lebanon: Eagle Eye Solutions, LLC, 2007. |
| [12] | LEISHMAN J G. Principles of helicopter aerodynamics[M].Cambridge: Cambridge University Press, 2006. |
| [13] | PETERS D A, NINH H. Dynamic inflow for practical applications[J].Journal of the American Helicopter Society, 1988, 33(4): 64–68.DOI:10.4050/JAHS.33.64 |
| [14] | 姜军, 齐俊桐, 韩建达. 面向机动飞行的多旋翼飞行器设计和建模与控制[J].科学通报, 2013, 58(S2): 135–144. JIANG J, QI J T, HAN J D. Aggressive maneuver oriented multi-rotor aircraft design, modeling and control[J].Chinese Science Bulletin, 2013, 58(S2): 135–144.(in Chinese) |
| [15] | WAN J. Ornicopter multidisciplinary analyses and conceptual design[D].Delft:Delft University of Technology, 2014. |
| [16] | PADFIELD G D. Helicopter flight dynamics[M].Chichester: John Wiley & Sons, Ltd., 2007. |
| [17] | MILNE R D. The analysis of weakly coupled dynamical systems[J].International Journal of Control, 1965, 2(2): 171–199.DOI:10.1080/00207176508905535 |
