CN103616816B - 一种高超声速飞行器升降舵故障控制方法 - Google Patents

Abstract

本发明提出了一种高超声速飞行器升降舵故障控制方法。所述方法首先基于高超声速飞行器巡航段的纵向模型建立其升降舵发生卡死故障的模型；然后设计状态反馈控制律并建立故障模式、大小、发生时刻已知情况下实现容错控制所需的匹配条件；进而考虑升降舵发生未知故障的情况，推导出控制器参数的自适应学习律。该发明能够通过快速有效地自适应调整正常冗余升降舵的控制输出，实现对发生故障的升降舵所损失的控制效力的补偿，保证高超声速飞行器在某一升降舵发生故障时的飞行安全，并且最大程度恢复飞行器的期望性能。

Description

一种高超声速飞行器升降舵故障控制方法

技术领域

本发明属于自动控制技术领域，尤其是涉及一种高超声速飞行器升降舵故障控制方法。

背景技术

高超声速飞行器是迄今为止飞行过程最为复杂的一种飞行器，与亚声速/超声速飞行器相比有许多不同的飞行特性。由于工作条件大范围变化，高超声速飞行器高低空的气动力特性的巨大差异和质量分布的快速变化，导致飞行器的动力学特征和模型参数在飞行中变化非常显著。发动机与机身一体化构形使得弹性机身、推进系统以及结构动力学之间的强耦合作用，导致高超声速飞行器的飞行动力学特征非常复杂，存在很大的不确定性。另外，由于重量、大小和成本的限制，高超声速飞行器通常装备最少组合的控制舵面，一个控制舵面发生故障往往对飞行器产生所需控制力矩的能力具有破坏性的影响。这就要求在出现控制面失效和各种异常故障状态时，高超声速飞行器的控制系统具有高度的自主容错能力。

目前针对高超声速飞行器的容错控制已有一些研究成果。美国Dryden飞行研究中心采用离线最优控制方法通过重新分配控制效力实现了对控制面故障的容错（ReconfigurablecontroldesignforthefullX-33flightenvelope.AIAA2001-4379,2001.）。Zhu等人提出了一种基于奇异摄动理论的直接容错控制方法（DirectfaulttolerantRLVattitudecontrolasingularperturbationapproach.AIAA2002-4778,2002）。Johnson等人基于直接自适应控制和神经网络动态逆设计终端能量管理段和进场着陆段的容错控制方案（Faulttolerancethroughdirectadaptivecontrolusingneuralnetworks.AIAA2006-6553,2006）。姜斌等人基于模糊控制、自适应控制、鲁棒控制和滑模控制等理论与方法设计了多种故障诊断与容错控制方案，增强了姿控系统的可靠性（近空间飞行器故障诊断与容错控制的研究进展.南京航空航天大学学报，2012）。目前针对高超声速飞行器的容错控制方法大多数需要检测和隔离发生故障的执行机构，需要获取发生故障的升降舵的故障信息，包括故障发生时间、具体哪个升降舵发生故障和故障的大小。

发明内容

本发明所要解决的技术问题在于克服现有技术的不足,提出了一种高超声速飞行器升降舵故障控制方法。所述方法不需要额外的故障诊断环节，即不需要检测和辨识升降舵的故障模式、故障大小和故障发生时间，能够通过在线自适应调节正常升降舵的偏转角度补偿故障升降舵损失的控制效力，保证飞行器的稳定性，并恢复其高度和速度跟踪性能。

为解决上述技术问题，本发明所采用的技术方案是：

一种高超声速飞行器升降舵故障控制方法，其特征在于，包括以下步骤：

步骤1，建立高超声速飞行器巡航段纵向运动学及动力学模型，所述模型表示为：

$\stackrel{\·}{V}=\frac{T\cos \mathrm{\α}-D}{m}-\frac{\mathrm{\μ}\sin \mathrm{\γ}}{r^2}$

$\stackrel{\·}{\mathrm{\γ}}=\frac{L+T\sin \mathrm{\α}}{\mathrm{mV}}-\frac{\left(\mathrm{\μ}{V}^{2}r\right)\cos \mathrm{\γ}}{{\mathrm{Vr}}^2}$

$\stackrel{\·}{h}=V\sin \mathrm{\γ}$

$\stackrel{\·}{\mathrm{\α}}=q-\stackrel{\·}{\mathrm{\γ}}$

$\stackrel{\·}{q}=M_{\mathrm{yy}}/I_{\mathrm{yy}}$

其中V、γ、h、α、q分别表示高超声速飞行器的速度、航迹角、高度、攻角和俯仰角速率；L、D、T分别为升力、阻力和推力；m、μ、r分别表示飞行器的质量、重力常数、飞行器到地心的距离，其中r=R_E+h，R_E为地球半径；I_yy、M_yy分别表示俯仰转动惯量和俯仰转动力矩；

步骤2，建立控制系统输入-输出模型

将所述高超声速飞行器巡航段纵向运动学及动力学模型转换为以高超声速飞行器的速度V和高度h为输出、油门开度β_c和升降舵舵偏角δ_e作为控制输入，建立输入-输出模型，其表达式如下：

$\stackrel{...}{V}=f_v+b_{11}{\mathrm{\β}}_c+b_{12}{\mathrm{\δ}}_e$

h⁽⁴⁾=f_h+b₂₁β_c+b₂₂δ_e

其中f_V、f_h、b₁₁、b₁₂、b₂₁、b₂₂分别是关于V、γ、h、α、q和飞行器气动参数的非线性函数；总的舵偏角δ_e与左、右两个升降舵的舵偏角δ_e1和δ_e2的关系为：δ_e=d₁δ_e1+d₂δ_e2，其中d₁和d₂为组合关系系数；

步骤3，考虑步骤2所述输入-输出模型升降舵发生卡死故障的情况，建立升降舵的故障模型，其表达式如下：

${\mathrm{\δ}}_{\mathrm{ej}}=v_j+{\mathrm{\σ}}_j({\stackrel{\‾}{\mathrm{\δ}}}_j-v_j),j=\mathrm{1,2}$

其中，表示升降舵卡死故障的大小；v_j表示第j个升降舵的控制输入；σ_j代表故障模式，当第j个升降舵发生故障时，σ_j=1，否则σ_j=0；

步骤4，基于速度追踪误差和高度追踪误差定义两个滑模面s₁和s₂：

s₁=λ₁₀e₁₀+λ₁₁e₁₁+λ₁₂e₁₂+e₁₃

s₂=λ₂₀e₂₀+λ₂₁e₂₁+λ₂₂e₂₂+λ₂₃e₂₃+e₂₄

其中，e₁₁=V-V_d，e₁₀=∫e₁₁，e₂₁=h-h_d，e₂₀=∫e₂₁， V_d和h_d为期望的速度和高度信号；λ₁₀,λ₁₁,λ₁₂,λ₂₀,λ₂₁,λ₂₂,λ₂₃为正常数，其选择满足多项式s³+λ₁₂s²+λ₁₁s+λ₁₀和s⁴+λ₂₃s³+λ₂₂s²+λ₂₁s+λ₂₀为Hurwitz多项式。

步骤5，确定标称控制输入，其表达式如下：

$\left[\begin{array}{c}{\mathrm{\β}}_c\\ v_0\end{array}\right]=B^{-1}\left[\begin{array}{cc}-f_1& -a_1\mathrm{sgn}\left(s_1\right)\\ -f_2& -a_2\mathrm{sgn}\left(s_2\right)\end{array}\right],B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]$

其中v₀为升降舵δ_e的标称偏转量；a₁和a₂为两个正常数，决定滑动模态的收敛速度；f₁和f₂分别由下式确定：

$f_1=-{\stackrel{...}{V}}_d+f_v+{\mathrm{\λ}}_{10}{e}_{11}+{\mathrm{\λ}}_{11}{e}_{12}+{\mathrm{\λ}}_{12}{e}_{13}$

$f_2=h_d^{\left(4\right)}+f_h+{\mathrm{\λ}}_{20}{e}_{21}+{\mathrm{\λ}}_{21}{e}_{22}+{\mathrm{\λ}}_{22}{e}_{23}+{\mathrm{\λ}}_{23}{e}_{24}$

sgn(·)为符号函数，其定义如下：

$\mathrm{sgn}\left(x\right)=\begin{array}{}\end{array}\left\{\begin{array}{cc}1& x>0\\ 0& x=0\\ -1& x<0\end{array}\right.$

或者为减小滑模控制产生的抖振，β_c、v₀分别由下式确定：

$\left[\begin{array}{c}{\mathrm{\&beta;}}_c\\ v_0\end{array}\right]=B^{-1}\left[\begin{array}{cc}-f_1& -a_1\mathrm{sat}\left(\frac{s_1}{\mathrm{\&epsiv;}}\right)\\ -f_2& -a_2\mathrm{sat}\left(\frac{s_2}{\mathrm{\&epsiv;}}\right)\end{array}\right],B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]$

其中，ε为一个小的正常数；sat(·)为饱和函数，其定义如下：

$\mathrm{sat}\left(x\right)=\left\{\begin{array}{cc}x& \mathrm{if}\left|x\right|\&le;1\\ \mathrm{sgn}\left(x\right)& \mathrm{if}\left|x\right|>1\end{array}\right.$

步骤6，基于步骤5给出的标称控制信号v₀，确定自适应容错控制律：

$v_j={\hat{k}}_{1j}{v}_0+{\hat{k}}_{2j},j=\mathrm{1,2}$

其中分别为容错控制器的参数估计值，分别由下式确定：

${\stackrel{\stackrel{.}{^}}{k}}_{1j}=-\mathrm{sgn}\left(d_j\right){\mathrm{\&Gamma;}}_{1j}(b_{12}{s}_1+b_{22}{s}_2)v_0$

${\stackrel{\stackrel{.}{^}}{k}}_{2j}=-\mathrm{sgn}\left(d_j\right){\mathrm{\&Gamma;}}_{2j}(b_{12}{s}_1+b_{22}{s}_2)$

其中Γ_1j和Γ_2j为两个正常数，决定了参数自适应律的学习速率。

步骤2中，所述非线性函数f_V、f_h、b₁₁、b₁₂、b₂₁、b₂₂的具体形式如下：

$f_V=\frac{{\mathrm{\&omega;}}_1{\stackrel{..}{x}}_0+{\stackrel{.}{x}}^T{\mathrm{\&Omega;}}_2\stackrel{.}{x}}{m}$

$f_h\begin{array}{c}=3\stackrel{..}{V}\stackrel{.}{\mathrm{\&gamma;}}\cos \mathrm{\&gamma;}-3\stackrel{.}{V}{\stackrel{.}{\mathrm{\&gamma;}}}^2\sin \mathrm{\&gamma;}+3\stackrel{.}{V}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}}\cos \mathrm{\&gamma;}-3V\stackrel{.}{\mathrm{\&gamma;}}\stackrel{..}{\mathrm{\&gamma;}}\sin \mathrm{\&gamma;}-V{\stackrel{.}{\mathrm{\&gamma;}}}^3\cos \mathrm{\&gamma;}+\frac{({\mathrm{\&omega;}}_1{\stackrel{..}{x}}_0+{\stackrel{.}{x}}^T{\mathrm{\&Omega;}}_2\stackrel{.}{x})\sin \mathrm{\&gamma;}}{m}\\ +V\cos \mathrm{\&gamma;}({\mathrm{\&pi;}}_1{\stackrel{..}{x}}_0+{\stackrel{.}{x}}^T{\mathrm{\&Pi;}}_2\stackrel{.}{x})\end{array}$

$b_{11}=\left(\frac{\mathrm{\&rho;}{V}^{2}S{c}_{\mathrm{\&beta;}}{\mathrm{\&omega;}}_n^2}{2m}\right)\cos \mathrm{\&alpha;},b_{12}=-\left(\frac{c_e\mathrm{\&rho;}{V}^{2}S\stackrel{\&OverBar;}{c}}{2m{I}_{\mathrm{yy}}}\right)(T\sin \mathrm{\&alpha;}+\frac{\&PartialD;D}{\&PartialD;\mathrm{\&alpha;}})$

$b_{21}=\left(\frac{\mathrm{\&rho;}{V}^{2}S{c}_{\mathrm{\&beta;}}{\mathrm{\&omega;}}_n^2}{2m}\right)\sin (\mathrm{\&alpha;}+\mathrm{\&gamma;}),b_{22}=\left(\frac{c_e\mathrm{\&rho;}{V}^{2}S\stackrel{\&OverBar;}{c}}{2m{I}_{\mathrm{yy}}}\right)[T\cos (\mathrm{\&alpha;}+\mathrm{\&gamma;})+\left(\frac{\&PartialD;L}{\&PartialD;\mathrm{\&alpha;}}\right)\cos \mathrm{\&gamma;}-\left(\frac{\&PartialD;D}{\&PartialD;\mathrm{\&alpha;}}\right)\sin \mathrm{\&gamma;}]$

其中

$x={\left[\mathrm{V\&gamma;\&alpha;\&beta;h}\right]}^T,{\stackrel{..}{x}}_0={\left[\stackrel{..}{V}\stackrel{..}{\mathrm{\&gamma;}}{\stackrel{..}{\mathrm{\&alpha;}}}_0{\stackrel{..}{\mathrm{\&beta;}}}_0\stackrel{..}{h}\right]}^T$

${\stackrel{..}{\mathrm{\&alpha;}}}_0=\frac{1}{2}\mathrm{\&rho;}{V}^{2}S\stackrel{\&OverBar;}{c}[C_M\left(\mathrm{\&alpha;}\right)+C_M\left(q\right)-c_e\mathrm{\&alpha;}]/I_{\mathrm{yy}}-\stackrel{..}{\mathrm{\&gamma;}},{\stackrel{..}{\mathrm{\&beta;}}}_0=-2\mathrm{\&xi;}{\mathrm{\&omega;}}_n\stackrel{.}{\mathrm{\&beta;}}-{\mathrm{\&omega;}}_n^2\mathrm{\&beta;}$

C_M(α)=-0.035α²+0.036617α+5.3261×10^-6

$C_M\left(q\right)=(\stackrel{\&OverBar;}{c}/2V)q(-6.796{\mathrm{\&alpha;}}^2+0.3015\mathrm{\&alpha;}-0.2289)$

C_M(δ_e)=c_e(δ_e-α)

其中C_M(α)、C_M(δ_e)和为C_M(q)为俯仰力矩系数；I_yy为俯仰转动惯量；S为参考气动面积；ρ为大气密度；ω_n、ζ、β为发动机动态模型的自然频率、阻尼比和输出；α₀和β₀为转换模型引入的辅助变量；为平均气动翼弦；c_e和c_β为计算俯仰力矩系数的相关参数；ω₁、π₁、Ω₂、Π₂为推导V和h的高阶微分所产生的向量及矩阵：

${\mathrm{\&omega;}}_1=[\left(\frac{\&PartialD;T}{\&PartialD;V}\right)\cos \mathrm{\&alpha;}-\frac{\&PartialD;D}{\&PartialD;V}\frac{-\mathrm{m\&mu;}\cos \mathrm{\&gamma;}}{r^2}-T\sin \mathrm{\&alpha;}\frac{\&PartialD;D}{\&PartialD;\mathrm{\&alpha;}}\left(\frac{\&PartialD;T}{\&PartialD;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\frac{2\mathrm{m\&mu;}\sin \mathrm{\&gamma;}}{r^3}]$

${\mathrm{\&pi;}}_1^T=\left[\begin{array}{c}\frac{\&PartialD;L/\&PartialD;V+(\&PartialD;T/\&PartialD;V)\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{L+T\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}+\frac{\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^2}+\frac{\cos \mathrm{\&gamma;}}{r}\\ \frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}-\frac{V\sin \mathrm{\&gamma;}}{r}\\ \frac{\&PartialD;L/\&PartialD;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{\mathrm{mV}}\\ \frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{\mathrm{mV}}\\ \frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V{r}^3}-\frac{V\cos \mathrm{\&gamma;}}{r^2}\end{array}\right]$

Ω₂=[ω₂₁ω₂₂ω₂₃ω₂₄ω₂₅],Π₂=[π₂₁π₂₂π₂₃π₂₄π₅]

${\mathrm{\&omega;}}_{21}=\left[\begin{array}{c}\left(\frac{{\&PartialD;}^{2}T}{{\&PartialD;V}^2}\right)\cos \mathrm{\&alpha;}-\frac{{\&PartialD;}^{2}D}{{\&PartialD;V}^2}\\ 0\\ -\left(\frac{\&PartialD;T}{\&PartialD;V}\right)\sin \mathrm{\&alpha;}-\frac{{\&PartialD;}^{2}D}{\&PartialD;V\&PartialD;\mathrm{\&alpha;}}\\ \left(\frac{{\&PartialD;}^{2}T}{\&PartialD;V\&PartialD;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\\ 0\end{array}\right],{\mathrm{\&omega;}}_{22}=\left[\begin{array}{c}0\\ \frac{\mathrm{m\&mu;}\sin \mathrm{\&gamma;}}{r^2}\\ 0\\ 0\\ \frac{2\mathrm{m\&mu;}\cos \mathrm{\&gamma;}}{r^3}\end{array}\right],{\mathrm{\&omega;}}_{23}=\left[\begin{array}{c}-\left(\frac{\&PartialD;T}{\&PartialD;V}\right)\sin \mathrm{\&alpha;}-\left(\frac{{\&PartialD;}^{2}D}{\&PartialD;V\&PartialD;\mathrm{\&alpha;}}\right)\\ 0\\ -T\cos \mathrm{\&alpha;}\left(\frac{{\&PartialD;}^{2}D}{{\&PartialD;\mathrm{\&alpha;}}^2}\right)\\ -\left(\frac{\&PartialD;T}{\&PartialD;\mathrm{\&beta;}}\right)\sin \mathrm{\&alpha;}\\ 0\end{array}\right]$

${\mathrm{\&omega;}}_{24}=\left[\begin{array}{c}\left(\frac{{\&PartialD;}^{2}T}{\&PartialD;V\&PartialD;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\\ 0\\ -\left(\frac{\&PartialD;T}{\&PartialD;\mathrm{\&beta;}}\right)\sin \mathrm{\&alpha;}\\ 0\\ 0\end{array}\right],{\mathrm{\&omega;}}_{25}=\left[\begin{array}{c}0\\ \frac{2\mathrm{m\&mu;}\cos \mathrm{\&gamma;}}{r^3}\\ 0\\ 0\\ \frac{-6\mathrm{m\&mu;}\sin \mathrm{\&gamma;}}{r^4}\end{array}\right]$

${\mathrm{\&pi;}}_{21}=\left[\begin{array}{c}\frac{{\&PartialD;}^{2}L/{\&PartialD;V}^2+({\&PartialD;}^{2}T/{\&PartialD;V}^2)\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{2[\&PartialD;L/\&PartialD;V+(\&PartialD;T/\&PartialD;V)\sin \mathrm{\&alpha;}]}{{\mathrm{mV}}^2}+\frac{2(L+T\sin \mathrm{\&alpha;})}{{\mathrm{mV}}^3}-\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^3{r}^2}\\ -\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{V^2{r}^2}-\frac{\sin \mathrm{\&gamma;}}{r}\\ \frac{({\&PartialD;}^{2}L/\&PartialD;\mathrm{\&alpha;}\&PartialD;V)+(\&PartialD;T/\&PartialD;V)\cos \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{\&PartialD;L/\&PartialD;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ \frac{({\&PartialD;}^{2}T/\&PartialD;\mathrm{\&beta;}\&PartialD;V)\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ -\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^3}-\frac{\cos \mathrm{\&gamma;}}{r^2}\end{array}\right]$

${\mathrm{\&pi;}}_{22}=\left[\begin{array}{c}-\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{V^2{r}^2}-\frac{\sin \mathrm{\&gamma;}}{r}\\ \frac{\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}-\frac{V\cos \mathrm{\&gamma;}}{r}\\ 0\\ 0\\ -\frac{2\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}+\frac{V\sin \mathrm{\&gamma;}}{r^2}\end{array}\right],{\mathrm{\&pi;}}_{23}=\left[\begin{array}{c}\frac{({\&PartialD;}^{2}L/\&PartialD;V\&PartialD;\mathrm{\&alpha;})+(\&PartialD;T/\&PartialD;V)\cos \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{\&PartialD;L/\&PartialD;\mathrm{\&alpha;}+T\cos }{{\mathrm{mV}}^2}\\ 0\\ \frac{{\&PartialD;}^{2}L/\&PartialD;{\mathrm{\&alpha;}}^2-T\sin }{\mathrm{mV}}\\ \frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\cos \mathrm{\&alpha;}}{\mathrm{mV}}\\ 0\end{array}\right]$

${\mathrm{\&pi;}}_{24}=\left[\begin{array}{c}\frac{({\&PartialD;}^{2}T/\&PartialD;V\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ 0\\ \frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\cos \mathrm{\&alpha;}}{\mathrm{mV}}\\ 0\\ 0\end{array}\right],{\mathrm{\&pi;}}_{25}=\left[\begin{array}{c}-\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^3}-\frac{\cos \mathrm{\&gamma;}}{r^2}\\ -\frac{2\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^3}+\frac{V\sin \mathrm{\&gamma;}}{r^2}\\ 0\\ 0\\ -\frac{6\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^4}+\frac{2V\cos \mathrm{\&gamma;}}{r^3}\end{array}\right]$

ω₁、π₁、Ω₂、Π₂的表达式中关于L、T、D的一阶、二阶偏导数可根据下列关系计算得到：

L=ρV²SC_L2,C_L=0.6203α

D=ρV²SC_D2,C_D=0.6450α²+0.0043378α+0.003772

T=ρV²SC_T2, $C_T=\left\{\begin{array}{cc}0.02576\mathrm{\&beta;}& \mathrm{if\&beta;}\&le;1\\ 0.0224+0.00336\mathrm{\&beta;}& \mathrm{if\&beta;}>1\end{array}\right..$

本发明的有益效果是：一种高超声速飞行器升降舵故障控制方法，所述方法能够在两个升降舵中任何一个发生故障的情况下，通过自适应调整正常升降舵的偏转角度，补偿故障升降舵损失的控制效力。在此过程中，不需要检测是哪一个升降舵发生了故障，也不需要知道故障发生的时间及大小，参数自适应律会自动更新容错控制器的参数，保证飞行器稳定并恢复其高度和速度追踪性能。该发明不需额外的故障诊断环节，不依赖故障诊断的及时性和准确性，能够充分利用飞行器的控制冗余，不仅能够保证飞行器的安全飞行，而且能够维持其期望的渐近追踪性能。

附图说明

图1是本发明控制方法的示意图。

图2是输出和输入响应图；其中，图2（a）是速度追踪响应图；图2（b）是高度追踪响应图；图2（c）是推力响应图；图2（d）是合成舵偏角响应图。

图3是状态响应图；其中，图3（a）是航迹角响应图；图3（d）是俯仰角响应图；图3（c）是攻角响应图；图3（d）是俯仰角速率响应图。

图4是两个升降舵的偏转角图；其中，图4（a）是左升降舵的偏转角；图4（b）是右升降舵的偏转角。

图5是控制器参数的自适应调整图；其中，图5（a）是参数的自适应调整图；图5（b）是参数的自适应调整图；图5（c）是参数的自适应调整图；图5（d）是参数的自适应调整图。

具体实施方式

下面结合附图，对本发明提出的一种高超声速飞行器升降舵故障控制方法进行详细说明：

本发明提供了一种针对高超声速飞行器升降舵故障控制方法，包括以下步骤：

步骤1，建立高超声速飞行器巡航段纵向运动学及动力学模型，所述模型表示为：

$\stackrel{\&CenterDot;}{V}=\frac{T\cos \mathrm{\&alpha;}-D}{m}-\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{r^2}$

$\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{L+T\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{\left(\mathrm{\&mu;}{V}^{2}r\right)\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}$

$\stackrel{\&CenterDot;}{h}=V\sin \mathrm{\&gamma;}$

$\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}=q-\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}$

$\stackrel{\&CenterDot;}{q}=M_{\mathrm{yy}}/I_{\mathrm{yy}}$

其中V、γ、h、α、q分别表示高超声速飞行器的速度、航迹角、高度、攻角和俯仰角速率；L、D、T分别为升力、阻力和推力；m、μ、r分别表示飞行器的质量、重力常数、飞行器到地心的距离（r=R_E+h，R_E为地球半径）；I_yy、M_yy分别表示俯仰转动惯量和俯仰转动力矩。

步骤2，将纵向模型转换为以V和h为输出、油门开度β_c和升降舵舵偏角δ_e为控制输入的输入-输出模型：

$\stackrel{...}{V}=f_v+b_{11}{\mathrm{\&beta;}}_c+b_{12}{\mathrm{\&delta;}}_e$

h⁽⁴⁾=f_h+b₂₁β_c+b₂₂δ_e

其中f_V、f_h、b₁₁、b₁₂、b₂₁、b₂₂是有关V、γ、h、α、q和飞行器气动参数等物理量的非线性函数：

$f_V=\frac{{\mathrm{\&omega;}}_1{\stackrel{..}{x}}_0+{\stackrel{.}{x}}^T{\mathrm{\&Omega;}}_2\stackrel{.}{x}}{m}$

$f_h\begin{array}{c}=3\stackrel{..}{V}\stackrel{.}{\mathrm{\&gamma;}}\cos \mathrm{\&gamma;}-3\stackrel{.}{V}{\stackrel{.}{\mathrm{\&gamma;}}}^2\sin \mathrm{\&gamma;}+3\stackrel{.}{V}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}}\cos \mathrm{\&gamma;}-3V\stackrel{.}{\mathrm{\&gamma;}}\stackrel{..}{\mathrm{\&gamma;}}\sin \mathrm{\&gamma;}-V{\stackrel{.}{\mathrm{\&gamma;}}}^3\cos \mathrm{\&gamma;}+\frac{({\mathrm{\&omega;}}_1{\stackrel{..}{x}}_0+{\stackrel{.}{x}}^T{\mathrm{\&Omega;}}_2\stackrel{.}{x})\sin \mathrm{\&gamma;}}{m}\\ +V\cos \mathrm{\&gamma;}({\mathrm{\&pi;}}_1{\stackrel{..}{x}}_0+{\stackrel{.}{x}}^T{\mathrm{\&Pi;}}_2\stackrel{.}{x})\end{array}$

$b_{11}=\left(\frac{\mathrm{\&rho;}{V}^{2}S{c}_{\mathrm{\&beta;}}{\mathrm{\&omega;}}_n^2}{2m}\right)\cos \mathrm{\&alpha;},b_{12}=-\left(\frac{c_e\mathrm{\&rho;}{V}^{2}S\stackrel{\&OverBar;}{c}}{2m{I}_{\mathrm{yy}}}\right)(T\sin \mathrm{\&alpha;}+\frac{\&PartialD;D}{\&PartialD;\mathrm{\&alpha;}})$

$b_{21}=\left(\frac{\mathrm{\&rho;}{V}^{2}S{c}_{\mathrm{\&beta;}}{\mathrm{\&omega;}}_n^2}{2m}\right)\sin (\mathrm{\&alpha;}+\mathrm{\&gamma;}),b_{22}=\left(\frac{c_e\mathrm{\&rho;}{V}^{2}S\stackrel{\&OverBar;}{c}}{2m{I}_{\mathrm{yy}}}\right)[T\cos (\mathrm{\&alpha;}+\mathrm{\&gamma;})+\left(\frac{\&PartialD;L}{\&PartialD;\mathrm{\&alpha;}}\right)\cos \mathrm{\&gamma;}-\left(\frac{\&PartialD;D}{\&PartialD;\mathrm{\&alpha;}}\right)\sin \mathrm{\&gamma;}]$

其中

$x={\left[\mathrm{V\&gamma;\&alpha;\&beta;h}\right]}^T,{\stackrel{..}{x}}_0={\left[\stackrel{..}{V}\stackrel{..}{\mathrm{\&gamma;}}{\stackrel{..}{\mathrm{\&alpha;}}}_0{\stackrel{..}{\mathrm{\&beta;}}}_0\stackrel{..}{h}\right]}^T$

${\stackrel{..}{\mathrm{\&alpha;}}}_0=\frac{1}{2}\mathrm{\&rho;}{V}^{2}S\stackrel{\&OverBar;}{c}[C_M\left(\mathrm{\&alpha;}\right)+C_M\left(q\right)-c_e\mathrm{\&alpha;}]/I_{\mathrm{yy}}-\stackrel{..}{\mathrm{\&gamma;}},{\stackrel{..}{\mathrm{\&beta;}}}_0=-2\mathrm{\&xi;}{\mathrm{\&omega;}}_n\stackrel{.}{\mathrm{\&beta;}}-{\mathrm{\&omega;}}_n^2\mathrm{\&beta;}$

C_M(α)=-0.035α²+0.036617α+5.3261×10^-6

$C_M\left(q\right)=(\stackrel{\&OverBar;}{c}/2V)q(-6.796{\mathrm{\&alpha;}}^2+0.3015\mathrm{\&alpha;}-0.2289)$

C_M(δ_e)=c_e(δ_e-α)

其中C_M(α)、C_M(δ_e)和为C_M(q)为俯仰力矩系数；I_yy为俯仰转动惯量；S为参考气动面积；ρ为大气密度；ω_n、ζ、β为发动机动态模型的自然频率、阻尼比和输出；α₀和β₀为转换模型引入的辅助变量；为平均气动翼弦；c_e和c_β为计算俯仰力矩系数的相关参数；ω₁、π₁、Ω₂、Π₂为推导V和h的高阶微分所产生的向量及矩阵：

${\mathrm{\&omega;}}_1=[\left(\frac{\&PartialD;T}{\&PartialD;V}\right)\cos \mathrm{\&alpha;}-\frac{\&PartialD;D}{\&PartialD;V}\frac{-\mathrm{m\&mu;}\cos \mathrm{\&gamma;}}{r^2}-T\sin \mathrm{\&alpha;}\frac{\&PartialD;D}{\&PartialD;\mathrm{\&alpha;}}\left(\frac{\&PartialD;T}{\&PartialD;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\frac{2\mathrm{m\&mu;}\sin \mathrm{\&gamma;}}{r^3}]$

${\mathrm{\&pi;}}_1^T=\left[\begin{array}{c}\frac{\&PartialD;L/\&PartialD;V+(\&PartialD;T/\&PartialD;V)\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{L+T\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}+\frac{\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^2}+\frac{\cos \mathrm{\&gamma;}}{r}\\ \frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}-\frac{V\sin \mathrm{\&gamma;}}{r}\\ \frac{\&PartialD;L/\&PartialD;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{\mathrm{mV}}\\ \frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{\mathrm{mV}}\\ \frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V{r}^3}-\frac{V\cos \mathrm{\&gamma;}}{r^2}\end{array}\right]$

Ω₂=[ω₂₁ω₂₂ω₂₃ω₂₄ω₂₅],Π₂=[π₂₁π₂₂π₂₃π₂₄π₅]

${\mathrm{\&omega;}}_{21}=\left[\begin{array}{c}\left(\frac{{\&PartialD;}^{2}T}{{\&PartialD;V}^2}\right)\cos \mathrm{\&alpha;}-\frac{{\&PartialD;}^{2}D}{{\&PartialD;V}^2}\\ 0\\ -\left(\frac{\&PartialD;T}{\&PartialD;V}\right)\sin \mathrm{\&alpha;}-\frac{{\&PartialD;}^{2}D}{\&PartialD;V\&PartialD;\mathrm{\&alpha;}}\\ \left(\frac{{\&PartialD;}^{2}T}{\&PartialD;V\&PartialD;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\\ 0\end{array}\right],{\mathrm{\&omega;}}_{22}=\left[\begin{array}{c}0\\ \frac{\mathrm{m\&mu;}\sin \mathrm{\&gamma;}}{r^2}\\ 0\\ 0\\ \frac{2\mathrm{m\&mu;}\cos \mathrm{\&gamma;}}{r^3}\end{array}\right],{\mathrm{\&omega;}}_{23}=\left[\begin{array}{c}-\left(\frac{\&PartialD;T}{\&PartialD;V}\right)\sin \mathrm{\&alpha;}-\left(\frac{{\&PartialD;}^{2}D}{\&PartialD;V\&PartialD;\mathrm{\&alpha;}}\right)\\ 0\\ -T\cos \mathrm{\&alpha;}\left(\frac{{\&PartialD;}^{2}D}{{\&PartialD;\mathrm{\&alpha;}}^2}\right)\\ -\left(\frac{\&PartialD;T}{\&PartialD;\mathrm{\&beta;}}\right)\sin \mathrm{\&alpha;}\\ 0\end{array}\right]$

${\mathrm{\&omega;}}_{24}=\left[\begin{array}{c}\left(\frac{{\&PartialD;}^{2}T}{\&PartialD;V\&PartialD;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\\ 0\\ -\left(\frac{\&PartialD;T}{\&PartialD;\mathrm{\&beta;}}\right)\sin \mathrm{\&alpha;}\\ 0\\ 0\end{array}\right],{\mathrm{\&omega;}}_{25}=\left[\begin{array}{c}0\\ \frac{2\mathrm{m\&mu;}\cos \mathrm{\&gamma;}}{r^3}\\ 0\\ 0\\ \frac{-6\mathrm{m\&mu;}\sin \mathrm{\&gamma;}}{r^4}\end{array}\right]$

${\mathrm{\&pi;}}_{21}=\left[\begin{array}{c}\frac{{\&PartialD;}^{2}L/{\&PartialD;V}^2+({\&PartialD;}^{2}T/{\&PartialD;V}^2)\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{2[\&PartialD;L/\&PartialD;V+(\&PartialD;T/\&PartialD;V)\sin \mathrm{\&alpha;}]}{{\mathrm{mV}}^2}+\frac{2(L+T\sin \mathrm{\&alpha;})}{{\mathrm{mV}}^3}-\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^3{r}^2}\\ -\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{V^2{r}^2}-\frac{\sin \mathrm{\&gamma;}}{r}\\ \frac{({\&PartialD;}^{2}L/\&PartialD;\mathrm{\&alpha;}\&PartialD;V)+(\&PartialD;T/\&PartialD;V)\cos \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{\&PartialD;L/\&PartialD;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ \frac{({\&PartialD;}^{2}T/\&PartialD;\mathrm{\&beta;}\&PartialD;V)\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ -\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^3}-\frac{\cos \mathrm{\&gamma;}}{r^2}\end{array}\right]$

${\mathrm{\&pi;}}_{22}=\left[\begin{array}{c}-\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{V^2{r}^2}-\frac{\sin \mathrm{\&gamma;}}{r}\\ \frac{\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}-\frac{V\cos \mathrm{\&gamma;}}{r}\\ 0\\ 0\\ -\frac{2\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}+\frac{V\sin \mathrm{\&gamma;}}{r^2}\end{array}\right],{\mathrm{\&pi;}}_{23}=\left[\begin{array}{c}\frac{({\&PartialD;}^{2}L/\&PartialD;V\&PartialD;\mathrm{\&alpha;})+(\&PartialD;T/\&PartialD;V)\cos \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{\&PartialD;L/\&PartialD;\mathrm{\&alpha;}+T\cos }{{\mathrm{mV}}^2}\\ 0\\ \frac{{\&PartialD;}^{2}L/\&PartialD;{\mathrm{\&alpha;}}^2-T\sin }{\mathrm{mV}}\\ \frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\cos \mathrm{\&alpha;}}{\mathrm{mV}}\\ 0\end{array}\right]$

${\mathrm{\&pi;}}_{24}=\left[\begin{array}{c}\frac{({\&PartialD;}^{2}T/\&PartialD;V\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{\mathrm{mV}}-\frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ 0\\ \frac{(\&PartialD;T/\&PartialD;\mathrm{\&beta;})\cos \mathrm{\&alpha;}}{\mathrm{mV}}\\ 0\\ 0\end{array}\right],{\mathrm{\&pi;}}_{25}=\left[\begin{array}{c}-\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^3}-\frac{\cos \mathrm{\&gamma;}}{r^2}\\ -\frac{2\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^3}+\frac{V\sin \mathrm{\&gamma;}}{r^2}\\ 0\\ 0\\ -\frac{6\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^4}+\frac{2V\cos \mathrm{\&gamma;}}{r^3}\end{array}\right]$

ω₁、π₁、Ω₂、Π₂的表达式中关于L、T、D的一阶、二阶偏导数可根据下列关系计算得到：

L=ρV²SC_L2,C_L=0.6203α

D=ρV²SC_D2,C_D=0.6450α²+0.0043378α+0.003772

T=ρV²SC_T2, $C_T=\left\{\begin{array}{cc}0.02576\mathrm{\&beta;}& \mathrm{if\&beta;}\&le;1\\ 0.0224+0.00336\mathrm{\&beta;}& \mathrm{if\&beta;}>1\end{array}\right.$

总的舵偏角δ_e与左、右两个升降舵的舵偏角δ_e1和δ_e2的关系为：δ_e=d₁δ_e1+d₂δ_e2，其中d₁和d₂为组合关系系数。

步骤3，考虑输入-输出模型中升降舵发生卡死故障的情况，建立升降舵的故障模型：

${\mathrm{\&delta;}}_{\mathrm{ej}}=v_j+{\mathrm{\&sigma;}}_j({\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j-v_j),j=\mathrm{1,2}$

其中表示升降舵卡死故障的大小；v_j是第j个升降舵的控制输入；σ_j代表故障模式，当第j个升降舵发生故障时，σ_j=1，否则σ_j=0。

步骤4，基于速度追踪和高度追踪误差选择两个滑模面s₁和s₂如下：

s₁=λ₁₀e₁₀+λ₁₁e₁₁+λ₁₂e₁₂+e₁₃

s₂=λ₂₀e₂₀+λ₂₁e₂₁+λ₂₂e₂₂+λ₂₃e₂₃+e₂₄

其中e₁₁=V-V_d，e₁₀=∫e₁₁，e₂₁=h-h_d，e₂₀=∫e₂₁， $e_{22}={\stackrel{.}{e}}_{21}=\stackrel{.}{h}-{\stackrel{.}{h}}_d,e_{23}={\stackrel{..}{e}}_{21}=\stackrel{..}{h}-{\stackrel{..}{h}}_d,e_{24}={\stackrel{...}{e}}_{21}=\stackrel{...}{h}-{\stackrel{...}{h}}_d;$ V_d和h_d为期望的速度和高度信号。λ₁₀,...,λ₁₂,λ₂₀,...,λ₂₃为正常数，其选择满足使多项式s³+λ₁₂s²+λ₁₁s+λ₁₀和s⁴+λ₂₃s³+λ₂₂s²+λ₂₁s+λ₂₀为Hurwitz多项式。

对两个滑模面求导得到：

${\stackrel{.}{s}}_1={\mathrm{\&lambda;}}_{10}{e}_{11}+{\mathrm{\&lambda;}}_{11}{e}_{12}+{\mathrm{\&lambda;}}_{12}{e}_{13}+\stackrel{...}{V}-{\stackrel{...}{V}}_d$

${\stackrel{.}{s}}_2={\mathrm{\&lambda;}}_{20}{e}_{21}+{\mathrm{\&lambda;}}_{21}{e}_{22}+{\mathrm{\&lambda;}}_{22}{e}_{23}+{\mathrm{\&lambda;}}_{23}{e}_{24}+h^{\left(4\right)}-h_d^{\left(4\right)}$

将和h⁽⁴⁾替换为前述动态输入-输出模型，得到：

${\stackrel{.}{s}}_1={\stackrel{...}{V}}_d+f_v+{\mathrm{\&lambda;}}_{10}{e}_{11}+{\mathrm{\&lambda;}}_{11}{e}_{12}+{\mathrm{\&lambda;}}_{12}{e}_{13}+b_{11}{\mathrm{\&beta;}}_c+b_{12}{\mathrm{\&delta;}}_e$

${\stackrel{.}{s}}_2=-h_d^{\left(4\right)}+f_h+{\mathrm{\&lambda;}}_{20}{e}_{21}+{\mathrm{\&lambda;}}_{21}{e}_{22}+{\mathrm{\&lambda;}}_{22}{e}_{23}+{\mathrm{\&lambda;}}_{23}{e}_{24}+b_{21}{\mathrm{\&beta;}}_c+b_{22}{\mathrm{\&delta;}}_e$

可进一步写为：

$\left[\begin{array}{c}{\stackrel{.}{s}}_1\\ {\stackrel{.}{s}}_2\end{array}\right]=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]+\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\left[\begin{array}{c}{\mathrm{\&beta;}}_c\\ {\mathrm{\&delta;}}_e\end{array}\right]$

其中，

$f_1=-{\stackrel{...}{V}}_d+f_v+{\mathrm{\&lambda;}}_{10}{e}_{11}+{\mathrm{\&lambda;}}_{11}{e}_{12}+{\mathrm{\&lambda;}}_{12}{e}_{13}$

$f_2-h_d^{\left(4\right)}+f_h+{\mathrm{\&lambda;}}_{20}{e}_{21}+{\mathrm{\&lambda;}}_{21}{e}_{22}+{\mathrm{\&lambda;}}_{22}{e}_{23}+{\mathrm{\&lambda;}}_{23}{e}_{24}$

步骤5，根据滑模控制思想，设计标称控制器如下：

$\left[\begin{array}{c}{\mathrm{\&beta;}}_c\\ v_0\end{array}\right]=B^{-1}\left[\begin{array}{cc}-f_1& -a_1\mathrm{sgn}\left(s_1\right)\\ -f_2& -a_2\mathrm{sgn}\left(s_2\right)\end{array}\right],B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]$

其中v₀表示升降舵正常时的标称偏转量；a₁和a₂为两个正常数，决定滑动模态的收敛速度；sgn()为符号函数，表达式如下：

$\mathrm{sgn}\left(x\right)=\begin{array}{}\end{array}\left\{\begin{array}{cc}1& x>0\\ 0& x=0\\ -1& x<0\end{array}\right.$

采用上述控制方法，当升降舵正常工作，即δ_e=v₀时，选择Lyapunov函数V=s^Ts2，对其求导得到说明闭环系统稳定且滑模面能够在有限时间内趋于零。

步骤6，设计容错控制器结构为：

v_j=k1_jv0+k2_j,j=1,2

其中k_1j和k_2j为标称控制器参数。

结合步骤3中的故障模型，将其代入步骤2中的速度和高度动态模型，得到：

$\stackrel{...}{V}=f_v+b_{11}{\mathrm{\&beta;}}_c+b_{12}\mathrm{d\&sigma;}\stackrel{-}{u}+b_{12}d(I-\mathrm{\&sigma;}){\left[k_{11}{k}_{12}\right]}^T{v}_0+b_{12}d(I-\mathrm{\&sigma;}){\left[k_{21}{k}_{22}\right]}^T$

$h^{\left(4\right)}=f_h+b_{21}{\mathrm{\&beta;}}_c+b_{22}\mathrm{d\&sigma;}\stackrel{-}{u}+b_{22}d(I-\mathrm{\&sigma;}){\left[k_{11}{k}_{12}\right]}^T{v}_0+b_{22}d(I-\mathrm{\&sigma;}){\left[k_{21}{k}_{22}\right]}^T$

其中d=[d₁,d₂]，σ=diag(σ₁,σ₂)，

为实现故障补偿，k_1j和k_2j应满足如下匹配条件：

d(I-σ)[k₁₁k₁₂]^T=1

$d(I-\mathrm{\&sigma;}){\left[k_{21}{k}_{22}\right]}^T=-\mathrm{d\&sigma;}\stackrel{-}{u}$

由于升降舵故障模式σ、大小和发生时间未知，我们无法从匹配条件中直接解出k₁₁、k₁₂、k₂₁和k₂₂，因此设计自适应容错控制器：

$v_j={\hat{k}}_{1j}{v}_0+{\hat{k}}_{2j},j=\mathrm{1,2}$

其中为k_1j和k_2j的估计值。

定义参数估计误差得到自适应容错控制器作用下的滑模面动态如下：

${\stackrel{.}{s}}_1=f_1+b_{11}{\mathrm{\&beta;}}_c+b_{12}{v}_0+b_{12}\underset{j\&NotEqual;j_p}{\overset{2}{\mathrm{\&Sigma;}}}{d}_j({\tilde{k}}_{1j}{v}_0+{\tilde{k}}_{2j})$

${\stackrel{.}{s}}_2=f_2+b_{21}{\mathrm{\&beta;}}_c+b_{22}{v}_0+b_{22}\underset{j\&NotEqual;j_p}{\overset{2}{\mathrm{\&Sigma;}}}{d}_j({\tilde{k}}_{1j}{v}_0+{\tilde{k}}_{2j})$

选择Lyapunov函数：

$V=\frac{1}{2}{s}_1^2+\frac{1}{2}{s}_2^2+\frac{\left|d_j\right|}{2}\underset{j\&NotEqual;j_p}{\overset{2}{\mathrm{\&Sigma;}}}{\mathrm{\&Gamma;}}_{1j}^{-1}{\tilde{k}}_{1j}^2+\frac{\left|d_j\right|}{2}\underset{j\&NotEqual;j_p}{\overset{2}{\mathrm{\&Sigma;}}}{\mathrm{\&Gamma;}}_{2j}^{-1}{\tilde{k}}_{2j}^2$

其导数计算如下：

$\stackrel{.}{V}\begin{array}{c}=a_1\left|s_1\right|-a_2\left|s_2\right|+\underset{j\&NotEqual;j_p}{\overset{2}{\mathrm{\&Sigma;}}}{\tilde{k}}_{1j}(b_{12}{d}_j{s}_1{v}_0+b_{22}{d}_j{s}_2{v}_0+|d_j\left|{\mathrm{\&Gamma;}}_{1j}^{-1}{\stackrel{\stackrel{.}{~}}{k}}_{1j}\right)\\ +\underset{j\&NotEqual;j_p}{\overset{2}{\mathrm{\&Sigma;}}}{\tilde{k}}_{2j}(b_{12}{d}_j{s}_1{v}_0+b_{22}{d}_j{s}_2{v}_0+|d_j\left|{\mathrm{\&Gamma;}}_{2j}^{-1}{\stackrel{\stackrel{.}{~}}{k}}_{2j}\right)\end{array}$

选择的参数自适应律如下：

${\stackrel{\stackrel{.}{^}}{k}}_{1j}=-\mathrm{sgn}\left(d_j\right){\mathrm{\&Gamma;}}_{1j}(b_{12}{s}_1+b_{22}{s}_2)v_0$

${\stackrel{\stackrel{.}{^}}{k}}_{2j}=-\mathrm{sgn}\left(d_j\right){\mathrm{\&Gamma;}}_{2j}(b_{12}{s}_1+b_{22}{s}_2)$

其中Γ_1j和Γ_2j为两个正常数，决定了参数自适应律的学习速率，得到

$\stackrel{.}{V}=-a_1\left|s_1\right|-a_2\left|s_2\right|\&le;0$

闭环系统稳定且滑模面能够在有限时间内趋于零。

为减小滑模控制产生的抖振，β_c、v₀分别由下式确定：

$\left[\begin{array}{c}{\mathrm{\&beta;}}_c\\ v_0\end{array}\right]=B^{-1}\left[\begin{array}{cc}-f_1& -a_1\mathrm{sat}\left(\frac{s_1}{\mathrm{\&epsiv;}}\right)\\ -f_2& -a_2\mathrm{sat}\left(\frac{s_2}{\mathrm{\&epsiv;}}\right)\end{array}\right],B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]$

其中，ε为一个小的正常数；sat(·)为饱和函数，其定义如下：

$\mathrm{sat}\left(x\right)=\left\{\begin{array}{cc}x& \mathrm{if}\left|x\right|\&le;1\\ \mathrm{sgn}\left(x\right)& \mathrm{if}\left|x\right|>1\end{array}\right..$

下面通过仿真验证本发明的有效性。仿真参数如下：

飞行器参数：m=9375slug，I_yy=7×10⁶slug·ft²，S=3603ft²，μ=3.31929×10^-11，R_E=20902244ft，ρ=0.24325×10^-4slug/ft²，c_e=0.0292，d₁=d₂=0.5，

状态初始值：[V(0),γ(0),θ_p(0),q(0),h(0)]=[15060ft/s,0,0,0,110000ft]；高度变化信号：速度变化信号：h_c=200ft，V_c=100ft/s，h_f(0)=V_f(0)=0；滑模控制器参数：[λ₁₀,λ₁₁,λ₁₂]=[64,48,12]，[λ₂₀,λ₂₁,λ₂₂,λ₂₃]=[4,8,8,4]，ε=0.1；控制器参数初始值： 自适应律参数：Γ₁₁=Γ₁₂=Γ₂₁=Γ₂₂=1×10^-5；舵偏角组合系数d₁=d₂=0.5。

仿真中考虑第二个升降舵在第12秒发生卡死故障：

${\mathrm{\&delta;}}_{e2}\left(t\right)=\left\{\begin{array}{cc}{v}_0\left(t\right)& t<12\sec \\ 0.2\mathrm{rad}& t\&GreaterEqual;12\sec \end{array}\right.$

仿真结果说明：

如图2所示是输出和输入响应图；其中，图2（a）是速度追踪响应图；图2（b）是高度追踪响应图；图2（c）是推力响应图；图2（d）是合成舵偏角响应图。从图2（a）和（b）可见，无论是在第二个升降舵发生卡死故障之前还是之后，本发明所提出的自适应容错控制方法均能获得满意的速度和高度追踪效果。由图2（d）可见合成舵偏角δ_e在12秒处发生了一定的超调，这是正常升降舵补偿故障升降舵效力时产生的瞬态响应。

如图3所示是状态响应图；其中，图3（a）是航迹角响应图；图3（b）是俯仰角响应图；图3（c）是攻角响应图；图3（d）是俯仰角速率响应图。从图中可观测到采用本发明提出的容错控制方法，在故障发生后，这些系统关键状态值均在允许值范围内。

如图4所示是两个升降舵的偏转角图；其中，图4（a）是左升降舵的偏转角；图4（b）是右升降舵的偏转角。从图4(b)中可见第二个升降舵在第12秒时卡死在0.2rad的位置，从图4(a)可见第一个升降舵在12秒时开始自动调整其偏转，通过增大其偏转角度补偿发生故障的升降舵损失的控制效力。

如图5所示是控制器参数的自适应调整图；其中，图5（a）是参数的自适应调整图；图5（b）是参数的自适应调整图；图5（c）是参数的自适应调整图；图5（d）是参数的自适应调整图。从图5（c）和（d）中可观察到参数和在12秒处发生了明显变化，从图5（a）和（b）可见和也发生了调整，正是通过自适应调整控制参数改变了正常升降舵的舵偏角，补偿了故障升降舵损失的控制效力。

Claims (2)

1.一种高超声速飞行器升降舵故障的自适应容错控制方法，其特征在于，包括以下步骤：

步骤1，建立高超声速飞行器巡航段纵向运动学及动力学模型，所述模型表示为：

$\stackrel{\&CenterDot;}{V}=\frac{Tcos\mathrm{\&alpha;}-D}{m}-\frac{\mathrm{\&mu;}sin\mathrm{\&gamma;}}{r^2}$

$\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=\frac{L+Tsin\mathrm{\&alpha;}}{mV}-\frac{(\mathrm{\&mu;}-V^{2}r)cos\mathrm{\&gamma;}}{{\mathrm{Vr}}^2}$

$\stackrel{\&CenterDot;}{h}=Vsin\mathrm{\&gamma;}$

$\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}=q-\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}$

$\stackrel{\&CenterDot;}{q}=M_{yy}/I_{yy}$

其中V、γ、h、α、q分别表示高超声速飞行器的速度、航迹角、高度、攻角和俯仰角速率；L、D、T分别为升力、阻力和推力；m、μ、r分别表示飞行器的质量、重力常数、飞行器到地心的距离，其中r＝R_E+h，R_E为地球半径；I_yy、M_yy分别表示俯仰转动惯量和俯仰转动力矩；

步骤2，建立控制系统输入-输出模型；

针对所述高超声速飞行器巡航段纵向运动学及动力学模型，以高超声速飞行器的速度V和高度h为输出、油门开度β_c和总的升降舵舵偏角δ_e作为控制输入，建立输入-输出模型，其表达式如下：

$\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{V}=f_V+b_{11}{\mathrm{\&beta;}}_c+b_{12}{\mathrm{\&delta;}}_e$

h⁽⁴⁾＝f_h+b₂₁β_c+b₂₂δ_e

其中f_V、f_h、b₁₁、b₁₂、b₂₁、b₂₂分别是关于V、γ、h、α、q和飞行器气动参数的非线性函数；总的升降舵舵偏角δ_e与左、右两个升降舵的舵偏角δ_e1和δ_e2的关系为：

δ_e＝d₁δ_e1+d₂δ_e2，其中d₁和d₂为组合关系系数；

步骤3，考虑步骤2所述输入-输出模型升降舵发生卡死故障的情况，建立升降舵的故障模型，其表达式如下：

${\mathrm{\&delta;}}_{ej}=v_j+{\mathrm{\&sigma;}}_j({\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j-v_j),j=1,2$

其中，δ_ej表示第j个升降舵的舵偏角，当j＝1时，δ_e1表示左升降舵的舵偏角；

当j＝2时，δ_e2表示右升降舵的舵偏角；表示升降舵卡死故障的大小；v_j表示第j个升降舵的控制输入；σ_j代表故障模式，当第j个升降舵发生故障时，σ_j＝1，否则σ_j＝0；

步骤4，基于速度追踪误差和高度追踪误差定义两个滑模面s₁和s₂：

s₁＝λ₁₀e₁₀+λ₁₁e₁₁+λ₁₂e₁₂+e₁₃

s₂＝λ₂₀e₂₀+λ₂₁e₂₁+λ₂₂e₂₂+λ₂₃e₂₃+e₂₄

其中，e₁₁＝V-V_d，e₁₀＝∫e₁₁，e₂₁＝h-h_d，e₂₀＝∫e₂₁， V_d和h_d为期望的速度和高度信号；λ₁₀,λ₁₁,λ₁₂,λ₂₀,λ₂₁,λ₂₂,λ₂₃为正常数，λ₁₀,λ₁₁,λ₁₂,λ₂₀,λ₂₁,λ₂₂,λ₂₃的取值满足：

s³+λ₁₂s²+λ₁₁s+λ₁₀和s⁴+λ₂₃s³+λ₂₂s²+λ₂₁s+λ₂₀

为Hurwitz多项式；

步骤5，确定标称控制输入，其表达式如下：

$\begin{array}{cc}\left[\begin{array}{c}{\mathrm{\&beta;}}_c\\ v_0\end{array}\right]=B^{-1}\left[\begin{array}{cc}-f_1& -a_1\mathrm{sgn}\left(s_1\right)\\ -f_2& -a_2\mathrm{sgn}\left(s_2\right)\end{array}\right],& B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\end{array}$

其中v₀为总的升降舵舵偏角δ_e的标称偏转量；a₁和a₂为两个正常数，决定滑动模态的收敛速度；f₁和f₂分别由下式确定：

$f_1=-{\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{V}}_d+f_V+{\mathrm{\&lambda;}}_{10}{e}_{11}+{\mathrm{\&lambda;}}_{11}{e}_{12}+{\mathrm{\&lambda;}}_{12}{e}_{13}$

$f_2=h_d^{\left(4\right)}+f_h+{\mathrm{\&lambda;}}_{20}{e}_{21}+{\mathrm{\&lambda;}}_{21}{e}_{22}+{\mathrm{\&lambda;}}_{22}{e}_{23}+{\mathrm{\&lambda;}}_{23}{e}_{24}$

sgn(·)为符号函数，其定义如下：

$\mathrm{sgn}\left(x\right)=\left\{\begin{array}{cc}1& x>0\\ 0& x=0\\ -1& x<0\end{array}\right.$

或者为减小滑模控制产生的抖振，β_c、v₀分别由下式确定：

$\begin{array}{cc}\left[\begin{array}{c}{\mathrm{\&beta;}}_c\\ v_0\end{array}\right]=B^{-1}\left[\begin{array}{cc}-f_1& -a_{1}sat\left(\frac{s_1}{\mathrm{\&epsiv;}}\right)\\ -f_2& -a_{2}sat\left(\frac{s_2}{\mathrm{\&epsiv;}}\right)\end{array}\right],& B=\left[\begin{array}{cc}{b}_{11}& b_{12}\\ b_{21}& b_{22}\end{array}\right]\end{array}$

其中，ε为正常数，其取值决定了滑模控制中状态轨迹边界层的厚度；sat(·)为饱和函数，其定义如下：

$sat\left(x\right)=\left\{\begin{array}{ccc}x& if& \left|x\right|\&le;1\\ \mathrm{sgn}\left(x\right)& if& \left|x\right|>1\end{array}\right.$

步骤6，基于步骤5给出的总的升降舵舵偏角δ_e的标称偏转量v₀，确定自适应容错控制律：

$v_j={\hat{k}}_{1j}{v}_0+{\hat{k}}_{2j},j=1,2$

其中分别为容错控制器的参数估计值，分别由下式确定：

${\stackrel{\&CenterDot;}{\hat{k}}}_{1j}=-\mathrm{sgn}\left(d_j\right){\mathrm{\&Gamma;}}_{1j}(b_{12}{s}_1+b_{22}{s}_2)v_0$

${\stackrel{\&CenterDot;}{\hat{k}}}_{2j}=-\mathrm{sgn}\left(d_j\right){\mathrm{\&Gamma;}}_{2j}(b_{12}{s}_1+b_{22}{s}_2)$

其中Γ_1j和Γ_2j为两个正常数，决定了参数自适应律的学习速率；d_j,j＝1,2分别表示左、右两个升降舵的舵偏角δ_e1和δ_e2的组合关系系数，其取值为d₁＝d₂＝0.5。

2.根据权利要求1所述的一种高超声速飞行器升降舵故障的自适应容错控制方法，其特征在于，步骤2中，所述非线性函数f_V、f_h、b₁₁、b₁₂、b₂₁、b₂₂的具体形式分别如下：

$f_V=\frac{{\mathrm{\&omega;}}_1{\stackrel{\&CenterDot;\&CenterDot;}{x}}_0+{\stackrel{\&CenterDot;}{x}}^T{\mathrm{\&Omega;}}_2\stackrel{\&CenterDot;}{x}}{m}$

$\begin{array}{c}{f}_h=3\stackrel{\&CenterDot;\&CenterDot;}{V}\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}cos\mathrm{\&gamma;}-3\stackrel{\&CenterDot;}{V}{\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}^2\sin \mathrm{\&gamma;}+3\stackrel{\&CenterDot;}{V}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}}\cos \mathrm{\&gamma;}-3V\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}}\sin \mathrm{\&gamma;}-V{\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}}^3\cos \mathrm{\&gamma;}+\frac{({\mathrm{\&omega;}}_1{\stackrel{\&CenterDot;\&CenterDot;}{x}}_0+{\stackrel{\&CenterDot;}{x}}^T{\mathrm{\&Omega;}}_2\stackrel{\&CenterDot;}{x})sin\mathrm{\&gamma;}}{m}\\ +V\cos \mathrm{\&gamma;}({\mathrm{\&pi;}}_1{\stackrel{\&CenterDot;\&CenterDot;}{x}}_0+{\stackrel{\&CenterDot;}{x}}^T{\mathrm{\&Pi;}}_2\stackrel{\&CenterDot;}{x})\end{array}$

$\begin{array}{cc}{b}_{11}=\left(\frac{{\mathrm{\&rho;V}}^2{\mathrm{Sc}}_{\mathrm{\&beta;}}{\mathrm{\&omega;}}_n^2}{2m}\right)cos\mathrm{\&alpha;},& b_{12}=-\left(\frac{c_e{\mathrm{\&rho;V}}^{2}S\stackrel{\&OverBar;}{c}}{2{\mathrm{mI}}_{yy}}\right)(Tsin\mathrm{\&alpha;}+\frac{\&part;D}{\&part;\mathrm{\&alpha;}})\end{array}$

$\begin{array}{cc}{b}_{21}=\left(\frac{{\mathrm{\&rho;V}}^2{\mathrm{Sc}}_{\mathrm{\&beta;}}{\mathrm{\&omega;}}_n^2}{2m}\right)\sin (\mathrm{\&alpha;}+\mathrm{\&gamma;}),& b_{22}=\left(\frac{c_e{\mathrm{\&rho;V}}^{2}S\stackrel{\&OverBar;}{c}}{2{\mathrm{mI}}_{yy}}\right)\&lsqb;T\cos (\mathrm{\&alpha;}+\mathrm{\&gamma;})+\left(\frac{\&part;L}{\&part;\mathrm{\&alpha;}}\right)\cos \mathrm{\&gamma;}-\left(\frac{\&part;D}{\&part;\mathrm{\&alpha;}}\right)\sin \mathrm{\&gamma;}\&rsqb;\end{array}$

其中，

x＝[Vγαβh]^T, ${\stackrel{\&CenterDot;\&CenterDot;}{x}}_0={\left[\begin{array}{ccccc}\stackrel{\&CenterDot;\&CenterDot;}{V}& \stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}}& {\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}}_0& {\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&beta;}}}_0& \stackrel{\&CenterDot;\&CenterDot;}{h}\end{array}\right]}^T$

$\begin{array}{cc}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&alpha;}}}_0=\frac{1}{2}{\mathrm{\&rho;V}}^{2}S\stackrel{\&OverBar;}{c}\&lsqb;C_M\left(\mathrm{\&alpha;}\right)+C_M\left(q\right)-c_e\mathrm{\&alpha;}\&rsqb;/I_{yy}-\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&gamma;}},& {\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&beta;}}}_0=-2{\mathrm{\&zeta;\&omega;}}_n\stackrel{\&CenterDot;}{\mathrm{\&beta;}}-{\mathrm{\&omega;}}_n^2\mathrm{\&beta;}\end{array}$

C_M(α)＝-0.035α²+0.036617α+5.3261×10^-6

$C_M\left(q\right)=(\stackrel{\&OverBar;}{c}/2V)q(-6.796{\mathrm{\&alpha;}}^2+0.3015\mathrm{\&alpha;}-0.2289)$

C_M(δ_e)＝c_e(δ_e-α)

其中C_M(α)、C_M(δ_e)和为C_M(q)为俯仰力矩系数；S为参考气动面积；ρ为大气密度；ω_n、ζ、β为发动机动态模型的自然频率、阻尼比和输出；α₀和β₀为转换模型引入的辅助变量；为平均气动翼弦；c_e和c_β为计算俯仰力矩系数的相关参数；ω₁、π₁、Ω₂、Π₂为推导V和h的高阶微分所产生的向量及矩阵：

${\mathrm{\&omega;}}_1=\left[\begin{array}{ccccc}\left(\frac{\&part;T}{\&part;V}\right)cos\mathrm{\&alpha;}-\frac{\&part;D}{\&part;V}& \frac{-m\mathrm{\&mu;}cos\mathrm{\&gamma;}}{r^2}& -Tsin\mathrm{\&alpha;}\frac{\&part;D}{\&part;\mathrm{\&alpha;}}& \left(\frac{\&part;T}{\&part;\mathrm{\&beta;}}\right)cos\mathrm{\&alpha;}& \frac{2m\mathrm{\&mu;}sin\mathrm{\&gamma;}}{r^3}\end{array}\right]$

${\mathrm{\&pi;}}_1^T=\left[\begin{array}{c}\frac{\&part;L/\&part;V+(\&part;T/\&part;V)sin\mathrm{\&alpha;}}{mV}-\frac{L+Tsin\mathrm{\&alpha;}}{{\mathrm{mV}}^2}+\frac{\mathrm{\&mu;}cos\mathrm{\&gamma;}}{V^2{r}^2}+\frac{cos\mathrm{\&gamma;}}{r}\\ \frac{\mathrm{\&mu;}sin\mathrm{\&gamma;}}{{\mathrm{Vr}}^2}-\frac{Vsin\mathrm{\&gamma;}}{r}\\ \frac{\&part;L/\&part;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{mV}\\ \frac{(\&part;T/\&part;\mathrm{\&beta;})sin\mathrm{\&alpha;}}{mV}\\ \frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^3}-\frac{V\cos \mathrm{\&gamma;}}{r^2}\end{array}\right]$

Ω₂＝[ω₂₁ω₂₂ω₂₃ω₂₄ω₂₅],Π₂＝[π₂₁π₂₂π₂₃π₂₄π₅]

$\begin{array}{ccc}{\mathrm{\&omega;}}_{21}=\left[\begin{array}{c}\left(\frac{{\&part;}^{2}T}{\&part;V^2}\right)\cos \mathrm{\&alpha;}-\frac{{\&part;}^{2}D}{\&part;V^2}\\ -\left(\frac{\&part;T}{\&part;V}\right)\sin \mathrm{\&alpha;}-\frac{{\&part;}^{2}D}{\&part;V\&part;\mathrm{\&alpha;}}\\ \left(\frac{{\&part;}^{2}T}{\&part;V\&part;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\\ 0\end{array}\right],& {\mathrm{\&omega;}}_{22}=\left[\begin{array}{c}0\\ \frac{m\mathrm{\&mu;}sin\mathrm{\&gamma;}}{r^2}\\ 0\\ 0\\ \frac{3m\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{r^3}\end{array}\right],& {\mathrm{\&omega;}}_{23}=\left[\begin{array}{c}-\left(\frac{\&part;T}{\&part;V}\right)\sin \mathrm{\&alpha;}-\left(\frac{{\&part;}^{2}D}{\&part;V\&part;\mathrm{\&alpha;}}\right)\\ -T\cos \mathrm{\&alpha;}-\left(\frac{{\&part;}^{2}D}{\&part;{\mathrm{\&alpha;}}^2}\right)\\ -\left(\frac{\&part;T}{\&part;\mathrm{\&beta;}}\right)\sin \mathrm{\&alpha;}\\ 0\end{array}\right]\end{array}$

$\begin{array}{cc}{\mathrm{\&omega;}}_{24}=\left[\begin{array}{c}\left(\frac{{\&part;}^{2}T}{\&part;V\&part;\mathrm{\&beta;}}\right)\cos \mathrm{\&alpha;}\\ 0\\ -\left(\frac{\&part;T}{\&part;\mathrm{\&beta;}}\right)\sin \mathrm{\&alpha;}\\ 0\\ 0\end{array}\right],& {\mathrm{\&omega;}}_{25}=\left[\begin{array}{c}0\\ \frac{2m\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{r^3}\\ 0\\ 0\\ \frac{-6m\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{r^4}\end{array}\right]\end{array}$

${\mathrm{\&pi;}}_{21}=\left[\begin{array}{c}\frac{{\&part;}^{2}L/\&part;V^2+({\&part;}^{2}T/\&part;V^2)\sin \mathrm{\&alpha;}}{mV}-\frac{2\&lsqb;\&part;L/\&part;V+(\&part;T/\&part;V)\sin \mathrm{\&alpha;}\&rsqb;}{{\mathrm{mV}}^2}+\frac{2(L+T\sin \mathrm{\&alpha;})}{{\mathrm{mV}}^3}-\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^3{r}^2}\\ -\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{V^2{r}^2}-\frac{\sin \mathrm{\&gamma;}}{r}\\ \frac{({\&part;}^{2}L/\&part;\mathrm{\&alpha;}\&part;V)+(\&part;T/\&part;V)\cos \mathrm{\&alpha;}}{mV}-\frac{\&part;L/\&part;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ \frac{({\&part;}^{2}T/\&part;\mathrm{\&beta;}\&part;V)sin\mathrm{\&alpha;}}{mV}-\frac{(\&part;T/\&part;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ -\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^3}-\frac{\cos \mathrm{\&gamma;}}{r^2}\end{array}\right]$

$\begin{array}{cc}{\mathrm{\&pi;}}_{22}=\left[\begin{array}{c}-\frac{\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{V^2{r}^2}-\frac{\sin \mathrm{\&gamma;}}{r}\\ \frac{\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^2}-\frac{V\cos \mathrm{\&gamma;}}{r}\\ 0\\ 0\\ -\frac{2\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^3}+\frac{V\sin \mathrm{\&gamma;}}{r^2}\end{array}\right],& {\mathrm{\&pi;}}_{23}=\left[\begin{array}{c}\frac{({\&part;}^{2}L/\&part;V\&part;\mathrm{\&alpha;})+(\&part;T/\&part;V)\cos \mathrm{\&alpha;}}{mV}-\frac{\&part;L/\&part;\mathrm{\&alpha;}+T\cos \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ 0\\ \frac{{\&part;}^{2}L/\&part;{\mathrm{\&alpha;}}^2-Tsin\mathrm{\&alpha;}}{mV}\\ \frac{(\&part;T/\&part;\mathrm{\&beta;})\cos \mathrm{\&alpha;}}{mV}\\ 0\end{array}\right]\end{array}$

$\begin{array}{cc}{\mathrm{\&pi;}}_{24}=\left[\begin{array}{c}\frac{({\&part;}^{2}T/\&part;V\&part;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{mV}-\frac{(\&part;L/\&part;\mathrm{\&beta;})\sin \mathrm{\&alpha;}}{{\mathrm{mV}}^2}\\ 0\\ \frac{(\&part;T/\&part;\mathrm{\&beta;})\cos \mathrm{\&alpha;}}{mV}\\ 0\\ 0\end{array}\right],& {\mathrm{\&pi;}}_{25}=\left[\begin{array}{c}-\frac{2\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{V^2{r}^3}-\frac{\cos \mathrm{\&gamma;}}{r^2}\\ -\frac{2\mathrm{\&mu;}\sin \mathrm{\&gamma;}}{{\mathrm{Vr}}^3}+\frac{V\sin \mathrm{\&gamma;}}{r^2}\\ 0\\ 0\\ -\frac{6\mathrm{\&mu;}\cos \mathrm{\&gamma;}}{{\mathrm{Vr}}^4}+\frac{2V\cos \mathrm{\&gamma;}}{r^3}\end{array}\right]\end{array}$

ω₁、π₁、Ω₂、Π₂的表达式中L、T、D的一阶、二阶偏导数计算式如下：

L＝ρV²SC_L/2,C_L＝0.6203α

D＝ρV²SC_D/2,C_D＝0.6450α²+0.0043378α+0.003772

T＝ρV²SC_T/2, $C_T=\left\{\begin{array}{ccc}0.02576\mathrm{\&beta;}& if& \mathrm{\&beta;}\&le;1\\ 0.0224+0.00336\mathrm{\&beta;}& if& \mathrm{\&beta;}>1\end{array}\right..$