CN102880182A - 一种存在网络随机延迟的微小型无人飞行器控制方法 - Google Patents

Abstract

本发明公开了一种存在网络随机延迟的微小型无人飞行器控制方法，属于飞行控制技术领域，包括建立无人飞行器系统模型并设计控制律，存在网络随机延迟的系统随机鲁棒性分析，确定网络更新周期、优化指标与指标权值，利用随机鲁棒设计方法设计鲁棒控制律，闭环六自由度非线性蒙特卡洛仿真验证。其中步骤二中的网络随机延迟包括传感器观测延迟、执行器控制延迟以及观测与控制叠加的混合延迟。本发明解决了传统的线性二次型调节器控制对通信质量要求较高和鲁棒性较差的缺点，将基于随机鲁棒分析与设计的控制方法引入到基于网络的无人飞行器的控制当中，降低了对无人机编队无线数据链更新率的要求，进而能够增强无人机群编队的鲁棒性。

Description

一种存在网络随机延迟的微小型无人飞行器控制方法

技术领域

本发明涉及存在网络随机延迟的微小型无人飞行器控制，属于飞行控制技术领域，具体涉及一种存在网络随机延迟的微小型无人飞行器控制方法。

背景技术

当前已经有多达三十多个国家投入大量人力和财力从事无人机的研究和生产。经过二十年的发展，该项技术已经比较成熟，在军民各个领域发挥着作用，尽管如此，单架无人机在遂行任务时存在着一些问题，例如单架无人机可能受到传感器的数量限制，不能从多角度全方位的对目标区域进行观察，面临大面积搜索任务时，不能有效的覆盖整个搜索区域；如果执行的是救援任务，单架无人机在载荷方面受到限制，往往影响整个救援的效能，带来更大损失，另外，一旦单架无人机出现故障，必须立即中断任务返回，可能会延误救援时机。

针对单架无人机的上述缺点，近些年提出了编队飞行控制的概念并取得了一定得研究成果，其中网络飞行控制系统设计是编队飞行控制是否成功的决定性因素，目前无人机群编队控制主要采用集中式控制结构：将编队成员分成一个长机和若干僚机，长机将自身的制导回路信息通过无线通信网络发送给僚机，僚机接收到长机的信息后与自身的制导回路信息比较，按照某种队形结构计算自身的制导指令进而与长机以一定队形结构飞行。由于无人机制导回路的惯性大，带宽较窄，而按照最优随机控制的观点，网络延迟往往可以等效为频率相对较高的白噪声，因此集中式控制体制下的网络延迟对编队队形控制的影响不大。近些年来，随着分布式编队控制的发展，要求编队成员不仅要在制导回路进行协调，甚至是控制回路也要进行协调，这就对网络飞行控制系统设计提出了较高的要求，由于控制回路的带宽相对制导回路宽，很有可能将网络延迟近似的噪声纳入到带宽之内，故网络延迟可能导致飞控系统控制品质恶化。

网络飞行控制系统设计的概念是：针对无人机编队飞行过程中的无线通信网络延迟、丢包和时序错乱情况，设计飞行控制系统，满足系统对超调量和调节时间等性能指标的要求。网络飞行控制系统设计方法的研究最早可追溯到上世纪九十年代，经过二十多年的发展形成了以下三个主要分支：

最优随机控制方法：该方法首先将网络延迟考虑为系统噪声，并假定该延迟小于一个采样周期并服从某种分布，采用线性随机系统模型来描述网络控制系统的随机延迟特性，将服从某种分布的网络随机延迟对系统的影响转化为LQG（Linear Quadratic Gaussian）问题，利用LQR（Linear Quadratic Regulator）和Kalman滤波的分离原理设计控制器。该方法的缺陷是实际系统控制中采用线性随机系统模型来描述网络控制系统的随机延迟特性，进而设计出的控制器鲁棒性较差。

模糊控制方法：1965年美国的L.Zadeh教授提出了模糊逻辑的概念，20世纪90年代在日本电气控制中得到了广泛应用。模糊控制方法一般应用在对象模型不确定、传统控制方法难以奏效的情况下，但是对模糊控制的严格数学分析方法并没有构建，其适用的问题也没有严格界定，另外在状态参数较多的情况下，模糊隶属度参数的选择存在困难。

鲁棒控制方法：鲁棒控制方法是为了解决现实工程中很难获得准确的被控对象模型的问题而发展起来的，该方法首先假设网络延迟存在上下界，从而将网络延迟看做是系统的一个不确定因素，同时综合考虑被控对象本身的不确定性，利用回路成型（Loop Shaping）、线性矩阵不等式（Linear Matrix Inequation）等传统鲁棒控制器设计方法设计控制律。由于没有对网络延迟分布进行描述，导致该方法设计的控制律往往过于保守。

除了上述的三种网络飞行控制系统设计方法以外，网络飞行控制系统设计方法还包括一些其他的控制方法，如网络摄动法、状态增广法和缓冲队列法等，这些方法计算复杂，设计结果保守，相对上面介绍的三种方法应用较少。

发明内容

本发明的目的是为了解决上述问题，针对无人飞行器编队飞行过程中不可避免的无线自组织网络随机延迟问题，以及现有网络飞行控制系统设计方法的过分保守性的缺点，引入了随机鲁棒分析与设计方法，解决了传统的线性二次型调节器（LQR）控制对通信质量要求较高和鲁棒性较差的缺点，降低了对编队无线数据链更新率的要求，使飞行器在存在较大网络随机延迟的情况下控制品质仍能满足要求。

一种存在网络随机延迟的微小型无人飞行器控制方法，包括以下几个步骤：

步骤一：建立无人飞行器系统模型并确定控制律；

步骤二：存在网络随机延迟的系统随机鲁棒性分析；

步骤三：确定网络更新周期、优化指标与指标权值；

步骤四：利用随机鲁棒设计方法设计鲁棒控制律；

步骤五：闭环六自由度非线性蒙特卡洛仿真验证；

本发明的优点在于：

（1）本发明在一定程度上缓解了网络延迟导致的飞行控制系统品质恶化，能够增强分布式系统的稳定性；其次控制系统的随机鲁棒分析与设计针对的网络传输模式是事件驱动的，这种传输模式有以下两个优点，一是降低了数据链设计的难度与成本，二是能较好的与现有民用网络协议兼容。

（2）将网络随机时延全部按照上界等效进而采用传统鲁棒控制方法设计的控制律的方法能够确保系统的鲁棒稳定性，但鲁棒性能一般不好，相对于时延按上界等效的鲁棒控制方法，本发明能够更好的对系统的鲁棒稳定性和鲁棒性能进行折中；

（3）在分析控制系统的同时，能够同时给出对分布式无人机编队自组织网络更新率的最低要求，对分布式无人机编队数据链的设计有一部分借鉴意义；

（4）本发明没有繁琐的数学推导，而是充分利用了计算机高速处理数据的能力，相对于以往的针对线性时滞系统依靠数学推导获得控制律的方法，该方法在线性系统LQR方法的基础上直接对非线性系统利用蒙特卡洛仿真与现代优化算法两种随机方法进行随机鲁棒分析与设计，工程性与可操作性较强。

附图说明

图1：网络飞行控制系统的结构示意图；

图2：基于事件驱动的均值μ_d为1000ms，方差σ_d为500ms的正态分布网络延迟效果图；

图3：随机鲁棒分析与设计方法的流程图；

图4：利用LQR方法设计的控制律对无网络延迟的非线性系统的控制效果图；

图5：均值μ_d为1000ms，方差σ_d为500ms的网络随机延迟蒙特卡洛仿真结果图；

图6：均值μ_d为500ms，方差σ_d为250ms的网络随机延迟蒙特卡洛仿真结果图；

图7：等效上界鲁棒控制设计保守性原因说明示意图；

图8：标准粒子群算法寻优代价收敛曲线；

图9：控制随机延迟、观测随机延迟同时存在的网络延迟效果图；

图10：利用随机鲁棒分析与设计方法设计后的控制律对无人飞行器控制的蒙特卡洛仿真验证图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种存在网络随机延迟的微小型无人飞行器控制方法，流程如图3所示，包括以下几个步骤：

步骤一：建立无人飞行器系统模型并确定控制律；

具体为：

（1）根据风洞吹风获得无人飞行器动力学参数和物理参数；

无人飞行器动力学参数和物理参数可以根据实际风洞吹风得到，采用英美坐标系，具体为：

①获取纵向力和力矩系数：包括升力系数C_L0、C_Lα、C_Lq、

阻力系数C_D0、C_Dα、C_Dq、俯仰力矩系数C_m0、C_mα、C_mq、

其中，C_L0为攻角为0度时的升力系数，C_Lα为升力关于攻角的升力系数，C_Lq为升力关于俯仰角速度的升力系数，为升力关于升降舵的升力系数，C_D0为攻角为0度时的阻力系数，C_Dα为阻力关于攻角的阻力系数，C_Dq为阻力关于俯仰角速度的阻力系数，

为阻力关于升降舵的阻力系数，C_m0为攻角为0度时的俯仰力矩系数，C_mα为俯仰力矩关于攻角的力矩系数，C_mq为俯仰力矩关于俯仰角速度的力矩系数、

为俯仰力矩关于升降舵的力矩系数；

②获取横侧向力和力矩系数：包括侧力系数C_Yβ、C_Yp、C_Yr、

滚转力矩系数C_lβ、C_lp、C_lr、

偏航力矩系数C_nβ、C_np、C_nr、

其中，C_Yβ为侧力关于侧滑角的侧力系数，C_Yp为侧力关于滚转角速度的侧力系数，C_Yr为侧力关于偏航角速度的侧力系数，为侧力关于副翼的侧力系数，

为侧力关于方向舵的侧力系数，C_lβ为滚转力矩关于侧滑角的力矩系数，C_lp为滚转力矩关于滚转角速度的力矩系数，C_lr为滚转力矩关于偏航角速度的力矩系数，

为滚转力矩关于副翼的力矩系数，

为滚转力矩关于方向舵的力矩系数，C_nβ为偏航力矩关于侧滑角的力矩系数，C_np为偏航力矩关于滚转角速度的力矩系数，C_nr为偏航力矩关于偏航角速度的力矩系数，

为偏航力矩关于副翼的力矩系数，

为偏航力矩关于方向舵的力矩系数；

③获取无人飞行器的质量m（kg）、平均几何弦长C_A（m）、翼展b（m）、飞行器参考面积S_w（m²）、X轴转动惯量I_X（kg·m²）、Y轴转动惯量I_Y（kg·m²）、Z轴转动惯量I_Z（kg·m²）、惯量积I_XZ（kg·m²）、舵机时间常数T_δ(s)、舵机放大系数K_δ、发动机时间常数T_t(s)、发动机放大系数K_t、怠速推力t_A（N）、最大推力t_max（N）和最小推力t_min（N）。

（2）建立非线性六自由度动力学与运动学方程；

选取无人机飞行状态向量为：

$\stackrel{\‾}{X}\left(t\right)={\left[\begin{array}{ccccccccccccccccc}V& \mathrm{\α}& \mathrm{\β}& p& q& r& \mathrm{\φ}& \mathrm{\θ}& \mathrm{\ψ}& x& y& h& m& {\mathrm{\δ}}_a& {\mathrm{\δ}}_e& {\mathrm{\δ}}_r& {\mathrm{\δ}}_t\end{array}\right]}^T$

其中，V表示速度、α表示攻角、β表示侧滑角、p表示滚转角速度、q表示俯仰角速度、r表示偏航角速度、φ表示滚转角、θ表示俯仰角、ψ表示偏航角、x表示东向位置、y表示南向位置、h表示高度、m表示质量、δ_a表示副翼偏角、δ_e表示升降舵偏角、δ_r表示方向舵偏角、δ_t表示发动机推力。

建立无人飞行器非线性六自由度运动学与动力学方程如下：

其中：

表示速度的导数、

表示攻角的导数、

表示侧滑角的导数、

表示滚转角速度的导数、

俯仰角速度的导数、

表示偏航角速度的导数、

表示滚转角的导数、

表示俯仰角的导数、

表示偏航角的导数、

表示东向位置导数、

表示南向位置的导数、

表示高度的导数、

表示质量消耗率、

表示副翼偏角导数、

表示升降舵偏角的导数、

表示方向舵偏角的导数、

表示发动机推力的导数、

表示副翼偏角指令、

表示升降舵偏角指令、

表示方向舵偏角指令、

表示油门指令、①为航迹倾角、为航迹偏角、K_m为发动机耗油率、δ_T为发动机推力、T_δ为舵机时间常数、T_t为发动机时间常数、K_δ为舵机放大系数、K_t为发动机放大系数，式（1）简记为：

其中

为飞行状态向量的导数， $u\left(t\right)=[{\mathrm{\δ}}_a^{*},{\mathrm{\δ}}_e^{*},{\mathrm{\δ}}_r^{*},{\mathrm{\δ}}_t^{*}{]}^T$ 为飞行控制向量。

$\left\{\begin{array}{c}{G}_x=m\·g\·(-\cos \left(\mathrm{\α}\right)\·\cos \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\θ}\right)+\sin \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right)+\sin \left(\mathrm{\α}\right)\·\cos \left(\mathrm{\β}\right)\·\cos \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right))\\ G_y=m\·g\·(\cos \left(\mathrm{\α}\right)\·\sin \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\θ}\right)+\cos \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right)-\sin \left(\mathrm{\α}\right)\·\sin \left(\mathrm{\β}\right)\·\cos \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right))\\ G_z=m\·g\·(\sin \left(\mathrm{\α}\right)\·\sin \left(\mathrm{\θ}\right)+\cos \left(\mathrm{\α}\right)\·\cos \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right))\\ \mathrm{\Σ}=I_X\·I_Z-I_{\mathrm{XZ}}\·I_{\mathrm{XZ}}\\ c_1=((I_Y-I_Z)\·I_Z-I_{\mathrm{XZ}}\·I_{\mathrm{XZ}})/\mathrm{\Σ}\\ c_2=(I_X-I_Y+I_Z)\·I_{\mathrm{XZ}}/\mathrm{\Σ}\\ c_3=I_Z/\mathrm{\Σ}\\ c_4=I_{\mathrm{XZ}}/\mathrm{\Σ}\\ c_5=(I_Z-I_X)/I_Y\\ c_6=I_{\mathrm{XZ}}/I_Y\\ c_7=1/I_Y\\ c_8=(I_X\·(I_X-I_Y)+I_{\mathrm{XZ}}\·I_{\mathrm{XZ}})/\mathrm{\Σ}\\ c_9=I_X/\mathrm{\Σ}\end{array}\right.---\left(2\right)$

其中：G_x为重力在气流坐标系x轴的投影，G_y为重力在气流坐标系y轴的投影，G_z为重力在气流坐标系z轴的投影，∑、c₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈为方便计算的中间量。I_X为X轴转动惯量、I_Z为Z轴转动惯量、I_Y为Y轴转动惯量、I_XZ为惯量积。

$\left\{\begin{array}{c}{C}_L=C_{L0}+C_{\mathrm{L\α}}\·\mathrm{\α}+C_{\mathrm{Lq}}\·q+C_{{\mathrm{L\δ}}_e}\·{\mathrm{\δ}}_e\\ C_D=C_{D0}+C_{\mathrm{D\α}}\·\mathrm{\α}+C_{\mathrm{Dq}}\·q+C_{D{\mathrm{\δ}}_e}\·{\mathrm{\δ}}_e\\ C_m=C_{m0}+C_{\mathrm{m\α}}\·\mathrm{\α}+C_{\mathrm{mq}}\·q+C_{m{\mathrm{\δ}}_e}\·{\mathrm{\δ}}_e\\ C_Y=C_{\mathrm{Y\β}}\·\mathrm{\β}+C_{\mathrm{Yp}}\·p+C_{\mathrm{Yr}}\·r+C_{{\mathrm{Y\δ}}_a}\·{\mathrm{\δ}}_a+C_{Y{\mathrm{\δ}}_r}\·{\mathrm{\δ}}_r\\ C_l=C_{\mathrm{l\β}}\·\mathrm{\β}+C_{\mathrm{lp}}\·p+C_{\mathrm{lr}}\·r+C_{{\mathrm{l\δ}}_a}\·{\mathrm{\δ}}_a+C_{l{\mathrm{\δ}}_r}\·{\mathrm{\δ}}_r\\ C_n=C_{\mathrm{n\β}}\·\mathrm{\β}+C_{\mathrm{np}}\·p+C_{\mathrm{nr}}\·r+C_{{\mathrm{n\δ}}_a}\·{\mathrm{\δ}}_a+C_{n{\mathrm{\δ}}_r}\·{\mathrm{\δ}}_r\\ L=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_L\\ D=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_D\\ M=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_Y\\ Y=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_Y\\ \stackrel{\‾}{L}=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_l\·b\\ N=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_n\·C_A\end{array}\right.---\left(3\right)$

其中：C_L为总升力系数、C_D为总阻力系数、C_Y为总侧力系数、C_m为总俯仰力矩系数、C_l为总滚转力矩系数、C_n为总偏航力矩系数、ρ为空气密度、L为总升力、D为总阻力、Y为总侧力、

为总滚转力矩、M为总俯仰力矩、N为总偏航力矩、S_w为飞行器参考面积、C_A为飞行器平均几何弦长、b为翼展。

（3）解耦线性化；

将状态向量[V α β p q r φ θ ψ x y h m δ_a δ_e δ_r δ_t]^T分成纵向和横侧向两组，分别为纵向[V α q θ h δ_e δ_t]^T，横侧向[β p r φ ψ δ_α δ_r]^T，无人机在飞行过程中的大部分航线属于直线段，横侧向控制的任务仅为沿航线飞行，所以本发明针对无人机直线段飞行时，假设无人飞行器处于定高稳定平飞状态，则必有升力-重力平衡、推力-阻力平衡以及俯仰力矩平衡，因此下面三个等式成立：

$\left\{\begin{array}{c}\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{L0}+C_{\mathrm{L\α}}\·\mathrm{\α}+C_{{\mathrm{L\δ}}_e}\·{\mathrm{\δ}}_e)=m\·g\\ \frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{D0}+C_{\mathrm{D\α}}\·\mathrm{\α}+C_{{\mathrm{D\δ}}_e}\·{\mathrm{\δ}}_e)=t_A\\ \frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{m0}+C_{\mathrm{m\α}}\·\mathrm{\α}+C_{m{\mathrm{\δ}}_e}\·{\mathrm{\δ}}_e)\·C_A=0\end{array}\right.---\left(4\right)$

解方程（4），得到无人飞行器的配平攻角α₀，配平升降舵偏角δ_e0，配平飞行速度V₀，在该状态点利用小扰动线性化原理可得无人飞行器纵向线性状态方程为：

$\left[\begin{array}{c}\stackrel{\·}{V}\\ \stackrel{\·}{\mathrm{\α}}\\ \stackrel{\·}{q}\\ \stackrel{\·}{\mathrm{\θ}}\\ \stackrel{\·}{h}\\ {\stackrel{\·}{\mathrm{\δ}}}_e\\ {\stackrel{\·}{\mathrm{\δ}}}_t\end{array}\right]=A_1\left[\begin{array}{c}V\\ \mathrm{\α}\\ q\\ \mathrm{\θ}\\ h\\ {\mathrm{\δ}}_e\\ {\mathrm{\δ}}_t\end{array}\right]+B_1\left[\begin{array}{c}{\mathrm{\δ}}_e^{*}\\ {\mathrm{\δ}}_t^{*}\end{array}\right]---\left(5\right)$

其中A₁为纵向线性系统状态矩阵，B₁为纵向线性系统控制矩阵；

无人飞行器横侧向线性状态方程为：

$\left[\begin{array}{c}\stackrel{\·}{\mathrm{\β}}\\ \stackrel{\·}{p}\\ \stackrel{\·}{r}\\ \stackrel{\·}{\mathrm{\φ}}\\ \stackrel{\·}{\mathrm{\ψ}}\\ {\stackrel{\·}{\mathrm{\δ}}}_a\\ {\stackrel{\·}{\mathrm{\δ}}}_r\end{array}\right]=A_2=\left[\begin{array}{c}\mathrm{\β}\\ p\\ r\\ \mathrm{\φ}\\ \mathrm{\ψ}\\ {\mathrm{\δ}}_a\\ {\mathrm{\δ}}_r\end{array}\right]+B_2\left[\begin{array}{c}{\mathrm{\δ}}_a^{*}\\ {\mathrm{\δ}}_r^{*}\end{array}\right]---\left(6\right)$

其中A₂为横侧向线性系统状态矩阵，B₂为横侧向线性系统控制矩阵；

（4）LQR设计控制律；

针对纵向：选取Q₁和R₁矩阵，采用线性二次型调节器设计纵向控制律K₁，K₁为2×7的矩阵。用k_ij表示K₁中的第i行第j列元素，为使控制系统结构尽量简化，速度反馈只加入发动机油门控制（k₂₂＝k₂₃＝k₂₄＝k₂₅＝0.0）中，而不加入升降舵控制（k₁₁＝0.0）中，无人飞行器舵机一般采用电动舵机，在飞行过程中无法准确测量舵面偏角，因此k₁₆＝k₁₇＝k₂₆＝k₂₇＝0.0，因为速度反馈不参与升降舵控制，而发动机油门调节只依赖速度，故应增加其他状态量对升降舵控制的反馈力度，最终选择以下控制器结构：

$\left\{\begin{array}{c}{\mathrm{\δ}}_e^{*}={\mathrm{\δ}}_{e0}-k_{12}\·(\mathrm{\α}-{\mathrm{\α}}_0)-k_{13}\·q-k_{14}\·(\mathrm{\θ}-{\mathrm{\α}}_0)-k_{15}\·(h-h^{*})\\ {\mathrm{\δ}}_t^{*}=t_A-k_{21}\·(V-V_0)\end{array}---\left(7\right)\right.$

其中：δ_e0为配平升降舵偏角，α₀为无人飞行器的配平攻角，h^*为高度指令，V₀为配平飞行速度。由于微小型飞行器由于结构较小无法安装精确的攻角传感器，故将攻角反馈换算为过载反馈，由于微小型飞行器在飞行过程中，机翼和水平安定面是主要升力面，升降舵面产生的升力相对较小，主要是提供俯仰力矩诱导主升力面发生变化，因此：

$\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{L0}+C_{\mathrm{L\α}}\·(\mathrm{\α}-{\mathrm{\α}}_0))+{\mathrm{\δ}}_t\·\sin \left(\mathrm{\α}\right)=(N_Z-1)\·m\·g---\left(8\right)$

其中：N_Z是无人飞行器的法向过载，C_L0相对C_Lα较小，可以忽略，δ_t·sin(α)相对

较小也可以忽略，代入数据得

$(\mathrm{\α}-{\mathrm{\α}}_0)=(N_Z-1)\·m\·g/(\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_{\mathrm{L\α}})---\left(9\right)$

联立式（7）、式（9）可得实际飞行中的纵向控制律：

$\left\{\begin{array}{c}{\mathrm{\δ}}_e^{*}={\mathrm{\δ}}_{e0}-k_1\·(N_z-1)-k_2\·q-k_3\·(\mathrm{\θ}-{\mathrm{\α}}_0)-k_4\·(h-h^{*})\\ {\mathrm{\δ}}_t^{*}=t_A-k_5\·(V-V_0)\end{array}---\left(10\right)\right.$

其中：k₁ k₂ k₃ k₄ k₅表示按照（7）和（9）换算后的最终应用的控制器反馈系数，

k₂＝k₁₃,k₃＝k₁₄,k₄＝k₁₅,k₅＝k₂₁。

针对横侧向：选取Q₂和R₂矩阵，采用线性二次型调节器设计纵向控制律K₂，K₂为2×7的矩阵。用

表示K₂中的第i行第j列元素，无人飞行器舵机一般采用电动舵机，在飞行过程中无法准确测量舵面偏角，

故K₂为2×5的矩阵。

实际飞行中的横侧向控制律：

$\left[\begin{array}{c}{\mathrm{\δ}}_a^{*}\\ {\mathrm{\δ}}_r^{*}\end{array}\right]=K_2\·\left[\begin{array}{c}\mathrm{\β}\\ p\\ r\\ \mathrm{\φ}\\ \mathrm{\ψ}\end{array}\right]---\left(11\right)$

步骤二：存在网络随机延迟的系统随机鲁棒性分析；

网络随机延迟的定义是：在传感器、执行器和控制器等多个节点通过网络交换数据时，由于网络带宽有限且网络中的数据流量变化不规则，不可避免的造成数据碰撞、多路传输、连接中断和网络拥塞等现象，因此出现的信息交换时间延迟。从延迟发生的位置可将网络随机延迟分为传感器观测延迟τ_o和执行器控制延迟τ_c两类，如图1所示，从延迟的类型上可将网络随机延迟分为常数延迟、相互独立的随机延迟和Markov特性的随机延迟。本发明主要采用相互独立的随机延迟作为网络飞行控制系统的主要延迟类型，假设传感器观测延迟和执行器控制延迟服从均值μ_d，方差σ_d的正态分布，即：

τ_o∈N(μ_d,σ_d),τ_c∈N(μ_d,σ_d)        （12）

如图2所示，图中圆圈位置代表延迟大小，横坐标代表延迟发生的时刻，由于控制系统采用事件驱动的方式处理数据，故延迟发生的时刻并不均匀，无人飞行器非线性系统方程由原来的 $\stackrel{\·}{\stackrel{\‾}{X}}\left(t\right)=f(\stackrel{\‾}{X}\left(t\right),u\left(t\right))$ 变为了：

$\stackrel{\·}{\stackrel{\‾}{X}}\left(t\right)=f(\stackrel{\‾}{X}(t-{\mathrm{\τ}}_o),u(t-{\mathrm{\τ}}_c))---\left(13\right)$

其中：

为飞行状态向量，为飞行状态向量的导数，

为加入传感器观测延迟的飞行状态向量，u(t-τ_c)为加入执行器控制延迟的飞行控制向量。网络延迟本身的随机性为随机鲁棒分析的应用提供了条件。对服从（12）式的网络随机延迟进行蒙特卡洛仿真考察（10）式控制律在网络延迟情况下对无人飞行器控制的鲁棒性，在超调量没有不稳定现象出现且不大于60%，此时有进一步优化的余地，将此时的μ_d作为步骤三设计网络更新率的依据，假设此时蒙特卡洛仿真确定的高度响应的平均超调量是σ，高度响应的平均调节时间是T_A。随机鲁棒分析过程本质是应用统计概率对由于无人飞行器编队网络随机延迟所引起的闭环不可接受行为（包括系统的不稳定特性和性能指标要求不满足）进行描述。

步骤三：确定网络更新周期、优化指标与指标权值；

根据步骤二的蒙特卡洛仿真结果可以确定，编队网络更新率为：

$r=\frac{1}{{\mathrm{\μ}}_d}---\left(14\right)$

其中：r为网络更新率，μ_d为传感器观测延迟和执行器控制延迟时间的均值，根据步骤二的蒙特卡洛仿真结果可以确定，超调量在网络存在不确定性的情况下超调量的方差相对调节时间的方差大，因此优化指标中超调量的权值归一化后应为0.8，调节时间的权值归一化后应为0.2，设定随机鲁棒设计的代价函数：

W＝ω_σ·σ+ω_T·T_A        （15）

其中σ是蒙特卡洛仿真高度响应的平均超调量，T_A是蒙特卡洛仿真高度响应的平均调节时间。

步骤四：利用随机鲁棒设计方法设计鲁棒控制律；

随机鲁棒设计过程包括控制器结构设计和现代优化算法两部分，控制器结构采用传统的PID控制，在步骤一中已叙述，现代优化算法采用粒子群算法，因为控制系统参数k₁,k₂,k₃,k₄,k₅变化是连续的，粒子群算法在处理多连续变量寻优相对遗传算法和蚁群算法具有优势。

采用标准粒子群算法求解该优化问题，问题的解对应于搜索空间中的一个粒子，每个粒子都有自己的位置和速度及一个由被优化函数决定的代价函数。各个粒子记忆、追随当前的最优粒子，每次迭代中，粒子通过跟随两个极值来更新自己的位置和速度：一个是粒子本身获得的最优解

一个是整个粒子群中全部粒子历代搜索获得的最优解

标准粒子群算法算法中，速度更新和位置更新如式(16)所示：

$\left\{\begin{array}{c}{v}_{\mathrm{id}}^{k+1}={\mathrm{\ωv}}_i^k+c_1{r}_1(p_{\mathrm{id}}^k-x_{\mathrm{id}}^k)+c_2{r}_2(g_{\mathrm{id}}^k-x_{\mathrm{id}}^k)\\ x_{\mathrm{id}}^{k+1}=x_{\mathrm{id}}^k+v_{\mathrm{id}}^{k+1}\end{array}---\left(16\right)\right.$

式中，i∈N(1,m),m为粒子群中粒子的个数，N代表整数；d∈N(1,n)，n为解向量的维数，k为迭代次数，c₁和c₂为学习因子，初始值为0.6，随迭代次数增加降为0.2；ω为惯量权重，初始值为0.9，随迭代次数增加降为0.4；r₁和r₂为[0,1]之间的随机数，判断粒子群中粒子优劣的标准是蒙特卡洛仿真获得的代价函数W，代价函数越小，在下一代该粒子存在的概率越大。

设定粒子最大的迭代次数，判断迭代是否到达最大迭代次数，如果到达，输出最好粒子代表的控制系统参数k₁,k₂,k₃,k₄,k₅，否则，继续进行迭代。

步骤五：闭环六自由度非线性蒙特卡洛仿真验证；

在同时存在传感器观测延迟τ_o和执行器控制延迟τ_c的条件下，且

τ_o∈N(μ_d,σ_d),τ_c∈N(μ_d,σ_d)

对步骤四中随机鲁棒设计出的控制律进行蒙特卡洛仿真验证，对比步骤二和步骤五的蒙特卡洛仿真结果可以看出，采用随机鲁棒设计后的飞行控制律较步骤一得出的控制律抗网络随机延迟能力增强，在相同传感器观测延迟τ_o和执行器控制延迟τ_c的条件下，最大超调量获得明显改善。所以采用本发明方法获得的控制律，相对于LQR设计出的控制律抗网络延迟的能力加强，把本发明方法设计出的k₁,k₂,k₃,k₄,k₅应用到无人机自动驾驶仪的纵向通道中，在无人机自主飞行过程中，能够在存在较大网络延迟的情况下维持无人机的正常飞行。

图3给出了本专利方法的总体流程图。本发明能够解决一类同时具有传感器观测延迟和执行器控制延迟的网络飞行控制律设计与优化问题，近些年来，随着分布式编队控制的发展，要求编队成员不仅要在制导回路进行协调，甚至是控制回路也要进行协调，这就对网络飞行控制系统设计提出了较高的要求，由于控制回路的带宽相对制导回路宽，很有可能将网络延迟近似的噪声纳入到带宽之内，故网络延迟可能导致飞控系统控制品质恶化，本发明在一定程度上缓解了网络延迟导致的飞行控制系统品质恶化，能够增强分布式系统的稳定性；其次控制系统的随机鲁棒分析与设计针对的网络传输模式是事件驱动的，这种传输模式有以下两个优点，一是降低了数据链设计的难度与成本，二是能较好的与现有民用网络协议兼容。

实施例：

步骤一：建立无人飞行器系统模型并确定控制律；

无人飞行器动力学参数和物理参数可以根据实际风洞吹风得到，采用英美坐标系，具体数值见下表：

（2）建立非线性六自由度动力学与运动学方程；

选取无人机飞行状态向量为：

$\stackrel{\‾}{X}\left(t\right)={\left[\begin{array}{ccccccccccccccccc}V& \mathrm{\α}& \mathrm{\β}& p& q& r& \mathrm{\φ}& \mathrm{\θ}& \mathrm{\ψ}& x& y& h& m& {\mathrm{\δ}}_a& {\mathrm{\δ}}_e& {\mathrm{\δ}}_r& {\mathrm{\δ}}_t\end{array}\right]}^T$

建立如公式（1）所示的无人机六自由度动力学与运动学方程；

（3）解耦线性化；

将状态向量[V α β p q r φ θ ψ x y h m δ_a δ_e δ_r δ_t]^T分成纵向和横侧向两组，分别为纵向[V α q θ h δ_e δ_t]^T，横侧向[β p  r φ ψ δ_a δ_r]^T解公式（4），得到无人飞行器的配平攻角α₀＝4.4821°，配平升降舵偏角δ_e0＝2.2736°，配平飞行速度V₀＝41.0214m/s，在该状态点利用小扰动线性化原理可得无人飞行器纵向线性系统状态方程为：

$\left[\begin{array}{c}\stackrel{\·}{V}\\ \stackrel{\·}{\mathrm{\α}}\\ \stackrel{\·}{q}\\ \stackrel{\·}{\mathrm{\θ}}\\ \stackrel{\·}{h}\\ {\stackrel{\·}{\mathrm{\δ}}}_e\\ {\stackrel{\·}{\mathrm{\δ}}}_t\end{array}\right]=\left[\begin{array}{ccccccc}-0.129607& -26.4232& 0& -9.81& 0& 2.40894& 0.0483686\\ -0.0109844& -5.57054& 1.00102& 0& 0& -0.321125& -0.000067\\ -0.0000225& -70.3391& -511.921& 0& 0& -54.0864& 0\\ 0& 0& 1& 0& 0& 0& 0\\ 0& -42& 0& 42& 0& 0& 0\\ 0& 0& 0& 0& 0& -20.0& 0\\ 0& 0& 0& 0& 0& 0& -5.0\end{array}\right]\left[\begin{array}{c}V\\ \mathrm{\α}\\ q\\ \mathrm{\θ}\\ h\\ {\mathrm{\δ}}_e\\ {\mathrm{\δ}}_t\end{array}\right]$

$+\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 20.0& 0\\ 0& 5.0\end{array}\right]\left[\begin{array}{c}{\mathrm{\δ}}_e^{*}\\ {\mathrm{\δ}}_t^{*}\end{array}\right]---\left(17\right)$

无人飞行器横侧向线性系统状态方程为：

$\left[\begin{array}{c}\stackrel{\·}{\mathrm{\β}}\\ \stackrel{\·}{p}\\ \stackrel{\·}{r}\\ \stackrel{\·}{\mathrm{\φ}}\\ \stackrel{\·}{\mathrm{\ψ}}\\ {\stackrel{\·}{\mathrm{\δ}}}_a\\ {\stackrel{\·}{\mathrm{\δ}}}_r\end{array}\right]=\left[\begin{array}{ccccccc}-0.522585& -1.99401& 0.967575& 0.233124& 0& -0.310635& 0.776925\\ -69.146& -380.061& 210.735& 0& 0& -99.6461& 26.6377\\ 7.91573& -10.785& -37.4416& 0& 0& -2.22708& -9.50665\\ 0& 1& 0.061951& 0& 0& 0& 0\\ 0& 0& 1.00192& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -20.0& 0\\ 0& 0& 0& 0& 0& 0& -20.0\end{array}\right]\left[\begin{array}{c}\mathrm{\β}\\ p\\ r\\ \mathrm{\φ}\\ \mathrm{\ψ}\\ {\mathrm{\δ}}_a\\ {\mathrm{\δ}}_r\end{array}\right]$

$+\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 0& 0\\ 20.0& 0\\ 0& 20.0\end{array}\right]\left[\begin{array}{c}{\mathrm{\δ}}_a^{*}\\ {\mathrm{\δ}}_r^{*}\end{array}\right]---\left(18\right)$

（4）LQR设计控制律；

针对纵向：选取Q₁和R₁矩阵如下：

$Q_1=\left[\begin{array}{ccccccc}0.01& & & & & & \\ & 131.3& & & & 0& \\ & & 820.82& & & & \\ & & & 131.3& & & \\ & & & & 0.0025& & \\ & 0& & & & 131.3& \\ & & & & & & 0.01\end{array}\right],$ $R_1=\left[\begin{array}{cc}131.3& 0\\ 0& 0.01\end{array}\right]---\left(19\right)$

采用线性二次型调节器设计纵向控制律：

$K_1=\left[\begin{array}{ccccccc}-0.00435& 0.133262& -0.016633& -2.4003& -0.00417& 0.443827& -0.00003\\ 0.46118& -4.6277& 0.100893& 56.3062& 0.147224& -0.0995& 0.417404\end{array}\right]---\left(20\right)$

根据公式（7）（8）（9）（10）的运算最终可得

$\left\{\begin{array}{c}{\mathrm{\δ}}_e^{*}={\mathrm{\δ}}_{e0}-k_1\·(N_z-1)-k_2\·q-k_3\·(\mathrm{\θ}-{\mathrm{\α}}_0)-k_4\·(h-h^{*})\\ {\mathrm{\δ}}_t^{*}=t_A-k_5\·(V-V_0)\end{array}\right.----\left(21\right)$

k₁＝-0.068,k₂＝0.02632,k₃＝2.4003,k₄＝0.006,k₅＝0.46305；

假设初始无人机飞行高度在300m，这时给高度指令h^*＝250m，控制计算机频率为50Hz，系统响应如图4所示，分别包括速度、攻角、俯仰角速度、俯仰角、高度、法向过载、升降舵偏角、油门、升降舵指令、油门指令的相应图，此时无人飞行器高度的超调量7.13%，定义调节时间为第一次到达h^*的时间，则调节时间为17.527s。

针对横侧向：选取Q₂和R₂矩阵如下：

$Q_2=\left[\begin{array}{ccccccc}131.3& & & & & & \\ & 820.82& & & & 0& \\ & & 820.82& & & & \\ & & & 131.3& & & \\ & & & & 131.3& & \\ & 0& & & & 131.3& \\ & & & & & & 131.3\end{array}\right],$ $R_2=\left[\begin{array}{cc}131.3& 0\\ 0& 131.3\end{array}\right]---\left(22\right)$

采用线性二次型调节器设计横侧向控制律：

$K_2=\left[\begin{array}{ccccccc}-0.172662& -0.0372412& -0.25231& -1.92896& -0.514471& 0.55969& -0.00405242\\ 1.29432& 0.0099507& -0.16617& 1.03541& 0.857508& -0.00405242& 0.511629\end{array}\right]---\left(23\right)$

无人飞行器舵机一般采用电动舵机，在飞行过程中无法准确测量舵面偏角，故：

$K_2=\left[\begin{array}{ccccc}-0.172662& -0.0372412& -0.26231& -1.92896& -0.514471\\ 1.29432& 0.0099507& -0.16617& 1.03541& 0.857508\end{array}\right]---\left(24\right)$

根据公式（11）可得实际飞行中的纵向控制律：

$\left[\begin{array}{c}{\mathrm{\δ}}_a^{*}\\ {\mathrm{\δ}}_r^{*}\end{array}\right]=\left[\begin{array}{ccccc}-0.172662& -0.0372412& -0.25231& -1.92896& -0.514471\\ 1.29432& 0.0099507& -0.16617& 1.03541& 0.857508\end{array}\right]\·\left[\begin{array}{c}\mathrm{\β}\\ p\\ r\\ \mathrm{\φ}\\ \mathrm{\ψ}\end{array}\right]---\left(25\right)$

步骤二：存在网络随机延迟的系统随机鲁棒性分析；

本发明主要采用相互独立的随机延迟作为网络飞行控制系统的主要延迟类型，假设传感器观测延迟和执行器控制延迟服从均值μ_d为1000ms，方差σ_d为500ms的正态分布，即：τ_o∈N(1000,500),τ_c∈N(1000,500)。对服从该正态分布的网络随机延迟进行蒙特卡洛仿真考察（21）式控制律在网络延迟情况下对无人飞行器控制的鲁棒性，仿真次数50次，仿真结果如图5所示。由图5可以得出如下结论：看出原有控制律的控制品质迅速下降到不可接受的程度，超调量均值大于80%，部分甚至超过100%，一般认为超调量超过20%是不满足无人机飞行品质的；而同时存在传感器观测延迟τ_o和执行器控制延迟τ_c情况下，不稳定的概率达到8%；

假设传感器观测延迟和执行器控制延迟服从均值μ_d为500ms，方差σ_d为250ms的正态分布，在该网络随机延迟情况下考察（21）式控制律对无人飞行器控制的鲁棒性，蒙特卡洛仿真次数为50，仿真结果如图6所示，对比图5和图6可以得出如下结论：

对于传感器观测延迟和执行器控制延迟同时存在，采用事件驱动方式的情况下，将网络延迟全部按照上界等效进而采用传统鲁棒控制方法设计的控制律太过保守。图6中在同时存在传感器观测延迟τ_o和执行器控制延迟τ_c的条件下，且τ_o∈N(500,250),τ_c∈N(500,250)，如果按照上界处理则等效的延迟分布τ_co∈N(1000,500)，由图5可以看出，这样的延迟分布τ_co无论等效在传感器一端还是在执行器一端，得到的控制超调量最大值均超过100%，而图6中相同控制律超调量最大值不超过60%。

原因如图7所示，图中空心圆圈代表执行器控制延迟时间，实心圆圈代表传感器观测延迟时间，在37.83s发生了一个执行器控制延迟670ms，38s发生了一个传感器观测延迟340ms，即在[37.83s,38.5s]执行器存在延迟，在[38s,38.34s]观测器存在延迟，而[38s,38.34s]是控制延迟遮掩观测延迟的部分，该部分属于可忽略部分，不应计入鲁棒控制的不确定性Δ进而增加不确定性的上界，控制延迟是造成飞行控制品质下降的直接原因，当控制延迟遮掩观测延迟的一部分时，可以认为这部分观测延迟对飞行控制品质的下降没有影响，所以网络延迟全部按照上界等效进而采用传统鲁棒控制方法设计的控制律太过保守。

图6传感器观测延迟和执行器控制延迟同时存在时相同控制律超调量最大值不超过60%，超调量平均值为40%，调节时间14.75s，没有不稳定现象出现，有进一步优化的余地，所以可为步骤三中代价函数选取提供参考。选择超调量平均值为40%，调节时间14.75s为计算代价函数权值的参考值。

步骤三：确定网络更新周期、优化指标与指标权值；

根据步骤二的蒙特卡洛仿真结果可以确定，编队网络更新率为：

$r=\frac{1}{{\mathrm{\μ}}_d}=2\mathrm{Hz}$

其中μ_d为传感器观测延迟和执行器控制延迟时间的均值取为500ms，根据步骤二的蒙特卡洛仿真结果可以确定，超调量在网络存在不确定性的情况下超调量的方差相对调节时间的方差大，因此优化指标中超调量的权值归一化后应为0.8，调节时间的权值归一化后应为0.2。传感器观测延迟和执行器控制延迟同时存在超调量均值40%，调节时间均值14.75s，因此公式（15）的权值为ω_σ＝0.8/0.4＝2,ω_T＝0.2/14.75＝0.0135。

步骤四：利用随机鲁棒设计方法设计鲁棒控制律；

控制器结构采用传统的PID控制，利用粒子群算法优化控制器参数，如式（16）所示,粒子群中粒子的个数m取为20，解向量的维数n为5，迭代次数k为20，对控制器参数k₁…k₅进行寻优，结果如图8所示，图8给出了标准粒子群算法的代价收敛曲线，进化8代后取得代价函数最小值，代价函数最小值W_min＝0.2344。对应的控制律为：

k₁＝-0.496721,k₂＝0.0370214,k₃＝2.79014,k₄＝0.00743004,k₅＝0.369074，此时最小超调量σ＝2.16091%，T_A＝17.039s。

步骤五：闭环六自由度非线性蒙特卡洛仿真验证；

对步骤四中随机鲁棒设计出的控制律在同时存在传感器观测延迟τ_o和执行器控制延迟τ_c的条件下，且τ_o∈N(500,250),τ_c∈N(500,250)时（如图9所示）进行蒙特卡洛仿真验证，仿真结果如图10所示，对比图6和图10可以看出，采用随机鲁棒设计后的飞行控制律较步骤一得出的控制律抗网络随机延迟能力增强，在相同传感器观测延迟τ_o和执行器控制延迟τ_c的条件下，最大超调量从60%下降到30%，超调量超过20%的概率为0.12。所以采用本发明方法获得的控制律，相对于LQR设计出的控制律抗网络延迟的能力加强。

Claims (2)

1.一种存在网络随机延迟的微小型无人飞行器控制方法，其特征在于，包括以下几个步骤：

步骤一：建立无人飞行器系统模型并确定控制律；

具体为：

（1）根据风洞吹风获得无人飞行器动力学参数和物理参数；

无人飞行器动力学参数和物理参数可以根据实际风洞吹风得到，采用英美坐标系，具体为：

①获取纵向力和力矩系数：包括升力系数C_L0、C_Lα、C_Lq、

阻力系数C_D0、C_Dα、C_Dq、

俯仰力矩系数C_m0、C_mα、C_mq、

其中，C_L0为攻角为0度时的升力系数，C_Lα为升力关于攻角的升力系数，C_Lq为升力关于俯仰角速度的升力系数，

为升力关于升降舵的升力系数，C_D0为攻角为0度时的阻力系数，C_Dα为阻力关于攻角的阻力系数，C_Dq为阻力关于俯仰角速度的阻力系数，

为阻力关于升降舵的阻力系数，C_m0为攻角为0度时的俯仰力矩系数，C_mα为俯仰力矩关于攻角的力矩系数，C_mq为俯仰力矩关于俯仰角速度的力矩系数、

为俯仰力矩关于升降舵的力矩系数；

②获取横侧向力和力矩系数：包括侧力系数C_Yβ、C_Yp、C_Yr、

滚转力矩系数C_1β、C_lp、C_lr、

偏航力矩系数C_nβ、C_np、C_nr、

其中，C_Yβ为侧力关于侧滑角的侧力系数，C_Yp为侧力关于滚转角速度的侧力系数，C_Yr为侧力关于偏航角速度的侧力系数，

为侧力关于副翼的侧力系数，为侧力关于方向舵的侧力系数，C_lβ为滚转力矩关于侧滑角的力矩系数，C_lp为滚转力矩关于滚转角速度的力矩系数，C_lr为滚转力矩关于偏航角速度的力矩系数，为滚转力矩关于副翼的力矩系数，为滚转力矩关于方向舵的力矩系数，C_nβ为偏航力矩关于侧滑角的力矩系数，C_np为偏航力矩关于滚转角速度的力矩系数，C_nr为偏航力矩关于偏航角速度的力矩系数，

为偏航力矩关于副翼的力矩系数，

为偏航力矩关于方向舵的力矩系数；

③获取无人飞行器的质量m、平均几何弦长C_A、翼展b、飞行器参考面积S_w、X轴转动惯量I_X、Y轴转动惯量I_Y、Z轴转动惯量I_Z、惯量积I_XZ、舵机时间常数T_δ(s)、舵机放大系数K_δ、发动机时间常数T_t(s)、发动机放大系数K_t、怠速推力t_A（N）、最大推力t_max（N）和最小推力t_min（N）；

（2）建立非线性六自由度动力学与运动学方程；

选取无人机飞行状态向量为：

$\stackrel{\‾}{X}\left(t\right)={\left[\begin{array}{ccccccccccccccccc}V& \mathrm{\α}& \mathrm{\β}& p& q& r& \mathrm{\φ}& \mathrm{\θ}& \mathrm{\ψ}& x& y& h& m& {\mathrm{\δ}}_a& {\mathrm{\δ}}_e& {\mathrm{\δ}}_r& {\mathrm{\δ}}_t\end{array}\right]}^T$

其中，V表示速度、α表示攻角、β表示侧滑角、p表示滚转角速度、q表示俯仰角速度、r表示偏航角速度、φ表示滚转角、θ表示俯仰角、ψ表示偏航角、x表示东向位置、y表示南向位置、h表示高度、m表示质量、δ_a表示副翼偏角、δ_e表示升降舵偏角、δ_r表示方向舵偏角、δ_t表示发动机推力；

建立无人飞行器非线性六自由度运动学与动力学方程如下：

其中：

表示速度的导数、

表示攻角的导数、

表示侧滑角的导数、

表示滚转角速度的导数、

俯仰角速度的导数、

表示偏航角速度的导数、

表示滚转角的导数、表示俯仰角的导数、

表示偏航角的导数、

表示东向位置导数、

表示南向位置的导数、

表示高度的导数、

表示质量消耗率、表示副翼偏角导数、

表示升降舵偏角的导数、

表示方向舵偏角的导数、

表示发动机推力的导数、

表示副翼偏角指令、

表示升降舵偏角指令、

表示方向舵偏角指令、

表示油门指令、μ为航迹倾角、

为航迹偏角、K_m为发动机耗油率、δ_T为发动机推力、T_δ为舵机时间常数、T_t为发动机时间常数、K_δ为舵机放大系数、K_t为发动机放大系数，式（1）简记为：

其中

为飞行状态向量的导数， $u\left(t\right)={[{\mathrm{\δ}}_a^{*},{\mathrm{\δ}}_e^{*},{\mathrm{\δ}}_r^{*},{\mathrm{\δ}}_t^{*}]}^T$ 为飞行控制向量；

$\left\{\begin{array}{c}{G}_x=m\·g\·(-\cos \left(\mathrm{\α}\right)\·\cos \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\θ}\right)+\sin \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right)+\sin \left(\mathrm{\α}\right)\·\cos \left(\mathrm{\β}\right)\·\cos \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right))\\ G_y=m\·g\·(\cos \left(\mathrm{\α}\right)\·\sin \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\θ}\right)+\cos \left(\mathrm{\β}\right)\·\sin \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right)-\sin \left(\mathrm{\α}\right)\·\sin \left(\mathrm{\β}\right)\·\cos \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right))\\ G_z=m\·g\·(\sin \left(\mathrm{\α}\right)\·\sin \left(\mathrm{\θ}\right)+\cos \left(\mathrm{\α}\right)\·\cos \left(\mathrm{\φ}\right)\·\cos \left(\mathrm{\θ}\right))\\ \mathrm{\Σ}=I_X\·I_Z-I_{\mathrm{XZ}}\·I_{\mathrm{XZ}}\\ c_1=((I_Y-I_Z)\·I_Z-I_{\mathrm{XZ}}\·I_{\mathrm{XZ}})/\mathrm{\Σ}\\ c_2=(I_X-I_Y+I_Z)\·I_{\mathrm{XZ}}/\mathrm{\Σ}\\ c_3=I_Z/\mathrm{\Σ}\\ c_4=I_{\mathrm{XZ}}/\mathrm{\Σ}\\ c_5=(I_Z-I_X)/I_Y\\ c_6=I_{\mathrm{XZ}}/I_Y\\ c_7=1/I_Y\\ c_8=(I_X\·(I_X-I_Y)+I_{\mathrm{XZ}}\·I_{\mathrm{XZ}})/\mathrm{\Σ}\\ c_9=I_X/\mathrm{\Σ}\end{array}\right.---\left(2\right)$

其中：G_x为重力在气流坐标系x轴的投影，G_y为重力在气流坐标系y轴的投影，G_z为重力在气流坐标系z轴的投影，∑、c₁、c₂、c₃、c₄、c₅、c₆、c₇、c₈为中间量；I_X为X轴转动惯量、I_Z为Z轴转动惯量、I_Y为Y轴转动惯量、I_XZ为惯量积；

$\left\{\begin{array}{c}{C}_L=C_{L0}+C_{\mathrm{L\α}}\·\mathrm{\α}+C_{\mathrm{Lq}}\·q+C_{{\mathrm{L\δ}}_e}\·{\mathrm{\δ}}_e\\ C_D=C_{D0}+C_{\mathrm{D\α}}\·\mathrm{\α}+C_{\mathrm{Dq}}\·q+C_{D{\mathrm{\δ}}_e}\·{\mathrm{\δ}}_e\\ C_m=C_{m0}+C_{\mathrm{m\α}}\·\mathrm{\α}+C_{\mathrm{mq}}\·q+C_{m{\mathrm{\δ}}_e}\·{\mathrm{\δ}}_e\\ C_Y=C_{\mathrm{Y\β}}\·\mathrm{\β}+C_{\mathrm{Yp}}\·p+C_{\mathrm{Yr}}\·r+C_{{\mathrm{Y\δ}}_a}\·{\mathrm{\δ}}_a+C_{Y{\mathrm{\δ}}_r}\·{\mathrm{\δ}}_r\\ C_l=C_{\mathrm{l\β}}\·\mathrm{\β}+C_{\mathrm{lp}}\·p+C_{\mathrm{lr}}\·r+C_{{\mathrm{l\δ}}_a}\·{\mathrm{\δ}}_a+C_{l{\mathrm{\δ}}_r}\·{\mathrm{\δ}}_r\\ C_n=C_{\mathrm{n\β}}\·\mathrm{\β}+C_{\mathrm{np}}\·p+C_{\mathrm{nr}}\·r+C_{{\mathrm{n\δ}}_a}\·{\mathrm{\δ}}_a+C_{n{\mathrm{\δ}}_r}\·{\mathrm{\δ}}_r\\ L=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_L\\ D=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_D\\ M=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_Y\\ Y=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_Y\\ \stackrel{\‾}{L}=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_l\·b\\ N=\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_n\·C_A\end{array}\right.---\left(3\right)$

其中：C_L为总升力系数、C_D为总阻力系数、C_Y为总侧力系数、C_m为总俯仰力矩系数、C_l为总滚转力矩系数、C_n为总偏航力矩系数、ρ为空气密度、L为总升力、D为总阻力、Y为总侧力、

为总滚转力矩、M为总俯仰力矩、N为总偏航力矩、S_w为飞行器参考面积、C_A为飞行器平均几何弦长、b为翼展；

（3）解耦线性化；

将状态向量[V α β p q r φ θ ψ x y h m δ_a δ_e δ_r δ_t]^T分成纵向和横侧向两组，分别为纵向[V α q θ h δ_e δ_t]^T，横侧向[β p r φ ψ δ_a δ_r]^T，当无人机的航线为直线段，则有升力-重力平衡、推力-阻力平衡以及俯仰力矩平衡，因此下面三个等式成立：

$\left\{\begin{array}{c}\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{L0}+C_{\mathrm{L\α}}\·\mathrm{\α}+C_{{\mathrm{L\δ}}_e}\·{\mathrm{\δ}}_e)=m\·g\\ \frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{D0}+C_{\mathrm{D\α}}\·\mathrm{\α}+C_{{\mathrm{D\δ}}_e}\·{\mathrm{\δ}}_e)=t_A\\ \frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{m0}+C_{\mathrm{m\α}}\·\mathrm{\α}+C_{m{\mathrm{\δ}}_e}\·{\mathrm{\δ}}_e)\·C_A=0\end{array}\right.---\left(4\right)$

解方程（4），得到无人飞行器的配平攻角α₀，配平升降舵偏角δ_e0，配平飞行速度V₀，在该状态点利用小扰动线性化原理可得无人飞行器纵向线性状态方程为：

$\left[\begin{array}{c}\stackrel{\·}{V}\\ \stackrel{\·}{\mathrm{\α}}\\ \stackrel{\·}{q}\\ \stackrel{\·}{\mathrm{\θ}}\\ \stackrel{\·}{h}\\ {\stackrel{\·}{\mathrm{\δ}}}_e\\ {\stackrel{\·}{\mathrm{\δ}}}_t\end{array}\right]=A_1\left[\begin{array}{c}V\\ \mathrm{\α}\\ q\\ \mathrm{\θ}\\ h\\ {\mathrm{\δ}}_e\\ {\mathrm{\δ}}_t\end{array}\right]+B_1\left[\begin{array}{c}{\mathrm{\δ}}_e^{*}\\ {\mathrm{\δ}}_t^{*}\end{array}\right]---\left(5\right)$

其中：A₁为纵向线性系统状态矩阵，B₁为纵向线性系统控制矩阵；

无人飞行器横侧向线性状态方程为：

$\left[\begin{array}{c}\stackrel{\·}{\mathrm{\β}}\\ \stackrel{\·}{p}\\ \stackrel{\·}{r}\\ \stackrel{\·}{\mathrm{\φ}}\\ \stackrel{\·}{\mathrm{\ψ}}\\ {\stackrel{\·}{\mathrm{\δ}}}_a\\ {\stackrel{\·}{\mathrm{\δ}}}_r\end{array}\right]=A_2=\left[\begin{array}{c}\mathrm{\β}\\ p\\ r\\ \mathrm{\φ}\\ \mathrm{\ψ}\\ {\mathrm{\δ}}_a\\ {\mathrm{\δ}}_r\end{array}\right]+B_2\left[\begin{array}{c}{\mathrm{\δ}}_a^{*}\\ {\mathrm{\δ}}_r^{*}\end{array}\right]---\left(6\right)$

其中A₂为横侧向线性系统状态矩阵，B₂为横侧向线性系统控制矩阵；

（4）LQR设计控制律；

针对纵向：选取Q₁和R₁矩阵，采用线性二次型调节器设计纵向控制律K₁，K₁为2×7的矩阵；用k_ij表示K₁中的第i行第j列元素，令k₂₂＝k₂₃＝k₂₄＝k₂₅＝0.0，k₁₁＝0.0，k₁₆＝k₁₇＝k₂₆＝k₂₇＝0.0，最终选择以下控制器结构：

$\left\{\begin{array}{c}{\mathrm{\δ}}_e^{*}={\mathrm{\δ}}_{e0}+k_{12}\·(\mathrm{\α}-{\mathrm{\α}}_0)-k_{13}\·q-k_{14}\·(\mathrm{\θ}-{\mathrm{\α}}_0)-k_{15}\·(h-h^{*})\\ {\mathrm{\δ}}_t^{*}=t_A-k_{21}\·(V-V_0)\end{array}\right.---\left(7\right)$

其中：δ_e0为配平升降舵偏角，α₀为无人飞行器的配平攻角，h^*为高度指令，V₀为配平飞行速度；

微小型飞行器在飞行过程中有：

$\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·(C_{L0}+C_{\mathrm{L\α}}\·(\mathrm{\α}-{\mathrm{\α}}_0))+{\mathrm{\δ}}_t\·\sin \left(\mathrm{\α}\right)=(N_Z-1)\·m\·g---\left(8\right)$

其中：N_Z是无人飞行器的法向过载，忽略C_L0、δ_t·sin(α)，代入数据得：

$(\mathrm{\α}-{\mathrm{\α}}_0)=(N_Z-1)\·m\·g/(\frac{1}{2}\mathrm{\ρ}{V}^2\·S_w\·C_{\mathrm{L\α}})---\left(9\right)$

联立式（7）、式（9）得到实际飞行中的纵向控制律：

$\left\{\begin{array}{c}{\mathrm{\δ}}_e^{*}={\mathrm{\δ}}_{e0}-k_1\·(N_z-1)-k_2\·q-k_3\·(\mathrm{\θ}-{\mathrm{\α}}_0)-k_4\·(h-h^{*})\\ {\mathrm{\δ}}_t^{*}=t_A-k_5\·(V-V_0)\end{array}\right.---\left(10\right)$

其中：k₁ k₂ k₃ k₄ k₅表示按照（7）和（9）换算后的最终应用的控制器反馈系数；

针对横侧向：选取Q₂和R₂矩阵，采用线性二次型调节器设计纵向控制律K₂，K₂为2×7的矩阵；用

表示K₂中的第i行第j列元素，令

故K₂为2×5的矩阵；

实际飞行中的横侧向控制律：

$\left[\begin{array}{c}{\mathrm{\δ}}_a^{*}\\ {\mathrm{\δ}}_r^{*}\end{array}\right]=K_2\·\left[\begin{array}{c}\mathrm{\β}\\ p\\ r\\ \mathrm{\φ}\\ \mathrm{\ψ}\end{array}\right]---\left(11\right)$

步骤二：存在网络随机延迟的系统随机鲁棒性分析；

采用相互独立的随机延迟作为网络飞行控制系统的延迟类型，假设传感器观测延迟τ_o和执行器控制延迟τ_c服从均值μ_d，方差σ_d的正态分布，即：

τ_o∈N(μ_d,σ_d),τ_c∈N(μ_d,σ_d)        （12）

无人飞行器非线性系统方程由原来的

变为了：

$\stackrel{\·}{\stackrel{\‾}{X}}\left(t\right)=f(\stackrel{\‾}{X}(t-{\mathrm{\τ}}_o),u(t-{\mathrm{\τ}}_c))---\left(13\right)$

其中：为飞行状态向量，

为飞行状态向量的导数，

为加入传感器观测延迟的飞行状态向量，u(t-τ_c)为加入执行器控制延迟的飞行控制向量；

假设此时蒙特卡洛仿真确定的高度响应的平均超调量是σ，高度响应的平均调节时间是T_A；

步骤三：确定网络更新周期、优化指标与指标权值；

根据步骤二的蒙特卡洛仿真结果可以确定，编队网络更新率为：

$r=\frac{1}{{\mathrm{\μ}}_d}---\left(14\right)$

其中：r为网络更新率，μ_d为传感器观测延迟和执行器控制延迟时间的均值，设定随机鲁棒设计的代价函数：

W＝ω_σ·σ+ω_T·T_A        （15）

其中：σ是蒙特卡洛仿真高度响应的平均超调量，T_A是蒙特卡洛仿真高度响应的平均调节时间；

步骤四：利用随机鲁棒设计方法设计鲁棒控制律；

采用标准粒子群算法求解代价函数的优化问题，问题的解对应于搜索空间中的一个粒子，每个粒子都有自己的位置和速度及一个由被优化函数决定的代价函数；各个粒子记忆、追随当前的最优粒子，每次迭代中，粒子通过跟随两个极值来更新自己的位置和速度：一个是粒子本身获得的最优解

一个是整个粒子群中全部粒子历代搜索获得的最优解

标准粒子群算法算法中，速度更新和位置更新如式(16)所示：

$\left\{\begin{array}{c}{v}_{\mathrm{id}}^{k+1}=\mathrm{\ω}{v}_i^k+c_1{r}_1(p_{\mathrm{id}}^k-x_{\mathrm{id}}^k)+c_2{r}_2(g_{\mathrm{id}}^k-x_{\mathrm{id}}^k)\\ x_{\mathrm{id}}^{k+1}=x_{\mathrm{id}}^k+v_{\mathrm{id}}^{k+1}\end{array}\right.---\left(16\right)$

式中，i∈N(1,m),m为粒子群中粒子的个数，N代表整数；d∈N(1,n)，n为解向量的维数，k为迭代次数，c₁和c₂为学习因子，初始值为0.6，随迭代次数增加降为0.2；ω为惯量权重，初始值为0.9，随迭代次数增加降为0.4；r₁和r₂为[0,1]之间的随机数，判断粒子群中粒子优劣的标准是蒙特卡洛仿真获得的代价函数W，代价函数越小，在下一代该粒子存在的概率越大；

设定粒子最大的迭代次数，判断迭代是否到达最大迭代次数，如果到达，输出最优粒子代表的控制系统参数k₁,k₂,k₃,k₄,k₅，进而对微小型无人飞行器的飞行进行控制，否则，继续进行迭代。

2.根据权利要求1所述的一种存在网络随机延迟的微小型无人飞行器控制方法，其特征在于，还包括步骤五，闭环六自由度非线性蒙特卡洛仿真验证，具体为：

在同时存在传感器观测延迟τ_o和执行器控制延迟τ_c的条件下，且：

τ_o∈N(μ_d,σ_d),τ_c∈N(μ_d,σ_d)

对步骤四中随机鲁棒设计出的控制律进行蒙特卡洛仿真验证。

Images