基于滚动时域优化的变后掠角飞机修正控制律设计

随着一些新型航空材料的应用，各式各样的变体飞机得以出现，如变后掠翼、折叠翼和扑翼等无人机。快速变体是新一代变体飞机应当具有的能力。例如变体战斗机在巡航的过程中，发现空中敌对目标后应当迅速地从巡航构型变化到机动构型，以更快的速度对目标实时跟踪打击。变体技术并不能改变飞行品质，反而会给飞机的动态特性带来不利的影响[1]。变体过程中，飞机的动态特性会出现较大变化，为了满足稳定飞行的要求，所设计的飞行控制系统需具有较高的鲁棒性，不受变体过程的影响。 Abdulrahim等基于H∞方法对一种仿生变体无人机进行控制器设计，仿真结果显示飞行稳定性会随着变体速率的提高而降低[2]。Hurst分别采用了多级补偿器和线性二次型调节器对一种仿生变体无人机进行着陆控制器设计，仿真结果表明基于线性二次型调节器的控制系统的轨迹跟踪误差较小，抗干扰能力较强[3]。Baldelli采用线性参数时变 (Linear parameter-varying，LPV) 方法对一种折叠翼无人机进行控制器设计，由于控制器增益较大，需要进行适当的降阶处理[4]。乐挺针对一种Z型翼无人机，提出一种多回路控制器，采用线性二次方法进行内回路控制器设计，采用增益自调度H∞方法进行外回路控制器设计，仿真结果表明所设计的控制器具有较强的鲁棒性，可以确保在机翼折叠过程中的飞行稳定性[5]。 上述控制方法虽然针对的控制对象不同，但都要求具有较高的建模精度。由于变体飞机在变体飞行过程中，很难对其动力学模型进行精确的数学描述，因此设计一种对建模精度要求不高而且能够确保变体过程中稳定飞行的控制系统是需要解决的一大难题。滚动时域优化 (Receding horizon optimal，RHO) 方法是一种在滚动时域内实时计算最优控制的预测控制方法，由于其对模型精度要求不高和具有强鲁棒性等特点被广泛地应用于各种领域[6]。Page基于串行修正体系结构[7-8]，采用RHO进行重构控制器设计，在不改变原有控制器结构的基础上，对原有控制器进行在线修正，保证了安全飞行[9]。鉴于RHO方法在重构控制中的成功应用，本文在文献[9]的基础上将RHO方法运用到变后掠角飞机控制律设计中，基于重构控制中的并行修正体系结构[10-11]提出了一种基于RHO的修正控制律设计方案，即：(1) 采用反步方法进行标称控制律设计，当变后掠角速率较小时，标称控制律能够提供基本的稳定性和跟踪性能；(2) 采用RHO方法进行并行修正控制律设计，在有限时域区间内实时解算修正控制量[12]，对标称控制律进行补偿，保证变后掠角过程中的飞行稳定性，使其基本不受建模误差和外界干扰的影响。 1 飞机变后掠角过程中的动态特性 飞机在变后掠角的过程中，质心位置会发生改变，而质心速度不变。研究飞机上某固定点的速度要比质心的速度更有意义。本文建立与机身固定的机体坐标轴系[13]，对图 1所示的小型变后掠角无人机进行研究。 图 1中给出了机翼后掠角为0°，30 °和45°时的构型。机翼后掠角可以在0°~45°之间任意变化，与MFX-2变后掠角飞机的机翼变形机构类似[14]。该变后掠角飞机在机翼后掠角变化过程中，机翼宽度会按照一定的比例随着机翼后掠角一起变化，即机翼后掠角确定，机翼气动弦长和机翼面积确定，可认为变后掠角自由度为1。在机翼后掠角最大时，机翼宽度最小。为了方便叙述，在后文中用机翼后掠角来表示飞机构型，如机翼后掠角为0°时，简称0°构型。表 1给出了变后掠角飞机在0°构型、30°构型和45°构型时的外形参数。飞机质量m=25 kg。 表 1 变后掠角飞机不同构型参数 Table 1 Different configuration parameters of variable sweep wing aircraft 表 1中:Sw为机翼参考面积；cA为机翼的平均几何弦长；b为机翼展长。变后掠角飞机从0°构型向45°构型变化过程中，飞机的机翼参考面积和展长逐渐减小，机翼的参考面积最大可改变36%，机翼展长最大可改变25.7%，机翼平均几何弦长最大可改变15%。 设左机翼、右机翼和机尾的质量分别为m1，m2和m3，其在机体坐标轴系中的坐标分别为r1=r1xi+r1yj+0k，r2=r2xi +r2y j+0k和r3=r3xi+0j+0k。飞机在机翼后掠角变化过程中保持面对称，则有m1=m2，r1x=r2x，r1y=-r2y，飞机的静矩$ \mathit{\boldsymbol{S}} = \sum\limits_{k = 1}^3 {{m_k}{\mathit{\boldsymbol{r}}_k}} = (2{m_1}{r_{1x}} + {m_3}{r_{3x}})\mathit{\boldsymbol{i}} + 0\mathit{\boldsymbol{j}} + 0\mathit{\boldsymbol{k}} $，沿着机体x轴的分量Sx=2m1r1x+m3r3x。变后掠角飞机的纵向运动方程为 $ \left\{ \begin{array}{l} \dot V = \frac{1}{m}\left( { - D + T\cos \alpha - mg\sin \gamma + {F_{Ix}}} \right)\\ \dot \gamma = \frac{1}{{mV}}\left( {L + T\sin \alpha - mg\cos \gamma + {F_{Iz}}} \right)\\ \dot \alpha = \frac{1}{{mV}}\left( { - L - T\sin \alpha + mg\cos \gamma + {F_{Iz}}} \right) + q\\ \dot q = - \frac{{{{\dot I}_y}}}{{{I_y}}}q + \frac{1}{{{I_y}}}\left( { - {S_x}g\cos \theta + {M_A} + T{Z_T} + {M_{Iy}}} \right)\\ \dot \theta = q \end{array} \right. $ (1) 式中：V为飞行速度；γ为航迹倾斜角；α为迎角；q为俯仰角速度；θ为俯仰角；m为飞机质量；g为重力加速度；Iy为飞机绕机体y轴的转动惯量；T为推力，与机体x轴方向平行；ZT为动力位置；FIx，FIz和MIy分别为变后掠角过程引起的惯性力和惯性力矩；L为升力；D为阻力；MA为空气动力产生的俯仰力矩。具体表达式为 $ \left\{ \begin{array}{l} {F_{Ix}} = {S_x}\left( {\dot q\sin \alpha + {q^2}\cos \alpha } \right) + 2{{\dot S}_x}q\sin \alpha - {{\ddot S}_x}\cos \alpha \\ {F_{Iz}} = {S_x}\left( {\dot q\cos \alpha + {q^2}\sin \alpha } \right) + 2{{\dot S}_x}q\cos \alpha - {{\ddot S}_x}\sin \alpha \\ {M_{Iy}} = {S_x}\left( {\dot V\sin \alpha + \dot V\alpha \cos \alpha - Vq\cos \alpha } \right) \end{array} \right. $ (2) $ \left\{ \begin{array}{l} L = Q{S_w}\left( {{C_{L0}} + {C_{L\alpha }}\alpha + {C_{L\delta e}}{\delta _e}} \right)\\ D = Q{S_w}\left( {{C_{D0}} + {C_{D\alpha }}\alpha + {C_{D\alpha 2}}{\alpha ^2}} \right)\\ {M_A} = Q{S_w}{c_A}\left( {{C_{m0}} + {C_{m\alpha }}\alpha + {C_{m\delta e}}{\delta _e} + {C_{mq}}\frac{{q{c_A}}}{{2V}}} \right) \end{array} \right. $ (3) 式中：CDα和CDα2分别为阻力系数对α和α2的导数；CL0为零迎角升力系数；CLα和CLδe分别为升力系数对α和δe的导数；δe为升降舵面偏角；cA为机翼的平均几何弦长；Sw为机翼参考面积；Q=1/2 ρV2，ρ为空气密度。 2 标称反步控制律设计 考虑式 (1)，采用推力T对飞行速度V进行直接控制。将迎角α作为航迹倾斜角γ的虚拟控制量，在式 (1) 航迹倾斜角方程中存在$\frac{1}{{mV}}Q{S_w}{C_{L\alpha }}\alpha $，$ \frac{1}{{mV}}T{\rm{sin}}\alpha $，Sxq2sinα, $ -{\ddot S_x}{\rm{sin}}\alpha $多个带迎角α的项，但$ \frac{1}{{mV}}Q{S_w}{C_{L\alpha }}\alpha $对航迹倾斜角γ可起到主导控制作用，因此选取$ \frac{1}{{mV}}Q{S_w}{C_{L\alpha }}\alpha $作为控制项。同理选取q和$ \frac{1}{{{I_y}}}Q{S_w}{c_A}{C_{m\delta e}}{\delta _e} $作为迎角α和俯仰角速度q的控制项。考虑外界干扰，将式 (1) 写成如下仿射非线性系统形式 $ \left\{ \begin{array}{l} \dot V = {f_V}\left( \mathit{\boldsymbol{x}} \right) + {g_V}\left( \mathit{\boldsymbol{x}} \right)T + {d_V}\left( t \right)\\ \dot \gamma = {f_\gamma }\left( \mathit{\boldsymbol{x}} \right) + {g_\gamma }\left( \mathit{\boldsymbol{x}} \right)\alpha + {d_\gamma }\left( t \right)\\ \dot \alpha = {f_\alpha }\left( \mathit{\boldsymbol{x}} \right) + q + {d_\alpha }\left( t \right)\\ \dot q = {f_q}\left( \mathit{\boldsymbol{x}} \right) + {g_q}\left( \mathit{\boldsymbol{x}} \right){\delta _e} + {d_q}\left( t \right)\\ \dot \theta = q \end{array} \right. $ (4) 式中：dV(t)，dγ(t)，d α(t) 和dq(t) 为未知的外界干扰；fV(x)，fγ(x)，fα(x)，fq(x)，gV(x)，gγ(x) 和gq (x) 为系统参数，表达式为 $ \left\{ \begin{array}{l} {f_V}\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{m}\left( { - D - mg\sin \gamma + {F_{Ix}}} \right)\\ {f_\gamma }\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{{mV}}\left( {Q{S_w}\left( {{C_{L0}} + {C_{L\delta e}}{\delta _e}} \right) + T\sin \alpha - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {mg\cos \gamma - {F_{Iz}}} \right)\\ {f_\alpha }\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{{mV}}\left( { - L - T\sin \alpha + mg\cos \gamma + {F_{Iz}}} \right)\\ {f_q}\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{{{I_y}}}\left( {Q{S_w}{c_A}\left( {{C_{m0}} + {C_{m\alpha }}\alpha + {C_{mq}}\frac{{q{c_A}}}{{2V}}} \right) - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;{S_x}g\cos \theta - {{\dot I}_y}q + T{Z_T} + {M_{Iy}}\\ {g_V}\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{m}\cos \alpha \\ {g_\gamma }\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{{mV}}Q{S_w}{C_{L\alpha }}\\ {g_q}\left( \mathit{\boldsymbol{x}} \right) = \frac{1}{{{I_y}}}Q{S_w}{c_A}{C_{m\delta e}} \end{array} \right. $ (5) 式中x=[V γ α q θ]T为状态向量。由于在实际飞行过程中，气动参数存在一定的摄动，假设系统参数的精确值未知，可以写成 $ \left\{ \begin{array}{l} {f_V}\left( \mathit{\boldsymbol{x}} \right) = {f_{V0}}\left( \mathit{\boldsymbol{x}} \right) + \Delta {f_V}\left( \mathit{\boldsymbol{x}} \right)\\ {f_\gamma }\left( \mathit{\boldsymbol{x}} \right) = {f_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right) + \Delta {f_\gamma }\left( \mathit{\boldsymbol{x}} \right)\\ {f_\alpha }\left( \mathit{\boldsymbol{x}} \right) = {f_{\alpha 0}}\left( \mathit{\boldsymbol{x}} \right) + \Delta {f_\alpha }\left( \mathit{\boldsymbol{x}} \right)\\ {f_q}\left( \mathit{\boldsymbol{x}} \right) = {f_{q0}}\left( \mathit{\boldsymbol{x}} \right) + \Delta {f_q}\left( \mathit{\boldsymbol{x}} \right)\\ {g_\gamma }\left( \mathit{\boldsymbol{x}} \right) = {g_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right) + \Delta {g_\gamma }\left( \mathit{\boldsymbol{x}} \right)\\ {g_q}\left( \mathit{\boldsymbol{x}} \right) = {g_{q0}}\left( \mathit{\boldsymbol{x}} \right) + \Delta {g_q}\left( \mathit{\boldsymbol{x}} \right) \end{array} \right. $ (6) 式中：fV0(x)，fγ0(x)，fα0(x)，fq0(x)，gγ0(x) 和gq0(x) 为已知的标称系统参数；ΔfV(x)，Δfγ(x)，Δfα(x)，Δfq(x)，Δgγ(x) 和Δgq(x) 为未知的建模误差。则式 (4) 表示的仿射非线性系统为[16] $ \left\{ \begin{array}{l} \dot V = {f_{V0}}\left( \mathit{\boldsymbol{x}} \right) + {g_V}\left( \mathit{\boldsymbol{x}} \right)T + {\Delta _V}\left( \mathit{\boldsymbol{x}} \right)\\ \dot \gamma = {f_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right) + {g_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right)\alpha + {\Delta _\gamma }\left( \mathit{\boldsymbol{x}} \right)\\ \dot \alpha = {f_{\alpha 0}}\left( \mathit{\boldsymbol{x}} \right) + q + {\Delta _\alpha }\left( \mathit{\boldsymbol{x}} \right)\\ \dot q = {f_{q0}}\left( \mathit{\boldsymbol{x}} \right) + {g_{q0}}\left( \mathit{\boldsymbol{x}} \right){\delta _e} + {\Delta _q}\left( \mathit{\boldsymbol{x}} \right)\\ \dot \theta = q \end{array} \right. $ (7) 式中：ΔV(x)，Δγ(x)，Δα(x) 和Δq(x) 为模型不确定项，表达式为 $ \left\{ \begin{array}{l} {\Delta _V}\left( \mathit{\boldsymbol{x}} \right) = \Delta {f_V}\left( \mathit{\boldsymbol{x}} \right) + {d_V}\left( t \right)\\ {\Delta _\gamma }\left( \mathit{\boldsymbol{x}} \right) = \Delta {f_\gamma }\left( \mathit{\boldsymbol{x}} \right) + \Delta {g_\gamma }\left( \mathit{\boldsymbol{x}} \right)\alpha + {d_\gamma }\left( t \right)\\ {\Delta _\alpha }\left( \mathit{\boldsymbol{x}} \right) = \Delta {f_\alpha }\left( \mathit{\boldsymbol{x}} \right) + {d_\alpha }\left( t \right)\\ {\Delta _q}\left( \mathit{\boldsymbol{x}} \right) = \Delta {f_q}\left( \mathit{\boldsymbol{x}} \right) + {g_q}\left( \mathit{\boldsymbol{x}} \right){\delta _e} + {d_q}\left( t \right) \end{array} \right. $ (8) 基于已知的标称系统参数进行标称反步控制器设计，保证基本的飞行稳定性和目标跟踪性能。设计相应虚拟控制信号、升降舵控制信号和推力控制信号为 $ \left\{ \begin{array}{l} {\alpha _d} = \frac{1}{{{g_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right)}}\left( { - {k_\gamma }\tilde \gamma - {f_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right) + {{\dot \gamma }_d}} \right)\\ {q_d} = - {k_\alpha }\tilde \alpha - {g_{\gamma 0}}\left( \mathit{\boldsymbol{x}} \right)\tilde \gamma - {f_{\alpha 0}}\left( \mathit{\boldsymbol{x}} \right) + {{\dot \alpha }_d}\\ {\delta _e} = \frac{1}{{{g_{q0}}\left( \mathit{\boldsymbol{x}} \right)}}\left( { - {k_q}\tilde q - \tilde a - {f_{q0}}\left( \mathit{\boldsymbol{x}} \right) + {{\dot q}_d}} \right)\\ T = \frac{1}{{{g_V}\left( \mathit{\boldsymbol{x}} \right)}}\left( { - {k_V}\tilde V - {f_{V0}}\left( \mathit{\boldsymbol{x}} \right) + {{\dot V}_d}} \right) \end{array} \right. $ (9) 式中：kγ>0，kα>0，kq>0，kV>0为设计参数；γd和Vd分别为航迹倾斜角和飞行速度指令；$\tilde \gamma $和$\tilde V $分别为航迹倾斜角和飞行速度指令跟踪误差；αd和qd为虚拟控制指令；$ \tilde q $和$ \tilde \alpha $为虚拟控制指令跟踪误差。采用以下二阶指令滤波器[17]得到Vd，γd及一阶微分信号$ {\dot V_d}, {\dot \gamma _d} $。 $ \left\{ \begin{array}{l} \frac{{{V_d}}}{{{V_{d0}}}} = \frac{{\omega _V^2}}{{{s^2} + 2{\xi _V}{\omega _V}s + \omega _V^2}}\\ \frac{{{\gamma _d}}}{{{\gamma _{d0}}}} = \frac{{\omega _\gamma ^2}}{{{s^2} + 2{\xi _\gamma }{\omega _\gamma }s + \omega _\gamma ^2}} \end{array} \right. $ (10) 式中：Vd0，γd0为指令滤波器输入；ωV，ωγ为带宽；ξV，ξγ为阻尼比。 式 (9) 的具体推导过程和控制器的稳定性证明可参照文献[16]，这里不再赘述。在变体过程中，飞机的动态特性处于相当不稳定的过渡状态，标称反步控制不能满足所需飞行品质要求。下面进行修正控制律设计，保证飞机在变后掠角过程中能够稳定飞行。 3 运动模型线性化 采用小扰动线性化方法对运动方程 (1) 进行线性化可得 $ \mathit{\boldsymbol{E}}\Delta \mathit{\boldsymbol{\dot x}} = \mathit{\boldsymbol{\bar A}}\Delta \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{\bar B}}\Delta \mathit{\boldsymbol{u}} $ (11) 式中：Δx=[ΔV Δγ Δα Δq Δθ] T为全状态向量x=[V γ α q θ]T相对于基准值x0的增量，x0=[V0 0 α0 0 θ0]T，V0为基准飞行速度，α0和θ0分别为配平迎角和俯仰角，α0=θ0；Δu=[Δδe ΔT]T为控制量u=[δe T]T相对于配平控制量u 0=[δe0 T0]T的改变量，Δδe为升降舵面偏角相对配平偏角δe0的改变量，ΔT为推力相对配平推力T0的改变量；E，A和B的表达式为 $ \begin{array}{l} \mathit{\boldsymbol{E = }}\\ \left[ {\begin{array}{*{20}{c}} 1&0&0&{ - \frac{1}{m}{S_x}\sin {\alpha _0}}&0\\ 0&1&0&{\frac{1}{{m{V_0}}}{S_x}\cos {\alpha _0}}&0\\ 0&0&1&{ - \frac{1}{{m{V_0}}}{S_x}\cos {\alpha _0}}&0\\ { - \frac{1}{{{I_y}}}{S_x}\sin {\alpha _0}}&0&{ - \frac{1}{{{I_y}}}{S_x}{V_0}\cos {\alpha _0}}&1&0\\ 0&0&0&0&1 \end{array}} \right] \end{array} $ (12) $ \mathit{\boldsymbol{\bar A}} = \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}&{{a_{13}}}&{{a_{14}}}&{{a_{15}}}\\ {{a_{21}}}&{{a_{22}}}&{{a_{23}}}&{{a_{24}}}&{{a_{25}}}\\ {{a_{31}}}&{{a_{32}}}&{{a_{33}}}&{{a_{34}}}&{{a_{35}}}\\ {{a_{41}}}&{{a_{42}}}&{{a_{43}}}&{{a_{44}}}&{{a_{45}}}\\ {{a_{51}}}&{{a_{52}}}&{{a_{53}}}&{{a_{54}}}&{{a_{55}}} \end{array}} \right] $ (13) $ \begin{array}{l} \mathit{\boldsymbol{\bar B}} = \\ {\left[ {\begin{array}{*{20}{c}} 0&{\frac{{{Q_0}{S_w}{C_{L\delta e}}}}{{m{V_0}}}}&{ - \frac{{{Q_0}{S_w}{C_{L\delta e}}}}{{m{V_0}}}}&{\frac{{{Q_0}{S_w}{c_A}{C_{m\delta e}}}}{{{I_y}}}}&0\\ {\frac{{\cos {\alpha _0}}}{m}}&{\frac{{\sin {\alpha _0}}}{{m{V_0}}}}&{ - \frac{{\sin {\alpha _0}}}{{m{V_0}}}}&{\frac{{{Z_T}}}{{{I_y}}}}&0 \end{array}} \right]^{\rm{T}}} \end{array} $ (14) 式中：$ {Q_0} = \frac{1}{2}\rho V_0^2 $；A中元素具体表达式为 $ \begin{array}{l} \left\{ \begin{array}{l} {a_{11}} = \frac{{ - \rho {V_0}{S_w}\left( {{C_{D0}} + {C_{{D_\alpha }{\alpha _0}}} + {C_{{D_\alpha }2}}\alpha _0^2} \right)}}{m},{a_{12}} = - g\\ {a_{13}} = \frac{{ - {Q_0}{S_w}\left( {{C_{{D_\alpha }}} + 2{C_{{D_\alpha }2}}{\alpha _0}} \right) - {T_0}\sin {\alpha _0} + {{\ddot S}_x}\sin {\alpha _0}}}{m}\\ {a_{14}} = \frac{{2{{\dot S}_x}\sin {\alpha _0}}}{m},{a_{15}} = 0 \end{array} \right.\\ \left\{ \begin{array}{l} {a_{21}} = \frac{{\rho \left( {{C_{L0}} + {C_{L\alpha }}{\alpha _0} + {C_{L\delta e}}{\delta _{e0}}} \right){S_w}}}{{2m}} - \\ \;\;\;\;\;\;\;\;\frac{{{T_0}\sin {\alpha _0} - mg - {{\ddot S}_x}\sin {\alpha _0}}}{{mV_0^2}}\\ {a_{22}} = 0,{a_{23}} = \frac{{{Q_0}{S_w}{C_{L\alpha }} + {T_0}\cos {\alpha _0} - {{\ddot S}_x}\cos {\alpha _0}}}{{m{V_0}}}\\ {a_{24}} = \frac{{ - 2{{\dot S}_x}\cos {\alpha _0}}}{{m{V_0}}},{a_{25}} = 0 \end{array} \right. \end{array} $ $ \begin{array}{l} \left\{ \begin{array}{l} {a_{31}} = - \frac{{\rho \left( {{C_{L0}} + {C_{L\alpha }}{\alpha _0} + {C_{L\delta e}}{\delta _{e0}}} \right){S_w}}}{{2m}} + \\ \;\;\;\;\;\;\;\;\frac{{{T_0}\sin {\alpha _0} - mg - {{\ddot S}_x}\sin {\alpha _0}}}{{mV_0^2}}\\ {a_{32}} = 0,{a_{32}} = - \frac{{{Q_0}{S_w}{C_{L\alpha }} + {T_0}\cos {\alpha _0} - {{\ddot S}_x}\cos {\alpha _0}}}{{m{V_0}}}\\ {a_{34}} = \frac{{2{{\dot S}_x}\cos {\alpha _0}}}{{m{V_0}}} + 1,{a_{25}} = 0 \end{array} \right.\\ \left\{ \begin{array}{l} {a_{41}} = \frac{{\rho {V_0}{S_w}{c_A}\left( {{C_{m0}} + {C_{m\alpha }} + {C_{m\delta e}}{\delta _{e0}}} \right)}}{{{I_y}}},{a_{42}} = 0\\ {a_{43}} = \frac{{{Q_0}{S_w}{c_A}{C_{m\alpha }}}}{{{I_y}}}\\ {a_{44}} = - \frac{{{{\dot I}_y} + {S_x}{V_0}\cos {\alpha _0}}}{{{I_y}}} + \frac{{\rho {V_0}{S_w}c_A^2{C_{mq}}}}{{4{I_y}}}\\ {a_{45}} = \frac{{{S_x}g\sin {\theta _0}}}{{{I_y}}} \end{array} \right.\\ {a_{51}} = 0,{a_{52}} = 0,{a_{53}} = 0,{a_{54}} = 1,{a_{52}} = 0 \end{array} $ 飞行速度和航迹倾斜角为系统输出，考虑建模误差和外界干扰，将式 (11) 写成标准线性状态方程的形式，即 $ \left\{ \begin{array}{l} \Delta \mathit{\boldsymbol{\dot x = A}}\Delta \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{B}}\Delta \mathit{\boldsymbol{u + }}{\mathit{\boldsymbol{b}}_0}\\ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{C}}\Delta \mathit{\boldsymbol{x}} \end{array} \right. $ (15) 式中：A=E-1A；B=E-1B；b0为模型不确定项；y=[ΔV Δγ]T为系统输出；输出矩阵C为 $ \mathit{\boldsymbol{C = }}\left[ {\begin{array}{*{20}{c}} 1&0&0&0&0\\ 0&1&0&0&0 \end{array}} \right] $ (16) 4 基于RHO的修正控制律设计 4.1 增广状态方程设计 基于重构控制中的并行修正体系结构进行RHO修正控制律设计，并行修正体系结构如图 2所示[10-11]。图中：un为标称控制律输出；ur为修正控制律输出；u=un+ur为控制量。修正控制律根据飞机的输出、指令信号和标称控制律输出实时计算修正控制量ur。 式 (15) 中的Δu可表达为 $ \Delta \mathit{\boldsymbol{u}}{\rm{ = }}\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}_0} = {\mathit{\boldsymbol{u}}_r} + {\mathit{\boldsymbol{u}}_n} - {\mathit{\boldsymbol{u}}_0} $ (17) 将式 (17) 代入式 (15) 得 $ \left\{ \begin{array}{l} \Delta \mathit{\boldsymbol{\dot x = A}}\Delta \mathit{\boldsymbol{x}} + \mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{u}}_r} + {\mathit{\boldsymbol{u}}_n} - {\mathit{\boldsymbol{u}}_0}} \right)\mathit{\boldsymbol{ + }}{\mathit{\boldsymbol{b}}_0}\\ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{C}}\Delta \mathit{\boldsymbol{x}} \end{array} \right. $ (18) 将二阶指令滤波器 (10) 写成如下状态空间的形式，即 $ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_r} = {\mathit{\boldsymbol{A}}_r}{\mathit{\boldsymbol{x}}_r} + {\mathit{\boldsymbol{B}}_r}{\mathit{\boldsymbol{r}}_{d0}}\\ {\mathit{\boldsymbol{y}}_r} = {\mathit{\boldsymbol{C}}_r}{\mathit{\boldsymbol{x}}_r} \end{array} \right. $ (19) 式中：${{\boldsymbol{x}}_r} = {[\Delta {V_d}\;\;\Delta {\dot V_d}\;\;\;\Delta {\gamma _d}\;\;\;\Delta {\dot \gamma _d}]^{\rm{T}}} $等于飞行速度指令信号减去基准飞行速度，即ΔVd=Vd-V0，同样Δγd=γd-γ0，由于γ0=0，有Δγd=γd；$ {\mathit{\boldsymbol{r}}_{d0}} = \left[{\begin{array}{*{20}{c}} {{V_{d0}}-{V_0}}\\ {{\gamma _{d0}}} \end{array}} \right] $；Ar，Br，Cr的表达式为 $ \begin{array}{l} {\mathit{\boldsymbol{A}}_r} = \left[ {\begin{array}{*{20}{c}} { - 2{\xi _V}{\omega _V}}&1&0&0\\ { - \omega _V^2}&0&0&0\\ 0&0&{ - 2{\xi _\gamma }{\omega _\gamma }}&1\\ 0&0&{ - \omega _\gamma ^2}&0 \end{array}} \right]\\ {\mathit{\boldsymbol{B}}_r} = \left[ {\begin{array}{*{20}{c}} 0&0\\ {\omega _V^2}&0\\ 0&0\\ 0&{\omega _\gamma ^2} \end{array}} \right],{\mathit{\boldsymbol{C}}_r} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&0&1&0 \end{array}} \right] \end{array} $ 基于式 (18，19) 在有限时域区间内实时计算修正控制器增益，使得最小性能指标为 $ \begin{array}{l} \mathop {\min }\limits_{{\delta _{er}}} {J_R} = \frac{1}{2}\int_{{t_0}}^{{t_f}} {\left[ {{{\left( {{\mathit{\boldsymbol{y}}_r} - \mathit{\boldsymbol{y}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_p}\left( {{\mathit{\boldsymbol{y}}_r} - \mathit{\boldsymbol{y}}} \right) + \mathit{\boldsymbol{y}}_I^{\rm{T}}{\mathit{\boldsymbol{Q}}_I}{\mathit{\boldsymbol{y}}_I} + } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\mathit{\boldsymbol{u}}_r^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{r}}}} \right]{\rm{d}}t \end{array} $ (20) 式中：tf和t0分别为有限时域区间的上下界；Qp和QI分别为跟踪误差和跟踪误差积分的加权阵；R为修正控制量的加权阵；yI为积分误差，可由式 (21) 得到 $ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_I} = {\mathit{\boldsymbol{y}}_r} - \mathit{\boldsymbol{y}}\\ {\mathit{\boldsymbol{y}}_I} = {\mathit{\boldsymbol{x}}_I} \end{array} \right. $ (21) 综合飞机线性状态方程、指令滤波器状态方程和积分误差方程得到增广状态方程为 $ \begin{array}{*{20}{c}} {\left[ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_r}\\ \Delta \mathit{\boldsymbol{\dot x}}\\ {{\mathit{\boldsymbol{\dot x}}}_I} \end{array} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_r}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{A}}&\mathit{\boldsymbol{0}}\\ {{\mathit{\boldsymbol{C}}_r}}&{ - \mathit{\boldsymbol{C}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} {\mathit{\boldsymbol{x}}_r}\\ \Delta \mathit{\boldsymbol{x}}\\ {\mathit{\boldsymbol{x}}_I} \end{array} \right] + \left[ \begin{array}{l} \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{B}}\\ \mathit{\boldsymbol{0}} \end{array} \right]{\mathit{\boldsymbol{u}}_r} + }\\ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}_r}{\mathit{\boldsymbol{r}}_{d0}}}\\ {\mathit{\boldsymbol{B}}\left( {{\mathit{\boldsymbol{u}}_n} - {\mathit{\boldsymbol{u}}_0}} \right)}\\ \mathit{\boldsymbol{0}} \end{array}} \right] + \left[ \begin{array}{l} \mathit{\boldsymbol{0}}\\ {\mathit{\boldsymbol{b}}_0}\\ \mathit{\boldsymbol{0}} \end{array} \right]} \end{array} $ (22) 将式 (22) 写成对应的如下形式 $ {{\mathit{\boldsymbol{\dot x}}}_R} = {\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{u}}_r} + {\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}} $ (23) 式中：$ {{\boldsymbol{A}}_R} = \left[{\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_r}}&{\boldsymbol{0}}&{\boldsymbol{0}}\\ {\boldsymbol{0}}&{\boldsymbol{A}}&{\boldsymbol{0}}\\ {{{\boldsymbol{C}}_r}}&{\boldsymbol{-C}}&{\boldsymbol{0}} \end{array}} \right];{{\boldsymbol{B}}_R} = \left[{\begin{array}{*{20}{c}} {\boldsymbol{0}}\\ {\boldsymbol{B}}\\ {\boldsymbol{0}} \end{array}} \right];{\boldsymbol{b}_R} = \left[{\begin{array}{*{20}{c}} {{{\boldsymbol{B}}_r}{\boldsymbol{r}_{d0}}}\\ {{\boldsymbol{B}}({\rm{ }}{\boldsymbol{u}_n}-{\boldsymbol{u}_0})}\\ {\boldsymbol{0}} \end{array}} \right];{\boldsymbol{b}} = \left[{\begin{array}{*{20}{c}} {\boldsymbol{0}}\\ {{\boldsymbol{b}_0}}\\ {\boldsymbol{0}} \end{array}} \right] $。 4.2 修正控制量设计 式 (20) 可以重新写成 $ \mathop {\min }\limits_{{\delta _R}} {J_R} = \frac{1}{2}\int_{{t_0}}^{{t_f}} {\left( {\mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{Q}}_R}{\mathit{\boldsymbol{x}}_R} + \mathit{\boldsymbol{u}}_r^T\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{u}}_r}} \right){\rm{d}}t} $ (24) 式中 $ {\mathit{\boldsymbol{Q}}_R} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_r^{\rm{T}}{\mathit{\boldsymbol{Q}}_p}{\mathit{\boldsymbol{C}}_r}}&{ - \mathit{\boldsymbol{C}}_r^{\rm{T}}{\mathit{\boldsymbol{Q}}_p}\mathit{\boldsymbol{C}}}&\mathit{\boldsymbol{0}}\\ { - {\mathit{\boldsymbol{C}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_p}{\mathit{\boldsymbol{C}}_r}}&{{\mathit{\boldsymbol{C}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_p}{\mathit{\boldsymbol{C}}_r}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_I}} \end{array}} \right] $ (25) 对修正控制量进行求解，先构造如下Hamilton函数[18] $ \begin{array}{l} \mathit{\boldsymbol{H}}\left( {{\mathit{\boldsymbol{x}}_R},{\mathit{\boldsymbol{u}}_r},{\mathit{\boldsymbol{\lambda }}_R},t} \right) = \frac{1}{2}\mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{Q}}_R}{\mathit{\boldsymbol{x}}_R} + \frac{1}{2}\mathit{\boldsymbol{u}}_r^T\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{u}}_r} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {{\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{u}}_r} + {\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{\lambda }}_R} \end{array} $ (26) 式中λR为时变拉格朗日乘子向量。 ur使得Hamilton函数值最小的必要条件为∂H/∂ur=0，有 $ \frac{{\partial \mathit{\boldsymbol{H}}}}{{\partial {\mathit{\boldsymbol{u}}_r}}}\left( {{\mathit{\boldsymbol{x}}_R},{\mathit{\boldsymbol{u}}_r},{\mathit{\boldsymbol{\lambda }}_R},t} \right) = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{u}}_r} + \mathit{\boldsymbol{B}}_R^{\rm{T}}{\mathit{\boldsymbol{\lambda }}_R} = 0 $ (27) 由于 $ \frac{{{\partial ^2}\mathit{\boldsymbol{H}}}}{{\partial \mathit{\boldsymbol{u}}_r^2}}\left( {{\mathit{\boldsymbol{x}}_R},{\mathit{\boldsymbol{u}}_r},{\mathit{\boldsymbol{\lambda }}_R},t} \right) = \mathit{\boldsymbol{R}} $ (28) 式中：R为正定矩阵；H为ur的二次型，则满足式 (27) 的u r可使H值最小。求式 (27) 的解为 $ {\mathit{\boldsymbol{u}}_r} = - {\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{B}}_R^{\rm{T}}{\mathit{\boldsymbol{\lambda }}_R} $ (29) 将时变拉格朗日乘子向量λR设计成状态xR的线性组合，有 $ {\mathit{\boldsymbol{\lambda }}_R} = {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{k}}_R} $ (30) 式中KR和kR为参数变量。 对于时变拉格朗日乘子向量λR有如下等式成立[19] $ \begin{array}{l} - {{\mathit{\boldsymbol{\dot \lambda }}}_R} = \frac{\partial }{{\partial {\mathit{\boldsymbol{x}}_R}}}\left( {\frac{1}{2}\mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{Q}}_R}{\mathit{\boldsymbol{x}}_R} + \frac{1}{2}\mathit{\boldsymbol{u}}_r^T\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{u}}_r}} \right) + \\ \;\;\;\;\;\;\;\;\;\frac{\partial }{{\partial {\mathit{x}_R}}}{\left( {{\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{B}}_\mathit{\boldsymbol{R}}}{\mathit{\boldsymbol{u}}_r} + {\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{\lambda }}_R} \end{array} $ (31) 整理式 (31) 得 $ {{\mathit{\boldsymbol{\dot \lambda }}}_R} = - {\mathit{\boldsymbol{Q}}_R}{\mathit{\boldsymbol{x}}_R} - \mathit{\boldsymbol{A}}_R^{\rm{T}}{\mathit{\boldsymbol{\lambda }}_R} $ (32) 将式 (30) 两边对时间t求导得 $ \begin{array}{l} {{\mathit{\boldsymbol{\dot \lambda }}}_R} = {{\mathit{\boldsymbol{\dot K}}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{K}}_R}{{\mathit{\boldsymbol{\dot x}}}_R} + {{\mathit{\boldsymbol{\dot k}}}_R} = {{\mathit{\boldsymbol{\dot K}}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{K}}_R}\left( {{\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} + } \right.\\ \;\;\;\;\;\;\;\left. {{\mathit{\boldsymbol{B}}_\mathit{\boldsymbol{R}}}{\mathit{\boldsymbol{u}}_r} + {\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right) + {{\mathit{\boldsymbol{\dot k}}}_R} \end{array} $ (33) 将式 (29) 代入式 (33) 得 $ \begin{array}{l} {{\mathit{\boldsymbol{\dot \lambda }}}_R} = {{\mathit{\boldsymbol{\dot K}}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{K}}_R}\left( {{\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} - {\mathit{\boldsymbol{B}}_\mathit{\boldsymbol{R}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{B}}_R^{\rm{T}}{\mathit{\boldsymbol{\lambda }}_R} + {\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right) + \\ \;\;\;\;\;\;\;{{\mathit{\boldsymbol{\dot k}}}_R} = {{\mathit{\boldsymbol{\dot K}}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} - {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{x}}_R} - \\ \;\;\;\;\;\;\;{\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} + {\mathit{\boldsymbol{K}}_R}\left( {{\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right) + {{\mathit{\boldsymbol{\dot k}}}_R} \end{array} $ (34) 根据式 (32，34) 得 $ \begin{array}{l} - {\mathit{\boldsymbol{Q}}_R}{\mathit{\boldsymbol{x}}_R} - \mathit{\boldsymbol{A}}_R^{\rm{T}}\left( {{\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{k}}_R}} \right) = {{\mathit{\boldsymbol{\dot K}}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{A}}_R}{\mathit{\boldsymbol{x}}_R} - \\ {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{x}}_R} - {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} + {\mathit{\boldsymbol{K}}_R}\left( {{\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right) + \\ {{\mathit{\boldsymbol{\dot k}}}_R} \end{array} $ (35) 整理式 (35) 得 $ \begin{array}{l} \mathit{\boldsymbol{0 = }}\left( {{{\mathit{\boldsymbol{\dot K}}}_R} + {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{A}}_R} + \mathit{\boldsymbol{A}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R} - {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R} + } \right.\\ \left. {{\mathit{\boldsymbol{Q}}_R}} \right){\mathit{\boldsymbol{x}}_R} + {{\mathit{\boldsymbol{\dot k}}}_R} + \mathit{\boldsymbol{A}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R} - {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} + \\ {\mathit{\boldsymbol{K}}_R}\left( {{\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right) \end{array} $ (36) 由式 (36) 对于任意的状态向量xR都成立，则可以得到如下两个黎卡提微分方程 $ {{\mathit{\boldsymbol{\dot K}}}_R} = - {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{A}}_R} - \mathit{\boldsymbol{A}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R} + {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R} - {\mathit{\boldsymbol{Q}}_R} $ (37) $ {{\mathit{\boldsymbol{\dot k}}}_R} = - {\mathit{\boldsymbol{K}}_R}\left( {{\mathit{\boldsymbol{b}}_R} + \mathit{\boldsymbol{b}}} \right) - \mathit{\boldsymbol{A}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} + {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} $ (38) 通过求解式 (37，38) 所示的黎卡提微分方程得到参数变量KR和kR。由于式 (38) 中模型不确定项b未知，且滚动时域优化方法对模型精度要求不高，可用如下黎卡提微分方程近似计算kR[12] $ {{\mathit{\boldsymbol{\dot k}}}_R} = - {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{b}}_R} - \mathit{\boldsymbol{A}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} + {\mathit{\boldsymbol{K}}_R}{\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} $ (39) 由于KR，kR在有限时域区间[t0, tf]边界t f可以稳定在任意一个常值上，认为KR和kR在时域边界tf为0。在有限时域内对式 (37，39) 进行反向积分，得到KR(t0) 和kR(t0)。通过式 (40) 计算得到修正控制律增益 $ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{F}}_R}}&{{\mathit{\boldsymbol{f}}_R}} \end{array}} \right] = - {\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{K}}_R}\left( {{t_0}} \right)}&{{\mathit{\boldsymbol{k}}_R}\left( {{t_0}} \right)} \end{array}} \right] $ (40) 式中：FR=[Fr Fx FI]为修正控制律增益；Fr，Fx和FI分别为对应指令滤波器状态、飞机状态和积分误差状态的修正控制律增益。 通过式 (41) 计算得到修正控制量 $ {\mathit{\boldsymbol{u}}_r} = {\mathit{\boldsymbol{F}}_R}{\mathit{\boldsymbol{x}}_R} + {\mathit{\boldsymbol{f}}_R} $ (41) 4.3 稳定性分析 将式 (41) 代入式 (23) 得 $ {{\mathit{\boldsymbol{\dot x}}}_R} = \left( {{\mathit{\boldsymbol{A}}_R} - {\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{K}}_R}} \right){\mathit{\boldsymbol{x}}_R} + \mathit{\boldsymbol{b}} + {\mathit{\boldsymbol{b}}_R} - {\mathit{\boldsymbol{B}}_R}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{R}}_R^{\rm{T}}{\mathit{\boldsymbol{k}}_R} $ (42) 令Ac=AR-BRR-1BRTKR，Ac为严格Hurwitz矩阵。对于任意的正定矩阵Qc存在正定对称矩阵Pc使式 (43) 成立，则 $ \mathit{\boldsymbol{A}}_c^{\rm{T}}{\mathit{\boldsymbol{P}}_c} + {\mathit{\boldsymbol{P}}_c}{\mathit{\boldsymbol{A}}_c} = - {\mathit{\boldsymbol{Q}}_c} $ (43) 令bc=b+bR-BRR-1BRTk R，由黎卡提微分方程 (39) 的收敛性可得 $ \left\| {{\mathit{\boldsymbol{b}}_c}} \right\| \le {\delta _b} $ (44) 式中δb>0为已知常数。设计李雅普诺夫函数VR为 $ {V_R} = \mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{P}}_c}{\mathit{\boldsymbol{x}}_R} $ (45) VR对时间t求导得 $ \begin{array}{l} {{\dot V}_R} = \mathit{\boldsymbol{\dot x}}_R^{\rm{T}}{\mathit{\boldsymbol{P}}_c}{\mathit{\boldsymbol{x}}_R} + \mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{P}}_c}{{\mathit{\boldsymbol{\dot x}}}_R} = \mathit{\boldsymbol{x}}_R^{\rm{T}}\left( {\mathit{\boldsymbol{A}}_c^{\rm{T}}{\mathit{\boldsymbol{P}}_c} + {\mathit{\boldsymbol{P}}_c}{\mathit{\boldsymbol{A}}_c}} \right){\mathit{\boldsymbol{x}}_R} + \\ \;\;\;\;\;\;\;\mathit{\boldsymbol{b}}_c^{\rm{T}}{\mathit{\boldsymbol{P}}_c}{\mathit{\boldsymbol{x}}_R} + \mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{P}}_c}{\mathit{\boldsymbol{b}}_c} \le - {\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_c}} \right){\left\| {{\mathit{\boldsymbol{x}}_R}} \right\|^2} + \\ \;\;\;\;\;\;\;{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_c}} \right)\left( {\mathit{\boldsymbol{b}}_c^{\rm{T}}{\mathit{\boldsymbol{x}}_R} + \mathit{\boldsymbol{x}}_R^{\rm{T}}{\mathit{\boldsymbol{b}}_c}} \right) \le - {\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_c}} \right){\left\| {{\mathit{\boldsymbol{x}}_R}} \right\|^2} + \\ \;\;\;\;\;\;\;{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_c}} \right)\left( {\mathit{\boldsymbol{v}}{{\left\| {{\mathit{\boldsymbol{x}}_R}} \right\|}^2} + {\mathit{\boldsymbol{v}}^{ - 1}}\delta _b^2} \right) \le \\ \;\;\;\;\;\;\; - \frac{{\left( {{\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_c}} \right) - \mathit{\boldsymbol{v}}{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_c}} \right)} \right)}}{{{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_c}} \right)}}{V_c} + {\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_c}} \right){\mathit{\boldsymbol{v}}^{ - 1}}\delta _b^2 \end{array} $ (46) 式中：ν>0；λmin(Qc) 为Qc的最小特征值；λmax(Pc) 为Pc的最大特征值。通过调节有限时域区间[t0, tf]的长度，可以使xR较快地收敛到式 (47) 所示范围内，同时具有较好的实时性。 $ \mathop {\lim }\limits_{t \to \infty } \left\| {{\mathit{\boldsymbol{x}}_R}} \right\| \le \sqrt {\frac{{\lambda _{\max }^2\left( {{\mathit{\boldsymbol{P}}_c}} \right){\mathit{\boldsymbol{v}}^{ - 1}}\delta _b^2}}{{{\lambda _{\min }}\left( {{\mathit{\boldsymbol{P}}_c}} \right)\left[ {{\lambda _{\min }}\left( {{\mathit{\boldsymbol{Q}}_c}} \right) - \mathit{\boldsymbol{v}}{\lambda _{\max }}\left( {{\mathit{\boldsymbol{P}}_c}} \right)} \right]}}} $ (47) 式中λmin(Pc) 为Pc的最小特征值。 5 仿真分析 为了验证控制系统的鲁棒性，在标称气动参数上加上标称气动参数的30%作为建模误差，并施加干扰信号dq(t)=6((°)/s2)。指令滤波器状态空间的形式为 $ {\mathit{\boldsymbol{A}}_r} = \left[ {\begin{array}{*{20}{c}} { - 2}&1&0&0\\ { - 1}&0&0&0\\ 0&0&{ - 2}&1\\ 0&0&{ - 1}&0 \end{array}} \right] $ (48) $ {\mathit{\boldsymbol{B}}_r} = \left[ {\begin{array}{*{20}{c}} 0&0\\ 1&0\\ 0&0\\ 0&1 \end{array}} \right] $ (49) $ {\mathit{\boldsymbol{C}}_r} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&0&1&0 \end{array}} \right] $ (50) 滚动时域优化修正控制律设计参数为$ {\boldsymbol{Q}_p} = \left[ {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right] $, ${{\boldsymbol{Q}}_{1}}=\left[ \begin{matrix} 20\ 000 & 0 \\ 0 & 2\ 000 \\ \end{matrix} \right]$,$\boldsymbol{R} = \left[ {\begin{array}{*{20}{c}} {0.0001}&0\\ 0&{0.2} \end{array}} \right] $。 仿真初始条件为：飞行速度V=30 m/s、飞行高度h=1 000 m，变后掠角飞机在0°构型、30°构型和45°构型时的配平值如表 2所示。γd0为幅值20°的方波信号，γd0作为指令滤波器状态方程 (19) 的输入，输出得到航迹倾斜角指令信号γd，飞行速度保持在30 m/s。总的仿真时间为70 s，仿真开始飞机处于0°构型，在20 s时让飞机以变后掠角速度ωΛ=9°/s从0°构型变化到45 °构型，40 s时以变后掠角速度ωΛ=9°/s从45°构型变化到0°构型。 表 2 配平状态 Table 2 Trim states 图 3给出了分别在文献[19]中最优控制、本文标称控制和加入修正控制时的变后掠角飞机纵向响应仿真结果。图 3中蓝线表示在文献[19]中最优控制下的仿真结果，红线表示在本文标称控制下的仿真结果，黑线表示在加入修正控制时的仿真结果。可以看出在只有标称控制时，由于建模误差和外界干扰的影响，飞行速度和航迹倾斜角跟踪误差最大，而且在20 s飞机从0°构型变化到45°构型的过程中，航迹倾斜角和飞行速度跟踪误差较大，采用标称控制加修正控制时，航迹倾斜角和飞行速度的跟踪误差较小，基本不受变后掠角过程的影响。在文献[19]中最优控制下的指令跟踪误差介于两者之间。 所设计的修正控制系统是一种混合控制系统，即基于两种控制方法的并行控制系统。反步标称控制基于非线性系统进行设计，能够有效利用非线性系统本身固有的非线性特性，对建模精度要求较高。RHO方法能够预测被控对象的期望响应特性，对模型的精度要求不高，相对文献[19]中的最优控制方法，区别主要在于所选取的优化区间。文献[19]中最优控制使用一个对全局相同的性能指标函数，而RHO方法的指标函数只涉及到当前时刻之后的有限时域，当进入下一计算周期，优化区间相应地向前推进。基于RHO的变后掠角飞机修正控制系统根据变后掠角飞机状态输出与指令信号之间的差值实时调节修正控制律增益，如图 3(g)所示，在整个仿真过程中修正控制律增益实时调整，可有效地抑制外界干扰和建模误差的影响。从仿真结果可以看出飞机在变后掠角过程中能够较好地跟踪指令信号，基本不受变后掠角过程、建模误差和外界干扰的影响。这表明所设计的修正控制系统具有较高的鲁棒性，可以保证变后掠角过程中的飞行安全性。 6 结论 (1) 反步控制通过Lyapunov函数进行控制信号递推设计，根据控制指令逐级递推得到最终的真正控制信号。由于标称反步控制量是基于已知模型参数进行解算得到的，因此，模型参数的准确性对控制效果的影响较大。 (2) RH O属于模型预测控制的一种，由预测模型、滚动优化和反馈校正3部分组成，对模型精度要求不高。RHO的优化时域空间随着时间的推移不断变化，实时对控制器参数进行优化更新，由其解算得到的修正控制量可对标称控制量进行在线补偿，可确保飞机在变后掠角过程中稳定飞行。