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

摘要

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

主权项

一种存在网络随机延迟的微小型无人飞行器控制方法，其特征在于，包括以下几个步骤：步骤一：建立无人飞行器系统模型并确定控制律；具体为：(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、平均几何弦长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)建立非线性六自由度动力学与运动学方程；选取无人机飞行状态向量为：$\bar{X}\left(t\right)={\left[\begin{array}{ccccccccccccccccc}V& \alpha & \beta & p& q& r& \phi & \theta & \psi & x& y& h& m& {\delta }_a& {\delta }_e& {\delta }_r& {\delta }_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)={[{\delta }_a^{*},{\delta }_e^{*},{\delta }_r^{*},{\delta }_t^{*}]}^T$为飞行控制向量；$\left\{\begin{array}{c}{G}_x=m·g·(-\cos \left(\alpha \right)·\cos \left(\beta \right)·\sin \left(\theta \right)+\sin \left(\beta \right)·\sin \left(\phi \right)·\cos \left(\theta \right)+\sin \left(\alpha \right)·\cos \left(\beta \right)·\cos \left(\phi \right)·\cos \left(\theta \right))\\ G_y=m·g·(\cos \left(\alpha \right)·\sin \left(\beta \right)·\sin \left(\theta \right)+\cos \left(\beta \right)·\sin \left(\phi \right)·\cos \left(\theta \right)-\sin \left(\alpha \right)·\sin \left(\beta \right)·\cos \left(\phi \right)·\cos \left(\theta \right))\\ G_z=m·g·(\sin \left(\alpha \right)·\sin \left(\theta \right)+\cos \left(\alpha \right)·\cos \left(\phi \right)·\cos \left(\theta \right))\\ \Sigma =I_X·I_Z-I_{\mathrm{XZ}}·I_{\mathrm{XZ}}\\ c_1=((I_Y-I_Z)·I_Z-I_{\mathrm{XZ}}·I_{\mathrm{XZ}})/\Sigma \\ c_2=(I_X{-I}_Y+I_Z)·I_{\mathrm{XZ}}/\Sigma \\ c_3=I_Z/\Sigma \\ c_4=I_{\mathrm{XZ}}/\Sigma \\ 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}})/\Sigma \\ c_9=I_X/\Sigma \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\alpha }}·\alpha +C_{\mathrm{Lq}}·q+C_{{\mathrm{L\delta }}_e}·{\delta }_e\\ C_D=C_{D0}+C_{\mathrm{D\alpha }}·\alpha +C_{\mathrm{Dq}}·q+C_{{\mathrm{D\delta }}_e}·{\delta }_e\\ C_m=C_{m0}+C_{\mathrm{m\alpha }}·\alpha +C_{\mathrm{mq}}·q+C_{{\mathrm{m\delta }}_e}·{\delta }_e\\ C_Y=C_{\mathrm{Y\beta }}·\beta +C_{\mathrm{Yp}}·p+C_{\mathrm{Yr}}·r+C_{{\mathrm{Y\delta }}_a}·{\delta }_a+C_{{\mathrm{Y\delta }}_r}·{\delta }_r\\ C_l=C_{\mathrm{l\beta }}·\beta +C_{\mathrm{lp}}·p+C_{\mathrm{lr}}·r+C_{{\mathrm{l\delta }}_a}·{\delta }_a+C_{{\mathrm{l\delta }}_r}·{\delta }_r\\ C_n=C_{\mathrm{n\beta }}·\beta +C_{\mathrm{np}}·p+C_{\mathrm{nr}}·r+C_{{\mathrm{n\delta }}_a}·{\delta }_a+C_{{\mathrm{n\delta }}_r}·{\delta }_r\\ L=\frac{1}{2}{\mathrm{\rho V}}^2·S_w·C_L\\ D=\frac{1}{2}{\mathrm{\rho V}}^2·S_w·C_D\\ M=\frac{1}{2}{\mathrm{\rho V}}^2·S_w·C_m·C_A\\ Y=\frac{1}{2}{\mathrm{\rho V}}^2·S_w·C_Y\\ \bar{L}=\frac{1}{2}{\mathrm{\rho V}}^2·S_w·C_l·b\\ N=\frac{1}{2}{\mathrm{\rho 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{\rho V}}^2·S_w·(C_{L0}+C_{\mathrm{L\alpha }}·\alpha +C_{{\mathrm{L\delta }}_e}·{\delta }_e)=m·g\\ \frac{1}{2}{\mathrm{\rho V}}^2·S_w·(C_{D0}+C_{\mathrm{D\alpha }}·\alpha +C_{{\mathrm{D\delta }}_e}·{\delta }_e)=t_A\\ \frac{1}{2}{\mathrm{\rho V}}^2·S_w·(C_{m0}+C_{\mathrm{m\alpha }}·\alpha +C_{{\mathrm{m\delta }}_e}·{\delta }_e)·C_A=0\end{array}\right.---\left(4\right)$解方程(4)，得到无人飞行器的配平攻角α₀，配平升降舵偏角δ_e0，配平飞行速度V₀，在该状态点利用小扰动线性化原理可得无人飞行器纵向线性状态方程为：$\left[\begin{array}{c}\stackrel{·}{V}\\ \stackrel{·}{\alpha }\\ \stackrel{·}{q}\\ \stackrel{·}{\theta }\\ \stackrel{·}{h}\\ {\stackrel{·}{\delta }}_e\\ {\stackrel{·}{\delta }}_t\end{array}\right]=A_1\left[\begin{array}{c}V\\ \alpha \\ q\\ \theta \\ h\\ {\delta }_e\\ {\delta }_t\end{array}\right]+B_1\left[\begin{array}{c}{\delta }_e^{*}\\ {\delta }_t^{*}\end{array}\right]---\left(5\right)$其中：A₁为纵向线性系统状态矩阵，B₁为纵向线性系统控制矩阵；无人飞行器横侧向线性状态方程为：$\left[\begin{array}{c}\stackrel{·}{\beta }\\ \stackrel{·}{p}\\ \stackrel{·}{r}\\ \stackrel{·}{\phi }\\ \stackrel{·}{\psi }\\ {\stackrel{·}{\delta }}_a\\ {\stackrel{·}{\delta }}_r\end{array}\right]=A_2\left[\begin{array}{c}\beta \\ p\\ r\\ \phi \\ \psi \\ {\delta }_a\\ {\delta }_r\end{array}\right]+B_2\left[\begin{array}{c}{\delta }_a^{*}\\ {\delta }_r^{*}\end{array}\right]---\left(6\right)$其中A₂为横侧向线性系统状态矩阵，B₂为横侧向线性系统控制矩阵；(4)LQR设计控制律；针对纵向：选取Q₁和R₁矩阵，采用线性二次型调节器设计纵向控制律K₁，K₁为2×7的矩阵；用ki_j表示K₁中的第i行第j列元素，令k₂₂＝k₂₃＝k₂₄＝k₂₅＝0.0，k₁₁＝0.0，k₁₆＝k₁₇＝k₂₆＝k₂₇＝0.0，最终选择以下控制器结构：$\left\{\begin{array}{c}{\delta }_e^{*}={\delta }_{e0}+k_{12}·\left({\alpha -\alpha }_0\right)-k_{13}·q-k_{14}·(\theta -{\alpha }_0)-k_{15}·(h-h^{*})\\ {\delta }_t^{*}=t_A-k_{21}·(V-V_0)\end{array}\right.---\left(7\right)$其中：δ_e0为配平升降舵偏角，α₀为无人飞行器的配平攻角，h^*为高度指令，V₀为配平飞行速度；微小型飞行器在飞行过程中有：$\frac{1}{2}{\mathrm{\rho V}}^2·S_w·(C_{L0}+C_{\mathrm{L\alpha }}·(\alpha -{\alpha }_0))+{\delta }_t·\sin \left(\alpha \right)=(N_Z-1)·m·g---\left(8\right)$其中：N_Z是无人飞行器的法向过载，忽略C_L0、δ_t·sin(α)，代入数据得：$(\alpha -{\alpha }_0)=(N_Z-1)·m·g/(\frac{1}{2}{\mathrm{\rho V}}^2·S_w·C_{\mathrm{L\alpha }})---\left(9\right)$联立式(7)、式(9)得到实际飞行中的纵向控制律：$\left\{\begin{array}{c}{\delta }_e^{*}={\delta }_{e0}-k_1·(N_Z-1)-k_2·q-k_3·(\theta -{\alpha }_0)-k_4·(h-h^{*})\\ {\delta }_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}{\delta }_a^{*}\\ {\delta }_r^{*}\end{array}\right]=K_2·\left[\begin{array}{c}\beta \\ p\\ r\\ \phi \\ \psi \end{array}\right]---\left(11\right)$步骤二：存在网络随机延迟的系统随机鲁棒性分析；采用相互独立的随机延迟作为网络飞行控制系统的延迟类型，假设传感器观测延迟τ_o和执行器控制延迟τ_c服从均值μ_d，方差σ_d的正态分布，即：τ_o∈N(μ_d,σ_d),τ_c∈N(μ_d,σ_d)    (12)无人飞行器非线性系统方程由原来的u(t))变为了：$\stackrel{·}{\stackrel{‾}{X}}\left(t\right)=f(\bar{X}(t-{\tau }_o),u(t-{\tau }_c))---\left(13\right)$其中：为飞行状态向量，为飞行状态向量的导数，为加入传感器观测延迟的飞行状态向量，u(t‑τ_c)为加入执行器控制延迟的飞行控制向量；假设此时蒙特卡洛仿真确定的高度响应的平均超调量是σ，高度响应的平均调节时间是T_A；步骤三：确定网络更新周期、优化指标与指标权值；根据步骤二的蒙特卡洛仿真结果可以确定，编队网络更新率为：$r=\frac{1}{{\mu }_d}---\left(14\right)$其中：r为网络更新率，μ_d为传感器观测延迟和执行器控制延迟时间的均值，设定随机鲁棒设计的代价函数：W＝ω_σ·σ+ω_T·T_A    (15)其中：σ是蒙特卡洛仿真高度响应的平均超调量，T_A是蒙特卡洛仿真高度响应的平均调节时间；步骤四：利用随机鲁棒设计方法设计鲁棒控制律；采用标准粒子群算法求解代价函数的优化问题，问题的解对应于搜索空间中的一个粒子，每个粒子都有自己的位置和速度及一个由被优化函数决定的代价函数；各个粒子记忆、追随当前的最优粒子，每次迭代中，粒子通过跟随两个极值来更新自己的位置和速度：一个是粒子本身获得的最优解一个是整个粒子群中全部粒子历代搜索获得的最优解标准粒子群算法算法中，速度更新和位置更新如式(16)所示：$\left\{\begin{array}{c}{v}_{\mathrm{id}}^{k+1}={\mathrm{\omega 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₅，进而对微小型无人飞行器的飞行进行控制，否则，继续进行迭代。