CN103235597A - 一种航天器的快速姿态机动快速稳定联合控制方法 - Google Patents

Abstract

本发明涉及一种航天器的快速姿态机动快速稳定联合控制方法，属于航天器姿态控制和振动控制领域。（1）吸取航天器姿态机动轨迹规划技术和输入成形技术的优势，发明一种航天器快速姿态机动快速稳定联合控制方法，能够使得航天器在任务要求时间内完成机动，并且机动后能够保证航天器姿态快速稳定到指标要求值以内；（2）本发明还综合考虑了航天器姿态控制执行机构的力矩输出能力和航天器的最大角速度机动能力，使航天器姿态控制执行机构的输出力矩能够易于实现。

Description

一种航天器的快速姿态机动快速稳定联合控制方法

技术领域

本发明涉及一种航天器的快速姿态机动快速稳定联合控制方法，属于航天器姿态控制和振动控制领域。

背景技术

近年来，高分辨率遥感卫星、激光通讯卫星等大型复杂航天器的发展得到了广泛的关注。该类航天器除了具有姿态高精度和高稳定度的控制能力以外，还需具备快速机动快速稳定的能力，以完成敏捷成像和快速跟踪瞄准等空间任务。

以上类型航天器在姿态快速机动过程中，难免会激起所携带的大型柔性附件的结构振动。这类振动靠自身阻尼特性很难衰减，如不采取有效的振动抑制措施，会导致有效载荷难以正常工作，使空间观测任务无法完成。如美国初期的哈勃望远镜（HST），由于没有对冷热交变环境引起的柔性帆板振动进行处理，导致成像质量很差，后经两次维修，才使得哈勃望远镜成为了太空望远镜的典范。

近些年很多学者围绕大型复杂航天器的快速姿态机动快速稳定控制问题进行了深入研究。如通过智能材料采用分布式进行柔性结构的振动控制。Gopinath等研究了表面粘贴式和内部嵌入式两种压电元件联合使用的振动抑制效果（Gopinath T,Raja S,Tadashige I.Finite element formulation of laminated platewith flexible piezoelectric actuators and vibration control analysis[C]，Proceedings ofthe SPIE-The International Society for Optical Engineering,San Diego,CA,USA,2011）；Orszulik等研究了正位置反馈算法与比例微分控制器相结合的振动抑制策略（Orszulik R R,Jinjun S.Vibration control using input shaping and adaptivepositive position feedback[J].Journal of Guidance,Control,and Dynamics,2011,34(4):1031-1044）；Mahmoodi等在一种改进正位置反馈控制器的基础上又设计出了一种具有自适应能力的振动抑制控制器（Mahmoodi N S,Ahmadian M,Inman D J.Adaptive modified positive position feedback for active vibration controlof structures[J].Journal of Intelligent Material Systems and Structures,April2010.21(6):571-580）。除分布式振动控制方法以外，研究较多的还有集中式振动控制，如HST上使用6阶双凹陷滤波器与PID控制器串联，并在SA3帆板根部安装一种被动阻尼器（Anandakrishnan S M,Connor C T,Lee S,et al.Hubble spacetelescope solar damper for improving control system stability.Aerospace ConferenceProceedings,2000IEEE,4:261-276）；Wie等提出了广义结构滤波器结构，并对结构滤波器的可用性进行了COFS-I(Control of flexible structures mast flightsystem)实验，不仅验证了共位配置下的最小相位结构滤波器的有效性，还验证了非共位配置下的非最小相位结构滤波器的有效性（Wie,B.Experimentaldemonstration of a classical approach for flexible structure control[J].Journal ofGuidance,Control,and Dynamics,1992,15(6):1327-1333）。另一种集中振动控制的方式是通过使用输入成形实现。输入成形思想的雏形被认为是Smith提出的Posicast控制。Singer和Seering将这一思想得到升华（Singer N C,Seering W P.Preshaping command inputs to reduce system vibration[J].Transactions of the ASME.Journal of Dynamic Systems,Measurement and Control.1990,112(1):76-82.）。在实际应用中，输入成形作为开环控制器，或者作为反馈控制系统的前置滤波器，与控制律的联合设计可提高系统的响应特性。

对于分布式振动控制而言，由于需要在大型挠性附件上布置智能元件，造成了工程实现性较差。对于集中式振动控制而言，也因为需要完全已知航天器系统的模态信息才能进行有效振动抑制。因此综合而言，这些振动控制技术都还存在着一些弊端：如在进行振动抑制时并没有能够很好的结合航天器的运动特性，也没有将航天器姿态控制执行机构的输出能力及响应能力得到合理利用和规划。这将会给这些振动控制技术在工程实施上带来困难，造成航天器姿态控制执行机构输出不了所期望的力矩，导致无法对振动进行抑制，甚至还会出现对振动放大的反效果。

发明内容

本发明的目的是为了解决现有传统姿态机动控制中面临的快速性和稳定性很难折衷的问题，提供一种航天器的快速姿态机动快速稳定联合控制方法，通过结合轨迹规划技术和输入成形技术的优势，使得航天器既能保证航天器快速机动快速稳定的能力，又能够避免姿态控制执行机构的饱和限制问题。

本发明的目的是通过下述技术方案实现的：

本发明公开的一种航天器的快速姿态机动快速稳定联合控制方法，具体步骤如下：

步骤一：根据欧拉旋转定理，刚体的姿态从给定方位到任一其它方位的改变可通过绕欧拉特征轴的旋转而成，其间经历的角行程是最短的，欧拉特征轴在刚体旋转过程中是固联于刚体的，相对惯性空间也是不动的。因此，为了实现航天器的快速机动快速稳定控制，让航天器按照欧拉特征轴-角方式，沿最短路径进行机动。这就需要首先规划出绕特征轴的机动路径，再对其进行跟踪控制。本步骤则是对航天器从一个姿态转变到另一个姿态进行了特征轴的计算及旋转角度的确定。具体如下：

首先定义航天器相对参考坐标系的姿态四元数为

$Q=q_0+q=\cos \frac{\mathrm{\σ}}{2}+\sin \frac{\mathrm{\σ}}{2}n---\left(1\right)$

其中，q₀是四元数中的标量参数，q是四元数中的矢量参数，而四元数矢量部分的n就代表了欧拉旋转轴的方向，标量部分的σ就代表了绕欧拉轴的旋转角度。因而通过初末姿态的四元数来计算特征轴和旋转角度。

设航天器的初始姿态四元数为Q₁=q₁₀+q₁，目标姿态四元数为Q_t=q_t0+q_t，其中q₁₀和q_t0分别是初始姿态四元数和目标姿态四元数中的标量参数，q₁和q_t分别是初始姿态四元数和目标姿态四元数中的矢量参数。航天器姿态机动的特征主轴四元数Q_e可表达为

Q_e=Q₁ ^-1Q_t=(q₁₀-q₁)(q_t0+q_t)=q_e0+q_e   (2)

根据公式(2)得到沿特征轴旋转的角度表达式为

σ=2arccos(q_e0)   (3)

同样根据公式(2)得到特征轴表达式为

$n_e=\frac{1}{\sin \frac{\mathrm{\σ}}{2}}{q}_e---\left(4\right)$

将公式(4)写成标量形式，即表示在航天器本体坐标系下为

$\left\{\begin{array}{c}{n}_{\mathrm{ex}}=q_{e1}/\sin \frac{\mathrm{\σ}}{2}\\ n_{\mathrm{ey}}=q_{e2}/\sin \frac{\mathrm{\σ}}{2}\\ n_{\mathrm{ez}}=q_{e3}/\sin \frac{\mathrm{\σ}}{2}\end{array}\right.---\left(\text{5}\right)$

其中n_ex，n_ey和n_ez分别是矢量n_e三个方向上的分量；q_e1，q_e2和q_e3分别是矢量q_e的三个分量。如令规划出的角速度ω为绕特征轴n_e的旋转角速度。则航天器的期望角速度轨迹表示为

ω_r=ωn_e   (6)

在航天器本体系下描述为

$\left\{\begin{array}{c}{\mathrm{\ω}}_{\mathrm{rx}}={\mathrm{\ωn}}_{\mathrm{ex}}\\ {\mathrm{\ω}}_{\mathrm{ry}}={\mathrm{\ωn}}_{\mathrm{ey}}\\ {\mathrm{\ω}}_{\mathrm{rz}}={\mathrm{\ωn}}_{\mathrm{ez}}\end{array}\right.---\left(7\right)$

经过上式换算之后，就得到了各坐标轴的期望角速度运动规律，大角度姿态机动就变成了姿态跟踪问题。

步骤二：把航天器的角加速度曲线划分为匀加速、匀速、匀减速三段，根据航天器姿态机动角加速度a_max和最大角速度ω_max的限制，完成对航天器姿态机动的轨迹规划，即完成对步骤一中航天器绕欧拉轴的旋转角度的规划。具体如下：

根据机动角度的大小判断是否需要匀速段。当机动角度较小时，机动过程可不包括匀速段，而只由匀加速和匀减速两段构成。以角加速度a_max从0加速到最大角速度ω_max所需的时间t_a0=ω_max/a_max为判定依据。若机动角度

有匀速段，匀加速段和匀减速段时间等长t₁=t_a0，匀速段时长t_y=Δθω_max-t₁，总机动时间t_f=2t₁+t_y。

若机动角度

无匀速段，匀加速段时长和匀减速时长相等，即

匀速段时长t_y=0，总机动时间t_f=2t₁+t_y。

根据以上对机动时长的判定，可得到航天器机动过程的角加速度如下

$a=\left\{\begin{array}{cc}{a}_{\max },& 0\&le;t<t_1\\ 0,& t_1\&le;t<t_1+t_y\\ -a_{\max },& t_1+t_y\&le;t<t_f\\ 0,& t\&GreaterEqual;t_f\end{array}\right.---\left(8\right)$

对方程(8)求一次微分可得航天器机动过程的角速度，求二次微分可得航天器机动过程中的角度，即完成了轨迹规划的设计。

步骤三：建立带有挠性附件的航天器姿态动力学模型，为以后的航天器姿态控制器设计和输入成形器设计提供模型，并为数值仿真验证提供基础。该步骤具体操作如下：

认为航天器上带有N个挠性附件，并且认为航天器的中心刚体角速度、挠性附件相对中心体的角速度和挠性附件弹性振动速度很小，由此引起的高阶非线性耦合项可忽略。则可得如下航天器姿态动力学方程：

$I_s{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{bak}}{I}_{\mathrm{ak}}^b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{bak}}{H}_{\mathrm{bak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}+{\mathrm{\&omega;}}_b^{\&times;}\left(I_s{\mathrm{\&omega;}}_b\right)=T_c{+T}_d---\left(\text{9}\right)$

${I_{\mathrm{ak}}^b}^T{A}_{\mathrm{akb}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+I_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+H_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}=T_{\mathrm{ak}}---\left(10\right)$

$H_{\mathrm{bak}}^T{A}_{\mathrm{akb}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+M_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}+K_{\mathrm{ak}}{q}_{\mathrm{ak}}+C_{\mathrm{ak}}{\stackrel{\&CenterDot;}{q}}_{\mathrm{ak}}=0,(k=\mathrm{1,2},...,N)---\left(11\right)$

其中方程(9)为航天器中心刚体的转动方程，方程(10)和方程(11)分别为第k个挠性附件的转动方程和振动方程。以上方程中，ω_b为航天器中心刚体（星体）的角速度；ω_ak为挠性附件k相对于中心刚体的角速度；q_ak为第k个附件的前l阶模态坐标组成的l×1列阵；I_s为整个航天器的转动惯量矩阵；M_ak为附件k的模态质量阵；I_ak为附件k对其体坐标系的转动惯量矩阵；H_ak为附件k对其体坐标系的模态角动量系数矩阵；为附件k对星本体坐标系的转动惯量矩阵；H_bak为附件k对星本体坐标系的耦合转动惯量矩阵；K_ak和C_ak分别为第k个挠性附件的模态刚度阵和模态阻尼阵；A_bak为第k个挠性附件的体坐标系到星本体坐标系的坐标转换矩阵；T_c为控制力矩；T_d为环境干扰力矩；T_ak为附件在铰接处受到的驱动力矩。

进一步对模态坐标进行归一化处理，令

q_ak=Φ_akη_ak   (12)

其中，Φ_ak为广义模态矩阵，则以上三式可重写为

$I_s{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{R}_{\mathrm{bak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_{\mathrm{bak}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+{\mathrm{\&omega;}}_b^{\&times;}\left(I_s{\mathrm{\&omega;}}_b\right)=T_c{+T}_d---\left(\text{13}\right)$

$I_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+F_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+R_{\mathrm{bak}}^T{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b=T_{\mathrm{ak}}---\left(14\right)$

${\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+{\mathrm{\&Lambda;}}_{\mathrm{ak}}^2{\mathrm{\&eta;}}_{\mathrm{ak}}+F_{\mathrm{bak}}^T{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b=0,(k=\mathrm{1,2},...,N)---\left(15\right)$

其中，Λ_ak为附件的模态频率对角阵，ξ_ak为附件的模态阻尼矩阵，R_bak为挠性附件转动对中心刚体转动的刚性耦合系数矩阵，表示为 $R_{\mathrm{bak}}=A_{\mathrm{bak}}{I}_{\mathrm{ak}}^b---\left(16\right)$

F_bak为挠性附件振动对中心刚体转动的柔性耦合系数矩阵，表示为

F_bak=A_bakH_bakΦ_ak   (17)

F_ak为挠性附件振动对自身转动的柔性耦合系数矩阵，表示为

F_ak=H_akΦ_ak   (18)

方程(13)、(14)和(15)组成了带有多个挠性附件的航天器姿态动力学模型，模型中保留了主要的非线性项、干扰力矩项。

步骤四：根据步骤二中得到的航天器姿态动力学模型，认为航天器在执行姿态机动任务时，大型挠性附件处于锁死状态，即在航天器姿态动力学方程中暂且忽略挠性附件的转动方程及其转动角速度。并且忽略环境干扰力矩。可得到简化的动力学模型，并将其写成状态方程的形式，用于计算航天器系统的模态振动参数，模态振动参数包含模态频率和阻尼比，并完成输入成形器的设计。具体操作如下：

将方程(13)和方程(15)写成如下形式

$\left[\begin{array}{cc}{I}_s{J}_{\mathrm{vb}}& F_{\mathrm{bak}}\\ F_{\mathrm{bak}}^T{J}_{\mathrm{vb}}& E\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{cc}{\mathrm{\&omega;}}_b^{\&times;}{I}_s{J}_{\mathrm{vb}}& 0\\ 0& 2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& {\mathrm{\&Lambda;}}_{\mathrm{ak}}^2\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_b\\ {\mathrm{\&eta;}}_{\mathrm{ak}}\end{array}\right]=\left[\begin{array}{c}{T}_c\\ 0\end{array}\right]---\left(19\right)$

其中，

T_c为控制力矩的表达式，本文选用工程易于实现的反馈控制器作为姿态控制器。如下所示：

$\left[\begin{array}{c}{T}_c\\ 0\end{array}\right]=\left[\begin{array}{cc}-K_{\mathrm{db}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{c}{K}_{\mathrm{db}}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_r\\ 0\end{array}\right]+\left[\begin{array}{cc}-K_{\mathrm{pb}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_b\\ {\mathrm{\&eta;}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{c}{K}_{\mathrm{pb}}{\mathrm{\&theta;}}_r\\ 0\end{array}\right]---\left(20\right)$

其中，θ_r为期望姿态角度，K_db和K_pb为控制参数。

令 $x={\left[\begin{array}{cc}{\mathrm{\&theta;}}_b^T& {\mathrm{\&eta;}}_{\mathrm{ak}}^T\end{array}\right]}^T,$ $X={\left[\begin{array}{cc}{x}^T& {\stackrel{\&CenterDot;}{x}}^T\end{array}\right]}^T,$ 则方程(19)和(20)可写成以下状态方程的形式

$\stackrel{\&CenterDot;}{X}=\mathrm{AX}+\mathrm{BU}---\left(21\right)$

其中，

$A=\left[\begin{array}{cc}0& E\\ -M^{-1}K& -M^{-1}C\end{array}\right];$ $B=\left[\begin{array}{c}0\\ M^{-1}T\end{array}\right];$ $M=\left[\begin{array}{cc}{I}_s{J}_{\mathrm{vb}}& F_{\mathrm{bak}}\\ F_{\mathrm{bak}}^T{J}_{\mathrm{vb}}& E\end{array}\right];$ $C=\left[\begin{array}{cc}{K}_{\mathrm{db}}& 0\\ 0& 2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}\end{array}\right]$

$K=\left[\begin{array}{cc}{K}_{\mathrm{pb}}& 0\\ 0& {\mathrm{\&Lambda;}}_{\mathrm{ak}}^2\end{array}\right];$ $T=\left[\begin{array}{c}{K}_{\mathrm{db}}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_r+K_{\mathrm{pb}}{\mathrm{\&theta;}}_r\\ 0\end{array}\right]$

根据闭环系统方程(21)，得到系统矩阵A的特征值。那么系统模态的振动频率和阻尼比的数值解由下式给出：

${\mathrm{\&lambda;}}_{\mathrm{sys}}=-{\mathrm{\&xi;}}_{\mathrm{sys}}{\mathrm{\&omega;}}_{\mathrm{sys}}\&PlusMinus;{\mathrm{j\&omega;}}_{\mathrm{sys}}\sqrt{1-{\mathrm{\&xi;}}_{\mathrm{sys}}^2}---\left(22\right)$

其中λ_sys为矩阵A的特征值，ω_sys和ξ_sys为系统模态振动频率和阻尼比。

为提高系统对参数变化的鲁棒性，并且从机动时间上折衷考虑，所述的输入成形器为ZVD(Zero Vibration and Derivative)输入成形器。针对任意一组模态振动参数的ZVD输入成形器表达式如下所示：

$\left[\begin{array}{ccc}0& t_1& t_2\\ A_1& A_2& A_3\end{array}\right]=\left[\begin{array}{ccc}0& \frac{\mathrm{\&pi;}}{{\mathrm{\&omega;}}_d}& \frac{2\mathrm{\&pi;}}{{\mathrm{\&omega;}}_d}\\ \frac{1}{1+2K+K^2}& \frac{2K}{1+2K+K^2}& \frac{K^2}{1+2K+K^2}\end{array}\right]---\left(23\right)$

其中，第一行表示各脉冲的作用时刻（其中定义第一个脉冲发生在0时刻），第二行表示各脉冲的幅值； $K=e^{-{\mathrm{\&xi;}}_{\mathrm{eys}}\mathrm{\&pi;}/\sqrt{1-{{\mathrm{\&xi;}}_{\mathrm{sys}}}^2}};$ ${\mathrm{\&omega;}}_d={\mathrm{\&omega;}}_{\mathrm{sys}}\sqrt{1-{{\mathrm{\&xi;}}_{\mathrm{sys}}}^2}.$

在使用过程中，可根据需要将多个输入成形器做卷积，合成所需要的输入成形器，即ZVD=ZVD1*ZVD2*...*ZVDn。

步骤五：将步骤一中计算得到的航天器姿态机动所绕的特征轴和特征转角、步骤二中设计得到的姿态角度轨迹规划的曲线、步骤四中设计得到的输入成形器以及反馈控制器和步骤三中推导得出的航天器姿态动力学共同构成航天器姿态控制系统完整回路，并且需要注意的是所述快速姿态机动快速稳定联合控制方法要将轨迹规划模块放置在输入成形器模块之前，并且属于姿态控制系统回路中的前馈环节。根据此完整回路，通过使用数值仿真软件搭建航天器姿态控制系统，以完成对本专利所述的航天器快速姿态机动快速稳定联合控制方法的验证。

有益效果：

1、本发明的一种航天器的快速姿态机动快速稳定联合控制方法，克服现有传统姿态机动控制中面临的快速性和稳定性很难折衷的问题，通过联合使用轨迹规划和输入成形器以使得航天器在快速机动中尽量不激起挠性附件的振动，并对已经激起的振动进行一定的振动控制。通过使用发明的航天器快速姿态机动快速稳定联合控制方法，能够实现航天器在30s内完成30°的机动，并且机动后的姿态精度能够达到0.002°，姿态稳定度能够达到1.3×10^-3°/s。

2、本发明的一种航天器的快速姿态机动快速稳定联合控制方法，为避免对姿态控制系统的稳定性造成影响，所设计的轨迹规划和输入成形器均是以前馈形式加在航天器姿态控制系统中。其次，本发明综合考虑了航天器姿态控制执行机构的力矩输出能力和航天器的最大角速度机动能力，使航天器30s内完成30°机动所使用的输出力矩不超过8Nm，该输出力矩对于目前航天器而言，易于实现。

3、本发明的一种航天器的快速姿态机动快速稳定联合控制方法，本专利对集中式振动控制方法进行改进，在综合考虑航天器执行机构输出能力和航天器最大角速度约束，合理规划航天器的运动轨迹，将运动轨迹规划和集中式振动控制方法中的输入成形技术联合起来，来完成带有大型挠性附件的航天器的快速姿态机动和快速稳定控制任务。

附图说明

图1为轨迹规划的航天器三轴角速度曲线图；

图2为轨迹规划的航天器三轴角度曲线图；

图3航天器结构简图；

图4为航天器快速姿态机动快速稳定联合控制系统回路；

图5为使用了本发明的联合控制方法后的航天器姿态角度随时间变化图；

图6为使用了本发明的联合控制方法后的航天器姿态角速度随时间变化图；

图7为使用了本发明的联合控制方法后的航天器期望力矩随时间变化图；

图8为航天器的传统姿态控制系统回路；

图9为传统姿态控制方法下的航天器姿态角度随时间变化图；

图10为传统姿态控制方法下的航天器姿态角速度随时间变化图；

图11为传统姿态控制方法下的航天器期望力矩随时间变化图。

具体实施方式

下面将结合附图对本发明作进一步的详细说明。

实施例1

一种航天器的快速姿态机动快速稳定联合控制方法，可通过下述步骤来完成：

步骤一：根据欧拉旋转定理，刚体的姿态从给定方位到任一其它方位的改变可通过绕欧拉特征轴的旋转而成，其间经历的角行程是最短的，欧拉特征轴在刚体旋转过程中是固联于刚体的，相对惯性空间也是不动的。因此，为了实现航天器的快速机动快速稳定控制，让航天器按照欧拉特征轴-角方式，沿最短路径进行机动。这就需要首先规划出绕特征轴的机动路径，再对其进行跟踪控制。本步骤则是对航天器从一个姿态转变到另一个姿态进行了特征轴的计算及旋转角度的确定。具体如下：

首先定义航天器相对参考坐标系的姿态四元数为

$Q=q_0+q=\cos \frac{\mathrm{\&sigma;}}{2}+\sin \frac{\mathrm{\&sigma;}}{2}n---\left(1\right)$

其中，q₀是四元数中的标量参数，q是四元数中的矢量参数，而四元数矢量部分的n就代表了欧拉旋转轴的方向，标量部分的σ就代表了绕欧拉轴的旋转角度。因而通过初末姿态的四元数来计算特征轴和旋转角度。

设航天器的初始姿态四元数为Q₁=q₁₀+q₁，目标姿态四元数为Q_t=q_t0+q_t，其中q₁₀和q_t0分别是初始姿态四元数和目标姿态四元数中的标量参数，q₁和q_t分别是初始姿态四元数和目标姿态四元数中的矢量参数。航天器姿态机动的特征主轴四元数Q_e可表达为

Q_e=Q₁ ^-1Q_t=(q₁₀-q₁)(q_t0+q_t)=q_e0+q_e   (2)

根据公式(2)得到沿特征轴旋转的角度表达式为

σ=2arccos(q_e0)   (3)

同样根据公式(2)得到特征轴表达式为

$n_e=\frac{1}{\sin \frac{\mathrm{\&sigma;}}{2}}{q}_e---\left(4\right)$

将公式(4)写成标量形式，即表示在航天器本体坐标系下为

$\left\{\begin{array}{c}{n}_{\mathrm{ex}}=q_{e1}/\sin \frac{\mathrm{\&sigma;}}{2}\\ n_{\mathrm{ey}}=q_{e2}/\sin \frac{\mathrm{\&sigma;}}{2}\\ n_{\mathrm{ez}}=q_{e3}/\sin \frac{\mathrm{\&sigma;}}{2}\end{array}\right.---\left(\text{5}\right)$

其中n_ex，n_ey和n_ez分别是矢量n_e三个方向上的分量；q_e1，q_e2和q_e3分别是矢量q_e的三个分量。如令规划出的角速度ω为绕特征轴n_e的旋转角速度。则航天器的期望角速度轨迹表示为

ω_r=ωn_e   (6)

在航天器本体系下描述为

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_{\mathrm{rx}}={\mathrm{\&omega;n}}_{\mathrm{ex}}\\ {\mathrm{\&omega;}}_{\mathrm{ry}}={\mathrm{\&omega;n}}_{\mathrm{ey}}\\ {\mathrm{\&omega;}}_{\mathrm{rz}}={\mathrm{\&omega;n}}_{\mathrm{ez}}\end{array}\right.---\left(7\right)$

经过上式换算之后，就得到了各坐标轴的期望角速度运动规律，大角度姿态机动就变成了姿态跟踪问题。

认为航天器姿态初始角度为0°，初始角速度为0°/s。三轴上的期望姿态角度分别为30°、-30°和0°。通过以上步骤计算得到特征轴和特征转角，如下所示：

σ=42.2°，n_e=[0.6974 -0.6974 -0.1862]^T

步骤二：把航天器的角加速度曲线划分为匀加速、匀速、匀减速三段，根据航天器姿态机动角加速度a_max和最大角速度ω_max的限制，完成对航天器姿态机动的轨迹规划，即完成对步骤一中航天器绕欧拉轴的旋转角度的规划。具体如下：

根据机动角度的大小判断是否需要匀速段。当机动角度较小时，机动过程可不包括匀速段，而只由匀加速和匀减速两段构成。以角加速度a_max从0加速到最大角速度ω_max所需的时间t_a0=ω_max/a_max为判定依据。若机动角度

有匀速段，匀加速段和匀减速段时间等长t₁=t_a0，匀速段时长t_y=Δθω_max-t₁，总机动时间t_f=2t₁+t_y。

若机动角度

无匀速段，匀加速段时长和匀减速时长相等，即

匀速段时长t_y=0，总机动时间t_f=2t₁+t_y。

根据以上对机动时长的判定，可得到航天器机动过程的角加速度如下

$a=\left\{\begin{array}{cc}{a}_{\max },& 0\&le;t<t_1\\ 0,& t_1\&le;t<t_1+t_y\\ -a_{\max },& t_1+t_y\&le;t<t_f\\ 0,& t\&GreaterEqual;t_f\end{array}\right.---\left(8\right)$

对方程(8)求一次微分可得航天器机动过程的角速度，求二次微分可得航天器机动过程中的角度，即完成了轨迹规划的设计。

认为卫星最大控制力矩T_max=25Nm，卫星绕旋转轴的转动惯量I_max=1500kg·m²，允许的最大机动角加速度a_max=T_maxI_max=0.95°/s²和允许的最大角速度V_max=2.5°/s。则通过以上步骤完成对步骤一中得到的特征转角42.2°的轨迹规划，将轨迹规划的一条曲线按照步骤一中得到的特征轴进行投影，即可得到航天器三个姿态方向上的轨迹规划角速度和角度的曲线，如图1和图2所示。

步骤三：建立带有挠性附件的航天器姿态动力学模型，为以后的航天器姿态控制器设计和输入成形器设计提供模型，并为数值仿真验证提供基础。该步骤具体操作如下：

认为航天器上带有N个挠性附件，并且认为航天器的中心刚体角速度、挠性附件相对中心体的角速度和挠性附件弹性振动速度很小，由此引起的高阶非线性耦合项可忽略。则可得如下航天器姿态动力学方程：

$I_s{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{bak}}{I}_{\mathrm{ak}}^b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{bak}}{H}_{\mathrm{bak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}+{\mathrm{\&omega;}}_b^{\&times;}\left(I_s{\mathrm{\&omega;}}_b\right)=T_c{+T}_d---\left(\text{9}\right)$

${I_{\mathrm{ak}}^b}^T{A}_{\mathrm{akb}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+I_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+H_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}=T_{\mathrm{ak}}---\left(10\right)$

$H_{\mathrm{bak}}^T{A}_{\mathrm{akb}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+M_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}+K_{\mathrm{ak}}{q}_{\mathrm{ak}}+C_{\mathrm{ak}}{\stackrel{\&CenterDot;}{q}}_{\mathrm{ak}}=0,(k=\mathrm{1,2},...,N)---\left(11\right)$

其中方程(9)为航天器中心刚体的转动方程，方程(10)和方程(11)分别为第k个挠性附件的转动方程和振动方程。以上方程中，ω_b为航天器中心刚体（星体）的角速度；ω_ak为挠性附件k相对于中心刚体的角速度；q_ak为第k个附件的前l阶模态坐标组成的l×1列阵；I_s为整个航天器的转动惯量矩阵；M_ak为附件k的模态质量阵；I_ak为附件k对其体坐标系的转动惯量矩阵；H_ak为附件k对其体坐标系的模态角动量系数矩阵；

为附件k对星本体坐标系的转动惯量矩阵；H_bak为附件k对星本体坐标系的耦合转动惯量矩阵；K_ak和C_ak分别为第k个挠性附件的模态刚度阵和模态阻尼阵；A_bak为第k个挠性附件的体坐标系到星本体坐标系的坐标转换矩阵；T_c为控制力矩；T_d为环境干扰力矩；T_ak为附件在铰接处受到的驱动力矩。

进一步对模态坐标进行归一化处理，令

q_ak=Φ_akη_ak   (12)

其中，Φ_ak为广义模态矩阵，则以上三式可重写为

$I_s{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{R}_{\mathrm{bak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_{\mathrm{bak}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+{\mathrm{\&omega;}}_b^{\&times;}\left(I_s{\mathrm{\&omega;}}_b\right)=T_c{+T}_d---\left(\text{13}\right)$

$I_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+F_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+R_{\mathrm{bak}}^T{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b=T_{\mathrm{ak}}---\left(14\right)$

${\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+{\mathrm{\&Lambda;}}_{\mathrm{ak}}^2{\mathrm{\&eta;}}_{\mathrm{ak}}+F_{\mathrm{bak}}^T{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b=0,(k=\mathrm{1,2},...,N)---\left(15\right)$

其中，Λ_ak为附件的模态频率对角阵，ξ_ak为附件的模态阻尼矩阵，R_bak为挠性附件转动对中心刚体转动的刚性耦合系数矩阵，表示为

$R_{\mathrm{bak}}=A_{\mathrm{bak}}{I}_{\mathrm{ak}}^b---\left(16\right)$

F_bak为挠性附件振动对中心刚体转动的柔性耦合系数矩阵，表示为

F_bak=A_bakH_bakΦ_ak   (17)

F_ak为挠性附件振动对自身转动的柔性耦合系数矩阵，表示为

F_ak=H_akΦ_ak   (18)

方程(13)、(14)和(15)组成了带有多个挠性附件的航天器姿态动力学模型，模型中保留了主要的非线性项、干扰力矩项。

在该实施例中认为航天器带有3块柔性太阳帆板，如图3所示。

航天器本体质量和惯量参数如下所示：

$m_b=1100.8\mathrm{kg},I_b=\left[\begin{array}{ccc}1099.7& -20& -10\\ -20& 899.7& -9.4\\ -10& -9.4& 822.8\end{array}\right]\mathrm{kg}\&CenterDot;m^2$

每块帆板的参数均相同，帆板的质量和惯量参数如下所示：

$m_{\mathrm{ak}}=40.25\mathrm{kg},I_{\mathrm{ak}}=\left[\begin{array}{ccc}71.0412& & \\ & 3.4213& \\ & & 74.4625\end{array}\right]\mathrm{kg}\&CenterDot;m^2$

帆板尺寸是2.3×1.0m。星本体固连坐标系原点到帆板安装点的矢量如下所示：

其中，r为星体柱体半径，r=0.75m。

星本体坐标系到各帆板固连坐标系的坐标转换矩阵：

帆板的阻尼比ξ=0.2%，前三阶模态频率分别如下所示：

${\mathrm{\&Lambda;}}_{a1}={\mathrm{\&Lambda;}}_{a2}={\mathrm{\&Lambda;}}_{a3}=\left[\begin{array}{ccc}0.102& & \\ & 0.486& \\ & & 0.653\end{array}\right]\mathrm{Hz}=\left[\begin{array}{ccc}0.6413& & \\ & 3.0522& \\ & & 3.9869\end{array}\right]\mathrm{rad}/s$

步骤四：根据步骤二中得到的航天器姿态动力学模型，认为航天器在执行姿态机动任务时，大型挠性附件处于锁死状态，即在航天器姿态动力学方程中暂且忽略挠性附件的转动方程及其转动角速度。并且忽略环境干扰力矩。可得到简化的动力学模型，并将其写成状态方程的形式，用于计算航天器系统的模态振动参数，模态振动参数包含模态频率和阻尼比，并完成输入成形器的设计。具体操作如下：

将方程(13)和方程(15)写成如下形式

$\left[\begin{array}{cc}{I}_s{J}_{\mathrm{vb}}& F_{\mathrm{bak}}\\ F_{\mathrm{bak}}^T{J}_{\mathrm{vb}}& E\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{cc}{\mathrm{\&omega;}}_b^{\&times;}{I}_s{J}_{\mathrm{vb}}& 0\\ 0& 2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& {\mathrm{\&Lambda;}}_{\mathrm{ak}}^2\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_b\\ {\mathrm{\&eta;}}_{\mathrm{ak}}\end{array}\right]=\left[\begin{array}{c}{T}_c\\ 0\end{array}\right]---\left(19\right)$

其中，

T_c为控制力矩的表达式，本文选用工程易于实现的反馈控制器作为姿态控制器。如下所示：

$\left[\begin{array}{c}{T}_c\\ 0\end{array}\right]=\left[\begin{array}{cc}-K_{\mathrm{db}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{c}{K}_{\mathrm{db}}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_r\\ 0\end{array}\right]+\left[\begin{array}{cc}-K_{\mathrm{pb}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_b\\ {\mathrm{\&eta;}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{c}{K}_{\mathrm{pb}}{\mathrm{\&theta;}}_r\\ 0\end{array}\right]---\left(20\right)$

其中，θ_r为期望姿态角度，K_db和K_pb为控制参数。本实施例中控制参数的选取如下所示：

K_db=diag(3000 3000 3000)，K_pb=diag(2500 2500 2500)

令 $x={\left[\begin{array}{cc}{\mathrm{\&theta;}}_b^T& {\mathrm{\&eta;}}_{\mathrm{ak}}^T\end{array}\right]}^T,$ $X={\left[\begin{array}{cc}{x}^T& {\stackrel{\&CenterDot;}{x}}^T\end{array}\right]}^T,$ 则方程(19)和(20)可写成以下状态方程的形式

$\stackrel{\&CenterDot;}{X}=\mathrm{AX}+\mathrm{BU}---\left(21\right)$

其中，

$A=\left[\begin{array}{cc}0& E\\ -M^{-1}K& -M^{-1}C\end{array}\right];$ $B=\left[\begin{array}{c}0\\ M^{-1}T\end{array}\right];$ $M=\left[\begin{array}{cc}{I}_s{J}_{\mathrm{vb}}& F_{\mathrm{bak}}\\ F_{\mathrm{bak}}^T{J}_{\mathrm{vb}}& E\end{array}\right];$ $C=\left[\begin{array}{cc}{K}_{\mathrm{db}}& 0\\ 0& 2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}\end{array}\right]$

$K=\left[\begin{array}{cc}{K}_{\mathrm{pb}}& 0\\ 0& {\mathrm{\&Lambda;}}_{\mathrm{ak}}^2\end{array}\right];$ $T=\left[\begin{array}{c}{K}_{\mathrm{db}}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_r+K_{\mathrm{pb}}{\mathrm{\&theta;}}_r\\ 0\end{array}\right]$

根据闭环系统方程(21)，得到系统矩阵A的特征值。那么系统模态的振动频率和阻尼比的数值解由下式给出：

${\mathrm{\&lambda;}}_{\mathrm{sys}}=-{\mathrm{\&xi;}}_{\mathrm{sys}}{\mathrm{\&omega;}}_{\mathrm{sys}}\&PlusMinus;{\mathrm{j\&omega;}}_{\mathrm{sys}}\sqrt{1-{\mathrm{\&xi;}}_{\mathrm{sys}}^2}---\left(22\right)$

其中λ_sys为矩阵A的特征值，ω_sys和ξ_sys为系统模态振动频率和阻尼比。

依据步骤三中航天器的参数和公式(21)和公式(22)，计算得到所关注的前三阶振动模态参数（包括0阶，1阶，2阶和3阶），如下所示：

根据以上对应的阻尼比和频率设计抑制挠性附件振动的输入成形器，再将多个输入成型器进行卷积计算，得到抑制柔性航天器多个系统模态振动的鲁棒输入成形器。

为提高系统对参数变化的鲁棒性，并且从机动时间上折衷考虑，所述的输入成形器为ZVD(Zero Vibration and Derivative)输入成形器。针对任意一组模态振动参数的ZVD输入成形器表达式如下所示：

$\left[\begin{array}{ccc}0& t_1& t_2\\ A_1& A_2& A_3\end{array}\right]=\left[\begin{array}{ccc}0& \frac{\mathrm{\&pi;}}{{\mathrm{\&omega;}}_d}& \frac{2\mathrm{\&pi;}}{{\mathrm{\&omega;}}_d}\\ \frac{1}{1+2K+K^2}& \frac{2K}{1+2K+K^2}& \frac{K^2}{1+2K+K^2}\end{array}\right]---\left(23\right)$

其中，第一行表示各脉冲的作用时刻（其中定义第一个脉冲发生在0时刻），第二行表示各脉冲的幅值； $K=e^{-{\mathrm{\&xi;}}_{\mathrm{eys}}\mathrm{\&pi;}/\sqrt{1-{{\mathrm{\&xi;}}_{\mathrm{sys}}}^2}};$ ${\mathrm{\&omega;}}_d={\mathrm{\&omega;}}_{\mathrm{sys}}\sqrt{1-{{\mathrm{\&xi;}}_{\mathrm{sys}}}^2}.$

将得到的阻尼比和频率参数分别代入到公式(13)中，可得到4个输入成形器，如下所示：

$\mathrm{ZVD}1=\left[\begin{array}{ccc}0& 4.9524& 9.9048\\ 0.2594& 0.4998& 0.2407\end{array}\right];$ $\mathrm{ZVD}2=\left[\begin{array}{ccc}0& 4.2248& 8.4496\\ 0.9912& 0.0087& 0.0001\end{array}\right]$

$\mathrm{ZVD}3=\left[\begin{array}{ccc}0& 1.0284& 2.0567\\ 0.2526& 0.5& 0.2474\end{array}\right];$ $\mathrm{ZVD}4=\left[\begin{array}{ccc}0& 0.7842& 1.5685\\ 0.2544& 0.5& 0.2456\end{array}\right]$

在使用过程中，需要将四个输入成形器做卷积，合成所需要的输入成形器，即ZVD=ZVD1*ZVD2*ZVD3*ZVD4。

步骤五：将步骤一中计算得到的航天器姿态机动所绕的特征轴和特征转角、步骤二中设计得到的姿态角度轨迹规划的曲线、步骤四中设计得到的输入成形器以及反馈控制器和步骤三中推导得出的航天器姿态动力学共同构成航天器姿态控制系统完整回路，并且需要注意的是所述快速姿态机动快速稳定联合控制方法要将轨迹规划模块放置在输入成形器模块之前，并且属于姿态控制系统回路中的前馈环节，如图4所示，其中干扰输入可根据需要进行施加外部扰动力矩。根据此完整回路，通过使用数值仿真软件搭建航天器姿态控制系统，可分别仿真得到航天器姿态角度曲线、航天器姿态角速度曲线和航天器姿态机动的期望力矩曲线，如图5、图6和图7所示。

由航天器姿态角度曲线和航天器姿态角速度曲线得知，航天器姿态角度由初始零值到期望值所用的时间为30秒，达到稳定后的最大姿态角度偏差为0.002°，姿态稳定度为1.3×10^-3°/s。由航天器姿态机动的期望力矩曲线可知，执行机构的力矩输出在8Nm，在工程中也是很容易实现的。

为能够体现本专利所述的航天器快速姿态机动快速稳定联合控制方法的优势，与不使用本专利提出的联合控制方法的航天器姿态机动控制效果进行比较。通过使用相同的仿真环境，搭建航天器姿态控制系统，如图8所示。可得航天器的姿态角度曲线、航天器姿态角速度曲线和航天器姿态机动的期望力矩曲线，如图9、图10和图11所示。

通过和步骤五中所得结果进行比较。可得：航天器姿态角度由初始零值到期望值所用的时间为10秒左右，但是由于存在挠性附件的振动，使得最后的姿态角度和角速度在稳定值附近产生振荡。最大角度偏差达到了0.6°，姿态稳定度也仅有0.4°/s。由于没有对航天器姿态控制执行机构输出力矩能力以及最大姿态机动角速度进行限制，导致力矩输出最大达到了107Nm，这样大的力矩输出在实际中很难实现。这也进一步说明了本发明方法的有效性和工程可行性。

Claims (1)

1.一种航天器的快速姿态机动快速稳定联合控制方法，其特征在于：具体步骤如下：

步骤一：根据欧拉旋转定理，刚体的姿态从给定方位到任一其它方位的改变可通过绕欧拉特征轴的旋转而成，其间经历的角行程是最短的，欧拉特征轴在刚体旋转过程中是固联于刚体的，相对惯性空间也是不动的。因此，为了实现航天器的快速机动快速稳定控制，让航天器按照欧拉特征轴-角方式，沿最短路径进行机动。这就需要首先规划出绕特征轴的机动路径，再对其进行跟踪控制。本步骤则是对航天器从一个姿态转变到另一个姿态进行了特征轴的计算及旋转角度的确定。具体如下：

首先定义航天器相对参考坐标系的姿态四元数为

$Q=q_0+q=\cos \frac{\mathrm{\&sigma;}}{2}+\sin \frac{\mathrm{\&sigma;}}{2}n---\left(1\right)$

其中，q₀是四元数中的标量参数，q是四元数中的矢量参数，而四元数矢量部分的n就代表了欧拉旋转轴的方向，标量部分的σ就代表了绕欧拉轴的旋转角度。因而通过初末姿态的四元数来计算特征轴和旋转角度。

设航天器的初始姿态四元数为Q₁=q₁₀+q₁，目标姿态四元数为Q_t=q_t0+q_t，其中q₁₀和q_t0分别是初始姿态四元数和目标姿态四元数中的标量参数，q₁和q_t分别是初始姿态四元数和目标姿态四元数中的矢量参数。航天器姿态机动的特征主轴四元数Q_e可表达为

Q_e=Q₁ ^-1Q_t=(q₁₀-q₁)(q_t0+q_t)=q_e0+q_e    (2)

根据公式(2)得到沿特征轴旋转的角度表达式为

σ=2arccos(q_e0)    (3)

同样根据公式(2)得到特征轴表达式为

$n_e=\frac{1}{\sin \frac{\mathrm{\&sigma;}}{2}}{q}_e---\left(4\right)$

将公式(4)写成标量形式，即表示在航天器本体坐标系下为

$\left\{\begin{array}{c}{n}_{\mathrm{ex}}=q_{e1}/\sin \frac{\mathrm{\&sigma;}}{2}\\ n_{\mathrm{ey}}=q_{e2}/\sin \frac{\mathrm{\&sigma;}}{2}\\ n_{\mathrm{ez}}=q_{e3}/\sin \frac{\mathrm{\&sigma;}}{2}\end{array}\right.---\left(\text{5}\right)$

其中n_ex，n_ey和n_ez分别是矢量n_e三个方向上的分量；q_e1，q_e2和q_e3分别是矢量q_e的三个分量。如令规划出的角速度ω为绕特征轴n_e的旋转角速度。则航天器的期望角速度轨迹表示为

ω_r=ωn_e    (6)

在航天器本体系下描述为

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_{\mathrm{rx}}={\mathrm{\&omega;n}}_{\mathrm{ex}}\\ {\mathrm{\&omega;}}_{\mathrm{ry}}={\mathrm{\&omega;n}}_{\mathrm{ey}}\\ {\mathrm{\&omega;}}_{\mathrm{rz}}={\mathrm{\&omega;n}}_{\mathrm{ez}}\end{array}\right.---\left(7\right)$

经过上式换算之后，就得到了各坐标轴的期望角速度运动规律，大角度姿态机动就变成了姿态跟踪问题。

步骤二：把航天器的角加速度曲线划分为匀加速、匀速、匀减速三段，根据航天器姿态机动角加速度a_max和最大角速度ω_max的限制，完成对航天器姿态机动的轨迹规划，即完成对步骤一中航天器绕欧拉轴的旋转角度的规划。具体如下：

根据机动角度的大小判断是否需要匀速段。当机动角度较小时，机动过程可不包括匀速段，而只由匀加速和匀减速两段构成。以角加速度a_max从0加速到最大角速度ω_max所需的时间t_a0=ω_max/a_max为判定依据。若机动角度

有匀速段，匀加速段和匀减速段时间等长t₁=t_a0，匀速段时长t_y=Δθω_max-t₁，总机动时间t_f=2t₁+t_y。

若机动角度

无匀速段，匀加速段时长和匀减速时长相等，即

匀速段时长t_y=0，总机动时间t_f=2t₁+t_y。

根据以上对机动时长的判定，可得到航天器机动过程的角加速度如下

$a=\left\{\begin{array}{cc}{a}_{\max },& 0\&le;t<t_1\\ 0,& t_1\&le;t<t_1+t_y\\ -a_{\max },& t_1+t_y\&le;t<t_f\\ 0,& t\&GreaterEqual;t_f\end{array}\right.---\left(8\right)$

对方程(8)求一次微分可得航天器机动过程的角速度，求二次微分可得航天器机动过程中的角度，即完成了轨迹规划的设计。

步骤三：建立带有挠性附件的航天器姿态动力学模型，为以后的航天器姿态控制器设计和输入成形器设计提供模型，并为数值仿真验证提供基础。该步骤具体操作如下：

认为航天器上带有N个挠性附件，并且认为航天器的中心刚体角速度、挠性附件相对中心体的角速度和挠性附件弹性振动速度很小，由此引起的高阶非线性耦合项可忽略。则可得如下航天器姿态动力学方程：

$I_s{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{bak}}{I}_{\mathrm{ak}}^b{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{\mathrm{bak}}{H}_{\mathrm{bak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}+{\mathrm{\&omega;}}_b^{\&times;}\left(I_s{\mathrm{\&omega;}}_b\right)=T_c{+T}_d---\left(\text{9}\right)$

${I_{\mathrm{ak}}^b}^T{A}_{\mathrm{akb}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+I_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+H_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}=T_{\mathrm{ak}}---\left(10\right)$

$H_{\mathrm{bak}}^T{A}_{\mathrm{akb}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+M_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{\mathrm{ak}}+K_{\mathrm{ak}}{q}_{\mathrm{ak}}+C_{\mathrm{ak}}{\stackrel{\&CenterDot;}{q}}_{\mathrm{ak}}=0,(k=\mathrm{1,2},...,N)---\left(11\right)$

其中方程(9)为航天器中心刚体的转动方程，方程(10)和方程(11)分别为第k个挠性附件的转动方程和振动方程。以上方程中，ω_b为航天器中心刚体（星体）的角速度；ω_ak为挠性附件k相对于中心刚体的角速度；q_ak为第k个附件的前l阶模态坐标组成的l×1列阵；I_s为整个航天器的转动惯量矩阵；M_ak为附件k的模态质量阵；I_ak为附件k对其体坐标系的转动惯量矩阵；H_ak为附件k对其体坐标系的模态角动量系数矩阵；

为附件k对星本体坐标系的转动惯量矩阵；H_bak为附件k对星本体坐标系的耦合转动惯量矩阵；K_ak和C_ak分别为第k个挠性附件的模态刚度阵和模态阻尼阵；A_bak为第k个挠性附件的体坐标系到星本体坐标系的坐标转换矩阵；T_c为控制力矩；T_d为环境干扰力矩；T_ak为附件在铰接处受到的驱动力矩。

进一步对模态坐标进行归一化处理，令

q_ak=Φ_akη_ak    (12)

其中，Φ_ak为广义模态矩阵，则以上三式可重写为

$I_s{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{R}_{\mathrm{bak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_{\mathrm{bak}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+{\mathrm{\&omega;}}_b^{\&times;}\left(I_s{\mathrm{\&omega;}}_b\right)=T_c{+T}_d---\left(\text{13}\right)$

$I_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_{\mathrm{ak}}+F_{\mathrm{ak}}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+R_{\mathrm{bak}}^T{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b=T_{\mathrm{ak}}---\left(14\right)$

${\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}{\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}+{\mathrm{\&Lambda;}}_{\mathrm{ak}}^2{\mathrm{\&eta;}}_{\mathrm{ak}}+F_{\mathrm{bak}}^T{\stackrel{\&CenterDot;}{\mathrm{\&omega;}}}_b=0,(k=\mathrm{1,2},...,N)---\left(15\right)$

其中，Λ_ak为附件的模态频率对角阵，ξ_ak为附件的模态阻尼矩阵，R_bak为挠性附件转动对中心刚体转动的刚性耦合系数矩阵，表示为

$R_{\mathrm{bak}}=A_{\mathrm{bak}}{I}_{\mathrm{ak}}^b---\left(16\right)$

F_bak为挠性附件振动对中心刚体转动的柔性耦合系数矩阵，表示为

F_bak=A_bakH_bakΦ_ak    (17)

F_ak为挠性附件振动对自身转动的柔性耦合系数矩阵，表示为

F_ak=H_akΦ_ak    (18)

方程(13)、(14)和(15)组成了带有多个挠性附件的航天器姿态动力学模型，模型中保留了主要的非线性项、干扰力矩项。

步骤四：根据步骤二中得到的航天器姿态动力学模型，认为航天器在执行姿态机动任务时，大型挠性附件处于锁死状态，即在航天器姿态动力学方程中暂且忽略挠性附件的转动方程及其转动角速度。并且忽略环境干扰力矩。可得到简化的动力学模型，并将其写成状态方程的形式，用于计算航天器系统的模态振动参数，模态振动参数包含模态频率和阻尼比，并完成输入成形器的设计。具体操作如下：

将方程(13)和方程(15)写成如下形式

$\left[\begin{array}{cc}{I}_s{J}_{\mathrm{vb}}& F_{\mathrm{bak}}\\ F_{\mathrm{bak}}^T{J}_{\mathrm{vb}}& E\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{cc}{\mathrm{\&omega;}}_b^{\&times;}{I}_s{J}_{\mathrm{vb}}& 0\\ 0& 2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& {\mathrm{\&Lambda;}}_{\mathrm{ak}}^2\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_b\\ {\mathrm{\&eta;}}_{\mathrm{ak}}\end{array}\right]=\left[\begin{array}{c}{T}_c\\ 0\end{array}\right]---\left(19\right)$

其中，

T_c为控制力矩的表达式，本文选用工程易于实现的反馈控制器作为姿态控制器。如下所示：

$\left[\begin{array}{c}{T}_c\\ 0\end{array}\right]=\left[\begin{array}{cc}-K_{\mathrm{db}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_b\\ {\stackrel{\&CenterDot;}{\mathrm{\&eta;}}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{c}{K}_{\mathrm{db}}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_r\\ 0\end{array}\right]+\left[\begin{array}{cc}-K_{\mathrm{pb}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{\&theta;}}_b\\ {\mathrm{\&eta;}}_{\mathrm{ak}}\end{array}\right]+\left[\begin{array}{c}{K}_{\mathrm{pb}}{\mathrm{\&theta;}}_r\\ 0\end{array}\right]---\left(20\right)$

其中，θ_r为期望姿态角度，K_db和K_pb为控制参数。

令 $x={\left[\begin{array}{cc}{\mathrm{\&theta;}}_b^T& {\mathrm{\&eta;}}_{\mathrm{ak}}^T\end{array}\right]}^T,$ $X={\left[\begin{array}{cc}{x}^T& {\stackrel{\&CenterDot;}{x}}^T\end{array}\right]}^T,$ 则方程(19)和(20)可写成以下状态方程的形式

$\stackrel{\&CenterDot;}{X}=\mathrm{AX}+\mathrm{BU}---\left(21\right)$

其中，

$A=\left[\begin{array}{cc}0& E\\ -M^{-1}K& -M^{-1}C\end{array}\right];$ $B=\left[\begin{array}{c}0\\ M^{-1}T\end{array}\right];$ $M=\left[\begin{array}{cc}{I}_s{J}_{\mathrm{vb}}& F_{\mathrm{bak}}\\ F_{\mathrm{bak}}^T{J}_{\mathrm{vb}}& E\end{array}\right];$ $C=\left[\begin{array}{cc}{K}_{\mathrm{db}}& 0\\ 0& 2{\mathrm{\&xi;}}_{\mathrm{ak}}{\mathrm{\&Lambda;}}_{\mathrm{ak}}\end{array}\right]$

$K=\left[\begin{array}{cc}{K}_{\mathrm{pb}}& 0\\ 0& {\mathrm{\&Lambda;}}_{\mathrm{ak}}^2\end{array}\right];$ $T=\left[\begin{array}{c}{K}_{\mathrm{db}}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_r+K_{\mathrm{pb}}{\mathrm{\&theta;}}_r\\ 0\end{array}\right]$

根据闭环系统方程(21)，得到系统矩阵A的特征值。那么系统模态的振动频率和阻尼比的数值解由下式给出：

${\mathrm{\&lambda;}}_{\mathrm{sys}}=-{\mathrm{\&xi;}}_{\mathrm{sys}}{\mathrm{\&omega;}}_{\mathrm{sys}}\&PlusMinus;{\mathrm{j\&omega;}}_{\mathrm{sys}}\sqrt{1-{\mathrm{\&xi;}}_{\mathrm{sys}}^2}---\left(22\right)$

其中λ_sys为矩阵A的特征值，ω_sys和ξ_sys为系统模态振动频率和阻尼比。

为提高系统对参数变化的鲁棒性，并且从机动时间上折衷考虑，所述的输入成形器为ZVD(Zero Vibration and Derivative)输入成形器。针对任意一组模态振动参数的ZVD输入成形器表达式如下所示：

$\left[\begin{array}{ccc}0& t_1& t_2\\ A_1& A_2& A_3\end{array}\right]=\left[\begin{array}{ccc}0& \frac{\mathrm{\&pi;}}{{\mathrm{\&omega;}}_d}& \frac{2\mathrm{\&pi;}}{{\mathrm{\&omega;}}_d}\\ \frac{1}{1+2K+K^2}& \frac{2K}{1+2K+K^2}& \frac{K^2}{1+2K+K^2}\end{array}\right]---\left(23\right)$

其中，第一行表示各脉冲的作用时刻（其中定义第一个脉冲发生在0时刻），第二行表示各脉冲的幅值； $K=e^{-{\mathrm{\&xi;}}_{\mathrm{eys}}\mathrm{\&pi;}/\sqrt{1-{{\mathrm{\&xi;}}_{\mathrm{sys}}}^2}};$ ${\mathrm{\&omega;}}_d={\mathrm{\&omega;}}_{\mathrm{sys}}\sqrt{1-{{\mathrm{\&xi;}}_{\mathrm{sys}}}^2}.$

在使用过程中，可根据需要将多个输入成形器做卷积，合成所需要的输入成形器，即ZVD=ZVD1*ZVD2*...*ZVDn。

步骤五：将步骤一中计算得到的航天器姿态机动所绕的特征轴和特征转角、步骤二中设计得到的姿态角度轨迹规划的曲线、步骤四中设计得到的输入成形器以及反馈控制器和步骤三中推导得出的航天器姿态动力学共同构成航天器姿态控制系统完整回路，并且需要注意的是所述快速姿态机动快速稳定联合控制方法要将轨迹规划模块放置在输入成形器模块之前，并且属于姿态控制系统回路中的前馈环节。根据此完整回路，通过使用数值仿真软件搭建航天器姿态控制系统，以完成对本专利所述的航天器快速姿态机动快速稳定联合控制方法的验证。

Images