CN102343985B - 带反作用飞轮的卫星时间最优姿态机动方法 - Google Patents

Abstract

本发明公开了一种带反作用飞轮的卫星时间最优姿态机动方法，包括以下几个步骤：第一步、建立考虑反作用轮动力学的航天器姿态运动模型，在此基础上建立卫星时间最优姿态机动模型；第二步、针对考虑了反作用轮动力学的航天器姿态运动模型获取开环最优控制；第三步、获取鲁棒反馈控制器，实现航天器重定向姿态机动；本发明实现航天器的最快姿态机动，而且机动控制的精度高、鲁棒性强，能满足执行机构的力矩饱和和动量饱和约束。

Description

带反作用飞轮的卫星时间最优姿态机动方法

技术领域

本发明涉及一种航天器快速，高精度机动的控制方法，具体涉及一种带反作用飞轮的卫星时间最优姿态机动方法，属于航天器控制技术领域。

背景技术

带反作用轮的航天器时间最短重定向机动问题就是要找到一种控制使得航天器在最短的时间内实现从某一稳定的姿态机动到另一稳定的姿态。一些最优控制律是通过采用直接方法或间接方法得到的。K.D.Bilimoria，and B.Wie，“Time-Optimal Three-AxisReorientation of a Rigid Spacecraft，”Journal of Guidance Control and Dynamics，Vol.16，No.3，1993，pp.446-452.公开了一个针对刚体航天器三轴重定向的bang-bang控制，并且证明了绕特征轴的旋转不是时间最优的。H.Shen，and P.Tsiotras，“Time-Optimal Control of Axisymmetric Rigid Spacecraft Using Two Controls，”Journal of Guidance Control and Dynamics，Vol.22，No.5，1999，pp.682-694.中，仅通过两个控制实现了轴对称刚体航天器的最优机动。L.C.Lai，C.C.Yang，and C.J.Wu，“Time-Optimal Maneuvering Control of a Rigid Spacecraft，”ACTA Astronautica，Vol.60，No.10，2007，pp.791-800.中，把时间最优机动控制问题转换成非线性规划问题，把控制参数作为优化设计变量，通过遗传算法求得了最优解。M.V.Levskii，“The Problemof the Time-Optimal Control of Spacecraft Reorientation，”Journal of AppliedMathematics and Mechanics，Vol.73，No.1，2009，pp.16-25.采用庞德里亚金极大值原理求解了考虑航天器角动量约束的最短时间旋转问题。A.Fleming，P.Sekhavat，and I.M.Ross，“Minimum-Time Reorientation of a Rigid Body，”Journal of GuidanceControl and Dynamics，Vol.33，No.1，2010，pp.160-170.采用了间接方法和伪谱法求得了考虑约束的最优重定向问题。S.Liu，and T.Singh，“Fuel/Time Optimal Control ofSpacecraft Maneuvers，”Journal of Guidance Control and Dynamics，Vol.20，No.2，1996，pp.394-397.开发了STO算法并解决了在三个独立有界的脉冲控制下的燃料最优和时间最优姿态机动问题。X.Bao，and J.L.Junkins，“New Results for Time-OptimalThree-Axis Reorientation of a Rigid Spacecraft”Journal of Guidance Control andDynamics，Vol.32，No.4，2009，pp.1071-1076.的研究显示航天器绕特征轴的机动在控制输入约束的前提下是时间最优的方法。R.G.Melton，“Boundary Points and Arcs inConstrained，Time-Optimal Satellite Reorientation Maneuvers，”AIAA/AASAstrodynamics Specialist Conference，2-5August 2010，Toronto，Ontario，Canada，pp.1-16中，求解卫星时间最优重定向机动问题时考虑边界点和边界弧作为约束。

在现有技术中的最优姿态重定向机动没有将执行机构的动力学考虑在卫星的姿态动力学中，精度有待提高。

发明内容

本发明的目的是为了解决上述问题，提出一种带反作用飞轮的卫星时间最优姿态机动方法。

本发明的带反作用飞轮的卫星时间最优姿态机动方法，包括以下几个步骤：

第一步、建立考虑反作用轮动力学的航天器姿态运动模型，在此基础上建立卫星时间最优姿态机动模型；

第二步、针对考虑了反作用轮动力学的航天器姿态运动模型获取开环最优控制；

第三步、获取鲁棒反馈控制器，实现航天器重定向姿态机动；

本发明的优点在于：

(1)实现航天器的最快姿态机动；

(2)机动控制的精度高；

(3)鲁棒性强；

(4)能满足执行机构的力矩饱和和动量饱和约束。

附图说明

图1是本发明的方法流程图；

图2是本发明的有3个反作用轮的卫星布局图；

图3是本发明的姿态重定向机动；

图4是本发明的实施例中开环最优四元数曲线；

图5是本发明的实施例中开环最优角速度曲线；

图6是本发明的实施例中反作用轮的开环最优角速度曲线；

图7是本发明的实施例中开环最优控制力矩曲线；

图8是本发明的实施例中在三种控制方案下的姿态重定向机动曲线。

图中：

具体实施方式

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

本发明是一种带反作用飞轮的卫星时间最优姿态机动方法，流程如图1所示，针对反作用轮安装在惯性轴上的航天器，从一种稳定的姿态机动到另一种稳定的姿态，包括以下几个步骤：

第一步、建立考虑反作用轮动力学的航天器姿态运动模型，在此基础上建立卫星时间最优姿态机动模型；

1、建立考虑反作用轮动力学的航天器姿态运动模型；

姿态运动模型包括姿态动力学模型和姿态运动学模型。

如图2所示，为一个有三个反作用轮安装在惯性轴上的刚体卫星的布局，图中Oxbybzb为飞行器坐标系，O是飞行器的质心，反作用轮主要用于吸收周期扰动力矩，偶尔用于卫星姿态重定向机动。

姿态是通过四元数来描述的，用四元数描述的姿态运动学模型如下所示。

$\stackrel{\·}{q}=\frac{1}{2}Q\left(\mathrm{\ω}\right)\·q=\frac{1}{2}\mathrm{\Ξ}\left(q\right)\mathrm{\ω}---\left(1\right)$

其中，q＝[q₁，q₂，q₃，q₄]^T是四元数向量，q₁，q₂，q₃，q₄分别为四元数的四个分量，ω＝[ω₁，ω₂，ω₃]^T是卫星的角速度向量，ω₁、ω₂、ω₃分别为卫星的角速度向量在飞行器坐标系三个轴上的分量，Q(ω)和Ξ(q)是如下矩阵：

$\begin{array}{cc}Q\left(\mathrm{\ω}\right)=\left[\begin{array}{cc}-{\mathrm{\ω}}^{\×}& \mathrm{\ω}\\ -{\mathrm{\ω}}^T& 0\end{array}\right]& \mathrm{\Ξ}\left(q\right)=\left[\begin{array}{c}{q}_4{I}_{3\×3}+q_{13}^{\×}\\ -q_{13}^T\end{array}\right]\end{array}---\left(2\right)$

其中，I_3×3表示3×3的单位矩阵，ω^×和是反对称矩阵，如下：

$\begin{array}{cc}{\mathrm{\ω}}^{\×}=\left[\begin{array}{ccc}0& -{\mathrm{\ω}}_3& {\mathrm{\ω}}_2\\ {\mathrm{\ω}}_3& 0& -{\mathrm{\ω}}_1\\ -{\mathrm{\ω}}_2& {\mathrm{\ω}}_1& 0\end{array}\right]& \end{array}{q}_{13}^{\×}=\left[\begin{array}{ccc}0& -q_3& q_2\\ q_3& 0& -q_1\\ -q_2& q_1& 0\end{array}\right]---\left(3\right)$

考虑反作用轮动力学的刚体卫星的姿态动力学模型如下：

$\stackrel{\·}{\mathrm{\ω}}=I_s^{-1}(-{\mathrm{\ω}}^{\×}{I}_s\mathrm{\ω}-{\mathrm{\ω}}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}-T_u+T_{\mathrm{ex}})---\left(4\right)$

其中，I_s和I_RW分别是卫星和反作用轮的惯性矩矩阵，Ω反作用轮的角速度向量，Ω＝[Ω₁，Ω₂，Ω₃]^T，Ω₁，Ω₂，Ω₃分别为安装在Oxb，Oyb和Ozb轴上的反作用轮的角速度，T_u是反作用轮的力矩向量，T_u＝[T_u1 T_u2 T_u3]^T，T_u1 T_u2 T_u3分别表示安装在Oxb，Oyb和Ozb轴上的反作用轮产生的力矩，T_ex是环境扰动力矩，不考虑此项。忽略扰动力矩的姿态动力学模型(4)为：

$\stackrel{\·}{\mathrm{\ω}}=I_s^{-1}(-{\mathrm{\ω}}^{\×}{I}_s\mathrm{\ω}-{\mathrm{\ω}}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}-T_u)---\left(5\right)$

反作用轮的姿态动力学模型如下：

$\stackrel{\·}{\mathrm{\Ω}}=I_{\mathrm{RW}}^{-1}{T}_u---\left(6\right)$

将考虑反作用轮的刚体卫星的姿态运动学、动力学模型(1)(5)(6)进行整合。相应的状态变量和控制变量描述如下：

x＝[q₁ q₂ q₃ q₄ ω₁ ω₂ ω₃ Ω₁ Ω₂ Ω₃]^T，u＝T_u＝[T_u1 T_u2 T_u3]^T    (7)

其中，x是状态变量，u是控制变量。状态变量包括卫星的姿态四元数和角速度以及反作用轮的角速度。控制变量为反作用轮的力矩。姿态运动学、动力学模型(1)，(5)和(6)可以描述成统一的形式即考虑反作用轮动力学的航天器姿态运动模型：

$\stackrel{\·}{x}=f(x,u)---\left(8\right)$

表示x的导数。

2、建立卫星时间最优姿态机动模型；

姿态重定向是指将卫星从一个稳定的姿态机动到另一个稳定的姿态。在图3中，姿态角是通过飞行器坐标系Oxbybzb和LVLH坐标系Oxyz的相对旋转定义的，图中A，B，C，D表示卫星飞行过程中经历的四种不同的姿态。最少时间重定向问题就是要设计一组控制力矩以最少的时间实现从两个不同的稳定姿态之间的机动。在最优姿态重定向机动过程中，执行机构的能力必须考虑。在此工作中，反作用轮是控制力矩的提供者，执行机构动力学以及反作用轮的力矩和动量饱和均考虑在卫星最优姿态重定向机动模型中。本发明的控制力矩能够保证在满足最大控制力矩和最大动量约束的前提下实现最短时间的姿态机动。

卫星姿态重定向机动模型的初始稳定状态如下：

$\begin{array}{cccc}{q}_1\left(t_0\right)=q_{{1t}_0}& q_2\left(t_0\right)=q_{{2t}_0}& q_3\left(t_0\right)=q_{{3t}_0}& q_4\left(t_0\right)=q_{{4t}_0}\end{array}$

$\begin{array}{cccc}{\mathrm{\ω}}_1\left(t_0\right)={\mathrm{\ω}}_{{1t}_0}& {\mathrm{\ω}}_2\left(t_0\right)={\mathrm{\ω}}_{{2t}_0}& {\mathrm{\ω}}_3\left(t_0\right)={\mathrm{\ω}}_{{3t}_0}& ---\left(9\right)\end{array}$

$\begin{array}{ccc}{\mathrm{\Ω}\mathrm{}}_1\left(t_0\right)={\mathrm{\Ω}}_{{1t}_0}& {\mathrm{\Ω}}_2\left(t_0\right)={\mathrm{\Ω}}_{{2t}_0}& {\mathrm{\Ω}}_3\left(t_0\right)={\mathrm{\Ω}}_{{3t}_0}\end{array}$

其中：

均为常数值；

q₁(t₀)、q₂(t₀)、q₃(t₀)、q₄(t₀)为初始时刻姿态四元数的值；ω₁(t₀)、ω₂(t₀)、ω₃(t₀)为卫星转动角速度三个分量在初始时刻的值；Ω₁(t₀)、Ω₂(t₀)、Ω₃(t₀)为三个反作用轮在初始时刻的转动角速度的值。

卫星姿态重定向机动模型的终端时刻的稳定状态如下：

$\begin{array}{cccc}{q}_1\left(t_f\right)=q_{{1t}_f}& q_2\left(t_f\right)=q_{{2t}_f}& q_3\left(t_f\right)=q_{{3t}_f}& q_4\left(t_f\right)=q_{{4t}_f}\end{array}$

$\begin{array}{cccc}{\mathrm{\ω}}_1\left(t_f\right)={\mathrm{\ω}}_{{1t}_f}& {\mathrm{\ω}}_2\left(t_f\right)={\mathrm{\ω}}_{{2t}_f}& {\mathrm{\ω}}_3\left(t_f\right)={\mathrm{\ω}}_{{3t}_f}& \end{array}---\left(10\right)$

$\begin{array}{ccc}{\mathrm{\Ω}}_1\left(t_f\right)={\mathrm{\Ω}}_{{1t}_f}& {\mathrm{\Ω}}_2\left(t_f\right)={\mathrm{\Ω}}_{{2t}_f}& {\mathrm{\Ω}}_3\left(t_f\right)={\mathrm{\Ω}}_{{3t}_f}\end{array}$

其中：

均为常数值；

q₁(t_f)、q₂(t_f)、q₃(t_f)、q₄(t_f)为终端时刻姿态四元数的值；ω₁(t_f)、ω₂(t_f)、ω₃(t_f)为卫星转动角速度三个分量在终端时刻的值；Ω₁(t_f)、Ω₂(t_f)、Ω₃(t_f)为三个反作用轮在终端时刻的转动角速度的值。

在重定向问题中，卫星在初始时刻和终端时刻的角速度都为0，则

${\mathrm{\ω}}_{1t_0}={\mathrm{\ω}}_{2t_0}={\mathrm{\ω}}_{3t_0}={\mathrm{\ω}}_{1t_f}={\mathrm{\ω}}_{2t_f}={\mathrm{\ω}}_{3t_f}=0.$

状态方程考虑了执行机构的动力学。反作用轮的最大力矩和最大动量分别是控制约束和状态约束。反作用轮的角动量的限制可以转换成角速度的限制，如下：

$\begin{array}{ccc}\left|{\mathrm{\Ω}}_1\right|\≤\stackrel{\‾}{\mathrm{\Ω}}& \left|{\mathrm{\Ω}}_2\right|\≤\stackrel{\‾}{\mathrm{\Ω}}& \left|{\mathrm{\Ω}}_3\right|\≤\stackrel{\‾}{\mathrm{\Ω}}\end{array}---\left(11\right)$

在此，反作用轮的最大角速度。最大控制力矩约束如下：

$\begin{array}{ccc}\left|T_{u1}\right|\≤{\stackrel{\‾}{T}}_u& \left|T_{u2}\right|\≤{\stackrel{\‾}{T}}_u& \left|T_{u3}\right|\≤{\stackrel{\‾}{T}}_u\end{array}---\left(12\right)$

在此，

反作用轮的最大力矩。卫星角速度也需要满足一定的约束，描述如下：

$\begin{array}{ccc}\left|{\mathrm{\ω}}_1\right|\≤\stackrel{\‾}{\mathrm{\ω}}& \left|{\mathrm{\ω}}_2\right|\≤\stackrel{\‾}{\mathrm{\ω}}& \left|{\mathrm{\ω}}_3\right|\≤\stackrel{\‾}{\mathrm{\ω}}\end{array}---\left(13\right)$

是最大角速度。根据四元数的定义，四元数分量必须满足一下的条件。

$q_1^2+q_2^2+q_3^2+q_4^2=1---\left(14\right)$

显然，|q_i|≤1(i＝1，2，3，4)。式(8)在优化问题中也是一项等式约束。为了得到时间最短重定向机动问题的最优控制，给出时间最短问题的最优性能指标如下：

$\underset{u}{\min }J=\underset{u}{\min }(t_f-t_0)---\left(15\right)$

t₀是初始时间，t_f是终端时间。

方程(8)-(15)就构成了卫星时间最优姿态机动模型。

第二步、针对考虑了反作用轮动力学的航天器姿态运动模型获取开环最优控制；

(1)对航天器姿态运动模型进行归一化处理

状态变量和控制变量归一化如下：

$\begin{array}{ccc}\widetilde{\mathrm{\ω}}=\frac{\mathrm{\ω}}{\stackrel{\‾}{\mathrm{\ω}}}& \widetilde{\mathrm{\Ω}}=\frac{\mathrm{\Ω}}{\stackrel{\‾}{\mathrm{\Ω}}}& {\tilde{T}}_u=\frac{T_u}{{\stackrel{\‾}{T}}_u}\end{array}---\left(16\right)$

或者

$\begin{array}{ccc}\mathrm{\ω}=\stackrel{\‾}{\mathrm{\ω}}\widetilde{\mathrm{\ω}}& \mathrm{\Ω}=\stackrel{\‾}{\mathrm{\Ω}}\widetilde{\mathrm{\Ω}}& T_u={\stackrel{\‾}{T}}_u{\tilde{T}}_u\end{array}---\left(17\right)$

用归一化变量描述的航天器姿态运动模型描述如下：

$\left\{\begin{array}{c}\stackrel{\·}{q}=\frac{1}{2}\stackrel{\‾}{\mathrm{\ω}}Q\left(\widetilde{\mathrm{\ω}}\right)\·q\\ \stackrel{\stackrel{\·}{~}}{\mathrm{\ω}}=I_s^{-1}\left(-\stackrel{\‾}{\mathrm{\ω}}{\widetilde{\mathrm{\ω}}}^{\×}{I}_s\widetilde{\mathrm{\ω}}-\stackrel{\‾}{\mathrm{\Ω}}{\widetilde{\mathrm{\ω}}}^{\×}{I}_{\mathrm{RW}}\widetilde{\mathrm{\Ω}}-\frac{{\stackrel{\‾}{T}}_u}{\stackrel{\‾}{\mathrm{\ω}}}{\widetilde{\stackrel{\‾}{T}}}_u\right)\\ \stackrel{\stackrel{\·}{~}}{\mathrm{\Ω}}=\frac{{\stackrel{\‾}{T}}_u}{\stackrel{\‾}{\mathrm{\Ω}}}{I}_{\mathrm{RW}}^{-1}{\widetilde{\stackrel{\‾}{T}}}_u\end{array}\right.---\left(18\right)$

四元数的最大数值为1，同时所有四元数的分量都位于区间[-1，1]。

(2)用归一化参数描述的最优控制问题

对应于最优控制问题，归一化变量描述的动力学方程(18)如下所示。

$\stackrel{\stackrel{\·}{~}}{x}=f(\tilde{x},\tilde{u})---\left(19\right)$

其中

$\tilde{x}={\left[\begin{array}{cccccccccc}{q}_1& q_2& q_3& q_4& {\widetilde{\mathrm{\ω}}}_1& {\widetilde{\mathrm{\ω}}}_2& {\widetilde{\mathrm{\ω}}}_3& {\widetilde{\mathrm{\Ω}}}_1& {\widetilde{\mathrm{\Ω}}}_2& {\widetilde{\mathrm{\Ω}}}_3\end{array}\right]}^T$

                    (20)

最优控制的性能指标描述如下：

J＝t_f-t₀                                (21)

不等式约束如下：

$\begin{array}{cc}{\tilde{x}}_{\min }\≤\tilde{x}\≤{\tilde{x}}_{\max }& {\tilde{u}}_{\min }\≤\stackrel{\text{~}}{u}\≤{\tilde{u}}_{\max }\end{array}---\left(22\right)$

其中， ${\tilde{x}}_{\max }=\left[\begin{array}{cccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1& 1\end{array}\right],{\tilde{x}}_{\min }=-{\tilde{x}}_{\max },{\tilde{u}}_{\max }=\left[\begin{array}{ccc}1& 1& 1\end{array}\right],{\tilde{u}}_{\min }=-{\tilde{u}}_{\max }.$

等式约束见方程(14)。优化的目标是要找到一个归一化的控制

使得卫星姿态重定向的时间最短。

(3)采用勒让德伪谱法将最优控制问题转化成非线性规划问题

采用勒让德伪谱法离散前面描述的最优控制问题，此方法是建立在用拉格朗日插值多项式估计状态变量和控制变量的基础上。未知的系数是插值节点变量的值，叫做求积点，或叫勒让德高斯兰伯特(LGL)点。因为LGL点位于区间[-1，1]，前面描述的最优控制问题描述在时间区间[t₀，t_f]，所以我们采用如下形式对LGL区间和物理时间区间进行转换：τ∈[τ₀，τ_N]＝[-1，1]

$t=\frac{(t_f-t_0)\mathrm{\τ}+(t_f+t_0)}{2}---\left(23\right)$

归一化的姿态运动模型如下：

$\stackrel{\stackrel{\·}{~}}{x}\left(\mathrm{\τ}\right)=\frac{t_f-t_0}{2}f(\tilde{x}\left(\mathrm{\τ}\right),\widetilde{u\left(\mathrm{\τ}\right)})---\left(24\right)$

$\begin{array}{cc}\tilde{x}(-1)={\tilde{x}}_{t_0}& \tilde{x}\left(1\right)={\tilde{x}}_{t_f}\end{array}$

其中：

表示τ对应时刻的归一化的状态变量的导数、归一化的状态变量和控制变量的值，

表示初始时刻和终端时刻归一化的状态变量的值，均为常值。经过转换后，通过N阶多项式的形式来估计连续状态变量和控制变量，如下所示。

$\begin{array}{cc}\tilde{x}\≈{\tilde{x}}^N\left(\mathrm{\τ}\right)=\underset{t=0}{\overset{N}{\mathrm{\Σ}}}{\tilde{x}}_l{\mathrm{\φ}}_1\left(\mathrm{\τ}\right)& \tilde{u}\≈{\tilde{u}}^N\left(\mathrm{\τ}\right)=\underset{t=0}{\overset{N}{\mathrm{\Σ}}}{\tilde{u}}_l{\mathrm{\φ}}_l\left(\mathrm{\τ}\right)\end{array}---\left(25\right)$

其中，l＝0，1，…，N，N表示一选定的正整数，

表示N个点拟合的状态变量和控制变量在τ对应时刻的值。

${\mathrm{\φ}}_l\left(\mathrm{\τ}\right)=\frac{1}{N(N+1)L_N\left({\mathrm{\τ}}_l\right)}\frac{({\mathrm{\τ}}^2-1){\stackrel{\·}{L}}_N\left(\mathrm{\τ}\right)}{\mathrm{\τ}-{\mathrm{\τ}}_l}=\left\{\begin{array}{ccc}1& \mathrm{if}& l=j\\ 0& \mathrm{if}& l\≠j\end{array}\right.---\left(26\right)$

上式是N阶的拉格朗日插值多项式，L_N(τ_l)勒让德多项式。

${\tilde{x}}_l={\tilde{x}}^N\left({\mathrm{\τ}}_l\right),{\tilde{u}}_l={\tilde{u}}^N\left({\mathrm{\τ}}_l\right),{\mathrm{\τ}}_l=\mathrm{\τ}\left(t_l\right)---\left(27\right)$

其中：t_l表示第l个节点对应的时刻，τ_l表示第l个节点对应的τ的值。表示第l个节点对应的归一化的状态变量和控制变量的值。为了根据在节点τ_l的值

来表达状态变量的导数

相应的状态方程(19)可描述成以下形式：

${\stackrel{\stackrel{\·}{~}}{x}}^N\left({\mathrm{\τ}}_k\right)=\underset{t=0}{\overset{N}{\mathrm{\Σ}}}{D}_{\mathrm{kl}}\tilde{x}\left({\mathrm{\τ}}_l\right)---\left(28\right)$

其中D_kl是(N+1)×(N+1)的差分矩阵D的分量：

$D=\left[D_{\mathrm{kl}}\right]=\left\{\begin{array}{cc}\frac{L_N\left({\mathrm{\τ}}_k\right)}{L_N\left({\mathrm{\τ}}_l\right)}\frac{1}{({\mathrm{\τ}}_k-{\mathrm{\τ}}_l)}& k\≠l\\ -\frac{N(N+1)}{4}& k=l=0\\ \frac{N(N+1)}{4}& k=l=N\\ 0& \mathrm{otherwise}\end{array}\right.---\left(29\right)$

因而，在最优控制问题中的状态方程的等式约束可以通过离散状态描述为：

$\underset{l=0}{\overset{N}{\mathrm{\Σ}}}{D}_{\mathrm{kl}}\tilde{x}\left({\mathrm{\τ}}_l\right)-\frac{t_f-t_0}{2}f\left(\underset{l=0}{\overset{N}{\mathrm{\Σ}}}\tilde{x}\left({\mathrm{\τ}}_l\right){\mathrm{\φ}}_l\left({\mathrm{\τ}}_k\right),\underset{l=0}{\overset{N}{\mathrm{\Σ}}}\tilde{u}\left({\mathrm{\τ}}_l\right){\mathrm{\φ}}_l\left({\mathrm{\τ}}_k\right)\right)=0---\left(30\right)$

其它的等式约束见方程(14)，不等式约束如下所示。

$\begin{array}{cc}{\tilde{x}}_{\min }\≤\tilde{x}\left({\mathrm{\τ}}_l\right)\≤{\tilde{x}}_{\max }& {\tilde{u}}_{\min }\≤\tilde{u}\left({\mathrm{\τ}}_l\right)\≤{\tilde{u}}_{\max }\end{array}---\left(31\right)$

最优控制问题的性能指标见方程(21)，优化变量是和

这样最优控制问题就转化成了一个非线性规划问题。

(4)优化计算求解非线性规划问题，得到开环最优姿态机动参数

通过利用一些数值优化软件例如SNOPT或matlab就可以求解此优化问题了，从而得到了最优姿态。

第三步，设计鲁棒反馈控制器，实现航天器重定向姿态机动；

上面所述步骤通过求解最优控制问题得到了最优姿态机动的开环控制。在实际应用中，由于存在不确定性和环境扰动，例如动力学建模的不确定性，空气阻力扰动力。所以需要一个具有较好鲁棒性的反馈控制律来对参考轨迹进行跟踪。根据姿态误差方程来推导鲁棒控制器。

1、建立姿态误差方程

对刚体卫星来说，姿态动力学方程在方程(4)中已给出，姿态运动学方程在方程(1)中给出。而期望的姿态运动轨迹也通过解前面的开环最优控制问题得到。q_d和ω_d定义为期望姿态四元素和转动角速度，q_e为坐标系F_b相对期望坐标系F_d的姿态四元数。相应的从F_d到F_b的转换矩阵C(q_e)如下所示。

$C=(q_{e4}^2-q_{e13}^T{q}_{e13})I_{3\×3}+2q_{e13}{q}_{e13}^T-{2q}_{e4}{q}_{e13}^{\×}---\left(32\right)$

其中：q_e1、q_e2、q_e3、q_e4为四元数q_e的四个分量，向量q_e13＝[q_e1、q_e2、q_e3]^T，

为q_e13的转置向量，

是q_e13的反对称矩阵，q与q_d，q_e的关系描述如下：

q＝mat(q_d)q_e                                         (33)

其中

$\mathrm{mat}\left(q_d\right)=\left[\begin{array}{cccc}{q}_{d4}& q_{d3}& -q_{d2}& q_{d1}\\ -q_{d3}& q_{d4}& -q_{d1}& q_{d2}\\ q_{d2}& q_{d1}& q_{d4}& q_{d3}\\ -q_{d1}& -q_{d2}& -q_{d3}& q_{d4}\end{array}\right]---\left(34\right)$

其中：q_d1、q_d2、q_d3、q_d4为q_d的四个分量，F_b相对F_d的角速度ω_e在坐标系F_b中描述如下：

ω_e＝ω-Cω_d                                      (35)

根据ω_e和q_e重写方程(4)和方程(1)，就得到了的误差方程。

$I_s{\stackrel{\·}{\mathrm{\ω}}}_e=-{({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)}^{\×}{I}_s({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)-{({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}+I_s({\mathrm{\ω}}_e^{\×}{\mathrm{C\ω}}_d-{C\stackrel{\·}{\mathrm{\ω}}}_d)-T_u$

                                  (36)

$\stackrel{\·}{q}=\frac{1}{2}\mathrm{\Ξ}\left(q_e\right){\mathrm{\ω}}_e$

在方程(36)中考虑了反作用轮的动力学影响。

2、获取鲁棒反馈控制器

我们选择李亚普洛夫函数如下：

$V=\frac{1}{2}{{\mathrm{\ω}}_e}^T{I}_s{\mathrm{\ω}}_e+k_1(q_{e1}^2+q_{e2}^2+q_{e3}^2+{(1-q_{e4})}^2)---\left(37\right)$

其中k₁＞0。可以推导出速度V的导数：

$\stackrel{\·}{V}={\mathrm{\ω}}_e^T\{-{({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)}^{\×}{I}_s({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)-{({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}+I_s({\mathrm{\ω}}_e^{\×}{\mathrm{C\ω}}_d-{C\stackrel{\·}{\mathrm{\ω}}}_d)-T_u\}+k_1{\mathrm{\ω}}_e^T{q}_{e13}^T$

$={\mathrm{\ω}}_e^T\{-{\mathrm{\ω}}_e^{\×}(I_s({\mathrm{\ω}}_e+{\mathrm{C\ω}}_d)+I_{\mathrm{RW}}\mathrm{\Ω})-({\left(C{\mathrm{\ω}}_d\right)}^{\×}{I}_s+I_s{\left({\mathrm{C\ω}}_d\right)}^{\×}){\mathrm{\ω}}_e-{\left({\mathrm{C\ω}}_d\right)}^{\×}{I}_s{\mathrm{C\ω}}_d-I_s{C\stackrel{\·}{\mathrm{\ω}}}_d$

                                         (38)

$-{\left({\mathrm{C\ω}}_d\right)}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}-T_u\}+k_1{\mathrm{\ω}}_e^T{q}_{e13}^T$

$={\mathrm{\ω}}_e^T(-{\left({\mathrm{C\ω}}_d\right)}^{\×}{I}_s{\mathrm{C\ω}}_d-I_s{C\stackrel{\·}{\mathrm{\ω}}}_d-{\left({\mathrm{C\ω}}_d\right)}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}-T_u)+k_1{\mathrm{\ω}}_e^T{q}_{e13}^T$

根据上式提出鲁棒反馈控制器为：

$u=T_u{=(-{\mathrm{C\ω}}_d)}^{\×}{I}_s{\mathrm{C\ω}}_d-I_s{C\stackrel{\·}{\mathrm{\ω}}}_d-{\left({\mathrm{C\ω}}_d\right)}^{\×}{I}_{\mathrm{RW}}\mathrm{\Ω}+k_1{q}_{e13}+k_2{\mathrm{\ω}}_e---\left(39\right)$

其中k₂＞0。将上面的方程代入方程(38)，可得到

$\stackrel{\·}{V}=-k_2{\mathrm{\ω}}_e^T{\mathrm{\ω}}_e\≤0---\left(40\right)$

因而，系统在方程(39)所示的控制律下是稳定的。闭环跟踪的主要任务是抵消在开环最优轨迹中没有考虑的扰动力矩的影响。方程(39)的控制律中前三项是用于补偿动力学模型的误差，为一个前馈补偿，最后两项实质为一抵消不确定扰动的PD反馈项。

通过说书的控制器去跟踪前面优化求解得到的开环最优姿态，就可以实现快速、高精度的航天器重定向姿态机动。并且本发明的方法作姿态机动具有很强的鲁棒性。

实施例：

实施例中，以拥有三个反作用轮的刚体卫星为例，三个反作用轮安装在卫星的惯性轴上。卫星动力学包含了反作用轮的动力学。刚体卫星的参数如下表1所示，初始和终端状态值见表2。

表1刚体卫星的相关参数

其中，I_xx、I_yy、I_zz分别为卫星在三个惯性主轴上的转动惯量，I_w为反作用轮的转动惯量。

表2此算例的初始和终端条件

q(t₀)、ω(t₀)、Ω(t₀)、q(t_f)、ω(t_f)、Ω(t_f)分别表示在初始时刻和终端时刻姿态四元数、姿态角速度和反作用轮转动角速度的值。

通过勒让德伪谱法和利用软件TOMLAB\SNOPT优化得到了开环最优控制。最短重定向机动时间为99s。相应的开环最优四元数曲线，卫星角速度曲线，反作用轮的角速度曲线和控制力矩曲线如图4-图7所示。

在现有技术中，最优重定向问题是去发现在体轴上的三个虚拟的控制力矩，执行机构动力学没有考虑在动力学方程中。其动力学模型如下：

$\stackrel{\·}{\mathrm{\ω}}=I_s^{-1}(-{\mathrm{\ω}}^{\×}{I}_s\mathrm{\ω}+T_u)---\left(41\right)$

和方程(5)比较，公式中的-ω^×I_RWΩ在方程(41)中没有考虑。

为了比较不同方法的结果，以方程(48)作为动力学方程的开环最短时间机动问题也通过勒让德伪谱法和利用软件SNOPT优化求解。在其他参数与前面一致的情况下，其相应的仿真结果结果显示最短机动时间为88秒。机动时间比本发明机动时间更短，主要是因为模型中忽略了反作用轮的影响，并且没有考虑反作用轮动量饱和的限制，所以由于模型的不够准确，此结果并非最优，可作次优跟踪用。接着采用非线性预测控制方法(NPC)对此次优轨迹进行跟踪。在跟踪控制时考虑了力矩和动量约束，相应的参数如图8中间栏所示，姿态四元数达到稳定的时间为130秒，卫星角速度达到稳定的时间为145秒，三个轮子都正常工作，仅一个轮子在姿态旋转过程中力矩和动量都达到了饱和。

同样为了对比，采用四元数反馈(QFC)的方法来求解同样的姿态重定向问题。模型中考虑了反作用轮的力矩和角动量饱和约束。相应的结果如图8左栏所示，此方法为特征轴旋转机动，姿态四元数达到稳定状态的时间为132秒，卫星角速度达到稳定状态的时间为146秒，在姿态机动过程中，仅一个反作用轮工作，其它两个轮子处于静止状态，所以卫星近绕一根惯性轴旋转。

本发明提出的最短时间重定向机动的结果图8右栏所示，采用的跟踪控制律见公式(39)。在姿态机动过程中，三个反作用轮均工作，有一个力矩达到过饱和，还有一个动量达到过饱和。姿态四元数达到稳定状态的时间为99秒，角速度达到稳定状态的时间为100秒，所以姿态重定向机动时间为100秒，比四元数反馈方法缩短了31.5％，比次优化和非线性预测控制方案缩短了31％。相应的对比数据如表3所示。

表3三种方法的性能比较

Claims (1)

1.一种带反作用轮的卫星时间最优姿态机动方法，针对反作用轮安装在惯性轴上的卫星，从一种稳定的姿态机动到另一种稳定的姿态，其特征在于，包括以下几个步骤： 

第一步、建立考虑反作用轮动力学的卫星姿态运动模型，在此基础上建立卫星时间最优姿态机动模型； 

（1）建立考虑反作用轮动力学的卫星姿态运动模型； 

卫星姿态运动模型包括卫星姿态动力学模型和卫星姿态运动学模型； 

卫星姿态运动学模型为： 

其中，q=[q₁,q₂,q₃,q₄]^T是四元数向量，q₁,q₂,q₃,q₄分别为四元数向量的四个分量，ω＝[ω₁,ω₂,ω₃]^T是卫星的角速度向量，ω₁、ω₂、ω₃分别为卫星的角速度向量在飞行器坐标系三个轴上的分量，Q(ω)和Ξ(q)为： 

其中，q₁₃=[q₁,q₂,q₃]^T，

为q₁₃的转置向量，

是q₁₃的反对称矩阵，I_3×3表示3×3的单位矩阵，ω^×是ω反对称矩阵，如下： 

卫星姿态动力学模型为： 

其中，I_s和I_RW分别是卫星和反作用轮的惯性矩矩阵，Ω反作用轮的角速度向量,Ω＝[Ω₁,Ω₂,Ω₃]^T，Ω₁,Ω₂,Ω₃为三个分别安装在Oxb,Oyb和Ozb轴上的反作用轮的角速度，T_u是反作用轮的力矩向量,T_u＝[T_u1 T_u2 T_u3]^T，T_u1T_u2T_u3表示三个分别安装在Oxb,Oyb和Ozb轴上的反作用轮产生的力矩，T_ex是环境扰动力矩，不考虑此项；忽略扰动力矩的卫星姿态动力学模型（4）为： 

反作用轮的姿态动力学模型为： 

通过式（1）、（5）、（6）得到考虑反作用轮动力学的卫星姿态运动模型为： 

其中：x是状态变量，u是控制变量；状态变量包括卫星的姿态四元数和角速度以及反 作用轮的角速度；控制变量为反作用轮的力矩；状态变量和控制变量为： 

x＝[q₁ q₂ q₃ q₄ ω₁ ω₂ ω₃ Ω₁ Ω₂ Ω₃]^T,u＝T_u＝[T_u1 T_u2 T_u3]^T    （8） 

（2）建立卫星时间最优姿态机动模型； 

卫星姿态重定向机动模型的初始稳定状态为： 

其中，

为初始时刻姿态四元数的值；

为卫星角速度三个分量在初始时刻的值；

为三个反作用轮在初始时刻的角速度的值； 

卫星姿态重定向机动模型的终端时刻的稳定状态为： 

其中，

为终端时刻姿态四元数的值；

为卫星角速度三个分量在终端时刻的值；

为三个反作用轮在终端时刻的角速度的值； 

卫星在初始时刻和终端时刻的角速度都为0，则 

反作用轮的角动量的限制转换成角速度的限制，如下： 

其中，

为反作用轮的最大角速度；最大控制力矩约束如下： 

其中，

为反作用轮的最大力矩；卫星角速度的约束，如下： 

其中，是最大角速度；四元数满足： 

|q_i|≤1，i＝1,2,3,4；式（7）也是一项等式约束；卫星时间最优姿态机动模型的最优性能指标为： 

其中：t₀是初始时间，t_f是终端时间；式（7）-（15）构成卫星时间最优姿态机动模型； 

第二步、针对考虑了反作用轮动力学的卫星姿态运动模型获取开环最优控制； 

（1）对卫星姿态运动模型进行归一化处理 

状态变量和控制变量归一化为： 

或者 

表示

的反对称矩阵，归一化处理后的卫星姿态运动模型为： 

四元数的最大数值为1，同时所有四元数都位于区间[-1,1]； 

（2）用归一化参数描述的最优控制问题； 

归一化变量描述的卫星姿态运动模型为； 

其中 

             （20） 

最优控制的性能指标为： 

J＝t_f-t₀             （21） 

不等式约束为： 

其中，

等式约束为

优化的目标为：找到一个归一化的控制

使得卫星姿态重定向的时间最短； 

（3）采用勒让德伪谱法将最优控制问题转化成非线性规划问题 

最优控制问题描述在时间区间[t₀,t_f]，采用如下形式对勒让德高斯兰伯特区间和物理时间区间进行转换:τ∈[τ₀,τ_N]＝[-1,1] 

归一化的卫星姿态运动模型如下： 

       （24） 

其中：

表示τ对应时刻的归一化的状态变量的导数、归一化的状态变量和控 制变量的值，

表示初始时刻和终端时刻归一化的状态变量的值，

均为常值，经过转换后，通过N阶多项式的形式来估计连续状态变量和控制变量，如下所示； 

其中，l＝0,1,...,N，N表示一选定的正整数，表示N个点拟合的状态变量和控制变量在τ对应时刻的值，

上式是N阶的拉格朗日插值多项式，L_N(τ_l)勒让德多项式； 

其中，t_l表示第l个节点对应的时刻，τ_l表示第l个节点对应的τ的值，表示第l个节点对应的归一化的状态变量和控制变量的值，为了根据在节点τ_l的值来表达状态变量的导数 

相应的卫星姿态运动模型式（19）为： 

其中：D_kl是(N+1)×(N+1)的差分矩阵D的分量： 

因而，在最优控制问题中的状态方程的等式约束通过离散状态表示： 

其它的等式约束为

不等式约束为： 

最优控制问题的性能指标为式（21），优化变量是

和

最优控制问题转化为非线性规划问题； 

（4）优化计算求解非线性规划问题，得到开环最优姿态机动参数； 

第三步、获取鲁棒反馈控制器，实现卫星重定向姿态机动； 

上面所述步骤通过求解最优控制问题得到了最优姿态机动模型的开环控制，根据姿态误差方程获取鲁棒反馈控制器； 

（1）建立姿态误差方程； 

q_d、ω_d为期望姿态四元数向量和角速度向量，q_e为坐标系F_b相对期望坐标系F_d的姿态 四元数向量；相应的从F_d到F_b的转换矩阵C(q_e)为； 

其中：q_e1、q_e2、q_e3、q_e4为四元数向量q_e的四个分量，向量q_e13=[q_e1、q_e2、q_e3]^T，

为q_e13的转置向量，

是q_e13的反对称矩阵，q与q_d,q_e的关系为： 

q＝mat(q_d)q_e          （33） 

其中 

其中：q_d1、q_d2、q_d3、q_d4为q_d的四个分量，F_b相对F_d的角速度ω_e在坐标系F_b中为： 

ω_e＝ω-Cω_d       （35） 

根据ω_e和q_e重写式（4）和式（1），则姿态误差方程为： 

（2）获取鲁棒反馈控制器 

李亚普洛夫函数为： 

其中：k₁＞0；可以得到速度V的导数为： 

    （38） 

根据式（38），则鲁棒反馈控制器为： 

通过鲁棒反馈控制器跟踪第二步（4）得到开环最优姿态机动参数，实现卫星重定向姿态机动。 

Images