CN105607485A - 基于故障特征模型的挠性充液卫星姿态自适应容错控制方法 - Google Patents

Abstract

本发明属于自动控制技术领域，涉及一种基于故障特征模型的挠性充液卫星姿态自适应容错控制方法。本发明在系统存在液体晃动、挠性振动和外部干扰的情况下，针对执行器故障设计了基于故障特征模型的容错控制器，步骤为：(1)根据挠性充液卫星姿态动力学方程建立其故障特征模型；(2)利用参数估计算法对故障特征模型的系数进行辨识，得到故障特征模型系数的估计值；(3)根据故障特征模型系数的估计值设计自适应容错控制器；(4)由自适应容错控制器得到控制力矩，并将其施加到卫星姿态控制系统以控制卫星的姿态角。使用该方法能保障故障发生前、后的控制精度均很高，且达到稳态的时间短。

Description

基于故障特征模型的挠性充液卫星姿态自适应容错控制方法

技术领域

本发明涉及自动控制技术领域，具体涉及一种基于故障特征模型的挠性充液卫星姿态自适应容错控制方法。

背景技术

随着航天技术的不断发展，卫星承担的任务越来越多，结构变得越来越复杂,同时卫星的在轨寿命要求大大延长，姿态控制精度要求不断提高。由于大型复杂卫星的液体燃料携带量加大，弹性附件变大，液体晃动和弹性附件振动对航天器姿态的影响越来越严重，液体晃动的抑制与挠性体振动的抑制问题一并成为姿态控制系统设计的难点之一。卫星长时间运行于真空、失重、高低温和强辐射的环境中，执行机构不可避免会发生故障，而执行器发生故障,轻则导致精度、性能降低,重则造成卫星失效。为保证卫星在执行机构故障情况下能够稳定可靠运行，设计一种容错控制方法来完成挠性充液卫星的高稳定度、高精度控制便显得尤为重要。现有技术尚未公开有有效的挠性充液卫星姿态容错控制的方法。

吴宏鑫院士从工程应用的角度出发，提出了特征建模的思想，为参数和阶数均未知的高阶复杂系统的建模提供了一种新思路，为一些高阶复杂系统实现低阶控制器设计打下了理论基础。基于特征模型的全系数自适应控制方法在许多工业控制领域中已取得成功应用，特别值得一提的是：该方法的理论思想和工程要义被创造性地应用于飞船返回再入控制，其开伞精度达到国际先进水平。尽管基于特征模型的自适应控制方法在工程上方便易行，但是，将这一思想方法应用于容错控制器的设计尚无文献记载与公开的技术方法。

发明内容

本发明的目的是提供一种实时性好，控制精度高，且能保证闭环系统稳定性的基于故障特征模型的挠性充液卫星姿态自适应容错控制方法，以提高挠性充液卫星姿态控制系统的可靠性，防止因执行器故障引起的卫星姿态失控。

实现本发明目的的技术方案是：一种基于故障特征模型的挠性充液卫星姿态自适应容错控制方法，其通过以下步骤实现：

第一步，根据(1)式所给出的挠性充液卫星姿态动力学方程，建立由(2)式所描述的故障特征模型：

$\left\{\begin{array}{c}I\stackrel{\·}{\mathrm{\ω}}+\mathrm{\ω}\×I\mathrm{\ω}+C_0\stackrel{\·\·}{\mathrm{\η}}+G\stackrel{\·\·}{q}=T_d+u+f,\\ \stackrel{\·\·}{\mathrm{\η}}+2\mathrm{\ξ}\mathrm{\Λ}\stackrel{\·}{\mathrm{\η}}+{\mathrm{\Λ}}^2\mathrm{\η}+{C_0}^T\stackrel{\·}{\mathrm{\ω}}=0,\\ \stackrel{\·\·}{q}+{\mathrm{\Λ}}_{2}q+C_1\stackrel{\·}{\mathrm{\ω}}+C_2\mathrm{\α}=0,\\ y=\mathrm{\α},\end{array}\right.---\left(1\right)$

(1)式中：为星体的姿态角矢量，为滚动角，θ为俯仰角，ψ为偏航角；u＝[u_1,u₂,u₃]^T为作用在星体上的控制力矩，是系统的输入；y为系统的输出；为星体的角速度矢量；η＝[η₁,η₂,…,η_m]^T为挠性体模态坐标矢量：q＝[η₁,η₂,…,η_3n]^T为等效摆晃动模态坐标矢量，n为等效摆的个数；T_d为作用在星体上的干扰力矩，为3维列矢量；I为卫星的转动惯量矩阵，为3阶方阵；C0为挠性体振动运动与星体转动运动耦合系数矩阵，其行数为3、列数为m；ξ为挠性体振动模态阻尼矩阵，为m阶对角阵；Λ为挠性体振动模态频率矩阵，为m阶方阵；Λ₂为等效摆振动模态频率矩阵，为3n阶对角阵；G，C₁，C₂为等效摆振动运动与星体转动运动耦合系数矩阵，G的行数为3、列数为3n，C₁，C₂的行数均为3n、列数均为3；f为故障项，其表达式为：

为关于α的矢量函数，为故障函数，diag表示对角矩阵，β_i(t-T_i)是故障开关函数，定义为

${\mathrm{\&beta;}}_i(t-T_i)=\{\begin{array}{cc}0,& t<T_i,\\ 1,& t\&GreaterEqual;T_i,\end{array},i=1,2,3,$

T_i为第i个执行器故障发生的时刻；

$\begin{array}{c}{y}_i(k+1)\\ =a_{i1}\left(k\right)y_i\left(k\right)+a_{i2}\left(k\right)y_i(k-1)+b_{i1}\left(k\right){\stackrel{\&OverBar;}{u}}_1\left(k\right)+b_{i2}\left(k\right){\stackrel{\&OverBar;}{u}}_2\left(k\right)+b_{i3}\left(k\right){\stackrel{\&OverBar;}{u}}_3\left(k\right)+E_i\left(k\right),\end{array}---\left(2\right)$

(2)式中：i＝1,2,3；y₁(k)为第k个采样时刻的滚动角y₂(k)为第k个采样时刻的俯仰角θ，y₃(k)为第k个采样时刻的偏航角ψ；为第k个采样时刻故障特征模型的总输入，包括：第k个采样时刻作用在星体上的控制力矩u、干扰力矩T_d、第i个执行器发生故障时故障函数的量值以及由轨道坐标系至卫星本体坐标系进行坐标变换所涉及的常量；E_i(k)为建模误差；a_i1(k)，a_i2(k)，b_i1(k)，b_i2(k)，b_i3(k)为故障特征模型的系数；

第二步，利用参数估计算法对由第一步得到的故障特征模型的系数a_i1(k)，a_i2(k)，b_ij(k)(i,j＝1,2,3)进行辨识得到a_i1(k)，a_i2(k)，b_ij(k)的估计值，分别为：

第三步，根据第二步得到的a_i1(k)，a_i2(k)，b_ij(k)的估计值设计容错控制器：

u_i(k)＝u_gi(k)+u_Ii(k)，i＝1,2,3，(3)

(3)式中：u_i(k)为第k个采样时刻控制力矩u的第i个分量，u_gi(k)为黄金分割自适应控制律，设计为

$u_{gi}\left(k\right)=\frac{-b_{0i}\&lsqb;l_1{\hat{a}}_{i1}{e}_i\left(k\right)+l_2{\hat{a}}_{i2}\left(k\right)e_i(k-1)\&rsqb;}{{\hat{b}}_{ii}\left(k\right)+{\mathrm{\&lambda;}}_{0i}},i=1,2,3,---\left(4\right)$

(4)式中：b_0i，λ_0i均为大于0的控制增益调节参数，视控制效果可适当增减；l₁＝0.382，l₂＝0.618，e_i(k)＝y_i(k)-y_ri，y_ri为姿态角的期望值，u_Ii(k)为逻辑积分控制律；

第四步，将第三步设计的u_i(k)施加到由(1)式所描述的系统中，以控制卫星的姿态角。

本发明的有益效果是：

①首次对含有故障的控制系统建立了特征模型，即故障特征模型，为挠性充液卫星姿态的容错控制打下了基础。

②基于故障特征模型设计的黄金分割自适应容错控制器不依赖于特征模型的控制增益矩阵，从而避免了对多输入多输出耦合系统通过参数辨识可能使控制增益矩阵奇异的问题。

③由于在黄金分割自适应控制律的分子中引入了控制增益调节参数，当由特征模型较小的辨识系数得到的较小控制量，达不到控制要求时，上述控制增益调节参数起到了放大作用。

④该方法能有效克服挠性振动、液体晃动和外干扰对卫星姿态控制的影响。

⑤使用该方法能保障故障发生前、后的控制精度均很高，且达到稳态的时间短。

⑥本发明提出的控制方法保证了闭环系统的稳定性。

附图说明

图1为本发明实施例中对滚动角期望值y_r1的跟踪曲线图；

图2为本发明实施例中对俯仰角期望值y_r2的跟踪曲线图；

图3为本发明实施例中对偏航角期望值y_r3的跟踪曲线图。

具体实施方式

为使本发明的内容和技术方案更加清楚明白，以下结合实施例对本发明进一步详细说明。

实施例：

本基于故障特征模型的挠性充液卫星姿态自适应容错控制方法的具体实施步骤：

第一步，根据(1)式所给出的挠性充液卫星姿态动力学方程，建立由(2)式所描述的故障特征模型：

$\left\{\begin{array}{c}I\stackrel{\&CenterDot;}{\mathrm{\&omega;}}+\mathrm{\&omega;}\&times;I\mathrm{\&omega;}+C_0\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}+G\stackrel{\&CenterDot;\&CenterDot;}{q}=T_d+u+f,\\ \stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}+2\mathrm{\&xi;}\mathrm{\&Lambda;}\stackrel{\&CenterDot;}{\mathrm{\&eta;}}+{\mathrm{\&Lambda;}}^2\mathrm{\&eta;}+{C_0}^T\stackrel{\&CenterDot;}{\mathrm{\&omega;}}=0,\\ \stackrel{\&CenterDot;\&CenterDot;}{q}+{\mathrm{\&Lambda;}}_{2}q+C_1\stackrel{\&CenterDot;}{\mathrm{\&omega;}}+C_2\mathrm{\&alpha;}=0,\\ y=\mathrm{\&alpha;},\end{array}\right.---\left(1\right)$

(1)式中：为星体的姿态角矢量，为滚动角，θ为俯仰角，ψ为偏航角；u＝[u_1,u₂,u₃]^T为作用在星体上的控制力矩，是系统的输入；y为系统的输出；为星体的角速度矢量；η＝[η₁,η₂,…,η_m]^T为挠性体模态坐标矢量：q＝[η₁,η₂,…,η_3n]^T为等效摆晃动模态坐标矢量，n为等效摆的个数；T_d为作用在星体上的干扰力矩，为3维列矢量；I为卫星的转动惯量矩阵，表达式为：

$I=\left[\begin{array}{ccc}{I}_x& -I_{xy}& -I_{xz}\\ -I_{xy}& I_y& -I_{yz}\\ -I_{xz}& -I_{yz}& I_z\end{array}\right],---(1-1)$

(1—1)式中：矩阵I的主对角元素I_x，I_y，I_z为卫星的三个主惯量，其它元素为惯量积；C0为挠性体振动运动与星体转动运动耦合系数矩阵，其行数为3、列数为m；ξ为挠性体振动模态阻尼矩阵，为m阶对角阵；Λ为挠性体振动模态频率矩阵，为m阶方阵；Λ₂为等效摆振动模态频率矩阵，为3n阶对角阵；G，C₁，C₂为等效摆振动运动与星体转动运动耦合系数矩阵，G的行数为3、列数为3n，C₁，C₂的行数均为3n、列数均为3；f为故障项，其表达式为：

为关于α的矢量函数，为故障函数，diag表示对角矩阵，β_i(t-T_i)是故障开关函数，定义为

${\mathrm{\&beta;}}_i(t-T_i)=\{\begin{array}{cc}0,& t<T_i,\\ 1,& t\&GreaterEqual;T_i,\end{array},i=1,2,3,$

T_i为第i个执行器故障发生的时刻；

$\begin{array}{c}{y}_i(k+1)\\ =a_{i1}\left(k\right)y_i\left(k\right)+a_{i2}\left(k\right)y_i(k-1)+b_{i1}\left(k\right){\stackrel{\&OverBar;}{u}}_1\left(k\right)+b_{i2}\left(k\right){\stackrel{\&OverBar;}{u}}_2\left(k\right)+b_{i3}\left(k\right){\stackrel{\&OverBar;}{u}}_3\left(k\right)+E_i\left(k\right),\end{array}---\left(2\right)$

(2)式中：i＝1,2,3；y₁(k)为第k个采样时刻的滚动角y₂(k)为第k个采样时刻的俯仰角θ，y₃(k)为第k个采样时刻的偏航角ψ；为第k个采样时刻故障特征模型的总输入，包括：第k个采样时刻作用在星体上的控制力矩u、干扰力矩T_d、第i个执行器发生故障时故障函数的量值以及由轨道坐标系至卫星本体坐标系进行坐标变换所涉及的常量；E_i(k)为建模误差；a_i1(k)，a_i2(k)，b_i1(k)，b_i2(k)，b_i3(k)为故障特征模型的系数。

(2)式所给出的特征模型的建模步骤如下：

先利用轨道坐标系至卫星本体坐标系的转换关系将(1)式化为(2—1)式：

$\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&Theta;}}+C^1\stackrel{\&CenterDot;}{\mathrm{\&Theta;}}+K^1\mathrm{\&Theta;}=T^1,---(2-1)$

(2—1)式中：Θ＝[α^т,η^т,q^т]^T,C¹＝M^-1L,K¹＝M^-1K,T¹＝M^-1T，M、L、K、T分别为：

$\begin{array}{cc}M=\left[\begin{array}{ccc}I& C_0& G\\ {C_0}^T& E_m& O_{m\&times;3n}\\ C_1& O_{3n\&times;m}& E_{3n}\end{array}\right],& L=\left[\begin{array}{ccc}F& O_{3\&times;m}& O_{3\&times;3n}\\ C_3& 2\mathrm{\&xi;}\mathrm{\&Lambda;}& O_{m\&times;3n}\\ C_4& O_{3n\&times;m}& O_{3n\&times;3n}\end{array}\right]\end{array},$

$\begin{array}{cc}K=\left[\begin{array}{ccc}W& O_{3\&times;m}& O_{3\&times;3n}\\ O_{m\&times;3}& {\mathrm{\&Lambda;}}^2& O_{m\&times;3n}\\ C_2& O_{3n\&times;m}& {\mathrm{\&Lambda;}}_2\end{array}\right],& T=\left[\begin{array}{c}{T}_d+u+f+T_{d0}\\ O_{m\&times;1}\\ O_{3n\&times;1}\end{array}\right]\end{array},$

其中，E_m、E_3n分别为m阶、n阶单位矩阵，O_i×j为i行、j列的零矩阵，

$\begin{array}{ccc}{C}_3={C_0}^T\left[\begin{array}{ccc}0& 0& -1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right]{\mathrm{\&omega;}}_0,& C_4=C_1\left[\begin{array}{ccc}0& 0& -1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right]{\mathrm{\&omega;}}_0,& T_{d0}=\left[\begin{array}{c}4{{\mathrm{\&omega;}}_0}^2{I}_{yz}\\ -3{{\mathrm{\&omega;}}_0}^2{I}_{xz}\\ -{{\mathrm{\&omega;}}_0}^2{I}_{xy}\end{array}\right]\end{array},$

$\begin{array}{cc}F={\mathrm{\&omega;}}_0\left[\begin{array}{ccc}0& 2I_{yz}& I_y-I_z-I_x\\ -2I_{yz}& 0& 2I_{xy}\\ -I_y+I_z+I_x& -2I_{xy}& 0\end{array}\right],& W={{\mathrm{\&omega;}}_0}^2\left[\begin{array}{ccc}4(I_y-I_z)& 3I_{xy}& -I_{xz}\\ 4I_{xy}& 3(I_x-I_z)& I_{yz}\\ -4I_{xz}& -3I_{yz}& I_y-I_x\end{array}\right]\end{array},$

ω₀为轨道坐标系相对空间惯性坐标系的角速度。

再将(2—1)式化为状态方程(2—2)：

$\{\begin{array}{c}\stackrel{\&CenterDot;}{X}=AX+BU,\\ Y=CX,\end{array}---(2-2)$

(2—2)式中： $\begin{array}{ccc}A={\left[\begin{array}{cc}-C^1& -K^1\\ E_p& O_{p\&times;p}\end{array}\right]}_{2p\&times;2p},& B=\left[\begin{array}{c}{M}^{-1}\\ O_{p\&times;p}\end{array}\right],& U=\left[\begin{array}{c}\stackrel{\&OverBar;}{u}\\ O_{m\&times;1}\\ O_{3n\&times;1}\end{array}\right]\end{array},$ C＝[O_3×pE₃O_3×(p-3)]_3×2p， $X={\left[\begin{array}{cc}\stackrel{\&CenterDot;}{\mathrm{\&Theta;}}& \mathrm{\&Theta;}\end{array}\right]}^T,$ p＝3+m+3n，E_p、E₃分别为p阶、3阶单位矩阵，

作拉普拉斯变换，将状态方程(2—2)化为输入输出形式(2—3)：

${\hat{y}}_i\left(s\right)=g_{i1}\left(s\right){\widehat{\stackrel{\&OverBar;}{u}}}_1\left(s\right)+g_{i2}\left(s\right){\widehat{\stackrel{\&OverBar;}{u}}}_2\left(s\right)+g_{i3}\left(s\right){\widehat{\stackrel{\&OverBar;}{u}}}_3\left(s\right),i=1,2,3,---(2-3)$

(2—3)式中:分别为y_i,的拉普拉斯变换，为传递函数矩阵的第i行第j列元素，其中，N_ij(s)＝c_ij，2p-1s^2p-1+c_ij，2p-2s^2p-2+…+c_ij，1s+c_ij，0，i＝1,2,3，j＝1，2，3，M(s)＝s^2p+d_2p-1s^2p-1+…+d₁s+d₀，且c_ij，2p-1，c_ij，2p-2，…，c_ij，0，d_2p-1，…，d₁，d₀是由(2—2)式中矩阵A，B，C的元素确定的常数。

根据(2—3)式得到(2—4)式：

$M\left(s\right){\hat{y}}_i\left(s\right)=\underset{j=1}{\overset{3}{\&Sigma;}}{N}_{ij}\left(s\right){\widehat{\stackrel{\&OverBar;}{u}}}_j\left(s\right),---(2-4)$

对(2—4)式作拉普拉斯逆变换得(2—5)式：

$y_i^{\left(2p\right)}+...+d_1{\stackrel{\&CenterDot;}{y}}_i+d_0{y}_i=\underset{j=1}{\overset{3}{\&Sigma;}}\&lsqb;c_{ij,2p-1}{\stackrel{\&OverBar;}{u}}_j^{(2p-1)}+...+c_{ij,1}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j+c_{ij,0}{\stackrel{\&OverBar;}{u}}_j\&rsqb;,---(2-5)$

当d₁≠0时，将(2—5)式改写成(2—6)式：

$d_1{\stackrel{\&CenterDot;}{y}}_i=-d_0{y}_i+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&OverBar;}{u}}_j+F_i\left(t\right),---(2-6)$

(2—6)式中： $F_i\left(t\right)=-y_i^{\left(2p\right)}-...-d_2{\stackrel{\&CenterDot;\&CenterDot;}{y}}_i+\underset{j=1}{\overset{3}{\&Sigma;}}\&lsqb;c_{ij,2p-1}{\stackrel{\&OverBar;}{u}}_j^{(2p-1)}+...+c_{ij,1}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j\&rsqb;.$

将(2—6)式两边对时间t求导得：

$d_1{\stackrel{\&CenterDot;\&CenterDot;}{y}}_i=-d_0{\stackrel{\&CenterDot;}{y}}_i+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j+{\stackrel{\&CenterDot;}{F}}_i\left(t\right),---(2-7)$

以下用Δt表示采样周期；对(2—6)式、(2—7)式两边进行差分，分别得

$d_1\frac{y_i\left(k\right)-y_i(k-1)}{\mathrm{\&Delta;}t}=-d_0{y}_i\left(k\right)+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&OverBar;}{u}}_j\left(k\right)+F_i\left(k\right),---(2-8)$

$\begin{array}{c}{d}_1\frac{y_i(k+1)-2y_i\left(k\right)+y_i(k-1)}{{\mathrm{\&Delta;t}}^2}\\ =-d_0\frac{y_i\left(k\right)-y_i(k-1)}{\mathrm{\&Delta;}t}+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j\left(k\right)+\frac{F_i\left(k\right)-F_i(k-1)}{\mathrm{\&Delta;}t},\end{array}---(2-9)$

将(2—8)式、(2—9)式两边相加，经整理即得(2)式，且有

$a_{i1}\left(k\right)=2-\mathrm{\&Delta;}t-\frac{d_0}{d_1}\mathrm{\&Delta;}t-\frac{d_0}{d_1}{\mathrm{\&Delta;t}}^2,a_{i2}\left(k\right)=-1+\mathrm{\&Delta;}t+\frac{d_0}{d_1}\mathrm{\&Delta;}t,b_{ij}\left(k\right)=\frac{c_{ij,0}}{d_1}{\mathrm{\&Delta;t}}^2,---(2-10)$

建模误差： $E_i\left(k\right)=\frac{\mathrm{\&Delta;}t\&lsqb;F_i\left(k\right)-F_i(k-1)\&rsqb;+{\mathrm{\&Delta;t}}^2{F}_i\left(k\right)+{\mathrm{\&Delta;t}}^2\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j\left(k\right)}{d_1}.$

当d₁＝0时，将(2—5)式改写成(2—11)式：

${\stackrel{\&CenterDot;}{y}}_i=-d_0{y}_i+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&OverBar;}{u}}_j+F_i\left(t\right),---(2-11)$

(2—11)式中： $F_i\left(t\right)={\stackrel{\&CenterDot;}{y}}_i-y_i^{\left(2p\right)}-...-d_2{\stackrel{\&CenterDot;\&CenterDot;}{y}}_i+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,2p-1}{\stackrel{\&OverBar;}{u}}_j^{(2p-1)}+...+c_{ij,1}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j.$

将(2—11)式两边对时间t求导得：

${\stackrel{\&CenterDot;\&CenterDot;}{y}}_i=-d_0\stackrel{\&CenterDot;}{y}+\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}+{\stackrel{\&CenterDot;}{F}}_i\left(t\right),---(2-12)$

对(2—11)式、(2—12)式两边进行差分，然后两边再相加，经整理即得(2)式，且有

a_i1(k)＝2-Δt-d₀Δt-d₀Δt²，a_i2(k)＝-1+Δt+d₀Δt，b_ij(k)＝c_ij,0Δt²，(2—13)建模误差： $E_i\left(k\right)=\mathrm{\&Delta;}t\&lsqb;F_i\left(k\right)-F_i(k-1)\&rsqb;+{\mathrm{\&Delta;t}}^2{F}_i\left(k\right)+{\mathrm{\&Delta;t}}^2\underset{j=1}{\overset{3}{\&Sigma;}}{c}_{ij,0}{\stackrel{\&CenterDot;}{\stackrel{\&OverBar;}{u}}}_j\left(k\right).$

第二步，利用参数估计算法对第一步得到的故障特征模型的系数a_i1(k)，a_i2(k)，b_ij(k)(i,j＝1,2,3)进行辨识得到a_i1(k)，a_i2(k)，b_ij(k)的估计值，分别为： 一些具体的参数估计算法参见：《自适应滤波、预测与控制》，G.C.古德温，孙贵生著，科学出版社，1992年版，第35页至第86页。在实施例中，取改进的投影算法：

(2—14)式中： ${\widehat{\mathrm{\&theta;}}}_i\left(k\right)={\&lsqb;{\hat{a}}_{i1}\left(k\right),{\hat{a}}_{i2}\left(k\right),{\hat{b}}_{i1}\left(k\right),{\hat{b}}_{i2}\left(k\right),{\hat{b}}_{i3}\left(k\right)\&rsqb;}^T;$ 根据(2—10)和(2—13)式，的初值取为[1.85,-0.9,0.1,0.1,0.1]^T，λ_i1＝λ_i2＝0.1_。

第三步，根据第二步得到的a_i1(k)，a_i2(k)，b_ij(k)的估计值设计容错控制器：

u_i(k)＝u_gi(k)+u_Ii(k)，i＝1,2,3，(3)

(3)式中：u_i(k)为第k个采样时刻控制力矩u的第i个分量，u_gi(k)为黄金分割自适应控制律，设计为

$u_{gi}\left(k\right)=\frac{-b_{0i}\&lsqb;l_1{\hat{a}}_{i1}{e}_i\left(k\right)+l_2{\hat{a}}_{i2}\left(k\right)e_i(k-1)\&rsqb;}{{\hat{b}}_{ii}\left(k\right)+{\mathrm{\&lambda;}}_{0i}},i=1,2,3,---\left(4\right)$

(4)式中：b_0i，λ_0i均为大于0的控制增益调节参数，视控制效果可适当增减；l₁＝0.382，l₂＝0.618，e_i(k)＝y_i(k)-y_ri，y_ri为姿态角的期望值，u_Ii(k)为逻辑积分控制律。在实施例中，取b_0i＝5000，λ_0i＝0.1，y_ri＝0，u_Ii(k)设计为：

u_Ii(k)＝u_Ii(k-1)-k_Iie_i(k)，i＝1,2,3，(4—1)

(4—1)式中：当e_i(k)y_i(k)≥0时，k_Ii＝0.05，当e_i(k)y_i(k)＜0时，k_Ii＝0.01。

第四步，将第三步设计的u_i(k)施加到由(1)式所描述的系统中，以控制卫星的姿态角。在实施例中，(1)式中的各个参数值见文献“挠性充液卫星动力学分析与姿态控制研究”(李英波.上海交通大学博士学位论文,2001,第118页和第119页)，姿态角初始值也与此文献相同，即分别取θ(0)＝0.8°，ψ(0)＝0.5°；故障函数取为 故障发生的时刻T_i＝25秒。

仿真实验结果如图1、图2、图3中的实线所示，图1中的实线是对滚动角的期望值y_r1的跟踪曲线，图2中的实线是对俯仰角θ的期望值y_r2的跟踪曲线，图3中的实线是对偏航角ψ的期望值y_r3的跟踪曲线；图1、图2、图3中的虚线是期望的输出。仿真实验结果表明：执行器发生故障前，采用本发明提出的方法，卫星姿态在很短的时间内达到稳定状态，滚动角稳态控制精度小于等于1.3×10^–5度、俯仰角稳态控制精度小于等于4.2×10^–6度、偏航角稳态控制精度小于等于1.5×10^–5度，滚动角速率稳态控制精度小于1.1×10^–8度/秒、俯仰角速率稳态控制精度小于5.4×10^–8度/秒、偏航角速率稳态控制精度小于1.1×10^–8度/秒；当执行器发生故障时，采用本发明提出的方法，系统可在较短的时间内恢复正常工作，卫星姿态稳定后，滚动角稳态控制精度小于8.2×10^–3度、俯仰角稳态控制精度小于1.5×10^–2度、偏航角稳态控制精度小于1.2×10^–2度，滚动角速率稳态控制精度小于2.5×10^–6度/秒、俯仰角速率稳态控制精度小于8.2×10^–6度/秒、偏航角速率稳态控制精度小于5.0×10^–6度/秒。

综上，通过对挠性充液卫星姿态控制的仿真实验，表明采用本发明不仅能够克服干扰、有效抑制挠性振动和液体晃动，而且具有较强的容错控制能力。

本发明未详细说明部分属于本领域技术人员公知常识。

Claims (1)

1.一种基于故障特征模型的挠性充液卫星姿态自适应容错控制方法，其特征在于通过以下步骤实现：

第一步，根据(1)式所给出的挠性充液卫星姿态动力学方程，建立由(2)式所描述的故障特征模型：

$\left\{\begin{array}{c}I\stackrel{\&CenterDot;}{\mathrm{\&omega;}}+\mathrm{\&omega;}\&times;I\mathrm{\&omega;}+C_0\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}+G\stackrel{\&CenterDot;\&CenterDot;}{q}=T_d+u+f,\\ \stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}+2\mathrm{\&xi;}\mathrm{\&Lambda;}\stackrel{\&CenterDot;}{\mathrm{\&eta;}}+{\mathrm{\&Lambda;}}^2\mathrm{\&eta;}+{C_0}^T\stackrel{\&CenterDot;}{\mathrm{\&omega;}}=0,\\ \stackrel{\&CenterDot;\&CenterDot;}{q}+{\mathrm{\&Lambda;}}_{2}q+C_1\stackrel{\&CenterDot;}{\mathrm{\&omega;}}+C_2\mathrm{\&alpha;}=0,\\ y=\mathrm{\&alpha;},\end{array}\right.---\left(1\right)$

(1)式中：为星体的姿态角矢量，为滚动角，θ为俯仰角，ψ为偏航角；u＝[u₁,u₂,u₃]^T为作用在星体上的控制力矩，是系统的输入；y为系统的输出；为星体的角速度矢量；η＝[η₁,η₂,…,η_m]^T为挠性体模态坐标矢量：q＝[η₁,η₂,…,η_3n]^T为等效摆晃动模态坐标矢量，n为等效摆的个数；T_d为作用在星体上的干扰力矩，为3维列矢量；I为卫星的转动惯量矩阵，为3阶方阵；C₀为挠性体振动运动与星体转动运动耦合系数矩阵，其行数为3、列数为m；ξ为挠性体振动模态阻尼矩阵，为m阶对角阵；Λ为挠性体振动模态频率矩阵，为m阶方阵；Λ₂为等效摆振动模态频率矩阵，为3n阶对角阵；G，C₁，C₂为等效摆振动运动与星体转动运动耦合系数矩阵，G的行数为3、列数为3n，C₁，C₂的行数均为3n、列数均为3；f为故障项，其表达式为：

为关于α的矢量函数，为故障函数，diag表示对角矩阵，β_i(t-T_i)是故障开关函数，定义为

${\mathrm{\&beta;}}_i(t-T_i)=\left\{\begin{array}{c}0,t<T_i,\\ 1,t\&GreaterEqual;T_i,\end{array}\right.,i=1,2,3,$

T_i为第i个执行器故障发生的时刻；

$\begin{array}{c}{y}_i\left(k+1\right)\\ =a_{i1}\left(k\right)y_i\left(k\right)+a_{i2}\left(k\right)y_i(k-1)+b_{i1}\left(k\right){\stackrel{\&OverBar;}{u}}_1\left(k\right)+b_{i2}\left(k\right){\stackrel{\&OverBar;}{u}}_2\left(k\right)+b_{i3}\left(k\right){\stackrel{\&OverBar;}{u}}_3\left(k\right)+E_i\left(k\right),\end{array}---\left(2\right)$

(2)式中：i＝1,2,3；y₁(k)为第k个采样时刻的滚动角y₂(k)为第k个采样时刻的俯仰角θ，y₃(k)为第k个采样时刻的偏航角ψ；为第k个采样时刻故障特征模型的总输入，包括：第k个采样时刻作用在星体上的控制力矩u、干扰力矩T_d、第i个执行器发生故障时故障函数的量值以及由轨道坐标系至卫星本体坐标系进行坐标变换所涉及的常量；E_i(k)为建模误差；a_i1(k)，a_i2(k)，b_i1(k)，b_i2(k)，b_i3(k)为故障特征模型的系数；

第二步，利用参数估计算法对由第一步得到的故障特征模型的系数a_i1(k)，a_i2(k)，b_ij(k)(i,j＝1,2,3)进行辨识得到a_i1(k)，a_i2(k)，b_ij(k)的估计值，分别为：

第三步，根据第二步得到的a_i1(k)，a_i2(k)，b_ij(k)的估计值设计容错控制器：

u_i(k)＝u_gi(k)+u_Ii(k)，i＝1,2,3，(3)

(3)式中：u_i(k)为第k个采样时刻控制力矩u的第i个分量，u_gi(k)为黄金分割自适应控制律，设计为

$u_{gi}\left(k\right)=\frac{-b_{0i}\&lsqb;l_1{\hat{a}}_{i1}{e}_i\left(k\right)+l_2{\hat{a}}_{i2}\left(k\right)e_i\left(k-1\right)\&rsqb;}{{\hat{b}}_{ii}\left(k\right)+{\mathrm{\&lambda;}}_{0i}},i=1,2,3,---\left(4\right)$

(4)式中：b_0i，λ_0i均为大于0的控制增益调节参数，视控制效果可适当增减；l₁＝0.382，l₂＝0.618，e_i(k)＝y_i(k)-y_ri，y_ri为姿态角的期望值，u_Ii(k)为逻辑积分控制律；

第四步，将第三步设计的u_i(k)施加到由(1)式所描述的系统中，以控制卫星的姿态角。