可重复使用运载器(Reusable launch vehicles, RLV)的再入制导与控制一直以来是各国研究的重点. RLV飞行包线大、再入环境复杂、约束苛刻、非线性与不确定特征强，同时，再入段大气较为稀薄，气动舵面更易发生饱和、部分失效等故障，这都为制导与控制系统的设计带来极大挑战^[1].因此，有必要研究能够抑制不确定/干扰、舵面失效及饱和影响的具备较强鲁棒性能的容错控制策略，并引导RLV在满足再入约束的条件下到达目标位置，这对于再入飞行任务的成功来说至关重要.

针对RLV的制导律设计问题，数值预测校正制导对初始误差不敏感，落点精度较高，鲁棒性较好，目前越来越受到关注^[2].Lu^[3]针对不同升阻比的再入飞行器，提出统一的预测校正制导算法，但计算量较大. Xue等^[4]用准平衡滑翔条件(Quasi-equilibrium glide condition，QEGC)将再入过程约束转换为倾侧角幅值约束，但在QEGC飞行特性降低的情况下，难以保证过载不突破限制.张鹏等^[5]将哥氏加速度引入QEGC并设计二次函数倾侧角幅值模型，在满足过程约束条件下以提高再入制导精度和鲁棒性.针对RLV的不确定性/干扰下的控制问题，Wang等^[6]在反步法中引入鲁棒自适应律来估计未知不确定的上界. Hu等^[7]通过自适应模糊逼近和动态面控制方法来减小不确定性的影响，但学习参数较多将影响实时性.对于舵面失效及饱和的飞行器容错控制问题，鲁棒、自适应和观测器方法仍是较为有效的手段.于彦波等^[8]针对航天器执行器故障和控制受限问题设计了基于积分滑模面的自适应姿态容错控制律来提高系统的鲁棒性能.Lee等^[9]等使用自适应线性未知输入观测器进行模型匹配，虽然可以有效抑制扰动且能处理系统中的未知非线性和舵面故障，但飞行器模型的线性化过程中存在一定的失真.

综上所述，国内外学者对于RLV再入制导律和姿态控制器的研究通常是两个相对独立的领域.然而在实际飞行中，制导系统与姿控系统须进行整合以形成完整的六自由度再入系统，才能检验飞行器能否完成再入任务，满足要求的制导精度.由于目前公开的RLV和高超声速飞行器气动资料较少，缺乏包含大攻角以及完整横侧向气动参数的飞行器模型^[10]，故该方面验证结果较少. Tian等^[11]同时考虑制导环和姿态环，采用伪谱法进行制导律设计并利用二阶滑模控制器和干扰观测器对制导指令进行跟踪. Qian等^[12]采用预测制导法生成标称轨迹，将轨迹指令输入到姿态环，应用反步法实现了飞行器的姿态控制.严晗^[13]针对导弹控制系统，直接进行制导控制一体化设计，但RLV的再入飞行存在强非线性特性，且受到各种约束的限制，姿态控制相比导弹更加复杂，直接进行制导与姿态一体化控制设计的难度极大.

HORUS-2B飞行器^[14]是近年来公开数据的一种RLV模型，它具备完整的纵横向气动参数数据，迎角范围为0°~45°，马赫数在1.2~20.0之间，完全能够开展再入段制导与控制研究.本文针对该飞行器，设计了约束预测校正制导律，对于姿态系统中的不确定/干扰、作动器部分失效及饱和问题，提出了改进跟踪微分干扰观测器(Improved tracking differential disturbance observer, ITDDO)来估计不确定/干扰及舵面部分失效作用，并设计辅助抗饱和系统来解决舵面饱和问题.之后，对上述制导与控制方法进行了六自由度联合仿真.实验结果表明该方法能够引导RLV在满足再入约束的条件下到达终端落点区域，并具备良好的鲁棒容错控制性能.

1 问题描述

考虑圆球形大地和地球自转的情况，RLV经量纲一的转化后以能量e为自变量的三自由度制导方程以及姿态系统方程^[11]如下：

$ \frac{\mathrm{d} R}{\mathrm{~d} e}=\frac{\sin \gamma}{D}, \frac{\mathrm{d} \theta}{\mathrm{d} e}=\frac{\cos \gamma \sin \psi}{R D \cos \varphi}, $

(1)

$ \frac{\mathrm{d} \varphi}{\mathrm{d} e}=\frac{\cos \gamma \sin \psi}{R D}, \frac{\mathrm{d} s_{\mathrm{f}}}{\mathrm{d} e}=\frac{-\cos \gamma}{R D}, $

(2)

$ \begin{array}{l} \frac{\mathrm{d} \gamma}{\mathrm{d} e}=\frac{1}{V^{2} D}\left[L \cos \sigma+\left(V^{2}-\frac{1}{R}\right)\left(\frac{\cos \gamma}{R}\right)+2 \omega_{\mathrm{E}} V \cos \varphi \sin \psi\right]+\\ \frac{\omega_{\mathrm{E}}^{2} R}{V^{2} D} \cos \varphi(\cos \gamma \cos \varphi+\sin \gamma \cos \psi \sin \varphi), \end{array} $

(3)

$ \begin{array}{l} \frac{\mathrm{d} \psi}{\mathrm{d} e}=\frac{L \sin \sigma}{V^{2} D \cos \gamma}+\frac{1}{R D} \cos \gamma \sin \psi \tan \varphi-\frac{1}{V D} 2 \omega_{\mathrm{E}} \cdot \\ (\tan \gamma \cos \psi \cos \varphi-\sin \varphi)+\frac{\omega_{\mathrm{E}}^{2} R}{V^{2} D \cos \gamma} \sin \psi \sin \varphi \cos \varphi, \end{array} $

(4)

$ \begin{aligned} {\dot \alpha}=& q-(p \cos \alpha+r \sin \alpha) \tan \beta-\sin \sigma / \cos \beta \cdot \\ &\left[\dot{\psi} \cos \gamma-\dot{\varphi} \sin \psi \sin \gamma+\left(\dot{\theta}+\omega_{\mathrm{E}}\right)(\cos \varphi \cos \psi \cdot\right.\\ &\sin \gamma-\sin \varphi \cos \gamma)]-\cos \sigma / \cos \beta[\dot{\gamma}-\dot{\varphi} \cos \psi-\\ &\left.\left(\dot{\theta}+\omega_{\mathrm{E}}\right) \cos \varphi \sin \psi\right], \end{aligned} $

(5)

$ \begin{aligned} \dot{\beta}=& p \sin \alpha-r \cos \alpha-\sin \sigma\left[\dot{\gamma}-\dot{\varphi} \cos \psi-\left(\dot{\theta}+\omega_{\mathrm{E}}\right)\right.\\ &\cos \varphi \sin \psi]+\cos \sigma[\dot{\psi} \cos \gamma-\dot{\varphi} \sin \psi \sin \gamma+\\ &\left.\left(\dot{\theta}+\omega_{\mathrm{E}}\right)(\cos \varphi \cos \psi \sin \gamma-\sin \varphi \cos \gamma)\right], \end{aligned} $

(6)

$ \begin{array}{l} \dot{\sigma}=p \cos \alpha \cos \beta+q \sin \beta+r \sin \alpha \cos \beta-\dot{\alpha} \sin \beta+\\ \ \ \ \ \ \ \ \dot{\psi} \sin \gamma+\dot{\varphi} \sin \psi \cos \gamma-\left(\dot{\theta}+\omega_{\mathrm{E}}\right)(\cos \varphi \cos \psi \cos \gamma+\\ \ \ \ \ \ \ \ \sin \varphi \sin \gamma) \text { , } \end{array} $

(7)

$ \dot{p}=\left(I_{y}-I_{z}\right) / I_{x} \cdot q r+M_{x} / I_{x} , $

(8)

$ \dot{q}=\left(I_{z}-I_{x}\right) / I_{y} \cdot p r+M_{y} / I_{y} , $

(9)

$ \dot{r}=\left(I_{x}-I_{y}\right) / I_{z} \cdot p q+M_{z} / I_{z}. $

(10)

式中：R为飞行器质心距地球球心的距离; θ、φ为飞行器所在经、纬度; γ、ψ分别为航迹倾角和航迹方位角; s_f为剩余待飞航程; α、β和σ分别为迎角、侧滑角和倾侧角(本文设右倾侧为正); p、q和r分别为滚转、俯仰和偏航角速率.式(1)~式(4)由再入轨迹状态方程^[11]经量纲一的转化并改变自变量为能量e之后得到，详细转化过程可参照文献[3]，其中：D=qSC_D/(mg₀)为归一化阻力，L=qSC_L/(mg₀)为归一化升力; q=0.5ρ(VV_c)²为动压; ρ为大气密度; V为归一化后的飞行速度; $ {V_c}{\rm{ = }}\sqrt {{g_0}{r_0}} $为速度归一化因子; g₀=9.8 m/s²、r₀=6 378 135 m分别为海平面的重力加速度和地球半径; S为机翼参考面积; m为飞行器质量.另外，I_i(i=x, y, z)为惯性矩; M_x、M_y和M_z分别为飞行器所受合外力矩在机体轴的3个分量; ω_E为地球自转角速率.

为方便设计非线性姿态控制律并考虑不确定/干扰和作动器失效，可将姿态系统方程式(5)~式(10)转化为如下仿射非线性方程组：

$ \left\{\begin{array}{l} \dot{\boldsymbol{\varOmega}}=\boldsymbol{f}_{\mathrm{s}}(\boldsymbol{\varOmega})+\boldsymbol{g}_{\mathrm{s}}(\boldsymbol{\varOmega}) \boldsymbol{\omega}+\boldsymbol{D}_{\mathrm{s}} , \\ \dot{\boldsymbol{\omega}}=\boldsymbol{f}_{\mathrm{f}}(\boldsymbol{\omega})+\boldsymbol{g}_{\mathrm{f}}(\boldsymbol{\omega}) \operatorname{sat}\left[\boldsymbol{M}_{\mathrm{RCS}}+(\boldsymbol{I}-\boldsymbol{E}) \boldsymbol{M}_{\mathrm{d}}\right]+\boldsymbol{D}_{\mathrm{f}}. \end{array}\right. $

(11)

式中：Ω =[α β σ]^T为气流姿态角向量; ω =[p q r]^T为角速度向量; f_s、f_f为系统矩阵; g_s、g_f为控制矩阵，具体表达式可由式(5)~式(10)推导得出，其中$ \dot \theta , \dot \varphi , \dot \gamma 和\dot \psi $须代入相应轨迹状态方程^[11]. $ {D_{\rm{s}}} = \Delta {f_{\rm{s}}} + \Delta {g_{\rm{s}}}\omega + {d_{\rm{s}}}和{D_{\rm{f}}} = \Delta {f_{\rm{f}}} + \Delta {g_{\rm{f}}}u + {d_{\rm{f}}} $分别为姿态角回路和角速率回路中的未知复合干扰，其中：Δ f_s、Δf_f、Δg_s和Δg_f为模型不确定项; d_s、d_f为外部干扰; sat(·)为饱和函数; M_RCS、M_d分别为反作用控制系统(RCS)与气动舵面产生的控制力矩，由M_d可反解出舵面的偏转指令.两者与非控制力矩的气动力矩项之和构成式(8)~式(10)中的合外力矩M =[M_x, M_y, M_z]^T. I为单位阵，E =diag([e₁, e₂, e₃])为舵面部分失效，其作用可合并到复合干扰D_f当中.

为保证RLV成功完成再入飞行任务，再入轨迹必须满足热流密度、过载和动压3个硬约束条件^[3]：

$ \dot{Q}=k_{\rm{Q}} \rho^{0.5} V^{3.15} \leqslant \dot{Q}_{\max } , $

(12)

$ n=\sqrt{L^{2}+D^{2}} \leqslant n_{\max } , $

(13)

$ \bar{q}=0.5 \rho V^{2} g_{0} r_{0} \leqslant \bar{q}_{\max }. $

(14)

式中：k_Q=91 089 918.35为模型系数; $ {\dot Q_{{\rm{max}}}} $为驻点处最大允许热流密度; n_max、q_max分别为最大允许过载与动压.

本文的主要目的为设计六自由度制导控制系统及鲁棒容错控制律，以引导HORUS飞行器从再入初始状态到末端能量管理段的飞行，在考虑姿态不确定/干扰、舵面部分失效及饱和的影响下实现飞行器鲁棒容错控制，并满足再入任务要求.

2 制导控制系统总体设计

RLV制导控制系统如图 1所示，该系统整合了制导环与姿控环.在六自由度仿真中，姿控系统的参考输入量Ω_r=[α_rβ_rσ_r]^T由制导环的改进预测校正制导算法在线生成.本文提出一种鲁棒容错控制方法来设计姿态控制器，用以跟踪制导环计算的指令.需要说明的是在制导系统当中，部分状态变量是量纲一的，故当实际飞行状态量传入到制导系统当中时，须经过归一化处理^[3].

本文的再入段包含再入起始段和准平衡滑翔段，飞行过程中考虑了姿态回路不确定/干扰因素的影响，也对舵面部分失效及饱和进行了容错控制设计.提出的ITDDO是容错控制方案的重要组成部分，它加在姿态角和角速率回路当中，结合反步法标称控制律在舵面不发生失效时用以估计姿态回路的不确定和干扰，部分失效发生时，对失效部分进行估计，同时也估计不确定及干扰，输出复合干扰估计$ {\hat D_{\rm{s}}}和{\hat D_{\rm{f}}} $，并计算姿态控制器的补偿量.针对舵面的饱和问题，设计辅助抗饱和系统来向容错控制器中引入误差状态量λ₁和λ₂.当舵面未发生饱和时，反馈状态量为0.由于预测校正制导算法本身具备一定的鲁棒性，且对初始误差不敏感，因此该联合制导控制系统理论上具备鲁棒容错性能.控制分配方案采用的是根据动压和马赫数的力矩分配.

3 六自由度制导控制系统设计 3.1 改进预测校正制导律设计

根据大气指数模型ρ=ρ₀e^ah3+bh2+ch+d(ρ₀、a、b、c和d是已知参数)可知，ρ仅是关于高度h(等于R减地球半径)的函数.因此，式(12)~式(14)能够进一步转换为关于高度h和速度V的不等式，根据此不等式可以确定RLV在高度-速度空间内的再入走廊下边界(如图 4虚线所示).

依据文献[5]提出的约束预测校正制导方法，RLV的迎角α按照事先确定的标称迎角剖面来变化，飞行器纵向运动采用倾侧角幅值|σ|来控制，侧向运动采用倾侧角符号sgn(σ)控制.设计中采用二次函数参数化模型^[5]来对倾侧角幅值剖面进行参数化设计.与常规线性参数化模型相比，该模型得到的再入轨迹更加平滑，再入过程中的热流密度和过载更小.另外，在QEGC条件中还引入哥氏加速度的影响，从而将高度-速度再入走廊转换为倾侧角幅值约束. α与σ代入式(1)~式(4)的制导方程，积分后得到RLV的预测落点.当预测落点与期望落点之间存在偏差时，返回重新校正初始设计的倾侧角幅值参数化模型，直至预测与期望落点间的偏差满足要求.一个周期内的再入飞行制导流程如图 2所示，其中s_f为飞行器的待飞航程.

然而，文献[5]中的侧向制导逻辑中将航向角偏差走廊分成多段来进行设计.此方法在编程实现时较为繁琐，本文采用S型函数来进行设计，将其中的再入走廊简化成如下所示的两个阶段：

$ \left\{\begin{array}{l} \left|\Delta \psi_{\mathrm{th}}\right|=\frac{9}{1+\mathrm{e}^{-60\left(s_{\mathrm{f}}-0.9\right)}}+14.35, s_{\mathrm{f}} \leqslant 0.95° ; \\ \left|\Delta \psi_{\mathrm{th}}\right|=\frac{8}{1+\mathrm{e}^{0.25\left(s_{\mathrm{f}}-20\right)}}+15, s_{\mathrm{f}}>0.95°. \end{array}\right. $

(15)

式中：|Δψ_th|为ψ与视场角(飞行器当前位置与目标位置间的夹角)之差的阈值.

3.2 改进跟踪微分干扰观测器的设计

为提高控制系统的鲁棒性并对舵面部分失效进行容错控制，本文针对不确定仿射非线性方程描述的MIMO系统，借鉴文献[15]提出了一种ITDDO，系统如下式所示，观测器形式如定理1所述.

$ \left\{\begin{array}{l} \dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})+\boldsymbol{g}(\boldsymbol{x}) \boldsymbol{u}+\boldsymbol{D}, \boldsymbol{x}(0)=\boldsymbol{x}_{0} ；\\ \boldsymbol{y}=\boldsymbol{x}. \end{array}\right. $

(16)

式中：x为可测状态量; u为系统控制量; x₀为状态初值; D为有界的未知复合干扰，包括模型不确定、外界干扰及舵面部分失效作用.

定理1  针对如式(16)所表达的MIMO系统，设计如下式的ITDDO.

$ \left\{\begin{array}{l} \dot{\boldsymbol{\hat{x}}}=\boldsymbol{f}(\boldsymbol{x})+\boldsymbol{g}(\boldsymbol{x}) \boldsymbol{u}+\hat{\boldsymbol{D}}, \hat{\boldsymbol{x}}(0)=\hat{x}_{0} ; \\ \dot{\boldsymbol{z}}=R^{2} f\left(\hat{\boldsymbol{x}}_{1}(t)-\boldsymbol{r}(t), \hat{\boldsymbol{D}}(t) / R\right), \\ \hat{\boldsymbol{D}}=\boldsymbol{z}+\boldsymbol{k}_{\mathrm{r}}\left[\boldsymbol{r}(t)-\boldsymbol{r}_{0}\right]+\hat{\boldsymbol{D}}_{0}, \boldsymbol{r}(0)=\boldsymbol{r}_{0}, \\ \ \ \ \ \ \ \ \ \hat{\boldsymbol{D}}(0)=\hat{\boldsymbol{D}}_{0} ; \\ \hat{\boldsymbol{x}}_{1}(t)=\int \hat{\boldsymbol{D}} \mathrm{d} t+\hat{\boldsymbol{x}}_{10}, \boldsymbol{\boldsymbol { x }}_{1}(0)=\hat{\boldsymbol{x}}_{10} ; \\ \boldsymbol{r}(t)=\boldsymbol{x}(t)-\boldsymbol{x}_{0}-\int[\boldsymbol{f}(\boldsymbol{x})+\boldsymbol{g}(\boldsymbol{x}) \boldsymbol{u}] \mathrm{d} t+\\ \ \ \ \ \ \ \ \ \boldsymbol{r}_{0}, \boldsymbol{x}(0)=\boldsymbol{x}_{0} . \end{array}\right. $

(17)

式中：第1式为状态估计方程; $ \hat x、\hat D $分别为对状态x及复合干扰D的估计; $ z、{\hat x_1}\left( t \right) $和r (t)为观测器中间状态; R为设计参数; k_r为设计参数矩阵; $ {x_0}、{\hat x_0}、{\hat x_{10}}、{D_0} $和r₀为初始值矩阵; f(·)为局部李普希兹连续函数且f(0, 0)=0，若采用如下函数：

$ \begin{aligned} f\left(z_{1}, z_{2}\right)=&-a_{1}\left[z_{1}+\left|z_{1}\right|^{p / q} \arctan \left(b_{1} \cdot z_{1}\right)\right]-\\ & a_{2}\left[z_{2}+\left|z_{2}\right|^{p / q} \arctan \left(b_{2} \cdot z_{2}\right)\right], \end{aligned} $

(18)

式中：a₁、a₂、b₁和b₂为大于零的调节参数，终端吸引子p和q设为正奇数，且0 < p/q < 1.那么，对于任意a>0, 当R→+∞时, $ \hat x $在t∈[a, +∞)上一致收敛于x, $ \hat D $在t∈[a, +∞)上一致收敛于D.

引理1^[16]  假设平衡点为(0, 0)的系统(19)全局渐进稳定，系统给定初始值为z₁(0)=z₁₀和z₂(0)=z₂₀.

$ \left\{\begin{array}{l} \dot{z}_{1}(t)=z_{2}(t), z_{1}(0)=z_{10} ; \\ \dot{z}_{2}(t)=f\left(z_{1}(t), z_{2}(t)\right), z_{2}(0)=z_{20}. \end{array}\right. $

(19)

式中f(·)为局部李普希兹连续函数且f(0, 0)=0.若参考输入r(t)可微且$ \mathop {{\rm{sup}}}\limits_{t \in [0, + \infty )} \left| {\dot r\left( t \right)} \right| < + \infty $，则系统：

$ \left\{\begin{array}{l} \dot{x}_{1 \mathrm{R}}(t)=x_{2 \mathrm{R}}(t), x_{1 \mathrm{R}}(0)=x_{10}, x_{2 \mathrm{R}}(0)=x_{20}; \\ \dot{x}_{2 \mathrm{R}}(t)=R^{2} f\left(x_{1 \mathrm{R}}(t)-r(t), x_{2 \mathrm{R}}(t) / R\right)+\bar{\alpha} \cdot \dot{r}(t) . \end{array}\right. $

(20)

的解在以下意义上是收敛的：对于任意a>0, 当R→+∞时，x_1R在t∈[a, +∞)上一致收敛于r(t).其中：x₁₀、x₂₀为任意初值，α为大于零的常数.

证明  将式(17)中的$ \dot {\hat x} = f\left( x \right) + g\left( x \right)u + \hat D $转化为

$ \dot{\hat{\boldsymbol{x}}}_{1}=\hat{\boldsymbol{D}}, $

(21)

其中

$ \dot{\hat{\boldsymbol{x}}}_{1}=\dot{\hat{\boldsymbol{x}}}-\bar{\boldsymbol{F}}(\boldsymbol{x}), $

(22)

此处 F(x)= f (x)+ g (x) u，则$ {\hat x_1} $参考输入的导数为

$ \dot{\boldsymbol{r}}(t)=\dot{\boldsymbol{x}}-\overline{\boldsymbol{F}}(\boldsymbol{x}), $

(23)

由式(21)、(22)可知$ \left[ {{{\hat x}_1} - r\left( t \right)} \right]' = \left( {\hat x - x} \right)' $，但并不表示$ {\hat x_1} - r\left( t \right) = \hat x - x $.于是，根据引理1可构造形如式(20)表达的跟踪微分器如下:

$ \left\{\begin{array}{l} \dot{\hat{\boldsymbol{x}}}_{1}=\hat{\boldsymbol{D}}, \hat{\boldsymbol{x}}_{1}(0)=\hat{\boldsymbol{x}}_{10}, \hat{\boldsymbol{D}}(0)=\hat{\boldsymbol{D}}_{0} ;\\ \dot{\hat{\boldsymbol{D}}}=R^{2} f\left(\hat{\boldsymbol{x}}_{1}-\boldsymbol{r}(t), \frac{\hat{\boldsymbol{D}}}{R}\right)+\boldsymbol{k}_{\mathrm{r}} \cdot \dot{\boldsymbol{r}}(t). \end{array}\right. $

(24)

由式(22)可知，式(23)中出现了状态x的微分，这在实际控制器设计中应尽量避免，因此可对式(24)的第2式两边进行积分，得到式(17)中的第2式与第3式，其中$ {\hat x_1} $和r (t)的计算可由式(21)与式(23)积分求得，结果为式(17)的第4式与第5式.

如f(·)取为式(18)，则根据文献[17]定理1的证明，式(24)对应的形如式(19)的系统全局渐进稳定.因此，根据引理1，对于任意a>0, 当R→+∞时，式(24)中的$ {\hat x_1} $在t∈[a, +∞)上一致收敛于r，即$ {\dot {\hat x}_1} $一致收敛于ṙ.根据式(22)和式(23)有$ \dot {\hat x} $一致收敛于$ \hat x $.如设两者状态初值$ {\hat x_0} = {x_0} $，则可得$ \hat x $一致收敛于x.根据式(16)、(17)的第1式，有$ \hat D $一致收敛于D.

定理1得证.

注1  ITDDO在反正切函数中添加了b₁和b₂，用于调节函数变化的速率，扩大了函数的调整范围.而且，ITDDO的推导过程考虑了F(x)的影响，同时借鉴文献[16]，在设计中引入前馈微分项k_r·ṙ(t)，进一步提高了干扰估计的精度.

根据姿态方程式(11)和设计ITDDO的式(17)、(18)，有姿态角和角速率回路的ITDDO形式如下：

$ \left\{\begin{array}{l} \dot{\boldsymbol{z}}_{\mathrm{s}}=R_{\mathrm{s}}^{2}\left\{-a_{1}\left[\hat{\boldsymbol{e}}_{\mathrm{ds}}+\left|\hat{\boldsymbol{e}}_{\mathrm{ds}}\right|^{p_{1} / q_{1}} \arctan \left(b_{11} \hat{\boldsymbol{e}}_{\mathrm{ds}}\right)\right]-a_{2} \cdot\right. \\ \ \ \ \ \ \ \ \ \left.\left[\hat{\boldsymbol{D}}_{\mathrm{s}} / R_{\mathrm{s}}+\left|\hat{\boldsymbol{D}}_{\mathrm{s}} / R_{\mathrm{s}}\right|^{p_{1} / q_{1}} \arctan \left(b_{12} \hat{\boldsymbol{D}}_{\mathrm{s}} / R_{\mathrm{s}}\right)\right]\right\} , \\ \hat{\boldsymbol{D}}_{\mathrm{s}}=\boldsymbol{z}_{\mathrm{s}}+\boldsymbol{k}_{\mathrm{rs}}\left\{\boldsymbol{\varOmega}(t)-\boldsymbol{\varOmega}_{0}-\int\left[\boldsymbol{f}_{\mathrm{s}}(\boldsymbol{\varOmega})+\boldsymbol{g}_{\mathrm{s}}(\boldsymbol{\varOmega}) \boldsymbol{\omega}\right] \mathrm{d} t\right\}, \end{array}\right. $

(25)

$ \left\{\begin{array}{ll} \dot{\boldsymbol{z}}_{\mathrm{f}}= R_{\mathrm{f}}^{2}\left\{-a_{3}\left[\hat{\boldsymbol{e}}_{\mathrm{df}}+\left|\hat{\boldsymbol{e}}_{\mathrm{df}}\right|^{p_{2} / q_{2}} \arctan \left(b_{21} \hat{\boldsymbol{e}}_{\mathrm{df}}\right)\right]-a_{4} \cdot\right. \\ \ \ \ \ \ \ \ \ \left.\left[\hat{\boldsymbol{D}}_{\mathrm{f}} / R_{\mathrm{f}}+\left|\hat{\boldsymbol{D}}_{\mathrm{f}} / R_{\mathrm{f}}\right|^{p_{2} / q_{2}} \arctan \left(b_{22} \hat{\boldsymbol{D}}_{\mathrm{f}} / R_{\mathrm{f}}\right)\right]\right\} , \\ \hat{\boldsymbol{D}}_{\mathrm{f}}= \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{k}_{\mathrm{rf}}\left\{\boldsymbol{\omega}(t)-\int\left[\boldsymbol{f}_{\mathrm{f}}(\boldsymbol{\omega})+\boldsymbol{g}_{\mathrm{f}}(\boldsymbol{\omega}) \boldsymbol{M}_{\mathrm{ctr}}\right] \mathrm{d} t\right\}. \end{array}\right. $

(26)

式中：ê_ds和ê_df为式(17)中的$ {\hat x_1}\left( t \right) - r\left( t \right) $，可根据式(11)和式(17)对应计算; R_i(i=s, f)为待设计的增益量; b_1i和b_2i(i=1, 2)为反正切函数速度调节因子; 终端吸引子p_i和q_i(i=1, 2)设为正奇数; a_i(i=1, 2, 3, 4)为调节因子; k_rs、k_rf为设计参数矩阵.通常，ITDDO在飞行器再入姿态稳定状态下切入，故角速度初值ω₀设为零，令干扰估计初值$ {\hat D_{{\rm{s0}}}} = {\hat D_{{\rm{f0}}}} = 0 $.

3.3 鲁棒容错姿态控制器设计

姿态系统的非线性控制律基于反步法进行设计.定义姿态角跟踪误差和角速率误差为：

$ \boldsymbol{e}_{\mathrm{s}}=\boldsymbol{\varOmega}-\boldsymbol{\varOmega}_{\mathrm{r}}, $

(27)

$ \boldsymbol{e}_{\mathrm{f}}=\boldsymbol{\omega}-\boldsymbol{\omega}_{\mathrm{r}}. $

(28)

式中：e_s=[e_αe_βe_σ]^T，Ω_r=[α_r β_r σ_r]^T为制导系统计算出的参考姿态角指令; ω_r=[p_rq_r r_r]^T是姿态角回路得到的虚拟控制量.将式(27)微分后得到$ {\dot e_{\rm{s}}} = \dot \Omega - {\dot \Omega _{\rm{r}}} $，代入式(11)中的$ \dot \Omega $后可得

$ \dot{\boldsymbol{e}}_{\mathrm{s}}=\boldsymbol{f}_{\mathrm{s}}+\boldsymbol{g}_{\mathrm{s}} \boldsymbol{\omega}+\boldsymbol{D}_{\mathrm{s}}-\dot{\boldsymbol{\varOmega}}_{\mathrm{r}}. $

(29)

根据式(29)设计姿态角回路控制律ω_r如下：

$ \boldsymbol{\omega}_{\mathrm{r}}=\boldsymbol{g}_{\mathrm{s}}^{-1}\left(-\boldsymbol{K}_{\mathrm{s}} \boldsymbol{e}_{\mathrm{s}}-\boldsymbol{f}_{\mathrm{s}}+\dot{\boldsymbol{\varOmega}}_{\mathrm{r}}-\hat{\boldsymbol{D}}_{\mathrm{s}}\right). $

(30)

式中：K_s>0为待设计的3×3对角方阵，$ {\hat D_{\rm{s}}} $为姿态角回路对复合干扰D_s的估计.

在不考虑舵面饱和的情况下，可将式(11)中的舵面失效合并至复合干扰当中，即D_f2= D_f-g_fE· M_d, 于是角速率回路改写成如下形式：

$ \dot{\boldsymbol{\omega}}=\boldsymbol{f}_{\mathrm{f}}+\boldsymbol{g}_{\mathrm{f}} \boldsymbol{M}_{\mathrm{r}}+\boldsymbol{D}_{\mathrm{f} 2}, $

(31)

式中M_r=M_RCS+ M_d，这里将舵面部分失效时未能提供的力矩视为角速率回路中不确定项的组成部分.在不考虑饱和时，将式(28)微分后可得

$ \dot{\boldsymbol{e}}_{\mathrm{f}}=\dot{\boldsymbol{\omega}}-\dot{\boldsymbol{\omega}}_{\mathrm{r}}=\boldsymbol{f}_{\mathrm{f}}+\boldsymbol{g}_{\mathrm{f}} \boldsymbol{M}_{\mathrm{r}}+\boldsymbol{D}_{\mathrm{f} 2}-\dot{\boldsymbol{\omega}}_{\mathrm{r}}. $

(32)

设计控制力矩指令M_r如下^[18]:

$ \boldsymbol{M}_{\mathrm{r}}=\boldsymbol{g}_{\mathrm{f}}^{-1}\left(-\boldsymbol{K}_{\mathrm{f}} \boldsymbol{e}_{\mathrm{f}}-\boldsymbol{f}_{\mathrm{f}}-\hat{\boldsymbol{D}}_{\mathrm{f} 2}+\dot{\boldsymbol{\omega}}_{\mathrm{r}}-\boldsymbol{g}_{\mathrm{s}}^{-1} \boldsymbol{e}_{\mathrm{s}}\right). $

(33)

为进一步加强控制器的鲁棒性能，向式(33)中引入改进sigmoid函数构成的饱和项：

$ \left\{\begin{array}{l} \operatorname{sat}\left(x_{i}\right)=\frac{2}{1+\mathrm{e}^{-3 \boldsymbol{x}_{i}}}-1, i=1, 2, 3 ;\\ \operatorname{sat}(\boldsymbol{x})=\left[\operatorname{sat}\left(x_{1}\right), \operatorname{sat}\left(x_{2}\right), \operatorname{sat}\left(x_{3}\right)\right]^{\mathrm{T}}; \\ \boldsymbol{x}=-\boldsymbol{e}_{\mathrm{f}} / \varepsilon. \end{array}\right. $

(34)

式中：e为自然常数，ε>0为调节因子，因此可得实际控制量为

$ \begin{aligned} \boldsymbol{M}_{\mathrm{r}}=& \boldsymbol{g}_{\mathrm{f}}^{-1}\left(-\boldsymbol{K}_{\mathrm{f}} \boldsymbol{e}_{\mathrm{f}}-\boldsymbol{f}_{\mathrm{f}}-\hat{\boldsymbol{D}}_{\mathrm{f} 2}+\dot{\boldsymbol{\omega}}_{\mathrm{r}}-\boldsymbol{g}_{\mathrm{s}}^{-1} \boldsymbol{e}_{\mathrm{s}}+\right.\\ &\left.\boldsymbol{K}_{c} \operatorname{sat}(\boldsymbol{x})\right), \end{aligned} $

(35)

式中$ {\dot \omega _{\rm{r}}} $的表达式需要代入式(30).

针对RLV再入过程中的舵面饱和问题，假设饱和前、后的力矩差值为Δ δ =sat(M)- M_r，且‖Δδ ‖≤ M，‖ · ‖为L₂范数，M为正常数，则角速率回路方程式(31)中的M_r应由实际力矩sat(M)代替.本文借鉴文献[19]设计了一种抗舵面饱和辅助系统:

$ \left\{\begin{array}{l} \dot{\boldsymbol{\lambda}}_{1}=-\boldsymbol{c}_{1} \boldsymbol{\lambda}_{1}+\boldsymbol{g}_{\mathrm{s}}(\boldsymbol{\varOmega}) \boldsymbol{\lambda}_{2} , \\ \dot{\boldsymbol{\lambda}}_{2}=-\boldsymbol{c}_{2} \boldsymbol{\lambda}_{2}+\boldsymbol{g}_{\mathrm{f}}(\boldsymbol{\omega}) \Delta \boldsymbol{\delta}. \end{array}\right. $

(36)

式中：c₁和c₂为待整定的对角矩阵，λ₁和λ₂为辅助系统中的状态向量.因此，当前系统的姿态角和角速率误差为

$ \left\{\begin{array}{l} \boldsymbol{z}_{\mathrm{s}}=\boldsymbol{e}_{\mathrm{s}}-\boldsymbol{\lambda}_{1}, \\ \boldsymbol{z}_{\mathrm{f}}=\boldsymbol{e}_{\mathrm{f}}-\boldsymbol{\lambda}_{2}. \end{array}\right. $

(37)

将上式取微分，代入式(28)、(29)、(32)、(36)和(37)，可得:

$ \dot{\boldsymbol{z}}_{\mathrm{s}}=\boldsymbol{f}_{\mathrm{s}}+\boldsymbol{g}_{\mathrm{s}}\left(\boldsymbol{z}_{\mathrm{f}}+\boldsymbol{\omega}_{\mathrm{r}}\right)+\boldsymbol{D}_{\mathrm{s}}-\dot{\boldsymbol{\varOmega}}_{\mathrm{r}}+\boldsymbol{c}_{1} \boldsymbol{\lambda}_{1} , $

(38a)

$ \dot{\boldsymbol{z}}_{\mathrm{f}}=\boldsymbol{f}_{\mathrm{f}}+\boldsymbol{g}_{\mathrm{f}} \boldsymbol{M}_{\mathrm{r}}+\boldsymbol{D}_{\mathrm{f} 2}-\dot{\boldsymbol{\omega}}_{\mathrm{r}}+\boldsymbol{c}_{2} \boldsymbol{\lambda}_{2}, $

(38b)

注意式(32)中的M_r已由sat(M)代替.

设计姿态角控制律和角速率控制律如下：

$ \boldsymbol{\omega}_{\mathrm{r}}= \boldsymbol{g}_{\mathrm{s}}^{-1}\left(-\boldsymbol{K}_{\mathrm{s}} z_{\mathrm{s}}-\boldsymbol{f}_{\mathrm{s}}+\dot{\boldsymbol{\varOmega}}_{\mathrm{r}}-\hat{\boldsymbol{D}}_{\mathrm{s}}-\boldsymbol{c}_{1} \boldsymbol{\lambda}_{1}\right), $

(39a)

$ \begin{aligned} \boldsymbol{M}_{\mathrm{r}}=& \boldsymbol{g}_{\mathrm{f}}^{-1}\left(-\boldsymbol{K}_{\mathrm{f}} \boldsymbol{z}_{\mathrm{f}}-\boldsymbol{f}_{\mathrm{f}}-\hat{\boldsymbol{D}}_{\mathrm{f} 2}+\dot{\boldsymbol{\omega}}_{\mathrm{r}}-\boldsymbol{g}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{z}_{\mathrm{s}}-\right.\\ &\left.\boldsymbol{c}_{2} \boldsymbol{\lambda}_{2}+\boldsymbol{K}_{\mathrm{c}} \operatorname{sat}\left(-\boldsymbol{z}_{\mathrm{f}} / \varepsilon\right)\right). \end{aligned} $

(39b)

式中：$ {\hat D_{\rm{s}}}和{\hat D_{f2}} $由改进跟踪微分干扰观测器的设计中的式(25)、(26)计算得出，λ₁、λ₂、z_s和z_f由式(36)~(38)得出.

3.4 稳定性分析

定理2  针对RLV的姿态控制系统，采用式(17)的ITDDO、式(36)的辅助抗饱和系统和式(39)的姿态控制律，闭环系统局部一致渐进有界.

证明  定义ITDDO的估计误差为$ e = \hat D - D $，由定理1可知，存在有界常数χ∈R⁺，使得‖ e ‖≤χ.

将式(39a)代入式(38a)，可得

$ \dot{\boldsymbol{z}}_{\mathrm{s}}=-\boldsymbol{K}_{\mathrm{s}} \boldsymbol{z}_{\mathrm{s}}+\boldsymbol{g}_{\mathrm{s}} \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{D}_{\mathrm{s}}-\hat{\boldsymbol{D}}_{\mathrm{s}}. $

(40)

姿态角回路稳定性证明的李雅普诺夫函数设为V₁=0.5z_s^Tz_s，因此:

$ \dot{V}_{1}=\boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}} \dot{\boldsymbol{z}}_{\mathrm{s}}=\boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}}\left(-\boldsymbol{K}_{\mathrm{s}} \boldsymbol{z}_{\mathrm{s}}+\boldsymbol{g}_{\mathrm{s}} \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{D}_{\mathrm{s}}-\hat{\boldsymbol{D}}_{\mathrm{s}}\right). $

(41)

令$ {e_{{\rm{DOS}}}} = {D_{\rm{s}}} - {\hat D_{\rm{s}}} $.由于在定理1中证明了$ \hat D $一致收敛于D，即干扰的估计误差有界，因此‖ e_DOS‖≤ χ_s，χ_s为有界正常数.

定义角速率回路李雅普诺夫函数为

$ V_{2}=0.5 \boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{z}_{\mathrm{f}}+V_{1}, $

(42)

将式(42)微分并代入式(38a)和式(41)可得

$ \begin{aligned} \dot{V}_{2}=&-\boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{s}} \boldsymbol{z}_{\mathrm{s}}-\boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{f}} \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}}\left(\boldsymbol{D}_{\mathrm{f}}-\hat{\boldsymbol{D}}_{\mathrm{f}}\right)+\\ & \boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOS}}+\boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{c}} \operatorname{sat}\left(-\boldsymbol{z}_{\mathrm{f}} / \boldsymbol{\varepsilon}\right). \end{aligned} $

(43)

令$ {e_{{\rm{DOS}}}} = {D_{\rm{f}}} - {\hat D_{\rm{f}}} $，‖ e_DOF‖≤χ_f，χ_f为有界正常数.根据不等式定理z₁^Tz₂≤ z₁^Tz₁/2+ z₂^Tz₂/2以及干扰估计误差有界，从式(43)可以得出:

$ \begin{aligned} \dot{V}_{2} \leqslant &-\boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{s}}-0.5 \boldsymbol{I}\right) \boldsymbol{z}_{\mathrm{s}}+0.5 \boldsymbol{e}_{\mathrm{DOS}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOS}}+0.5 \boldsymbol{e}_{\mathrm{DOF}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOF}}-\\ & \boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{f}}-0.5 \boldsymbol{I}\right) \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{c}} \operatorname{sat}\left(-\boldsymbol{z}_{\mathrm{f}} / \boldsymbol{\varepsilon}\right) \leqslant \\ &-\boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{s}}-0.5 \boldsymbol{I}\right) \boldsymbol{z}_{\mathrm{s}}+0.5 \bar{\chi}_{\mathrm{s}}^{2}+0.5 \bar{\chi}_{\mathrm{f}}^{2}-\\ & \boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{f}}-0.5 \boldsymbol{I}\right) \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{c}} \operatorname{sat}\left(-\boldsymbol{z}_{\mathrm{f}} / \boldsymbol{\varepsilon}\right), \end{aligned} $

(44)

式中, I为单位方阵.选取合适的K_s、K_f、K_c和ε能够保证$ {\dot V_2} $ < 0.于是，系统Lyaponov函数设为

$ V=V_{2}+\frac{1}{2} \boldsymbol{e}_{\mathrm{DOS}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOS}}+\frac{1}{2} \boldsymbol{e}_{\mathrm{DOF}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOF}}, $

(45)

可知，V₂≤V，将式(45)微分可得

$ \dot{V}=\dot{V}_{2}+\boldsymbol{e}_{\mathrm{DOS}}^{\mathrm{T}} \dot{\boldsymbol{e}}_{\mathrm{DOS}}+\boldsymbol{e}_{\mathrm{DOF}}^{\mathrm{T}} \cdot \dot{\boldsymbol{e}}_{\mathrm{DOF}}. $

(46)

将式(41)展开可得

$ \begin{aligned} \dot{V} \leqslant &-\boldsymbol{z}_{\mathrm{s}}^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{s}}-0.5 \boldsymbol{I}\right) \boldsymbol{z}_{\mathrm{s}}+0.5 \boldsymbol{e}_{\mathrm{DOS}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOS}}+0.5 \boldsymbol{e}_{\mathrm{DOF}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOF}}-\\ & \boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}}\left(\boldsymbol{K}_{\mathrm{f}}-0.5 \boldsymbol{I}\right) \boldsymbol{z}_{\mathrm{f}}+\boldsymbol{z}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{K}_{\mathrm{c}} \operatorname{sat}\left(-\boldsymbol{z}_{\mathrm{f}} / \varepsilon\right) \boldsymbol{e}_{\mathrm{DOS}}^{\mathrm{T}} \dot{\boldsymbol{e}}_{\mathrm{DOS}}+\boldsymbol{e}_{\mathrm{DOF}}^{\mathrm{T}} \boldsymbol{e}_{\mathrm{DOF}} \end{aligned}. $

(47)

由于干扰的估计误差及其一阶导数有界，所以e_DOS^Tė_DOS和e_DOF^T ė_DOF有界，设这两项上界之和为β. e _DOS^Te_DOS/2和e_DOF^T e_DOF/2的上界分别为0.5χ_s²和0.5χ_f²，设这些上界之和为c，又因z_f^T K_csat(- z_f/ε)≤ 0，故

$ \dot{V} \leqslant-\kappa V+c, $

(48)

式中$ k = {\rm{min}}\left( \begin{array}{l} 2{\lambda _{{\rm{min}}}}\left( {{K_{\rm{s}}} - 0.5I} \right)\\ 2{\lambda _{\min }}\left( {{K_{\rm{f}}} - 0.5I} \right) \end{array} \right) $.

式(48)的两边同乘以e^-κt，且在[0, t]上积分可得

$ V(t) \leqslant\left[V(0)-\frac{c}{\kappa}\right] e^{-\kappa t}+\frac{c}{\kappa}, $

(49)

因此

$ \begin{array}{c} \frac{1}{2}\left\|\boldsymbol{\varOmega}-\boldsymbol{\varOmega}_{\mathrm{r}}-\boldsymbol{\lambda}_{1}\right\|^{2}=\frac{1}{2}\left\|\boldsymbol{z}_{\mathrm{s}}\right\|^{2} \leqslant V(t) \leqslant \\ {\left[V(0)-\frac{c}{\kappa}\right] e^{-\kappa t}+\frac{c}{\kappa}}, \end{array} $

(50)

当系数κ>0时，下式成立:

$ \begin{aligned} \left\|\boldsymbol{\varOmega}-\boldsymbol{\varOmega}_{\mathrm{r}}-\boldsymbol{\lambda}_{1}\right\| \leqslant & \sqrt{2\left[V(0)-\frac{c}{\kappa}\right] e^{-\kappa t}+2 \frac{c}{\kappa}} \leqslant \\ & \sqrt{2 V(0)+2 \frac{c}{\kappa}}. \end{aligned} $

(51)

假设辅助系统状态向量的范数‖ λ₁‖有界，则构造如下Lyapunov函数:

$ V_{\lambda}=\frac{1}{2} \boldsymbol{\lambda}_{1}^{\mathrm{T}} \boldsymbol{\lambda}_{1}+\frac{1}{2} \boldsymbol{\lambda}_{2}^{\mathrm{T}} \boldsymbol{\lambda}_{2}, $

(52)

式(52)微分后得到：

$ \begin{aligned} \dot{V}_{\lambda} \leqslant& -\boldsymbol{\lambda}_{1}^{\mathrm{T}}\left[\boldsymbol{c}_{1}-0.5 \boldsymbol{g}_{\mathrm{s}}\right] \boldsymbol{\lambda}_{1}-\boldsymbol{\lambda}_{2}^{\mathrm{T}}\left[\boldsymbol{c}_{2}-0.5 \boldsymbol{g}_{\mathrm{s}}-0.5 \boldsymbol{g}_{\mathrm{f}}\right] \boldsymbol{\lambda}_{2}+\\ & 0.5 \lambda_{\max }\left(\boldsymbol{g}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{g}_{\mathrm{f}}\right)\|\Delta \delta\|^{2} \leqslant-\kappa_{\lambda} V_{\lambda}+c_{\lambda}, \end{aligned} $

(53)

其中

$ \left\{\begin{array}{l} \kappa_{\lambda}=\min \left(\begin{array}{l} 2 \lambda_{\min }\left(\boldsymbol{c}_{1}-0.5 \boldsymbol{g}_{\mathrm{s}}\right) \\ 2 \lambda_{\min }\left(\boldsymbol{c}_{2}-0.5 \boldsymbol{g}_{\mathrm{s}}-0.5 \boldsymbol{g}_{\mathrm{f}}\right) \end{array}\right), \\ c_{\lambda}=0.5 \lambda_{\max }\left(\boldsymbol{g}_{\mathrm{f}}^{\mathrm{T}} \boldsymbol{g}_{\mathrm{f}}\right)\|\Delta \delta\|^{2} . \end{array}\right. $

同理可得$ \left\| {{\lambda _1}} \right\| \le \sqrt {2{V_\lambda }\left( 0 \right) + 2{c_\lambda }/{k_\lambda }} $，故式(51)可转换为

$ \left\|\boldsymbol{\varOmega}-\boldsymbol{\varOmega}_{\mathrm{r}}\right\| \leqslant \sqrt{2 V(0)+2 \frac{c}{\kappa}}+\sqrt{2 V_{\lambda}(0)+2 \frac{c_{\lambda}}{\kappa_{\lambda}}}. $

(54)

因此，通过选取合适的系数值，可以使误差收敛到任意小.定理2证毕.

4 仿真验证

研究对象HORUS-2B有左右升降副翼舵δ₁和δ₂，将其偏转区间设为较小的范围[-18°, 18°]以验证抗饱和控制算法，方向舵δ_r的偏转区间为[-40°, 40°]，用于平衡的体襟翼偏转角设为零.m=26 029 kg，$ {\dot Q_{{\rm{max}}}} $ =450 kW/m²，n_max=2.5，q_max=12.8 kPa.再入初始状态和末端状态见表 1.大气密度模型参数为ρ₀=1.225 kg/m³、a=1.394 5×10^-5、b=2.322×10^-3、c=-0.257 1、d=1.758 2.

姿态角和角速率的初值分别为Ω₀=[40°,0,0]^T，ω₀=0,,0]^T，抗饱和辅助系统的对角矩阵设为c₁=diag([0.3, 0.3, 0])和c₂=diag([15,15, 0]).由于再入飞行中主要由偏航RCS来控制偏航方向，方向舵作用时间短，故未考虑其饱和，即$ {\dot \lambda _1}\left( 3 \right) = {\dot \lambda _2}\left( 3 \right) = 0 $.姿态角和角速率ITDDO的设计参数为：p=1、q=15、R_s=10、a₁=1.8、a₂=0.1、b₁₁=b₁₂=10、R_f=8、a₃=1、a₄=1和b₂₁=b₂₂=10，k_rf=diag(9, 9, 33)，k_rs=diag(5, 5, 5).式(39)中的K_s、K_f和K _c的取值分别为0.8、1.5和0.6的对角方阵，ε设为0.05.

首先，仅针对姿控系统实验了ITDDO对于复合干扰的学习能力，仿真时间为100 s，采样周期为20 ms.施加复合干扰模拟函数如下：

$ \left\{\begin{array}{ll} D_{\mathrm{s1}}=0.005 \sin (t+1) \cos (2 t), & t>5 ; \\ D_{\mathrm{s2}}=0.003 \cos (t-1) \sin (2 t+2), & t>5 ; \\ D_{\mathrm{s3}}=0.005 \sin (t+1) \sin (2 t), & t>5 \end{array}\right. $

(55)

$ \left\{\begin{array}{l} D_{\mathrm{f1}}=0.02 \sin (t+1), \quad t>5 ;\\ D_{\mathrm{f2}}=0.01 \cos (2 t+2), \ t>5 ;\\ D_{\mathrm{f3}}=0.005 \sin (t+1), \ t>5. \end{array}\right. $

(56)

如图 3所示，ITDDO与文献[17]的方法以及非线性干扰观测器(NDO)进行了比较.由于姿态角回路处在角速率回路的外环，其干扰学习效果更优，故仅针对角速率回路通过叠加所有采样点的误差平方和来进行3种方法的比较.可以看出，ITDDO三通道干扰估计的误差平方和最小，同时3个姿态角的控制误差平方和也达到了最小.

然后，将ITDDO结合定理2中的方法，应用到了制导与控制系统联合仿真当中.复合干扰仍按式(55)、(56)的形式.另外，在400~850 s时，设左升降副翼舵发生部分失效，失效因子为0.3，850~1 000 s时，右升降副翼舵发生部分失效，失效因子为0.2.为验证方法的有效性，分4种情景进行仿真：case1包含ITDDO、反步法和抗饱和辅助系统; case2包含ITDDO、改进的反步法和抗饱和辅助系统; case3包含NDO、改进反步法和抗饱和辅助系统; case4只包含了ITDDO和改进反步法.图 4、5显示出飞行器质点的位置变化，可以看出所有情景下飞行器都能满足制导控制要求，使终端落点到达允许的误差圆范围内(半径10 km)，同时也满足再入过程约束的要求(见式(12)~(14)).然而，case2和case3的精度略高于case1.

图 6、7为迎角和倾侧角的跟踪效果.再入过程中控制侧滑角为零，因篇幅限制，故略去该图.虽然这3种情景中产生的参考姿态角因在线预测制导(虚线)而不完全重合，但是此细微差别并不影响实验的结果，故只选择了一条轨迹作为参照.从图 6、7可知，当倾侧角每次发生翻转时，3种情景中的姿态角存在较大震荡.为了更清楚地描述这种情况，图 6中分别附有700~750 s和1 255~1 270 s的局部放大图.从放大图中可知，在720 s处，攻角的波动幅度大幅增加，这是因为在倾侧角翻转的过程中需要更多的控制力矩来保证系统的平稳运行.此时，由于舵面发生饱和及角速率回路输出限幅这两个原因，无法提供所需控制力矩而导致姿态角暂时偏离参考值.改进反步法的姿态角跟踪精度有了显著的提升，从图 7的1 095~1 100 s的局部放大区间也可以看出改进的反步法能够使跟踪性能得到显著的增强.而case4中的辅助抗饱和系统对控制效果的改善作用也较为明显.

在图 6、7中，每当倾侧角翻转前后，姿态角的跟踪就会发生震荡，它们出现的时间点和控制舵面发生饱和的时间点一致.由于左升降副翼舵和右升降副翼舵在670~750 s和1 080~1 140 s这两个区间发生比较严重的饱和，故图 8、9给出了这两个区间舵面的偏转情况.

从图 8、9可知，辅助抗饱和系统可以使舵面尽快离开饱和状态，这样有利于改善舵面饱和时引起的控制系统性能下降或系统无法运行的情况.由于篇幅限制，本文只显示了攻角的干扰估计和滚转通道的干扰估计效果.从图 10可以看出，在制导与控制的联合系统中，case2中的ITDDO比case3中的普通NDO具有更好的复合干扰估计效果.从图 11可知，ITDDO也具有更好的响应能力.

5 结论

1) 本文针对控制系统不确定/干扰和舵面部分失效问题，提出一种改进跟踪微分干扰观测器ITDDO来估计姿态系统复合干扰和舵面失效作用.该观测器估计效果良好，且设计参数较少，能够满足姿控系统实时性的需要.

2) 针对舵面饱和控制问题，在基于ITDDO的反步控制律之上结合辅助抗饱和系统来解决问题，应用改进sigmoid函数来改善角速率回路中反步法的控制性能.

3) 本文开展了改进预测校正制导回路和姿态控制回路的六自由度联合仿真验证工作.由仿真结果可知，提出的方法能够较好解决再入飞行中的姿态不确定/干扰、舵面部分失效及饱和控制问题，能够引导RLV在满足再入约束的条件下到达终端落点区域，系统具备良好的鲁棒容错控制性能.