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

摘要

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

主权项

一种基于鲁棒自适应的带有飞轮的航天器容错控制方法，其特征在于：该方法具体步骤如下：涉及的坐标系有：惯性坐标系，这里取地心赤道惯性坐标系作为参考系，原点固联于地心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 }{\omega }_b+q_0{\omega }_b)\\ {\stackrel{.}{q}}_0=-\frac{1}{2}{\hat{q}}^T{\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)$期望坐标系相对于惯性系的角速度在期望坐标系下表示为则在本体系下表示的角速度误差为：${\omega }_e={\omega }_b-A_{\mathrm{bd}}{\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{\omega }_e\\ {\stackrel{\stackrel{.}{^}}{q}}_e=\frac{1}{2}({\hat{q}}_e^{\times }+q_{e0}{E}_3){\omega }_e\end{array}\right.---\left(6\right)$进一步获得以下等式：${\stackrel{.}{\omega }}_b={\stackrel{.}{\omega }}_e-{{\omega }_{\mathrm{bd}}^{\times }A}_{\mathrm{bd}}{\omega }_d^D---\left(7\right)$步骤二：航天器动力学方程的建立假设航天器是刚体航天器，不存在柔性附件，则其带有反作用飞轮的动力学模型表示如下：$M=-{\omega }_b^{\times }{I}_b{\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{·}{q}=\frac{1}{2}F\left(q\right){\omega }_b---\left(12\right)$这里q为四元数，Crassidis&Markley提出角速度ω_b是由输入力矩控制T_w的，与此同时ω_b控制姿态q；令${\omega }_b=\frac{1}{{\lambda }_1}{F}^{-1}\left(q\right)q_e---\left(13\right)$这里q_e是(3)式的表达形式，结合(12)给出角度位置：${\lambda }_1\stackrel{.}{q}=q_e,---\left(14\right)$由于(14)是一阶微分方程，式子中的λ₁＞0，因而，姿态四元数q收敛到期望的姿态四元数q_d；这就意味着式子(13)作为航天器期望的角速度，表达式如下：$A_{\mathrm{bd}}{\omega }_d^D=-\frac{1}{{\lambda }_1}{F}^{-1}\left(q\right)q_e---15\left(7\right)$步骤四：控制器输出力矩的设计考虑当反作用飞轮的输出力矩存在部分损失，则通过控制器设计得到的指令力矩具有一定的容错能力；为了得到在线的故障信息，通过式子(8)中发生故障时的动力学表达式得到带有在线估计故障因子的动力学方程，由于干扰力矩T_d相对于控制力矩而言是小量，这里进行忽略：$I_b{\stackrel{\stackrel{.}{^}}{\omega }}_b=\hat{M}-{\hat{\omega }}_b^{\times }\hat{E}H+\hat{E}{T}_v---\left(16\right)$这里${\hat{\omega }}_b={\left[\begin{array}{ccc}{\hat{\omega }}_{b,1}& {\hat{\omega }}_{b,2}& {\hat{\omega }}_{b,3}\end{array}\right]}^T\in 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]\in R^{3\times 3}$是三轴方向上的故障系数的估计矩阵，中含有估计量观测的故障动力学的输入T_v＝[T_v1 T_v2 T_v3]^T∈R³将会接下来给出；定义观测误差向量${\tilde{\omega }}_b={\left[\begin{array}{ccc}{\tilde{\omega }}_{b,1}& {\tilde{\omega }}_{b,2}& {\tilde{\omega }}_{b,3}\end{array}\right]}^T\in R^3,$这里${\tilde{\omega }}_{b,i}={\hat{\omega }}_{b,i}-{\omega }_{b,i},$i∈{1,2,3}，为了进行角速度误差的补偿，选择力矩然后估计误差动力学方程表示如下：$I_b{\stackrel{\stackrel{.}{~}}{\omega }}_b=\tilde{E}{T}_v-k_{u}E{\tilde{\omega }}_b-{\hat{\omega }}_b^{\times }\hat{E}H+{\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]\in R^{3\times 3}$是故障系数的估计误差矩阵，其中k_u是正常数；设计一个输入T_v，使得估计的航天器角速度ω_b达到期望的角速度而观测的角速度和期望角速度之间的误差向量定义为：${\hat{\omega }}_e={\left[\begin{array}{ccc}{\hat{\omega }}_{e,1}& {\hat{\omega }}_{e,2}& {\hat{\omega }}_{e,3}\end{array}\right]}^T\in R^3,$其中i∈{1,2,3}，使用式子(16)，估计与期望的误差动力学方程写成：$I_b{\stackrel{\stackrel{.}{~}}{\omega }}_e=\tilde{M}{-\hat{\omega }}_b^{\times }\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{.}{\omega }}_d^D---\left(18\right)$定理1：考虑等式(3)中带有故障因子的动力学，输出力矩设计为：$T_w=T_v+k_u{\tilde{\omega }}_b---\left(19\right)$这里是观测器方程的输入，k_v是一个正常数，故障观测器为：$\begin{array}{c}{\stackrel{\stackrel{.}{^}}{e}}_i={\mathrm{Proj}}_{[{ϵ}_i,1]}\{-{\alpha }_i{\tilde{\omega }}_{b,i}{T}_{v_i}\}\\ =\left\{\begin{array}{cc}0,& \mathrm{if}{\hat{e}}_i={ϵ}_i,-{\alpha }_i{\tilde{\omega }}_{b,i}{T}_{v_i}\le 0\mathrm{or}{\hat{e}}_i=1,-{\alpha }_i{\tilde{\omega }}_{b,i}{T}_{v_i}\ge 0\\ -{\alpha }_i{\tilde{\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}{\hat{\omega }}_e^T{I}_b{\hat{\omega }}_e+\frac{1}{2}{\tilde{\omega }}_b^T{I}_b{\tilde{\omega }}_b^T+\frac{1}{2}\underset{i=1}{\overset{4}{\Sigma }}\frac{{\tilde{e}}_i^2\left(t\right)}{{\alpha }_i},---\left(21\right)$从方程式(17)和(18)而言，李雅普诺夫方程的时间导数满足$\begin{array}{c}\stackrel{.}{V}={\hat{\omega }}_e^T(M-{\hat{\omega }}_b^{\times }\hat{E}H+\hat{E}{T}_v-I_b{\stackrel{.}{\omega }}_d^D)+{\tilde{\omega }}_b^T(\tilde{E}{T}_v-k_{u}E{\tilde{\omega }}_b-{\hat{\omega }}_b^{\times }\hat{E}H+{\omega }_b^{\times }\mathrm{EH})+\underset{i=1}{\overset{4}{\Sigma }}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\stackrel{.}{~}}{e}}_i\left(t\right)}{{\alpha }_i}\\ \le -k_v{\left|\right|{\hat{\omega }}_e\left|\right|}^2-k_u{\tilde{\omega }}_b^{T}E{\tilde{\omega }}_b+{\tilde{\omega }}_b^T\tilde{E}{T}_v+\underset{i=1}{\overset{4}{\Sigma }}\frac{{\tilde{e}}_i\left(t\right){\stackrel{\stackrel{.}{~}}{e}}_i\left(t\right)}{{\alpha }_i}\end{array}---\left(22\right)$由于e_i是未知的常数，基于投影算子，得到以下不等式$\frac{{\tilde{e}}_i\left(t\right){\stackrel{\stackrel{.}{~}}{e}}_i\left(t\right)}{{\alpha }_i}\le -{\tilde{e}}_i\left(t\right){\tilde{\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|{\hat{\omega }}_e\left|\right|}^2-k_u\theta {\left|\right|\tilde{\omega }\left|\right|}^2+\underset{i=1}{\overset{3}{\Sigma }}({\tilde{\omega }}_i{\tilde{e}}_i{v}_i+\frac{{\tilde{e}}_i\left(t\right)\stackrel{\stackrel{.}{~}}{e}\left(t\right)}{{\alpha }_i})\\ \le -k_v{\left|\right|{\hat{\omega }}_e\left|\right|}^2-k_u\theta {\left|\right|\tilde{\omega }\left|\right|}^2\end{array}---\left(24\right)$因此，是负半定的，得到和将从0积分到∞，获得$V\left(0\right)-V(\infty )\ge k_v\underset{0}{\overset{\infty }{\int }}{\left|\right|{\hat{\omega }}_e\left(\eta \right)\left|\right|}^2\mathrm{d\eta }+k_u\underset{0}{\overset{\infty }{\int }}{\left|\right|{\tilde{\omega }}_b\left(\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\ge 10\\ e_z=0.5& t\ge 20\end{array}\right.$这个表明在y轴方向10s之后控制能力损失了80％，z轴方向在20s之后损失了50％，控制方法剩余的控制参数没有改变。