CN104880948A - 一种基于鲁棒自适应的带有飞轮的航天器容错控制方法 - Google Patents

Abstract

一种基于鲁棒自适应的带有飞轮的航天器容错控制方法，该方法有五大步骤：步骤一：航天器运动学方程的建立；步骤二：航天器动力学方程的建立；步骤三：期望的角速度轨迹设计；步骤四：控制器输出力矩的设计；步骤五：数值仿真。本发明提出在不采用FDD装置的情况下，在线估计出三轴力矩故障因子并设计鲁棒自适应容错控制器。这一研究旨在丰富航天器PFTCS方法，并为将来的航天器姿态控制提供技术支持。

Description

一种基于鲁棒自适应的带有飞轮的航天器容错控制方法

技术领域

本发明涉及航天器姿态控制技术领域。具体地说，是涉及一种基于鲁棒自适应的带有飞轮的航天器容错控制方法。利用此方法能保证当航天器存在外部扰动，同时作为执行机构的反作用飞轮出现力矩损失的故障情况下，所设计的容错控制方法仍然可以有效地实现航天器的姿态控制。

背景技术

容错控制是指当控制系统中的某些部件发生故障时，系统仍能按期望的性能指标或性能指标略有降低(但可接受)的情况下，还能安全地完成控制任务。容错控制的研究，使得提高复杂系统的安全性和可靠性成为可能。

容错控制系统(Fault Tolerant Control System，FTCS)可以分为两类：主动容错控制系统(Active Fault Tolerant Control System，AFTCS)和被动容错控制系统(PassiveFault Tolerant Control System，PFTCS)。AFTCS通过利用故障检测与诊断系统(FaultDetection And Diagnosis，FDD)的实时信息来重构控制器以处理系统部件故障。作为AFTCS的一个子系统，FDD必须对故障、未建模动态特性、动力学参数不确定性和其它扰动有很高的敏感度。重构控制的多数研究都是基于FDD可以提供准确和及时的故障信息的前提下进行的。因此对于AFTC，FDD的微小误差不仅可能导致控制性能的下降，甚至导致整个系统的失稳。

PFTCS利用鲁棒控制方法确保在传感器、执行机构故障情况下闭环系统的稳定性，而无需故障检测诊断机制。由于PFTCS不依赖于任何FDD过程或者在线控制器切换，不会出现在故障发生和控制响应之间的延迟，以及控制器切换的延迟，此外计算量也相对较低。基于这些优点，PFTCS成为航天器容错控制研究中的一个热门领域。虽然在已发表文献中，存在一些航天器执行机构的PFTCS研究成果，但是能够在不使用FDD装置的情况下，获得故障信息，并将其应用于PFTCS中的成果并不多见。

本发明正是针对这一难点问题，提出在不采用FDD装置的情况下，在线估计出三轴力矩故障因子并设计鲁棒自适应容错控制器。这一研究旨在丰富航天器PFTCS方法，并为将来的航天器姿态控制提供技术支持。

发明内容

本发明提出的鲁棒容错控制方案是基于一种新的自适应控制方法，它的设计主要包括两个部分：第一，设计一个三轴力矩有效性故障因子观测器，这个故障观测器是2014年由Qiang S.等人在四旋翼研究领域提出来的，在实际四旋翼工程领域得到了验证，但是由于应用背景存在差异，本发明对于此观测器的设计做出了自己的改进，使得其对航天器三轴力矩故障因子做出较为准确的估计，从而可以应用于航天器领域；第二，基于观测器获得的有效性故障因子估计值，设计了一种自适应控制方法，以实现在执行机构出现故障下的航天器机动控制。在本方案中，利用在线故障因子观测器实现了对于航天器三轴力矩有效性的估计，基于此设计的自适应容错控制器能够应对三轴力矩部分损失带来的不利影响，实现航天器的机动控制。稳定性分析表明，航天器姿态角速度和姿态角都可以渐进收敛到期望值。最后通过数值仿真，与传统的PD控制方法进行对比，验证本方案提出的容错控制方法的有效性和控制效果。

本发明一种基于鲁棒自适应的带有飞轮的航天器容错控制方法，该方法具体步骤如下：坐标系定义的说明：

本发明中涉及的坐标系有：惯性坐标系，这里取地心赤道惯性坐标系作为参考系，原点固联于地心o_I，o_Ix_I轴在赤道平面内，指向春分点，o_Iz_I轴垂直于赤道平面，与地球自转角速度矢量一致，o_Iy_I轴在赤道平面内按右手定则与o_Ix_I，o_Iz_I组成正交坐标系，表示为f_I；对于航天器本身而言，定义一个本体坐标系，原点为航天器的质心。o_bx_b、o_by_b和o_bz_b三轴固定在航天器本体上，且构成右手坐标系。令o_bx_b、o_by_b和o_bz_b三轴为航天器的惯量主轴，表示为f_b；期望的坐标系，定义为f_d。

步骤一：航天器运动学方程的建立

采用四元数来描述航天器的姿态。定义航天器相对惯性坐标系的姿态四元数为 $q={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T={\left[\begin{array}{cc}{q}_0& {\hat{q}}^T\end{array}\right]}^T,$ 其中q₀为四元数的标部，为四元数的矢部。四元数的四个参数满足如下的约束方程

$q_0^2+{\hat{q}}^T\hat{q}=1---\left(1\right)$

系统姿态运动学方程可以写为如下的形式

$\left\{\begin{array}{c}\stackrel{\·}{\hat{q}}=\frac{1}{2}({\hat{q}}^{:\×}{\mathrm{\ω}}_b+q_0{\mathrm{\ω}}_b)\\ {\stackrel{\·}{q}}_0=-\frac{1}{2}{\hat{q}}^T{\mathrm{\ω}}_b\end{array}\right.---\left(2\right)$

定义期望四元数 $q_d={\left[\begin{array}{cccc}{q}_{d0}& q_{d1}& q_{d2}& q_{d3}\end{array}\right]}^T={\left[\begin{array}{cc}{q}_{d0}& {\hat{q}}_d^T\end{array}\right]}^T,$ 为期望坐标系相对于惯性系姿态四元数；本体四元数 $q={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T={\left[\begin{array}{cc}{q}_0& {\hat{q}}^T\end{array}\right]}^T,$ 为本体坐标系相对于惯性系的姿态四元数；ω_b＝[ω_bx ω_by ω_bz]^T为航天器相对于惯性系的角速度在本体坐标系下的表示；姿态四元数误差定义为q和期望的姿态q_d之间的误差，表达式为：

$\left\{\begin{array}{c}{q}_{e0}=q_0{q}_{d0}+{\hat{q}}^T{\hat{q}}_d\\ {\hat{q}}_e=q_{d0}\hat{q}-q_0{\hat{q}}_d+{\hat{q}}^{\×}{\hat{q}}_d\end{array}\right.---\left(3\right)$

期望坐标系相对于惯性系的角速度在期望坐标系下表示为则在本体系下表示的角速度误差为：

${\mathrm{\ω}}_e={\mathrm{\ω}}_b-A_{\mathrm{bd}}{\mathrm{\ω}}_d^D---\left(4\right)$

其中，转换矩阵A_bd，将期望的坐标系S_d转换到本体坐标系S_b：

$A_{\mathrm{bd}}=(q_{e0}^2-{\hat{q}}_e^T{\hat{q}}_e)E^3+2{\hat{q}}_e{\hat{q}}_e^T-2q_{e0}{\hat{q}}_e^{\×}---\left(5\right)$

这里E³∈R^3×3是单位矩阵。误差四元数满足以下形式的运动学等式：

$\left\{\begin{array}{c}{\stackrel{\·}{q}}_{e0}=-\frac{1}{2}{\hat{q}}_e^T{\mathrm{\ω}}_e\\ {\stackrel{\·}{\hat{q}}}_e=\frac{1}{2}({\hat{q}}_e^{\×}+q_{e0}{E}_3){\mathrm{\ω}}_e\end{array}\right.---\left(6\right)$

可以进一步获得以下等式：

${\stackrel{\·}{\mathrm{\ω}}}_b={\stackrel{\·}{\mathrm{\ω}}}_e-{\mathrm{\ω}}_{\mathrm{bd}}^{\×}{A}_{\mathrm{bd}}{\mathrm{\ω}}_d^D---\left(7\right)$

步骤二：航天器动力学方程的建立

假设航天器是刚体航天器，不存在柔性附件，则其带有反作用飞轮的动力学模型表示如下：

$M=-{\mathrm{\ω}}_b^{\×}{I}_b{\mathrm{\ω}}_b---\left(9\right)$

H＝CI_wΩ   (10)

这里考虑刚体航天器的主惯量矩阵为I_b＝diag[I_b1 I_b2 I_b3]；外部干扰力矩为T_d∈R^3×1；通过反作用飞轮产生控制力矩T_w∈R^3×1；i个飞轮的角速度组成的列向量表示为Ω＝[Ω₁ Ω₂ ... Ω_i]^T；i个飞轮的转子轴向惯量组成的惯量对角阵为I_w＝diag[I_w1 I_w2 ... I_wn]；故障系数矩阵表示为

$E=\mathrm{diag}\left[\begin{array}{ccc}{e}_{\mathrm{bx}}& e_{\mathrm{by}}& e_{\mathrm{bz}}\end{array}\right],0<{\mathrm{\&epsiv;}}_i\&le;e_{\mathrm{bi}}1,i=x,y,z---\left(11\right)$

这里e_bi＝1,(i＝x,y,z)意味着相对于体坐标系的三轴方向没有力矩输出故障，

e_bi＝0,(i＝x,y,z)意味着在第i轴的方向上完全故障，没有力矩输出，ε_i＞0意味着并不存在第i轴方向完全故障的情况。

步骤三：期望的角速度轨迹设计

由(2)式可以改写为：

$\stackrel{\&CenterDot;}{q}=\frac{1}{2}F\left(q\right){\mathrm{\&omega;}}_b---\left(12\right)$

这里q为四元数。Crassidis&Markley提出角速度ω_b是由输入力矩控制T_w的，与此同时ω_b控制姿态q。令

${\mathrm{\&omega;}}_b=-\frac{1}{{\mathrm{\&lambda;}}_1}{F}^{-1}\left(q\right)q_e---\left(13\right)$

这里q_e是(3)式的表达形式。结合(12)给出角度位置：

${\mathrm{\&lambda;}}_1\stackrel{\&CenterDot;}{q}=q_e,---\left(14\right)$

由于(14)是一阶微分方程，式子中的λ₁＞0，因而，姿态四元数q收敛到期望的姿态四元数q_d。这就意味着式子(13)可以作为航天器期望的角速度，表达式如下：

$A_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D=-\frac{1}{{\mathrm{\&lambda;}}_1}{F}^{-1}\left(q\right)q_e,---15\left(7\right)$

步骤四：控制器输出力矩的设计

考虑当反作用飞轮的输出力矩存在部分损失，则通过控制器设计得到的指令力矩具有一定的容错能力。为了得到在线的故障信息，通过式子(8)中发生故障时的动力学表达式可以得到带有在线估计故障因子的动力学方程，由于干扰力矩T_d相对于控制力矩而言是小量，在这里进行忽略：

$I_b{\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}}_b=\hat{M}-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\hat{E}T}_v---\left(16\right)$

这里 ${\widehat{\mathrm{\&omega;}}}_b={\left[\begin{array}{ccc}{\widehat{\mathrm{\&omega;}}}_{b,1}& {\widehat{\mathrm{\&omega;}}}_{b,2}& {\widehat{\mathrm{\&omega;}}}_{b,3}\end{array}\right]}^T\&Element;R^3$ 是估计的航天器本体角速度， $\hat{E}=\mathrm{diag}\left[\begin{array}{ccc}{\hat{e}}_1\left(t\right)& {\hat{e}}_w\left(t\right)& {\hat{e}}_3\left(t\right)\end{array}\right]\&Element;R^{3\&times;3}$ 是三轴方向上的故障系数的估计矩阵，中含有估计量观测的故障动力学的输入T_v＝[T_v1 T_v2 T_v3]^T∈R³将会接下来给出。

定义观测误差向量 ${\widetilde{\mathrm{\&omega;}}}_b={\left[\begin{array}{ccc}{\widetilde{\mathrm{\&omega;}}}_{b,1}& {\widetilde{\mathrm{\&omega;}}}_{b,2}& {\widetilde{\mathrm{\&omega;}}}_{b,3}\end{array}\right]}^T\&Element;R^3,$ 这里i∈{1,2,3}，为了可以进行角速度误差的补偿，选择力矩然后估计误差动力学方程表示如下：

$I_b{\stackrel{\&CenterDot;}{\widetilde{\mathrm{\&omega;}}}}_b=\tilde{E}{T}_v-k_{u}E{\widetilde{\mathrm{\&omega;}}}_b-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\mathrm{\&omega;}}_b^{\&times;}\mathrm{EH}---\left(17\right)$

这里 $\tilde{E}=\mathrm{diag}\left[\begin{array}{ccc}{\tilde{e}}_1& {\tilde{e}}_2& {\tilde{e}}_3\end{array}\right]\&Element;R^{3\&times;3}$ 是故障系数的估计误差矩阵，其中k_u是正常数。设计一个输入T_v，使得估计的航天器角速度ω_b可以达到期望的角速度而观测的角速度和期望角速度之间的误差向量定义为： ${\widehat{\mathrm{\&omega;}}}_e={\left[\begin{array}{ccc}{\widehat{\mathrm{\&omega;}}}_{e,1}& {\widehat{\mathrm{\&omega;}}}_{e,2}& {\widehat{\mathrm{\&omega;}}}_{e,3}\end{array}\right]}^T\&Element;R^3,$ 其中i∈{1,2,3}。使用式子(16)，估计与期望的误差动力学方程可写成：

$I_b{\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}}_e=\hat{M}-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_d^D---\left(18\right)$

定理1：考虑等式(3)中带有故障因子的动力学，输出力矩设计为：

$T_w=T_v+k_u{\widetilde{\mathrm{\&omega;}}}_b---\left(19\right)$

这里是观测器方程的输入，k_v是一个正常数。故障观测器为：

$\begin{array}{c}{\stackrel{\&CenterDot;}{\hat{e}}}_i={\mathrm{Proj}}_{[{\mathrm{\&epsiv;}}_i,1]}\{-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\}\\ =\left\{\begin{array}{ccc}0,& \mathrm{if}& {\hat{e}}_i={\mathrm{\&epsiv;}}_i,-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\&le;0\mathrm{or}{\hat{e}}_i=1,-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\&GreaterEqual;0\\ -{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i},& & \mathrm{otherwise}\end{array}\right.\end{array}---\left(20\right)$

这里的α_i是正自适应增益，并且这个投影算子Proj{·}用来在参数边界内保持参数的估计。然后，每个角速度ω_i都渐进地收敛到期望的角速度ω_d,i，尽管有效性部分损失，但是期望的输出力矩可以有效输出。

证明：如下构造李雅普诺夫方程

$V=\frac{1}{2}{\widehat{\mathrm{\&omega;}}}_e^T{I}_b{\widehat{\mathrm{\&omega;}}}_e+\frac{1}{2}{\widetilde{\mathrm{\&omega;}}}_b^T{I}_b{\widetilde{\mathrm{\&omega;}}}_b^T+\frac{1}{2}\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i^2\left(t\right)}{{\mathrm{\&alpha;}}_i},---\left(21\right)$

从方程式(17)和(18)而言，李雅普诺夫方程的时间导数满足

$\begin{array}{c}\stackrel{\&CenterDot;}{V}={\widehat{\mathrm{\&omega;}}}_e^T(M-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_d^D)+{\widetilde{\mathrm{\&omega;}}}_b^T(\tilde{E}{T}_v-k_{u}E{\widetilde{\mathrm{\&omega;}}}_b-{\widetilde{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\mathrm{\&omega;}}_b^{\&times;}\mathrm{EH})+\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\&CenterDot;}{\tilde{e}}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\\ \&le;{-k}_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u{\widetilde{\mathrm{\&omega;}}}_b^{T}E{\widetilde{\mathrm{\&omega;}}}_b+{\widetilde{\mathrm{\&omega;}}}_b^T\tilde{E}{T}_v+\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\&CenterDot;}{\tilde{e}}}_e\left(t\right)}{{\mathrm{\&alpha;}}_i}\end{array}---\left(22\right)$

由于e_i是未知的常数，基于投影算子，可以得到以下不等式

$\frac{{\tilde{e}}_i\left(t\right){\stackrel{\&CenterDot;}{\tilde{e}}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\&le;{\tilde{e}}_i\left(t\right){\widetilde{\mathrm{\&omega;}}}_i{v}_i---\left(23\right)$

令θ＝ζ_min{E}，这里ζ_min{.}算子代表一个矩阵的最小特征值。因为e_i满足0＜ε_i≤e_bi≤1，所以θ＞0。通过采用自适应律式子(18)和(22)可以进一步得到

$\begin{array}{c}\stackrel{\&CenterDot;}{V}\&le;-k_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u\mathrm{\&theta;}{\left|\right|\widetilde{\mathrm{\&omega;}}\left|\right|}^2+\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}({\widetilde{\mathrm{\&omega;}}}_i{\tilde{e}}_i{v}_i+\frac{{\tilde{e}}_i\left(t\right)\stackrel{\&CenterDot;}{\tilde{e}}\left(t\right)}{{\mathrm{\&alpha;}}_i})\\ \&le;{-k}_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u\mathrm{\&theta;}{\left|\right|\widetilde{\mathrm{\&omega;}\left|\right|}}^2\end{array}---\left(24\right)$

因此，是负半定的。可以得到和将从0积分到∞，可以获得

$V\left(0\right)-V(\&infin;)\&GreaterEqual;k_v\underset{0}{\overset{\&infin;}{\&Integral;}}{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2\mathrm{d\&eta;}+k_u\underset{0}{\overset{\&infin;}{\&Integral;}}{\left|\right|{\widetilde{\mathrm{\&omega;}}}_b\left|\right|}^2\mathrm{d\&eta;}---\left(25\right)$

因为以上不等式左边的项是有界的，所以可以得到和除此而外，通过式子(18)和(19)可以容易的证明和因此，通过Barbalat引理可以得到以下等式这个可以推导出尽管存在三轴力矩的部分损失，但是设计的姿态控制力矩得以保持，并使得航天器的角速度渐进收敛到期望的值。总之，在出现三轴输出力矩损失的情况下，可以获得以上设计的姿态控制力矩。

步骤五：数值仿真

为了证明上述方案的有效性，下面通过数值仿真，将上述容错控制方案与传统PD控制方案相比较。刚体航天器本体的惯量矩阵为I_b＝diag(295 130 210)(kg.m²),，假设航天器上装有四个飞轮，飞轮组的惯量阵为：

I_w＝diag(0.01044 0.01044 0.01044 0.01044)(kg·m²)

其安装采用四斜装构型，安装阵为

$C=\frac{1}{2}\left[\begin{array}{cccc}1& 1& -1& -1\\ -1& 1& 1& -1\\ \sqrt{2}& \sqrt{2}& \sqrt{2}& \sqrt{2}\end{array}\right]---\left(26\right)$

初始时刻本体坐标系相对于惯性坐标系的姿态四元数为

$q\left(0\right)=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right],---\left(27\right)$

期望坐标系相对于惯性坐标系的四元数为

$q_d=\left[\begin{array}{c}0.93\\ 0.1\\ -0.2\\ 0.3\end{array}\right]---\left(28\right)$

在控制器设计时，本发明忽略了外部干扰力矩，为了更加符合实际的空间环境，这里加入外部干扰力矩T_d，假设T_d为周期变化形式

$T_d=\left[\begin{array}{c}0.3\cos \left(0.01t\right)+0.1\\ 0.15\sin \left(0.02t\right)+0.3\cos \left(0.025t\right)\\ 0.3\sin \left(0.01t\right)+0.1\end{array}\right]\mathrm{Nm}---\left(29\right)$

a)无故障条件仿真:

在没有故障的条件下，采用传统的PD控制方法，控制参数设定如下：

k₁＝50,k₂＝50

本发明的容错控制方法，控制参数设定如下。有效性故障因子的初始估计自适应律参数ε_i＝0.4,{i＝(1,2,3)}，α_i＝0.5,{i＝(1,2,3)}，k_u＝30，k_v＝8。

图2到图7展示了没有故障的仿真结果。

b)故障条件的仿真:

在仿真中设定以下故障情况

$\left\{\begin{array}{cc}{e}_y=0.2& t\&GreaterEqual;10\\ e_z=0.5& t\&GreaterEqual;20\end{array}\right.$

这个表明在y轴方向10s之后控制能力损失了80％，z轴方向在20s之后损失了50％。控制方法剩余的控制参数没有改变。图8到图13展示了故障情况下的仿真结果。

本方案不需要基于故障信息的检测进行控制器的重构，属于PFTC。然而相比AFTC或其它PFTC方法，本方案提出的自适应容错控制方法有如下优点：

1)故障因子在线观测器，只需要很少数量的系统信息，便能估计出故障因子，为自适应容错控制方法提供参数信息；

2)方案不需要FDD机制和准确的故障信息。相应的，通过自适应容错控制本身对干扰力矩主动补偿和实时地对部分输出力矩损失进行快速估计的方法，该PFTC方案可以一定程度地减少环境干扰力矩和执行机构故障对控制性能和稳定性的不良影响；

3)该方法具有更好的灵活性和容错能力，可用于处理执行机构故障的情况；

4)它计算量小，并能更好地适用于工程实际应用。

附图说明

图1容错控制方案示意图。

图2航天器无故障情况下误差四元数响应曲线(虚部)(传统PD控制方法)。

图3航天器无故障情况下误差四元数响应曲线(虚部)(容错控制方法)。

图4航天器无故障情况下本体角速度响应曲线(传统PD控制方法)。

图5航天器无故障情况下本体角速度响应曲线(容错控制方法)。

图6航天器无故障情况下飞轮转速响应曲线(传统PD控制方法)。

图7航天器无故障情况下飞轮转速响应曲线(容错控制方法)。

图8航天器发生故障情况下误差四元数响应曲线(虚部)(传统PD控制方法)。

图9航天器发生故障情况下误差四元数响应曲线(虚部)(容错控制方法)。

图10航天器发生故障情况下角速度响应曲线(传统PD控制方法)。

图11航天器发生故障情况下角速度响应曲线(容错控制方法)。

图12航天器发生故障情况下飞轮转速响应曲线(传统PD控制方法)。

图13航天器发生故障情况下飞轮转速响应曲线(容错控制方法)。

图中符号说明如下：

图2-图13中横坐标，Time均代表时间，s代表以秒为单位；图2，图3，图8和图9纵坐标q_e代表误差四元数(虚部)，图例中q_e1,q_e2,q_e3分别代表四元数虚部的三个参数；图4，图5，图10，图11纵坐标代表本体角速度ω_b，图例中ω_b1,ω_b2,ω_b3分别代表航天器本体三正交轴的角速度，单位为弧度每秒(rad/s)；图6，图7，图12和图13纵坐标航天器飞轮转速Ω，图例中Ω₁,Ω₂,Ω₃,Ω₄分别代表航天器安装的四个飞轮的转速，单位为转每分(rpm)。

具体实施方式

下面结合附图对本方案作具体的说明。图1为本发明容错控制方案示意图。综上所述，本发明一种基于鲁棒自适应的带有飞轮的航天器容错控制方法，该方法具体步骤如下：

坐标系定义的说明：

本发明中涉及的坐标系有：惯性坐标系，这里取地心赤道惯性坐标系作为参考系，原点固联于地心o_I，o_Ix_I轴在赤道平面内，指向春分点，o_Iz_I轴垂直于赤道平面，与地球自转角速度矢量一致，o_Iy_I轴在赤道平面内按右手定则与o_Ix_I，o_Iz_I组成正交坐标系，表示为f_I；对于航天器本身而言，定义一个本体坐标系，原点为航天器的质心。o_bx_b、o_by_b和o_bz_b三轴固定在航天器本体上，且构成右手坐标系。令o_bx_b、o_by_b和o_bz_b三轴为航天器的惯量主轴，表示为f_b；期望的坐标系，定义为f_d。

步骤一：航天器运动学方程的建立

采用四元数来描述航天器的姿态。定义航天器相对惯性坐标系的姿态四元数为 $q={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T={\left[\begin{array}{cc}{q}_0& {\hat{q}}^T\end{array}\right]}^T,$ 其中q₀为四元数的标部，为四元数的矢部。四元数的四个参数满足如下的约束方程

$q_0^2+{\hat{q}}^T\hat{q}=1---\left(1\right)$

系统姿态运动学方程可以写为如下的形式

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{\hat{q}}=\frac{1}{2}({\hat{q}}^{:\&times;}{\mathrm{\&omega;}}_b+q_0{\mathrm{\&omega;}}_b)\\ {\stackrel{\&CenterDot;}{q}}_0=-\frac{1}{2}{\hat{q}}^T{\mathrm{\&omega;}}_b\end{array}\right.---\left(2\right)$

定义期望四元数 $q_d={\left[\begin{array}{cccc}{q}_{d0}& q_{d1}& q_{d2}& q_{d3}\end{array}\right]}^T={\left[\begin{array}{cc}{q}_{d0}& {\hat{q}}_d^T\end{array}\right]}^T,$ 为期望坐标系相对于惯性系姿态四元数；本体四元数 $q={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T={\left[\begin{array}{cc}{q}_0& {\hat{q}}^T\end{array}\right]}^T,$ 为本体坐标系相对于惯性系的姿态四元数；ω_b＝[ω_bx ω_by ω_bz]^T为航天器相对于惯性系的角速度在本体坐标系下的表示；姿态四元数误差定义为q和期望的姿态q_d之间的误差，表达式为：

$\left\{\begin{array}{c}{q}_{e0}=q_0{q}_{d0}+{\hat{q}}^T{\hat{q}}_d\\ {\hat{q}}_e=q_{d0}\hat{q}-q_0{\hat{q}}_d+{\hat{q}}^{\&times;}{\hat{q}}_d\end{array}\right.---\left(3\right)$

期望坐标系相对于惯性系的角速度在期望坐标系下表示为则在本体系下表示的角速度误差为：

${\mathrm{\&omega;}}_e={\mathrm{\&omega;}}_b-A_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D---\left(4\right)$

其中，转换矩阵A_bd，将期望的坐标系S_d转换到本体坐标系S_b：

$A_{\mathrm{bd}}=(q_{e0}^2-{\hat{q}}_e^T{\hat{q}}_e)E^3+2{\hat{q}}_e{\hat{q}}_e^T-2q_{e0}{\hat{q}}_e^{\&times;}---\left(5\right)$

这里E³∈R^3×3是单位矩阵。误差四元数满足以下形式的运动学等式：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{q}}_{e0}=-\frac{1}{2}{\hat{q}}_e^T{\mathrm{\&omega;}}_e\\ {\stackrel{\&CenterDot;}{\hat{q}}}_e=\frac{1}{2}({\hat{q}}_e^{\&times;}+q_{e0}{E}_3){\mathrm{\&omega;}}_e\end{array}\right.---\left(6\right)$

可以进一步获得以下等式：

${\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b={\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_e-{\mathrm{\&omega;}}_{\mathrm{bd}}^{\&times;}{A}_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D---\left(7\right)$

步骤二：航天器动力学方程的建立

假设航天器是刚体航天器，不存在柔性附件，则其带有反作用飞轮的动力学模型表示如下：

$M=-{\mathrm{\&omega;}}_b^{\&times;}{I}_b{\mathrm{\&omega;}}_b---\left(9\right)$

H＝CI_wΩ   (10)

这里考虑刚体航天器的主惯量矩阵为I_b＝diag[I_b1 I_b2 I_b3]；外部干扰力矩为T_d∈R^3×1；通过反作用飞轮产生控制力矩T_w∈R^3×1；i个飞轮的角速度组成的列向量表示为Ω＝[Ω₁ Ω₂ ... Ω_i]^T；i个飞轮的转子轴向惯量组成的惯量对角阵为I_w＝diag[I_w1 I_w2 ... I_wn]；故障系数矩阵表示为

$E=\mathrm{diag}\left[\begin{array}{ccc}{e}_{\mathrm{bx}}& e_{\mathrm{by}}& e_{\mathrm{bz}}\end{array}\right],0<{\mathrm{\&epsiv;}}_i\&le;e_{\mathrm{bi}}1,i=x,y,z---\left(11\right)$

这里e_bi＝1,(i＝x，y，z)意味着相对于体坐标系的三轴方向没有力矩输出故障，e_bi＝0,(i＝x,y,z)意味着在第i轴的方向上完全故障，没有力矩输出，ε_i＞0意味着并不存在第i轴方向完全故障的情况。

步骤三：期望的角速度轨迹设计

由(2)式可以改写为：

$\stackrel{\&CenterDot;}{q}=\frac{1}{2}F\left(q\right){\mathrm{\&omega;}}_b---\left(12\right)$

这里q为四元数。Crassidis&Markley提出角速度ω_b是由输入力矩控制T_w的，与此同时ω_b控制姿态q。令

${\mathrm{\&omega;}}_b=-\frac{1}{{\mathrm{\&lambda;}}_1}{F}^{-1}\left(q\right)q_e---\left(13\right)$

这里q_e是(3)式的表达形式。结合(12)给出角度位置：

${\mathrm{\&lambda;}}_1\stackrel{\&CenterDot;}{q}=q_e,---\left(14\right)$

由于(14)是一阶微分方程，式子中的λ₁＞0，因而，姿态四元数q收敛到期望的姿态四元数q_d。这就意味着式子(13)可以作为航天器期望的角速度，表达式如下：

$A_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D=-\frac{1}{{\mathrm{\&lambda;}}_1}{F}^{-1}\left(q\right)q_e,---15\left(7\right)$

步骤四：控制器输出力矩的设计

考虑当反作用飞轮的输出力矩存在部分损失，则通过控制器设计得到的指令力矩具有一定的容错能力。为了得到在线的故障信息，通过式子(8)中发生故障时的动力学表达式可以得到带有在线估计故障因子的动力学方程，由于干扰力矩T_d相对于控制力矩而言是小量，在这里进行忽略：

$I_b{\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}}_b=\hat{M}-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\hat{E}T}_v---\left(16\right)$

这里 ${\widehat{\mathrm{\&omega;}}}_b={\left[\begin{array}{ccc}{\widehat{\mathrm{\&omega;}}}_{b,1}& {\widehat{\mathrm{\&omega;}}}_{b,2}& {\widehat{\mathrm{\&omega;}}}_{b,3}\end{array}\right]}^T\&Element;R^3$ 是估计的航天器本体角速度， $\hat{E}=\mathrm{diag}\left[\begin{array}{ccc}{\hat{e}}_1\left(t\right)& {\hat{e}}_w\left(t\right)& {\hat{e}}_3\left(t\right)\end{array}\right]\&Element;R^{3\&times;3}$ 是三轴方向上的故障系数的估计矩阵，中含有估计量观测的故障动力学的输入T_v＝[T_v1 T_v2 T_v3]^T∈R³将会接下来给出。

定义观测误差向量 ${\widetilde{\mathrm{\&omega;}}}_b={\left[\begin{array}{ccc}{\widetilde{\mathrm{\&omega;}}}_{b,1}& {\widetilde{\mathrm{\&omega;}}}_{b,2}& {\widetilde{\mathrm{\&omega;}}}_{b,3}\end{array}\right]}^T\&Element;R^3,$ 这里i∈{1,2,3}，为了可以进行角速度误差的补偿，选择力矩然后估计误差动力学方程表示如下：

$I_b{\stackrel{\&CenterDot;}{\widetilde{\mathrm{\&omega;}}}}_b=\tilde{E}{T}_v-k_{u}E{\widetilde{\mathrm{\&omega;}}}_b-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\mathrm{\&omega;}}_b^{\&times;}\mathrm{EH}---\left(17\right)$

这里 $\tilde{E}=\mathrm{diag}\left[\begin{array}{ccc}{\tilde{e}}_1& {\tilde{e}}_2& {\tilde{e}}_3\end{array}\right]\&Element;R^{3\&times;3}$ 是故障系数的估计误差矩阵，其中k_u是正常数。设计一个输入T_v，使得估计的航天器角速度ω_b可以达到期望的角速度而观测的角速度和期望角速度之间的误差向量定义为： ${\widehat{\mathrm{\&omega;}}}_e={\left[\begin{array}{ccc}{\widehat{\mathrm{\&omega;}}}_{e,1}& {\widehat{\mathrm{\&omega;}}}_{e,2}& {\widehat{\mathrm{\&omega;}}}_{e,3}\end{array}\right]}^T\&Element;R^3,$ 其中i∈{1,2,3}。使用式子(16)，估计与期望的误差动力学方程可写成：

$I_b{\stackrel{\&CenterDot;}{\widehat{\mathrm{\&omega;}}}}_e=\hat{M}-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_d^D---\left(18\right)$

定理1：考虑等式(3)中带有故障因子的动力学，输出力矩设计为：

$T_w=T_v+k_u{\widetilde{\mathrm{\&omega;}}}_b---\left(19\right)$

这里是观测器方程的输入，k_v是一个正常数。故障观测器为：

$\begin{array}{c}{\stackrel{\&CenterDot;}{\hat{e}}}_i={\mathrm{Proj}}_{[{\mathrm{\&epsiv;}}_i,1]}\{-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\}\\ =\left\{\begin{array}{ccc}0,& \mathrm{if}& {\hat{e}}_i={\mathrm{\&epsiv;}}_i,-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\&le;0\mathrm{or}{\hat{e}}_i=1,-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\&GreaterEqual;0\\ -{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i},& & \mathrm{otherwise}\end{array}\right.\end{array}---\left(20\right)$

这里的α_i是正自适应增益，并且这个投影算子Proj{·}用来在参数边界内保持参数的估计。然后，每个角速度ω_i都渐进地收敛到期望的角速度ω_d,i，尽管有效性部分损失，但是期望的输出力矩可以有效输出。

证明：如下构造李雅普诺夫方程

$V=\frac{1}{2}{\widehat{\mathrm{\&omega;}}}_e^T{I}_b{\widehat{\mathrm{\&omega;}}}_e+\frac{1}{2}{\widetilde{\mathrm{\&omega;}}}_b^T{I}_b{\widetilde{\mathrm{\&omega;}}}_b^T+\frac{1}{2}\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i^2\left(t\right)}{{\mathrm{\&alpha;}}_i},---\left(21\right)$

从方程式(17)和(18)而言，李雅普诺夫方程的时间导数满足

$\begin{array}{c}\stackrel{\&CenterDot;}{V}={\widehat{\mathrm{\&omega;}}}_e^T(M-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_d^D)+{\widetilde{\mathrm{\&omega;}}}_b^T(\tilde{E}{T}_v-k_{u}E{\widetilde{\mathrm{\&omega;}}}_b-{\widetilde{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\mathrm{\&omega;}}_b^{\&times;}\mathrm{EH})+\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\&CenterDot;}{\tilde{e}}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\\ \&le;{-k}_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u{\widetilde{\mathrm{\&omega;}}}_b^{T}E{\widetilde{\mathrm{\&omega;}}}_b+{\widetilde{\mathrm{\&omega;}}}_b^T\tilde{E}{T}_v+\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\&CenterDot;}{\tilde{e}}}_e\left(t\right)}{{\mathrm{\&alpha;}}_i}\end{array}---\left(22\right)$

由于e_i是未知的常数，基于投影算子，可以得到以下不等式

$\frac{{\tilde{e}}_i\left(t\right){\stackrel{\&CenterDot;}{\tilde{e}}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\&le;{\tilde{e}}_i\left(t\right){\widetilde{\mathrm{\&omega;}}}_i{v}_i---\left(23\right)$

令θ＝ζ_min{E}，这里ζ_min{.}算子代表一个矩阵的最小特征值。因为e_i满足0＜ε_i≤e_bi≤1，所以θ＞0。通过采用自适应律式子(18)和(22)可以进一步得到

$\begin{array}{c}\stackrel{\&CenterDot;}{V}\&le;-k_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u\mathrm{\&theta;}{\left|\right|\widetilde{\mathrm{\&omega;}}\left|\right|}^2+\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}({\widetilde{\mathrm{\&omega;}}}_i{\tilde{e}}_i{v}_i+\frac{{\tilde{e}}_i\left(t\right)\stackrel{\&CenterDot;}{\tilde{e}}\left(t\right)}{{\mathrm{\&alpha;}}_i})\\ \&le;{-k}_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u\mathrm{\&theta;}{\left|\right|\widetilde{\mathrm{\&omega;}\left|\right|}}^2\end{array}---\left(24\right)$

因此，是负半定的。可以得到和将从0积分到∞，可以获得

$V\left(0\right)-V(\&infin;)\&GreaterEqual;k_v\underset{0}{\overset{\&infin;}{\&Integral;}}{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2\mathrm{d\&eta;}+k_u\underset{0}{\overset{\&infin;}{\&Integral;}}{\left|\right|{\widetilde{\mathrm{\&omega;}}}_b\left|\right|}^2\mathrm{d\&eta;}---\left(25\right)$

因为以上不等式左边的项是有界的，所以可以得到和除此而外，通过式子(18)和(19)可以容易的证明和因此，通过Barbalat引理可以得到以下等式这个可以推导出尽管存在三轴力矩的部分损失，但是设计的姿态控制力矩得以保持，并使得航天器的角速度渐进收敛到期望的值。总之，在出现三轴输出力矩损失的情况下，可以获得以上设计的姿态控制力矩。

步骤五：数值仿真

为了证明上述方案的有效性，下面通过数值仿真，将上述容错控制方案与传统PD控制方案相比较。刚体航天器本体的惯量矩阵为I_b＝diag(295 130 210)(kg.m²),，假设航天器上装有四个飞轮，飞轮组的惯量阵为：

I_w＝diag(0.01044 0.01044 0.01044 0.01044)(kg·m²)

其安装采用四斜装构型，安装阵为

$C=\frac{1}{2}\left[\begin{array}{cccc}1& 1& -1& -1\\ -1& 1& 1& -1\\ \sqrt{2}& \sqrt{2}& \sqrt{2}& \sqrt{2}\end{array}\right]---\left(26\right)$

初始时刻本体坐标系相对于惯性坐标系的姿态四元数为

$q\left(0\right)=\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right],---\left(27\right)$

期望坐标系相对于惯性坐标系的四元数为

$q_d=\left[\begin{array}{c}0.93\\ 0.1\\ -0.2\\ 0.3\end{array}\right]---\left(28\right)$

在控制器设计时，本发明忽略了外部干扰力矩，为了更加符合实际的空间环境，这里加入外部干扰力矩T_d，假设T_d为周期变化形式

$T_d=\left[\begin{array}{c}0.3\cos \left(0.01t\right)+0.1\\ 0.15\sin \left(0.02t\right)+0.3\cos \left(0.025t\right)\\ 0.3\sin \left(0.01t\right)+0.1\end{array}\right]\mathrm{Nm}---\left(29\right)$

a)无故障条件仿真:

在没有故障的条件下，采用传统的PD控制方法，控制参数设定如下：

k₁＝50,k₂＝50

本发明的容错控制方法，控制参数设定如下。有效性故障因子的初始估计自适应律参数ε_i＝0.4,{i＝(1,2,3)}，α_i＝0.5,{i＝(1,2,3)}，k_u＝30，k_v＝8。

图2到图7展示了没有故障的仿真结果。

b)故障条件的仿真:

在仿真中设定以下故障情况

$\left\{\begin{array}{cc}{e}_y=0.2& t\&GreaterEqual;10\\ e_z=0.5& t\&GreaterEqual;20\end{array}\right.$

这个表明在y轴方向10s之后控制能力损失了80％，z轴方向在20s之后损失了50％。控制方法剩余的控制参数没有改变。图8到图13展示了故障情况下的仿真结果。

图2展示传统PD控制方法可以使得无故障航天器的误差四元数(虚部)在50s中收敛到0。图3显示采用所提出的容错控制方法在110s内将航天器误差四元数(虚部)收敛到0。图4仿真结果显示，采用传统的PD控制方法可以使得航天器本体角速度在50s之后收敛到0；图5表明应用所提出的容错控制方法可以在110s内将本体角速度稳定0；图6表明采用传统PD方法使得无故障航天器飞轮转速在50s之后稳定变化；图7显示了所提出的容错控制方法在无故障情况下，在110s内飞轮转速达到稳定状态。两种方法都获得了期望的效果，和传统的控制方法相比，所提出的方法的调整时间要略慢一些，但是控制性能更加平稳。

图8,10,12显示采用传统的PD控制方法无法使航天器在三轴力矩有效性损失的故障情况下，将航天器机动到期望的姿态。作为对比，图9显示，采用所提出的容错控制方法，在120s内，航天器的误差四元数(虚部)收敛到了0；图11显示航天器的本体角速度在120s内到达0；图13显示在120s内航天器飞轮转速达到稳定转速。

综上所述，本发明阐述使用间接鲁棒自适应控制技术实现航天器存在三轴力矩部分损失时的姿态机动问题。利用四元数建立运动学模型，在控制器的设计中，利用一个在线的观测器，估计三轴力矩有效性因子，在没有FDD机构的情况下去弥补航天器故障以及外部干扰带来的不利影响，从而达到姿态机动控制目标。仿真结果可以看到，相比于传统的PD控制方法，在无故障情况下，所提出的控制方法，控制速度方面并没有优势，但是存在故障的时候，该容错控制方案可以处理有多个执行机构故障的情况，能最大程度上维持控制系统的稳定性和保证其良好的控制性能。此外，不依赖于FDD，计算量小，具有很好的实时性，更适合于实际航天器在轨应用。

Claims (1)

1.一种基于鲁棒自适应的带有飞轮的航天器容错控制方法，其特征在于：该方法具体步骤如下：

涉及的坐标系有：惯性坐标系，这里取地心赤道惯性坐标系作为参考系，原点固联于地心o_I，o_Ix_I轴在赤道平面内，指向春分点，o_Iz_I轴垂直于赤道平面，与地球自转角速度矢量一致，o_Iy_I轴在赤道平面内按右手定则与o_Ix_I，o_Iz_I组成正交坐标系，表示为f_I；对于航天器本身而言，定义一个本体坐标系，原点为航天器的质心；o_bx_b、o_by_b和o_bz_b三轴固定在航天器本体上，且构成右手坐标系；令o_bx_b、o_by_b和o_bz_b三轴为航天器的惯量主轴，表示为f_b；期望的坐标系，定义为f_d；

步骤一：航天器运动学方程的建立

采用四元数来描述航天器的姿态；定义航天器相对惯性坐标系的姿态四元数为 $q={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T={\left[\begin{array}{cc}{q}_0& {\hat{q}}^T\end{array}\right]}^T,$ 其中q₀为四元数的标部，为四元数的矢部，四元数的四个参数满足如下的约束方程

$q_0^2+{\hat{q}}^T\hat{q}=1---\left(1\right)$

系统姿态运动学方程写为如下的形式

$\left\{\begin{array}{c}\stackrel{\stackrel{.}{^}}{q}=\frac{1}{2}({\hat{q}}^{\&times;}{\mathrm{\&omega;}}_b+q_0{\mathrm{\&omega;}}_b)\\ {\stackrel{.}{q}}_0=-\frac{1}{2}{\hat{q}}^T{\mathrm{\&omega;}}_b\end{array}\right.---\left(2\right)$

定义期望四元数 $q_d={\left[\begin{array}{cccc}{q}_{d0}& q_{d1}& q_{d2}& q_{d3}\end{array}\right]}^T={\left[\begin{array}{cc}{q}_{d0}& {\hat{q}}_d^T\end{array}\right]}^T,$ 为期望坐标系相对于惯性系姿态四元数；本体四元数 $q={\left[\begin{array}{cccc}{q}_0& q_1& q_2& q_3\end{array}\right]}^T={\left[\begin{array}{cc}{q}_0& {\hat{q}}^T\end{array}\right]}^T,$ 为本体坐标系相对于惯性系的姿态四元数；ω_b＝[ω_bx ω_by ω_bz]^T为航天器相对于惯性系的角速度在本体坐标系下的表示；姿态四元数误差定义为q和期望的姿态q_d之间的误差，表达式为：

$\left\{\begin{array}{c}{q}_{e0}=q_0{q}_{d0}+{\hat{q}}^T{\hat{q}}_d\\ {\hat{q}}_e=q_{d0}\hat{q}-q_0{\hat{q}}_d+{\hat{q}}^{\&times;}{\hat{q}}_d\end{array}\right.---\left(3\right)$

期望坐标系相对于惯性系的角速度在期望坐标系下表示为则在本体系下表示的角速度误差为：

${\mathrm{\&omega;}}_e={\mathrm{\&omega;}}_b-A_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D---\left(4\right)$

其中，转换矩阵A_bd，将期望的坐标系S_d转换到本体坐标系S_b：

$A_{\mathrm{bd}}=(q_{e0}^2-{\hat{q}}_e^T{\hat{q}}_e)E^3+2{\hat{q}}_e{\hat{q}}_e^T-2q_{e0}{\hat{q}}_e^{\&times;}---\left(5\right)$

这里E³∈R^3×3是单位矩阵，误差四元数满足以下形式的运动学等式：

$\left\{\begin{array}{c}{\stackrel{.}{q}}_{e0}=-\frac{1}{2}{\stackrel{.}{q}}_e^T{\mathrm{\&omega;}}_e\\ {\stackrel{\stackrel{.}{^}}{q}}_e=\frac{1}{2}({\hat{q}}_e^{\&times;}+q_{e0}{E}_3){\mathrm{\&omega;}}_e\end{array}\right.---\left(6\right)$

进一步获得以下等式：

${\stackrel{.}{\mathrm{\&omega;}}}_b={\stackrel{.}{\mathrm{\&omega;}}}_e-{{\mathrm{\&omega;}}_{\mathrm{bd}}^{\&times;}A}_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D---\left(7\right)$

步骤二：航天器动力学方程的建立

假设航天器是刚体航天器，不存在柔性附件，则其带有反作用飞轮的动力学模型表示如下：

$M=-{\mathrm{\&omega;}}_b^{\&times;}{I}_b{\mathrm{\&omega;}}_b---\left(9\right)$

H＝CI_wΩ            (10)

这里考虑刚体航天器的主惯量矩阵为I_b＝diag[I_b1 I_b2 I_b3]；外部干扰力矩为T_d∈R^3×1；通过反作用飞轮产生控制力矩T_w∈R^3×1；i个飞轮的角速度组成的列向量表示为Ω＝[Ω₁ Ω₂ ... Ω_i]^T；i个飞轮的转子轴向惯量组成的惯量对角阵为I_w＝diag[I_w1 I_w2 ... I_wn]；故障系数矩阵表示为

E＝diag[e_bx e_by e_bz],0＜ε_i≤e_bi≤1,i＝x,y,z       (11)

这里e_bi＝1,(i＝x,y,z)意味着相对于体坐标系的三轴方向没有力矩输出故障，e_bi＝0,(i＝x,y,z)意味着在第i轴的方向上完全故障，没有力矩输出，ε_i＞0意味着并不存在第i轴方向完全故障的情况；

步骤三：期望的角速度轨迹设计

由(2)式改写为：

$\stackrel{\&CenterDot;}{q}=\frac{1}{2}F\left(q\right){\mathrm{\&omega;}}_b---\left(12\right)$

这里q为四元数，Crassidis&Markley提出角速度ω_b是由输入力矩控制T_w的，与此同时ω_b控制姿态q；令

${\mathrm{\&omega;}}_b=\frac{1}{{\mathrm{\&lambda;}}_1}{F}^{-1}\left(q\right)q_e---\left(13\right)$

这里q_e是(3)式的表达形式，结合(12)给出角度位置：

${\mathrm{\&lambda;}}_1\stackrel{.}{q}=q_e,---\left(14\right)$

由于(14)是一阶微分方程，式子中的λ₁＞0，因而，姿态四元数q收敛到期望的姿态四元数q_d；这就意味着式子(13)作为航天器期望的角速度，表达式如下：

$A_{\mathrm{bd}}{\mathrm{\&omega;}}_d^D=-\frac{1}{{\mathrm{\&lambda;}}_1}{F}^{-1}\left(q\right)q_e---15\left(7\right)$

步骤四：控制器输出力矩的设计

考虑当反作用飞轮的输出力矩存在部分损失，则通过控制器设计得到的指令力矩具有一定的容错能力；为了得到在线的故障信息，通过式子(8)中发生故障时的动力学表达式得到带有在线估计故障因子的动力学方程，由于干扰力矩T_d相对于控制力矩而言是小量，这里进行忽略：

$I_b{\stackrel{\stackrel{.}{^}}{\mathrm{\&omega;}}}_b=\hat{M}-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v---\left(16\right)$

这里 ${\widehat{\mathrm{\&omega;}}}_b={\left[\begin{array}{ccc}{\widehat{\mathrm{\&omega;}}}_{b,1}& {\widehat{\mathrm{\&omega;}}}_{b,2}& {\widehat{\mathrm{\&omega;}}}_{b,3}\end{array}\right]}^T\&Element;R^3$ 是估计的航天器本体角速度， $\hat{E}=\mathrm{diag}\left[\begin{array}{ccc}{\hat{e}}_1\left(t\right)& {\hat{e}}_2\left(t\right)& {\hat{e}}_3\left(t\right)\end{array}\right]\&Element;R^{3\&times;3}$ 是三轴方向上的故障系数的估计矩阵，中含有估计量观测的故障动力学的输入T_v＝[T_v1 T_v2 T_v3]^T∈R³将会接下来给出；

定义观测误差向量 ${\widetilde{\mathrm{\&omega;}}}_b={\left[\begin{array}{ccc}{\widetilde{\mathrm{\&omega;}}}_{b,1}& {\widetilde{\mathrm{\&omega;}}}_{b,2}& {\widetilde{\mathrm{\&omega;}}}_{b,3}\end{array}\right]}^T\&Element;R^3,$ 这里 ${\widetilde{\mathrm{\&omega;}}}_{b,i}={\widehat{\mathrm{\&omega;}}}_{b,i}-{\mathrm{\&omega;}}_{b,i},$ i∈{1,2,3}，为了进行角速度误差的补偿，选择力矩然后估计误差动力学方程表示如下：

$I_b{\stackrel{\stackrel{.}{~}}{\mathrm{\&omega;}}}_b=\tilde{E}{T}_v-k_{u}E{\widetilde{\mathrm{\&omega;}}}_b-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\mathrm{\&omega;}}_b^{\&times;}\mathrm{EH}---\left(17\right)$

这里 $\tilde{E}=\mathrm{diag}\left[\begin{array}{ccc}{\tilde{e}}_1& {\tilde{e}}_2& {\tilde{e}}_3\end{array}\right]\&Element;R^{3\&times;3}$ 是故障系数的估计误差矩阵，其中k_u是正常数；设计一个输入T_v，使得估计的航天器角速度ω_b达到期望的角速度而观测的角速度和期望角速度之间的误差向量定义为： ${\widehat{\mathrm{\&omega;}}}_e={\left[\begin{array}{ccc}{\widehat{\mathrm{\&omega;}}}_{e,1}& {\widehat{\mathrm{\&omega;}}}_{e,2}& {\widehat{\mathrm{\&omega;}}}_{e,3}\end{array}\right]}^T\&Element;R^3,$ 其中i∈{1,2,3}，使用式子(16)，估计与期望的误差动力学方程写成：

$I_b{\stackrel{\stackrel{.}{~}}{\mathrm{\&omega;}}}_e=\tilde{M}{-\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{.}{\mathrm{\&omega;}}}_d^D---\left(18\right)$

定理1：考虑等式(3)中带有故障因子的动力学，输出力矩设计为：

$T_w=T_v+k_u{\widetilde{\mathrm{\&omega;}}}_b---\left(19\right)$

这里是观测器方程的输入，k_v是一个正常数，故障观测器为：

$\begin{array}{c}{\stackrel{\stackrel{.}{^}}{e}}_i={\mathrm{Proj}}_{[{\mathrm{\&epsiv;}}_i,1]}\{-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\}\\ =\left\{\begin{array}{cc}0,& \mathrm{if}{\hat{e}}_i={\mathrm{\&epsiv;}}_i,-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\&le;0\mathrm{or}{\hat{e}}_i=1,-{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i}\&GreaterEqual;0\\ -{\mathrm{\&alpha;}}_i{\widetilde{\mathrm{\&omega;}}}_{b,i}{T}_{v_i},& \mathrm{otherwise}\end{array}\right.\end{array}---\left(20\right)$

这里的α_i是正自适应增益，并且这个投影算子Proj{·}用来在参数边界内保持参数的估计；然后，每个角速度ω_i都渐进地收敛到期望的角速度ω_d,i，尽管有效性部分损失，但是期望的输出力矩能有效输出；

证明：如下构造李雅普诺夫方程

$V=\frac{1}{2}{\widehat{\mathrm{\&omega;}}}_e^T{I}_b{\widehat{\mathrm{\&omega;}}}_e+\frac{1}{2}{\widetilde{\mathrm{\&omega;}}}_b^T{I}_b{\widetilde{\mathrm{\&omega;}}}_b^T+\frac{1}{2}\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i^2\left(t\right)}{{\mathrm{\&alpha;}}_i},---\left(21\right)$

从方程式(17)和(18)而言，李雅普诺夫方程的时间导数满足

$\begin{array}{c}\stackrel{.}{V}={\widehat{\mathrm{\&omega;}}}_e^T(M-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{.}{\mathrm{\&omega;}}}_d^D)+{\widetilde{\mathrm{\&omega;}}}_b^T(\tilde{E}{T}_v-k_{u}E{\widetilde{\mathrm{\&omega;}}}_b-{\widehat{\mathrm{\&omega;}}}_b^{\&times;}\hat{E}H+{\mathrm{\&omega;}}_b^{\&times;}\mathrm{EH})+\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\stackrel{.}{~}}{e}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\\ \&le;-k_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u{\widetilde{\mathrm{\&omega;}}}_b^{T}E{\widetilde{\mathrm{\&omega;}}}_b+{\widetilde{\mathrm{\&omega;}}}_b^T\tilde{E}{T}_v+\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\stackrel{.}{~}}{e}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\end{array}---\left(22\right)$

由于e_i是未知的常数，基于投影算子，得到以下不等式

$\frac{{\tilde{e}}_i\left(t\right){\stackrel{\stackrel{.}{~}}{e}}_i\left(t\right)}{{\mathrm{\&alpha;}}_i}\&le;-{\tilde{e}}_i\left(t\right){\widetilde{\mathrm{\&omega;}}}_i{v}_i---\left(23\right)$

令θ＝ζ_min{E}，这里ζ_min{.}算子代表一个矩阵的最小特征值，因为e_i满足0＜ε_i≤e_bi≤1，所以θ＞0；通过采用自适应律式子(18)和(22)进一步得到

$\begin{array}{c}\stackrel{.}{V}\&le;-k_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u\mathrm{\&theta;}{\left|\right|\widetilde{\mathrm{\&omega;}}\left|\right|}^2+\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}({\widetilde{\mathrm{\&omega;}}}_i{\tilde{e}}_i{v}_i+\frac{{\tilde{e}}_i\left(t\right)\stackrel{\stackrel{.}{~}}{e}\left(t\right)}{{\mathrm{\&alpha;}}_i})\\ \&le;-k_v{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left|\right|}^2-k_u\mathrm{\&theta;}{\left|\right|\widetilde{\mathrm{\&omega;}}\left|\right|}^2\end{array}---\left(24\right)$

因此，是负半定的，得到和将从0积分到∞，获得

$V\left(0\right)-V(\&infin;)\&GreaterEqual;k_v\underset{0}{\overset{\&infin;}{\&Integral;}}{\left|\right|{\widehat{\mathrm{\&omega;}}}_e\left(\mathrm{\&eta;}\right)\left|\right|}^2\mathrm{d\&eta;}+k_u\underset{0}{\overset{\&infin;}{\&Integral;}}{\left|\right|{\widetilde{\mathrm{\&omega;}}}_b\left(\mathrm{\&eta;}\right)\left|\right|}^2\mathrm{d\&eta;}---\left(25\right)$

因为以上不等式左边的项是有界的，所以得到和除此而外，通过式子(18)和(19)容易证明和因此，通过Barbalat引理得到以下等式这个能推导出尽管存在三轴力矩的部分损失，但是设计的姿态控制力矩得以保持，并使得航天器的角速度渐进收敛到期望的值；总之，在出现三轴输出力矩损失的情况下，能获得以上设计的姿态控制力矩；

步骤五：数值仿真

为了证明上述方案的有效性，下面通过数值仿真，将上述容错控制方案与传统PD控制方案相比较，刚体航天器本体的惯量矩阵为I_b＝diag(295 130 210)(kg.m²),，假设航天器上装有四个飞轮，飞轮组的惯量阵为：

I_w＝diag(0.01044 0.01044 0.01044 0.01044)(kg·m²)

其安装采用四斜装构型，安装阵为

$C=\frac{1}{2}\left[\begin{array}{cccc}1& 1& -1& -1\\ -1& 1& 1& -1\\ \sqrt{2}& \sqrt{2}& \sqrt{2}& \sqrt{2}\end{array}\right]---\left(26\right)$

初始时刻本体坐标系相对于惯性坐标系的姿态四元数为

$q\left(0\right)=\begin{array}{}\end{array}\left[\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right],---\left(27\right)$

期望坐标系相对于惯性坐标系的四元数为

$q_d=\left[\begin{array}{c}0.93\\ 0.1\\ -0.2\\ 0.3\end{array}\right]---\left(28\right)$

在控制器设计时，忽略了外部干扰力矩，为了更加符合实际的空间环境，这里加入外部干扰力矩T_d，假设T_d为周期变化形式

$T_d=\left[\begin{array}{c}\text{0.3cos}\left(0.01t\right)+0.1\\ \text{0.15sin}\left(0.02t\right)+0.3\cos \left(0.025t\right)\\ \text{0.3sin}\left(0.01t\right)+0.1\end{array}\right]\mathrm{Nm}---\left(29\right)$

a)无故障条件仿真:

在没有故障的条件下，采用传统的PD控制方法，控制参数设定如下：

k₁＝50,k₂＝50

其容错控制方法，控制参数设定如下：有效性故障因子的初始估计自适应律参数ε_i＝0.4,{i＝(1,2,3)}，α_i＝0.5,{i＝(1,2,3)}，k_u＝30，k_v＝8；

b)故障条件的仿真:

在仿真中设定以下故障情况

$\left\{\begin{array}{cc}{e}_y=0.2& t\&GreaterEqual;10\\ e_z=0.5& t\&GreaterEqual;20\end{array}\right.$

这个表明在y轴方向10s之后控制能力损失了80％，z轴方向在20s之后损失了50％，控制方法剩余的控制参数没有改变。