0 引言

模型飞行试验可以作为验证飞行器性能与技术的一种有效手段^[1]。飞机大迎角过失速机动存在气动特性明显非线性、惯性耦合等强非线性因素，传统的基于线性化小扰动理论设计的线性控制律已不再满足要求，需要发展非线性控制律^[2]。

非线性动态逆控制律是一种广为研究的方法。研究表明过失速机动中动态逆控制律比线性控制律具有更好的性能^[3]。动态逆方法采用时标分离原理将系统分成快、慢回路，快回路增益足够大才能在理论上保证系统的稳定性^[4]，但其受制于舵回路的性能，并不能任意取大。基于Lyapunov方法的反步控制律是得到广泛研究的另一种非线性控制律，其对时标分离尺度没有要求^[5]，但需要对中间虚拟控制信号进行求导，对于飞行控制律设计问题这是极为复杂繁琐的。对此，学者们进一步发展了动态面飞行控制律，通过引入一阶滤波器近似得到虚拟控制信号的导数^[6-7]。

控制律设计中需要用到飞机的气动力模型。过失速机动中飞机迎角会超过失速迎角，此时的空气绕流流场十分复杂，气动力具有明显非定常特征，由静态试验数据、动导数等数据构成的常规气动力模型不能很好地满足过失速机动控制的需要。控制律设计中考虑非定常效应有利于减少气动模型的误差，改善控制效果^[8]。另一方面，飞机舵回路作动器动态响应引起的控制滞后可能对控制效果产生不利影响，作动器偏转速率饱和甚至可能诱发飞机振荡^[9]，因此作动器响应能力是控制律设计中需要考虑的重要问题。Bates等^[10]、Choi等^[11]设计控制律时考虑了作动器的动力学响应。Hess等^[12]采用速率限制器技术有效消除了舵偏速率约束引起的控制时延，该方法实质是将舵偏速率作为控制量。

本文发展了一种改进的动态面飞行控制律，通过有效预测非定常气动力、考虑作动器动力学模型并采用舵偏速率作为实质控制量，改善了飞机过失速机动品质。

1 数学模型

首先定义主要的飞机相关坐标系。飞机机体系b与飞机固连，原点o_b位于飞机质心，o_bx_b轴在飞机对称面内并指向机头，o_by_b轴垂直于飞机对称平面指向机身右方，o_bz_b轴在飞机对称面内指向机身下方。航迹坐标系h与飞机固连，原点o_h位于飞机质心，o_hx_h轴与飞机速度矢量重合，o_hz_h轴位于包含飞机速度矢量的铅垂面内，指向下方，o_hy_h轴通过右手定则确定。地面系g固定于地面，原点o_g位于地面某点，o_gx_g轴在水平面内指向某一方向，o_gz_g轴垂直于水平面并指向地心，o_gy_g轴通过右手定则确定。

描述飞机迎角α、侧滑角β与绕速度矢滚转角μ的微分方程为^[13]:

$ \left[ {\begin{array}{*{20}{c}} {\dot \alpha }\\ {\mathit{\dot \beta }}\\ {\dot \mu } \end{array}} \right]{\rm{ = }}{\mathit{\boldsymbol{M}}_s}\left[ {\begin{array}{*{20}{c}} {\mathit{\dot \gamma }}\\ {\mathit{\dot \chi }{\rm{ }}} \end{array}} \right]{\rm{ + }}{\mathit{\boldsymbol{G}}_s}\left[ {\begin{array}{*{20}{c}} p\\ q\\ r \end{array}} \right] $

(1)

其中，p、q与r分别为飞机滚转、俯仰与偏航角速度；

$ {\mathit{\boldsymbol{M}}_s} = \left[ {\begin{array}{*{20}{c}} { - \frac{{\cos \mathit{\mu }}}{{\cos \beta }}}&{ - \frac{{\sin \mathit{\mu }{\rm{cos}}\gamma }}{{\cos \beta }}}\\ { - \sin \mathit{\mu }}&{\cos \mathit{\mu }{\rm{cos}}\gamma }\\ {\cos \mathit{\mu }{\rm{tan}}\beta }&{\sin \gamma {\rm{ + }}\sin \mathit{\mu }{\rm{cos}}\gamma \tan \beta } \end{array}} \right], $

$ {\mathit{\boldsymbol{G}}_\mathit{s}} = \left[ {\begin{array}{*{20}{c}} { - \cos \alpha \tan \beta }&1&{ - \sin \alpha \tan \beta }\\ {\sin \alpha }&0&{ - \cos \alpha }\\ {\frac{{\cos \alpha }}{{\cos \beta }}}&0&{\frac{{\sin \alpha }}{{\cos \beta }}} \end{array}} \right]; $

γ为航迹倾角，χ为航迹方位角，它们的微分方程为:

$ \begin{array}{*{20}{l}} {\mathit{\dot \gamma } = \frac{1}{{\mathit{mV}}}\left[ {\mathit{L}\cos \mu - \mathit{Y}\sin \mu - mg{\rm{cos}}\gamma {\rm{ + }}} \right.}\\ {{T_x}\left( {\sin \mu {\rm{cos}}\alpha {\rm{sin}}\mathit{\beta }{\rm{ + cos}}\mu {\rm{sin}}\alpha } \right) - }\\ {{T_y}\sin \mu \cos \mathit{\beta }{\rm{ + }}\left. {{T_\mathit{z}}\left( {\sin \mu {\rm{sin}}\beta \sin \alpha - \cos \mu \cos \alpha } \right)} \right]} \end{array} $

(2)

$ \begin{array}{*{20}{l}} {\mathop {\;\;\;\chi = }\limits^. \frac{1}{{\mathit{mVcos\gamma }}}\left[ {L\sin \mu {\rm{ + }}Y{\rm{cos}}\mu {\rm{ + }}} \right.}\\ {{T_x}\left( {\sin \mu {\rm{sin}}\alpha - {\rm{cos}}\mu {\rm{cos}}\alpha {\rm{sin}}\mathit{\beta }} \right){\rm{ + }}}\\ {{T_y}\cos \mu \cos \beta - {T_\mathit{z}}\left( {\cos \mu {\rm{sin}}\mathit{\beta }{\rm{sin}}\alpha + \sin \mu \cos \alpha } \right)} \end{array} $

(3)

其中，V为飞机速度幅值，D、L与Y分别为飞机受到的升力、阻力与侧向力，T_x、T_y与T_z分别为发动机推力在飞机体系x、y与z坐标方向的分量。

飞机角速度微分方程为：

$ \begin{array}{*{20}{l}} {\left[ {\begin{array}{*{20}{c}} {\dot p}\\ {\dot q}\\ {\dot r} \end{array}} \right] = {\mathit{\boldsymbol{I}}^{ - 1}}\left\{ { - \left[ {\begin{array}{*{20}{c}} 0&{ - r}&\mathit{q}\\ r&0&{ - p}\\ { - q}&p&0 \end{array}} \right]\mathit{\boldsymbol{I}}\left[ {\begin{array}{*{20}{c}} \mathit{p}\\ \mathit{q}\\ \mathit{r} \end{array}} \right] + } \right.}\\ {\left. {\;\;\;\;\;\;\;\;\;\left[ {\begin{array}{*{20}{c}} {{l_\alpha }}\\ {{m_\alpha }}\\ {{n_\alpha }} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{l_t}}\\ {{m_t}}\\ {{n_t}} \end{array}} \right]} \right\}} \end{array} $

(4)

其中，I为飞机惯量张量，l_a、m_a与n_a分别为飞机受到的滚转、俯仰与偏航力矩，l_t、m_t与n_t分别为飞机受到的滚转、俯仰与偏航推力矩。式(4) 中力矩项可进一步改写为：

$ \left[ {\begin{array}{*{20}{c}} {{l_\alpha }}\\ {{m_\alpha }}\\ {{n_\alpha }} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{l_t}}\\ {{m_t}}\\ {{n_t}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\hat l}\\ {\hat m}\\ {\hat n} \end{array}} \right] + {\mathit{\boldsymbol{G}}_f}\mathit{\boldsymbol{\delta }} + \Delta \mathit{\boldsymbol{M}} $

(5)

其中，δ=[δ_a  δ_r  δ_e  δ_tvy  δ_tvz]^T分别为副翼、方向舵、升降舵、偏航推力矢量偏角与俯仰推力矢量偏角，${{\left[ \hat{l}\quad \hat{m}\quad \hat{n} \right]}^{\rm{T}}}$为舵面偏角为零时飞机受到的力矩，ΔM表示气动力矩与推力矢量建模的不确定项，上标“T”为矢量转置符号，操纵导数矩阵G_f考虑主要项，为：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{G}}_f} = }\\ {\left[ {\begin{array}{*{20}{c}} {QSb{C_{l{\delta _a}}}}&{QSb{C_{l{\delta _r}}}}&0&{{{\left( {{l_t}} \right)}_{{\delta _{tvy}}}}}&0\\ 0&0&{QSc{C_{m{\delta _e}}}}&0&{{{\left( {{m_t}} \right)}_{{\delta _{tvz}}}}}\\ {QSb{C_{n{\delta _\alpha }}}}&{QSb{C_{n{\delta _r}}}}&0&{{{\left( {{n_t}} \right)}_{{\delta _{tvy}}}}}&0 \end{array}} \right]} \end{array} $

(6)

其中，Q为动压，S为飞机参考面积，b、c分别为飞机横向参考长度和纵向参考长度，C_lδa等代表相应的操纵导数，如C_lδa表示副翼在滚转方向的操纵导数，(m_t)_δtvz代表俯仰推力矩对δ_tvz的偏导数，其它操纵导数意义类推。

飞机姿态控制伺服机构包括气动舵面与推力矢量，它们的作动器模型可以采用一阶惯性环节表示为^[14]：

$ \mathit{\dot\delta = v = }\mathop {{\rm{sat}}}\limits_{{{\mathit{\dot \delta }}_{\max }}} \left( {{\omega _\mathit{\delta }}\left( {{\mathit{\delta }_{{\rm{cmd}}}} - \mathit{\delta }} \right)} \right) $

(7)

其中，δ为舵偏角(指代[δ_a  δ_r  δ_e  δ_tvy  δ_tvz]^T中任一量)，v为舵偏速率，δ_cmd为舵偏指令，ω_δ为作动器带宽，$\underset{{{{\dot{\delta }}}_{\rm{max}}}}{\mathop{\rm{sat}}}\,$()为表示舵偏速率约束的饱和函数。此外，实际舵偏角还存在位置约束，即δ∈[δ_min, δ_max]。

2 过失速机动控制律设计 2.1 控制律设计方法

考虑作动器动力学模型后被控对象为3阶系统，文章首先针对考虑控制约束的一般3阶级联仿射系统给出控制律设计流程。考虑如下系统：

$ \begin{array}{l} {{\mathit{\boldsymbol{\dot x}}}_1} = {\mathit{\boldsymbol{f}}_1}\left( {{\mathit{\boldsymbol{x}}_1}} \right) + {\mathit{\boldsymbol{g}}_1}\left( {{\mathit{\boldsymbol{x}}_1}} \right){\mathit{\boldsymbol{x}}_2} + {\mathit{\boldsymbol{d}}_1}\\ {{\mathit{\boldsymbol{\dot x}}}_2} = {\mathit{\boldsymbol{f}}_2}\left( {{\mathit{\boldsymbol{x}}_1},{\mathit{\boldsymbol{x}}_2}} \right) + {\mathit{\boldsymbol{g}}_2}\left( {{\mathit{\boldsymbol{x}}_1},{\mathit{\boldsymbol{x}}_2}} \right){\mathit{\boldsymbol{x}}_3} + {\mathit{\boldsymbol{d}}_2}\\ {{\mathit{\boldsymbol{\dot x}}}_3} = {\mathit{\boldsymbol{f}}_3}\left( {{\mathit{\boldsymbol{x}}_1},{\mathit{\boldsymbol{x}}_2},{\mathit{\boldsymbol{x}}_3}} \right) + {\mathit{\boldsymbol{g}}_3}\left( {{\mathit{\boldsymbol{x}}_1},{\mathit{\boldsymbol{x}}_2},{\mathit{\boldsymbol{x}}_3}} \right)\mathit{\boldsymbol{u + }}{\mathit{\boldsymbol{d}}_3} \end{array} $

(8)

其中，x₁∈R^n₁、x₂∈R^n₂、x₃∈R^n₃为状态变量，u∈R^m为控制输入，其受到约束|u_i|≤u_imax(i=1, 2, …, m)，变量维数满足m≥n₃≥n₂≥n₁，f₁、f₂与f₃为相应维度的矢量函数，g₁、g₂与g₃为相应维度的矩阵函数，满足rank(g₁)=n₁、rank(g₂)=n₂、rank(g₃)=n₃，d₁、d₂与d₃为相应维度的不确定项。

不妨设原点为平衡点，本文提出的控制律设计步骤可分为两步：

1) 针对x₁、x₂系统，将x₃视为控制(为与状态区别，用x_3d指代)，设计反步控制律^[15]:

$ {\mathit{\boldsymbol{x}}_{3d}} = \mathit{\boldsymbol{g}}_2^{ - 1}\left[ { - {\mathit{\boldsymbol{f}}_2} - {\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {\mathit{\boldsymbol{x}}_{2d}}} \right) - {\bf{g}}_1^{\rm{T}}{\mathit{\boldsymbol{x}}_1} + {{\mathit{\boldsymbol{\dot x}}}_{2d}}} \right] $

(9)

其中，x_2d为中间虚拟控制，其为:

$ {\mathit{\boldsymbol{x}}_{2d}} = \mathit{\boldsymbol{g}}_1^{ - 1}\left( { - \mathit{\boldsymbol{f}} - {\mathit{\boldsymbol{K}}_1}{\mathit{\boldsymbol{x}}_1}} \right) $

(10)

g₁^－1与g₂^－1分别表示矩阵g₁与g₂的(广义)逆，K₁、K₂为控制增益矩阵。引入一阶滤波器:

$ \mathit{\tau }{{\mathit{\boldsymbol{\dot {\tilde x}}}}_{2d}} + {{\mathit{\boldsymbol{\tilde x}}}_{2d}} = {\mathit{\boldsymbol{x}}_{2d}} $

(11)

其中，τ为滤波器时间参数。在式(9) 中用${\mathit{\boldsymbol{\tilde{x}}}}$_2d代替x_2d可得到动态面控制律为:

$ {\mathit{\boldsymbol{x}}_{3d}} = \mathit{\boldsymbol{g}}_2^{ - 1}\left[ { - {\mathit{\boldsymbol{f}}_2} - {\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) - \mathit{\boldsymbol{g}}_1^{\rm{T}}{\mathit{\boldsymbol{x}}_1} + {{\mathit{\boldsymbol{\dot {\tilde x}}}}_{2d}}} \right] $

(12)

注意相对于文献[6]中的动态面控制律，式(12) 多考虑了耦合项g₁^Tx₁。

2) 进一步考虑x₃状态，基于式(12) 定义滑模变量s为:

$ \mathit{\boldsymbol{s = }}{\mathit{\boldsymbol{g}}_2}{\mathit{\boldsymbol{x}}_3} + {\mathit{\boldsymbol{f}}_2} + {\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) + \mathit{\boldsymbol{g}}_1^{\rm{T}}{\mathit{\boldsymbol{x}}_1} - {{\mathit{\boldsymbol{\dot {\tilde x}}}}_{2d}} $

(13)

则通过使滑模变量幅值平方‖s‖²单调减小、并利用最大控制能力可取控制律u为^[16]:

$ \mathit{\boldsymbol{u = }} - {\rm{diag}}\left( {{u_{1\max }},{u_{2\max }}, \cdots ,{u_{m\max }}} \right){\mathop{\rm sgn}} \left( \mathit{\boldsymbol{z}} \right) $

(14)

其中，z=g₃^Tg₂^Ts，diag(a₁, a₂, …, a_m)表示对角元素为a₁, a₂, …, a_m的m维对角方阵，sgn()为符号函数，并定义$\rm{sgn}\left( \mathit{\boldsymbol{z}} \right)=\left[ \begin{matrix} \rm{sgn}({{\mathit{z}}_{1}}) \\ \rm{sgn}({{\mathit{z}}_{2}}) \\ \cdots \\ \rm{sgn}({{\mathit{z}}_{m}}) \\ \end{matrix} \right]$。

关于闭环控制系统的稳定性，在下述假设条件下，将建模不确定性作为扰动输入，可以证明系统为输入-状态稳定。

定义1^[15]：对于系统$\mathit{\boldsymbol{\dot{x}}}$=f(t, x, u)，其中函数f:[0, ∞)×Rⁿ×R^m→Rⁿ关于时间t分段连续，关于变量x与u满足局部Lipschitz条件，如果存在一个KL类函数β和一个K类函数α，对任意的初始状态x(t₀)和任意的有界输入u，对于任意t＞t₀，系统的解x(t)存在并且满足‖x(t)‖≤β(‖x(t₀)‖, t－t₀)+α($\underset{\tau \in \left[ {{t}_{0}},t \right]}{\mathop{\rm{sup}}}\,$‖u(τ)‖), 那么系统为输入-状态稳定。

假设1：矢量函数f₁、f₂、f₃，矩阵函数g₁、g₂、g₃均为Lipschitz函数。

假设2：在控制律(14) 作用下，存在一个正常数δ，使得滑模变量s导数满足s_i${{\dot{s}}_{i}}$≤－δs_isgn(s_i)，(i=1, 2, …, n₂)。

引理1^[17]：对于系统$\mathit{\boldsymbol{\dot{x}}}$=f(x)+g(x)d，如果存在某个Lyapunov函数V，使得‖x‖≥ρ(‖d‖)时，满足$\dot{V}=\frac{\partial V}{\partial \mathit{\boldsymbol{x}}}\mathit{\boldsymbol{\dot{x}}}＜0$，其中函数ρ为K_∞类函数，那么系统为输入-状态稳定。

定理1：在假设条件1、2下，存在一定的滤波器时间参数τ，使得采用控制律(14) 的闭环控制系统(8) 为输入-状态稳定。

证明：因为由Lipschitz函数复合与求和构成的函数也为Lipschitz函数，故存在常数L₁与L₂，使得：

$ {\left\| {{{\mathit{\boldsymbol{\dot x}}}_{2d}}} \right\|^2} \le {L_1}\left( {{{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|}^2} + {{\left\| {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right\|}^2}} \right) $

(15)

$ \begin{array}{l} - {\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^2} - {\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^2} - {\left\| {{{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right\|^2} - {\left\| \mathit{\boldsymbol{s}} \right\|^2} \le \\ - \frac{{\sigma _{\min }^2\left( {{\mathit{\boldsymbol{g}}_2}} \right)}}{{{L_2}}}{\left\| {{\mathit{\boldsymbol{x}}_3}} \right\|^2} \end{array} $

(16)

其中，σ_min(·)表示矩阵的最小奇异值。构造Lyapunov函数：

$ \begin{array}{l} \mathit{V = }{\rm{0}}{\rm{.5}}\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{x}}_1} + 0.5{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)^{\rm{T}}}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) + \\ 0.5a{\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)^{\rm{T}}}\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) + 0.5b{\mathit{\boldsymbol{s}}^{\rm{T}}}\mathit{\boldsymbol{s}} \end{array} $

(17)

其中，$a=\frac{2{{C}_{2}}}{{{L}_{1}}}c,b=\frac{\|{{\mathit{\boldsymbol{s}}}_{0}}\|}{\delta {{\sigma }_{\rm{min}}}({{\mathit{\boldsymbol{K}}}_{2}})},c=\frac{1}{\tau }+\sqrt{\frac{1}{{{\tau }^{2}}}-\frac{{{L}_{1}}{{C}_{1}}}{{{C}_{2}}}},$${{C}_{1}}=\frac{{{\sigma }_{\rm{max}}}{{({{\mathit{\boldsymbol{g}}}_{1}})}^{2}}}{{{\sigma }_{\rm{min}}}({{\mathit{\boldsymbol{K}}}_{1}})},{{C}_{2}}=\frac{1}{8}\rm{min}{{\sigma }_{\rm{min}}}({{\mathit{\boldsymbol{K}}}_{1}}),\rm{ }{{\sigma }_{\rm{min}}}({{\mathit{\boldsymbol{K}}}_{2}})$，s₀为初始时刻t₀的滑模变量值，σ_max(·)表示矩阵的最大奇异值，min(·, ·)表示取最小值。显然当τ满足下述关系式时，c存在实数解。

$ \mathit{\tau } \le \sqrt {\frac{{{C_2}}}{{{L_1}{C_1}}}} $

(18)

对式(17) 求导，有：

$ \begin{array}{*{20}{l}} {\mathit{ \dot {\;\;V =} }\mathit{\boldsymbol{x}}_1^{\rm{T}}{{\mathit{\boldsymbol{\dot x}}}_1} + {{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}\left( {{{\mathit{\boldsymbol{\dot x}}}_2} - {{\mathit{\boldsymbol{\dot {\tilde x}}}}_{2d}}} \right) + }\\ {\;\;\;\;\;\;a{{\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}\left( {{{\mathit{\boldsymbol{\dot x}}}_{2d}} - {{\mathit{\boldsymbol{\dot {\tilde x}}}}_{2d}}} \right) + b{\mathit{\boldsymbol{s}}^{\rm{T}}}\mathit{\boldsymbol{\dot s}}}\\ { = - \mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{K}}_1}{\mathit{\boldsymbol{x}}_1} - {{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) - }\\ {\;\;\;\frac{a}{\tau }{{\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) + }\\ {\;\;\;a{{\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}{{\mathit{\boldsymbol{\dot x}}}_{2d}} + \mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{g}}_1}\left( {{\mathit{\boldsymbol{x}}_1}} \right)\left( {{{\mathit{\boldsymbol{\tilde x}}}_{2d}} - {\mathit{\boldsymbol{x}}_{2d}}} \right) + }\\ {\;\;\;{{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}s + b{\mathit{\boldsymbol{s}}^{\rm{T}}}\mathit{\boldsymbol{\dot s}} + }\\ {\;\;\;\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{d}}_1} + {{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{d}}_2}} \end{array} $

(19)

利用Young’s不等式，可得：

$ {\left( {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)^{\rm{T}}}{{\mathit{\boldsymbol{\dot x}}}_{2d}} \le \frac{c}{2}{\left\| {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right\|^2} + \frac{1}{{2c}}{\left\| {{{\dot x}_{2d}}} \right\|^2} $

(20)

$ \mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{g}}_1}\left( {{\mathit{\boldsymbol{x}}_1}} \right)\left( {{{\mathit{\boldsymbol{\tilde x}}}_{2d}} - {\mathit{\boldsymbol{x}}_{2d}}} \right) \le \frac{1}{2}\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{K}}_1}{\mathit{\boldsymbol{x}}_1} + \frac{{{C_1}}}{2}{\left\| {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right\|^2} $

(21)

$ \begin{array}{*{20}{c}} {{{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{s}} \le \frac{1}{2}{{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) + }\\ {\frac{1}{{2{\sigma _{\min }}\left( {{\mathit{\boldsymbol{K}}_2}} \right)}}{{\left\| \mathit{\boldsymbol{s}} \right\|}^2}} \end{array} $

(22)

$ \mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{d}}_1} \le \frac{1}{4}\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{K}}_1}{\mathit{\boldsymbol{x}}_1} + \mathit{\boldsymbol{d}}_1^{\rm{T}}\mathit{\boldsymbol{K}}_1^{ - 1}{\mathit{\boldsymbol{d}}_1} $

(23)

$ {\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{d}}_2} \le \frac{1}{4}{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) + {\mathit{\boldsymbol{d}}_2}\mathit{\boldsymbol{K}}_2^{ - 1}{\mathit{\boldsymbol{d}}_2} $

(24)

将式(20~24)、(15) 代入式(19) 并化简，有

$ \mathit{\dot {\;\;\;\;V < q}} - \frac{1}{{2{\sigma _{\min }}\left( {{\mathit{\boldsymbol{K}}_2}} \right)}}{\left\| \mathit{\boldsymbol{s}} \right\|^2} + \frac{1}{{8{C_2}}}{\left\| \mathit{\boldsymbol{d}} \right\|^2}\; $

(25)

其中，d=[d₁^T  d₂^T  d₃^T]^T，

$ \begin{array}{l} q = - \frac{1}{8}\mathit{\boldsymbol{x}}_1^{\rm{T}}{\mathit{\boldsymbol{K}}_1}{\mathit{\boldsymbol{x}}_1} - \frac{1}{8}{\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{K}}_2}\left( {{\mathit{\boldsymbol{x}}_2} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right) - \\ \frac{{{C_1}}}{2}{\left\| {{\mathit{\boldsymbol{x}}_{2d}} - {{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right\|^2} \end{array} $

(26)

对于q，利用反证法可以证明存在一个常数e＞0，使得：

$ q \le - e{\left\| {{\mathit{\boldsymbol{x}}_1}} \right\|^2} - e{\left\| {{\mathit{\boldsymbol{x}}_2}} \right\|^2} - e{\left\| {{{\mathit{\boldsymbol{\tilde x}}}_{2d}}} \right\|^2} $

(27)

取$k=\rm{min}(\frac{e}{2},\rm{ }\frac{1}{2{{\sigma }_{\rm{min}}}({{\mathit{\boldsymbol{K}}}_{2}})})$，利用式(16) 并适当推导可得

$ \mathit{\dot {\;\;V <} } - {\rm{min}}\left( {\frac{e}{2},\frac{{k{\sigma _{\min }}^2\left( {{\mathit{\boldsymbol{g}}_2}} \right)}}{{{L_2}}}} \right){\left\| \mathit{\boldsymbol{x}} \right\|^2} + \frac{1}{{8{C_2}}}{\left\| \mathit{\boldsymbol{d}} \right\|^2} $

(28)

其中，$\mathit{\boldsymbol{x}}={{\left[ \mathit{\boldsymbol{x}}_{1}^{T}\quad \mathit{\boldsymbol{x}}_{2}^{T}\quad \mathit{\boldsymbol{x}}_{3}^{T}\quad \mathit{\boldsymbol{\tilde{x}}}_{2d}^{T} \right]}^{\rm{T}}}$。显然当$\|\mathit{\boldsymbol{x}}\|\ge \sqrt{\frac{1}{8{{C}_{2}}\rm{min}(\frac{e}{2},\rm{ }\frac{k{{\sigma }_{\rm{min}}}^{2}({{\mathit{\boldsymbol{g}}}_{2}})}{{{L}_{2}}})}}\left\| \mathit{\boldsymbol{d}} \right\|$时，${\dot{V}}$＜0，故根据引理1，闭环系统为输入-状态稳定。

2.2 动态面飞行控制律设计

依据2.1节中的控制律设计流程，本文首先针对气流角与角速度回路设计基本动态面飞行控制律。飞机前馈控制指令包括迎角指令α_c、侧滑角指令β_c与绕速度矢滚转角指令μ_c，为了提高迎角控制品质，设置迎角指令参考模型^[13]：

$ \frac{{{\alpha _r}\left( \mathit{s} \right)}}{{{\alpha _c}\left( \mathit{s} \right)}} = \frac{{\omega _\alpha ^2}}{{{s^2} + 2{\mathit{\xi }_\alpha }{\omega _\alpha }\mathit{s + }\omega _\alpha ^2}} $

(29)

其中，α_r为参考迎角，ξ_α与ω_α分别为模型阻尼参数与频率参数。为方便后续舵面控制分配，引入角动量h=I[p  q  r]^T以消除惯性积的交联影响。针对系统(1)、(4)，将${\mathit{\boldsymbol{\dot{h}}}}$视为控制，构造Lyapunov函数：

$ {\mathit{V}_2} = 0.5{\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mathit{\mu } - {\mathit{\mu }_c}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mathit{\mu } - {\mathit{\mu }_c}} \end{array}} \right]{\rm{ + }}0.5{\left( {\mathit{\boldsymbol{h}} - {\mathit{\boldsymbol{h}}_c}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{h}} - {\mathit{\boldsymbol{h}}_c}} \right) $

(30)

其中，h_c为中间虚拟控制，其为：

$ {\mathit{\boldsymbol{h}}_c} = \left[ {\begin{array}{*{20}{c}} {{h_{cx}}}\\ {{h_{cy}}}\\ {{h_{cz}}} \end{array}} \right] = \mathit{\boldsymbol{IG}}_S^{ - 1}\left\{ { - {\mathit{\boldsymbol{M}}_s}\left[ {\begin{array}{*{20}{l}} {\mathit{\dot \gamma }}\\ {\mathit{\dot \chi }} \end{array}} \right] + \left[ {\begin{array}{*{20}{l}} {{{\dot \alpha }_r}}\\ {{{\dot \beta }_c}}\\ {{{\dot \mu }_c}} \end{array}} \right] - {\mathit{\boldsymbol{K}}_1}\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mu - {\mu _c}} \end{array}} \right]} \right\} $

(31)

其中，K₁为气流角回路增益矩阵，$\mathit{\boldsymbol{G}}_{s}^{-1}=\left[ \begin{matrix} 0 & \rm{sin}\alpha & \rm{cos}\,\alpha \,\rm{cos}\beta \\ 1 & 0 & \rm{cos}\,\beta \,\rm{tan}\beta \\ 0 & -\rm{cos}\alpha & \rm{sin}\,\alpha \,\rm{cos}\beta \\ \end{matrix} \right]$。通过使${\dot{V}}$₂≤0可推导出控制${\mathit{\boldsymbol{\dot{h}}}}$_d为：

$ {{\mathit{\boldsymbol{\dot h}}}_d} = {{\mathit{\boldsymbol{\dot h}}}_c} - {\mathit{\boldsymbol{K}}_2}(\mathit{\boldsymbol{h}} - {\mathit{\boldsymbol{h}}_c}) - {\mathit{\boldsymbol{I}}^{ - 1}}\mathit{\boldsymbol{G}}_{\rm{s}}^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mathit{\mu } - {\mathit{\mu }_c}} \end{array}} \right] $

(32)

其中，K₂为角动量回路增益矩阵。引入一阶滤波器

$ {\mathit{\tau }_{hi}}{{\dot {\tilde h}}_{ci}} + {{\tilde h}_{ci}} = {h_{ci}}\left( {i = x,y,z} \right) $

(33)

其中，τ_hi(i=x, y, z)为时间参数。在式(32) 中用${{{\mathit{\boldsymbol{\tilde{h}}}}}_{c}}$代替h_c可得控制${\mathit{\boldsymbol{\dot{h}}}}$_d为：

$ {{\mathit{\boldsymbol{\dot h}}}_d}{{\mathit{\boldsymbol{\dot {\tilde h}}}}_c} - {\mathit{\boldsymbol{K}}_2}(\mathit{\boldsymbol{h}} - {{\mathit{\boldsymbol{\tilde h}}}_c}) - {\mathit{\boldsymbol{I}}^{ - 1}}\mathit{\boldsymbol{G}}_\mathit{s}^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mathit{\mu } - {\mathit{\mu }_c}} \end{array}} \right] $

(34)

利用式(4)、(5) 可以得到舵偏角指令δ_cmd为：

$ {\mathit{\boldsymbol{\delta }}_{{\rm{cmd}}}} = \mathit{\boldsymbol{G}}_f^{ - 1}\left( {{{\mathit{\boldsymbol{\dot {\tilde h}}}}_c} - {\mathit{\boldsymbol{K}}_2}(\mathit{\boldsymbol{h}} - {{\mathit{\boldsymbol{\tilde h}}}_c}) - {\mathit{\boldsymbol{I}}^{ - 1}}\mathit{\boldsymbol{G}}_\mathit{s}^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mathit{\mu } - {\mathit{\mu }_c}} \end{array}} \right] - \left[ \begin{array}{l} {f_{hx}}\\ {f_{hy}}\\ {f_{hz}} \end{array} \right]} \right) $

(35)

其中，G_f^－1为非方矩阵G_f的广义逆，

$ \left[ \begin{array}{l} {f_{hx}}\\ {f_{hy}}\\ {f_{hz}} \end{array} \right] = - \left[ {\begin{array}{*{20}{c}} 0&{ - r}&q\\ r&0&{ - p}\\ { - q}&p&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{h_x}}\\ {{h_y}}\\ {{h_\mathit{z}}} \end{array}} \right] + \left[ \begin{array}{l} {\hat l}\\ {\hat m}\\ {\hat n} \end{array} \right]。$

2.3 改进动态面飞行控制律设计

进一步考虑作动器动力学模型导出滑模控制指令，舵面运动模型为：

$ \left[ \begin{array}{l} {{\mathit{\dot \delta }}_a}\\ {{\mathit{\dot \delta }}_r}\\ {{\mathit{\dot \delta }}_e}\\ {{\mathit{\dot \delta }}_{tvy}}\\ {{\mathit{\dot \delta }}_{tvz}} \end{array} \right] = \left[ \begin{array}{l} {\mathit{v}_a}\\ {\mathit{v}_r}\\ {\mathit{v}_e}\\ {v_{tvy}}\\ {v_{tvz}} \end{array} \right] $

(36)

其中，v_a、v_r与v_e分别为副翼偏转速率、方向舵偏转速率与升降舵偏转速率，v_tvy与v_tvz分别为偏航推力矢量偏转速率与俯仰推力矢量偏转速率。

将式(34) 等号右边项移到左边，可以定义滑模变量s=[s_x  s_y  s_z]^T为：

$ \mathit{\boldsymbol{s}} = {{\mathit{\boldsymbol{\dot e}}}_\mathit{\boldsymbol{h}}} + {\mathit{\boldsymbol{K}}_{\rm{2}}}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{h}}} + {I^{ - 1}}\mathit{\boldsymbol{G}}_{\rm{s}}^{\rm{T}}\left[ {\begin{array}{*{20}{c}} {\alpha - {\alpha _r}}\\ {\beta - {\beta _c}}\\ {\mathit{\mu } - {\mathit{\mu }_c}} \end{array}} \right] $

(37)

其中，${{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{h}}}}=\mathit{\boldsymbol{h}}-{{{\mathit{\boldsymbol{\tilde{h}}}}}_{c}},{{{\mathit{\boldsymbol{\dot{e}}}}}_{\mathit{\boldsymbol{h}}}}=\mathit{\boldsymbol{\dot{h}}}-{{{\mathit{\boldsymbol{\dot{\tilde{h}}}}}}_{c}}$。依式(14)，可推导控制律为：

$ \left[ \begin{array}{l} {\mathit{v}_a}\\ {\mathit{v}_r}\\ {\mathit{v}_e}\\ {v_{tvy}}\\ {v_{tvz}} \end{array} \right] = - \left[ {\begin{array}{*{20}{c}} {{{\mathit{\dot \delta }}_{a\max }}}&{}&{}&{}&{}\\ {}&{{{\mathit{\dot \delta }}_{r\max }}}&{}&{}&{}\\ {}&{}&{{{\mathit{\dot \delta }}_{e\max }}}&{}&{}\\ {}&{}&{}&{{{\mathit{\dot \delta }}_{y\max }}}&{}\\ {}&{}&{}&{}&{{{\mathit{\dot \delta }}_{z\max }}} \end{array}} \right]{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{G}}_f}^{\rm{T}}\mathit{\boldsymbol{s}}} \right) $

(38)

其中，${\dot{\delta }}$_amax、${\dot{\delta }}$_rmax与${\dot{\delta }}$_emax分别为副翼、方向舵与升降舵最大偏转速率，${\dot{\delta }}$_ymax与${\dot{\delta }}$_zmax分别为俯仰推力矢量与偏航推力矢量最大偏转速率。

对于俯仰通道相关的升降舵δ_e与推力矢量δ_tvz，由于sgn(ab)=sgn(a)·sgn(b)对任意标量a、b成立，同时考虑到符号函数会引起抖振效应，取

$ {\mathit{v}_e} = - {{\mathit{\dot \delta }}_{e\max }}{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{G}}_f}\left( {2,3} \right)} \right){{\mathop{\rm sgn}} _m}\left( {{s_y}} \right) $

(39)

$ {v_{tvz}} = - {{\mathit{\dot \delta }}_{z\max }}{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{G}}_f}\left( {2,5} \right)} \right){{\mathop{\rm sgn}} _m}\left( {{s_y}} \right) $

(40)

其中，sgn_m()是为削弱滑模控制抖振引入的改进符号函数^[16]，其定义为：

$ {{\mathop{\rm sgn}} _m}\left( {{s_y}} \right) = \left\{ {\begin{array}{*{20}{c}} {{\mathop{\rm sgn}} \left( {{s_y}} \right),}&{\left| {{s_y}} \right| > \mathit{\varepsilon }}\\ {\frac{{{s_y}}}{\mathit{\varepsilon }},}&{\left| {{s_y}} \right| \le \mathit{\varepsilon }} \end{array}} \right. $

(41)

其中，ε为边界层厚度。滑模控制中引入边界层会影响s的控制精度，但其仍能保证有界误差，它是滑模控制应用中削弱抖振问题的有效手段。

对于与偏航-滚转通道相关的副翼δ_a、方向舵δ_r与推力矢量δ_tvy，考虑到副翼主要提供滚转力矩，方向舵主要提供偏航力矩，(l_t)_δtvy一般为0，故取:

$ {\mathit{v}_a} = - {{\mathit{\dot \delta }}_{a\max }}{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{G}}_f}\left( {1,1} \right)} \right){{\mathop{\rm sgn}} _m}\left( {{s_x}} \right) $

(42)

$ {\mathit{v}_r} = - {{\mathit{\dot \delta }}_{r\max }}{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{G}}_f}\left( {3,2} \right)} \right){{\mathop{\rm sgn}} _m}\left( {{s_z}} \right) $

(43)

$ {v_{tvy}} = - {{\mathit{\dot \delta }}_{y\max }}{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{G}}_f}\left( {3,4} \right)} \right){{\mathop{\rm sgn}} _m}\left( {{s_z}} \right) $

(44)

采用这种控制分配，副翼在偏航通道的干扰项QSbC_nδav_a由偏航推力矢量消除，方向舵在滚转通道的干扰项QSbC_lδrv_r由副翼消除。

2.4 舵偏角指令计算

实际作动器接受的是舵偏角指令，2.3节给出了舵偏速率指令${\dot{\delta }}$_cmd(指代[v_a  v_r  v_e  v_tvy  v_tvz]^T中任一量)，根据式(7)，同时考虑舵偏位置限制，可将舵偏速率指令转化为舵偏角指令:

$ {\mathit{\delta }_{{\rm{cmd}}}} = \mathop {{\rm{sat}}}\limits_{{\mathit{\delta }_{\max }}} \left( {\frac{{{{\mathit{\dot \delta }}_{{\rm{cmd}}}}}}{{{\mathit{\omega }_\mathit{\delta }}}} + \mathit{\delta }} \right) $

其中，$\underset{{{\delta }_{\rm{max}}}}{\mathop{\rm{sat}\left( {} \right)}}\,$为表示舵偏位置约束的饱和函数。由于舵偏角指令是通过作动器模型逆解算得到，控制律给出的控制量实质仍是舵偏速率，因此ω_δ基本不影响最终控制效果。

计算舵偏角指令时式(45) 中ω_δ的数值可能与实际作动器带宽存在误差，这对舵偏指令结果的影响程度是一个需要考虑的问题。分析表明ω_δ采用保守取值仅引起控制律中改进符号函数(41) 中边界层厚度ε的改变，对控制结果的影响不大，分析如下：

假设计算舵偏角指令时采用的名义作动器带宽参数为${\hat{\omega }}$_δ，它满足${\hat{\omega }}$_δ＜ω_δ，将式(45) 带入式(7)，显然按真实ω_δ计算的舵偏角指令$\left( {\frac{{{{\dot \delta }_{{\rm{cmd}}}}}}{{{\omega _\delta }}} + \delta } \right)$受到位置约束作用时，采用${\hat{\omega }}$_δ的计算结果并不影响最终指令；若$\left( \frac{{{{\dot{\delta }}}_{\rm{cmd}}}}{{{\omega }_{\delta }}}+\delta \right)\in \left[ {{\delta }_{\rm{min}}},{{\delta }_{\rm{max}}} \right]$，采用${\hat{\omega }}$_δ计算时的实际舵偏速率${\dot{\delta }}$为

$ \mathit{\dot {\;\;\delta =} }\mathop {{\rm{sat}}}\limits_{{{\mathit{\dot \delta }}_{\max }}} \left( {\frac{{{\mathit{\omega }_\mathit{\delta }}}}{{{{\mathit{\hat \omega }}_\mathit{\delta }}}} + {{\mathit{\dot \delta }}_{{\rm{cmd}}}}} \right) $

(46)

当$\frac{{{\omega }_{\delta }}}{{{{\hat{\omega }}}_{\delta }}}{{{\dot{\delta }}}_{\rm{cmd}}}\in \left[ -{{{\dot{\delta }}}_{\rm{max}}},{{{\dot{\delta }}}_{\rm{max}}} \right]$时，$\dot{\delta }=\frac{{{\omega }_{\delta }}}{{{{\hat{\omega }}}_{\delta }}}{{{\dot{\delta }}}_{\rm{cmd}}}$；当$\left| \frac{{{\omega }_{\delta }}}{{{{\hat{\omega }}}_{\delta }}}{{{\dot{\delta }}}_{\rm{cmd}}} \right|\ge {{{\dot{\delta }}}_{\rm{max}}}$时，$\left| {\dot{\delta }} \right|={{{\dot{\delta }}}_{\rm{max}}}$。综上，采用${\hat{\omega }}$_δ计算舵偏角指令等价于式(41) 中边界层厚度变为$\frac{{{{\hat{\omega }}}_{\delta }}}{{{\omega }_{\delta }}}\varepsilon $。

3 非线性、非定常气动力(矩)预测

非线性控制律设计中需要预测飞机气动力(矩)，准确的预测结果有利于提高过失速机动的控制品质。控制律中计算h_c时，${\dot{\gamma }}$与${\dot{\chi }}$与飞机受到的气动力有关，对于过失速机动中复杂的非线性、非定常气动力，文章不是通过气动力模型进行预测，而是利用飞机上的过载等测量量予以求解，包括气动力在内的外力在航迹坐标系下的表示F^h为:

$ {\mathit{\boldsymbol{F}}^h} = \left[ \begin{array}{l} F_x^h\\ F_y^h\\ F_z^h \end{array} \right] = mg\left( {\mathit{\boldsymbol{R}}_b^h\mathit{\boldsymbol{n + R}}_g^h\left[ \begin{array}{l} 0\\ 0\\ 1 \end{array} \right]} \right) $

(47)

其中，m为飞机质量，g为重力加速度，n为飞机过载，R_b^h=M_x(μ)M_z(－β)M_y(－α)为体系b到航迹坐标系h的方向余弦阵，R_g^h=M_y(γ)M_z(χ)为地面坐标系g到航迹坐标系h的方向余弦阵，M_y(α)表示绕y轴旋转角度α，其它类推。基于航迹动力学方程则可以求得${\dot{\gamma }}$与${\dot{\chi }}$为:

$ \mathit{\dot {\;\;\;\gamma =} } - \frac{{F_z^h}}{{m\mathit{V}}} $

(48)

$ \mathit{\dot \chi } = \frac{{\mathit{F}_\mathit{y}^\mathit{h}}}{{\mathit{mV}\cos \mathit{\gamma }}} $

(49)

在式(37) 中求解滑模变量时需要估计飞机受到的气动力矩，由于过失速机动中气动力非定常效应明显，控制律中内嵌的纵向气动力矩模型除考虑静态与动态气动力矩外，还考虑非定常气动力矩。非定常气动力矩增量ΔC_{m, unst}采用微分方程模型^[18]描述，具体形式为:

$ \begin{array}{*{20}{l}} {\frac{{{\rm{d}}\Delta {C_{m,{\rm{unst}}}}}}{{{\rm{d}}t}} = \frac{V}{c}\left\{ { - {k_1}\left( \alpha \right)\Delta {C_{m,{\rm{unst}}}} + {e_0}\left( \alpha \right) + {e_1}\left( \alpha \right){\mathit{\beta }^2}} \right\}{\rm{ + }}}\\ {\frac{{\dot \alpha }}{2}\left( {{\rm{sgn}}\left( {\dot \alpha } \right) + 1} \right)\left[ {{f_0}\left( \alpha \right) + {f_1}\left( \alpha \right){\mathit{\beta }^2}} \right] + }\\ {\frac{{\ddot \alpha }}{2}\left( {{\rm{sgn}}\left( {\dot \alpha } \right) + 1} \right)\left[ {{g_0}\left( \alpha \right) + {g_1}\left( \alpha \right){\mathit{\beta }^2}} \right]} \end{array} $

(50)

其中，k₁(α)，e₀(α)，e₁(α)，f₀(α)，f₁(α)，g₀(α)，g₁(α)为与迎角α有关的系数。

4 算例仿真与分析

针对缩比模型飞机进行控制律验证，文章同时还与动态面控制律进行对比，动态面控制律中也采用第3节的方式预测非定常气动力(矩)。仿真计算在Matlab/Simulink平台上进行，其中假设气动模型与推力矢量模型均存在40%的建模误差，同时考虑舵偏位置约束与速率约束。

4.1 小迎角跟踪控制

首先针对小迎角控制考察不同作动器带宽参数对控制效果的影响。迎角指令为周期3s、振幅1°的方波信号，飞机初始迎角为10°，图 1给出了飞机迎角对参考迎角的跟踪误差曲线，对于动态面控制律，图中给出了ω_δ=20rad/s、ω_δ=40rad/s与ω_δ=+∞(代表理想作动器)三种情形的仿真结果，可以看到随着作动器性能的提高(即带宽的增大)，迎角跟踪误差不断减小；对于本文控制律，由于实质控制量为舵偏速率，ω_δ=20rad/s与ω_δ=40rad/s时的仿真结果重合，最大跟踪误差优于采用理想作动器的动态面控制律约0.01°，虽然数值很小，但也说明本文控制律具有更好的控制效果。

4.2 “眼镜蛇”机动

“眼镜蛇”机动是典型的过失速机动，动作简单，便于考验控制律的控制品质。飞机初始迎角为10°，初始侧滑角为0°，机动时长为2.5s。在0~1s，α_c=70°，1~2.5s，α_c=10°；机动中β_c与μ_c一直为0°，ω_δ=40rad/s，发动机推力为T=200N。图 2给出了迎角仿真曲线，从图中看到采用动态面控制律时，作动器的动态响应过程对仿真结果存在较大影响，而采用本文控制律的迎角仿真结果与采用理想作动器的动态面控制律结果接近，能较好的克服作动器响应滞后对控制效果的影响。

“眼镜蛇”机动主要是纵平面的俯仰运动。图 3给出了升降舵与俯仰推力矢量偏转曲线，机动中升降舵与俯仰推力矢量偏转发生了位置饱和与速率饱和；动态面控制律的舵面运动在飞机迎角拉起中比本文控制律结果略有滞后；相较采用理想作动器的动态面控制律，本文控制律满足作动器的性能约束；此外本文控制律的舵面运动曲线总体比较光滑，这说明控制律中采用改进符号函数有效削弱了滑模控制的抖振问题。

考察非定常气动力效应对过失速机动的影响，图 4给出了机动过程中的非定常气动力矩示意图，伴随着迎角的改变其大致呈顺时针环形。图 5对比了控制器中考虑非定常气动力矩与否的迎角曲线，根据超调量等指标对比，可以看到控制律中考虑非定常气动力效应有利于改善过失速机动品质。

4.3 Herbst机动

Herbst机动综合了飞机动态拉起迎角进入过失速状态、大迎角下绕速度矢滚转等基本过失速机动动作，是检验飞机过失速机动能力的标准验证动作。文章进一步针对Herbst机动考察控制律效果，飞机初始迎角为10°，初始侧滑角为5°，机动时长为22s，机动中β_c=0°，α_c、μ_c与发动机推力指令通过跟踪离线优化的轨迹计算^[13]，ω_δ=40rad/s。图 6给出了三维空间中的Herbst机动轨迹，可以看到本文控制律与动态面控制律均较好地实现了180°快速转弯机动。

图 7给出了Herbst机动中控制律对迎角与侧滑角的控制误差，两种控制律在迎角深失速区域由于飞机绕速度矢量滚转都出现了一定的误差，不过相较于动态面控制律，本文控制律的误差更小。

图 8与图 9分别给出了Herbst机动中的气动舵面与推力矢量偏转曲线，可以看到两种控制律的舵面运动总体类似，机动中气动舵面与推力矢量的运动曲线没有明显的抖振现象。

同样考察Herbst机动中非定常气动力效应的影响，图 10对比了控制律中考虑非定常气动力矩与否的迎角控制曲线，结果再次表明考虑非定常气动力效应总体上能达到更高的机动控制精度。

5 结论

本文发展了一种考虑非定常气动力效应与舵回路作动器模型的改进动态面飞行控制律，考察了飞机小迎角控制、“眼镜蛇”机动与Herbst机动，主要结论为：

1) 控制律中对非定常气动力效应的有效预测有利于改善过失速机动品质。

2) 由于综合了作动器动力学模型，控制律对作动器带宽变化不敏感，可以有效消除作动器动态响应引起控制效果变差的问题。