CN106218922A - 挠性敏捷卫星的联合执行机构控制方法 - Google Patents

Abstract

挠性敏捷卫星的联合执行机构控制方法，涉及卫星姿态控制技术领域，解决现有针对挠性卫星姿态进行快速机动或高精度跟踪控制时，具有强非线性、受多种约束以及挠性附件易振动的特点，进而导致难以满足控制需求等问题，针对联合两种执行机构的挠性敏捷卫星，结合挠性卫星姿态动力学，运动学和挠性附件的振动方程，建立了面向卫星姿态最优轨迹规划及预测未来姿态信息的非线性状态空间方程；在建立了姿态机动快速性、挠性附件振动抑制性、CMG群奇异性等综合指标基础上，基于伪光谱方法，实现对原连续优化问题的离散化，进而计算出卫星姿态机动的最优轨迹及CMG群框架角速度的最优轨迹；基于非线性模型预测控制技术，设计了反作用飞轮的控制。

Description

挠性敏捷卫星的联合执行机构控制方法

技术领域

本发明涉及卫星姿态控制技术领域，具体涉及一种挠性敏捷卫星的联合执行机构控制方法，特别适用于具有挠性附件卫星姿态的大角度快速机动及高精度跟踪控制过程。

背景技术

目前，我国遥感光学卫星通常采用星下点及侧摆的推扫成像模式，其要求卫星姿态控制系统具有一定的姿态机动性及较高的机动后姿态稳定性。为解决此问题，现有技术主要是通过预留足够长的机动及稳定时间，以克服挠性附件或干扰力矩等造成的姿态波动。但这极大地限制了遥感卫星的应用效能。为了缩短与国外遥感卫星技术水平的差距，使我国遥感卫星具备如多点目标成像、凝视成像、非沿轨灵巧成像等能力，设计高效、先进的姿态控制算法至关重要！

新颖的成像模式需要卫星姿态控制系统在规定时间内完成卫星姿态的快速机动并保持较高精度的稳定性，或实现对期望姿态及姿态角速度的高精度跟踪控制，且对各种建模不确定性、刚柔耦外部干扰力合作用及矩等具有较强的鲁棒性，这对卫星姿态控制系统的设计提出了新的挑战。当前，以控制力矩陀螺(简记CMG)为主要执行机构的卫星姿态机动控制受到广泛研究，但其高效的操纵律设计存在难点，限制了CMG的应用。反作用飞轮虽然力矩输出精度较高，但相对于卫星大惯量来讲，不足以实现快速机动和跟踪控制。而同时联合这两种执行机构的控制策略通常是分时的，其目的是保护反作用飞轮角动量不饱和，这种策略并没有实现协调优化控制，效率较低。

针对以上问题，结合CMG的大力矩输出特性及反作用飞轮的高精度力矩输出特性，设计先进的协同控制策略和算法，是实现挠性卫星姿态快速机动和高精度跟踪控制的有效途径之一。该技术值得深入研究。

发明内容

本发明为解决现有针对挠性卫星姿态进行快速机动或高精度跟踪控制时，具有强非线性、受多种约束以及挠性附件易振动的特点，进而导致难以满足控制需求等问题，提供一种挠性敏捷卫星的联合执行机构控制方法。

挠性敏捷卫星的联合执行机构控制方法，该方法由以下步骤实现：

步骤一、选择惯性坐标系为参考坐标系，建立以金字塔构型CMG群和反作用飞轮为执行机构的挠性卫星姿态动力学及运动学方程；以挠性附件的模态坐标系为基准，建立挠性附件振动的动力学方程，通过定义新的状态变量，基于建立的三个方程获得用于CMG群框架角速度最优轨迹规划及预测卫星姿态未来信息的非线性连续状态空间方程，通过离散化方法，获得离散非线性状态空间方程；

步骤二、对步骤一建立的非线性连续状态空间方程，建立兼顾卫星姿态机动快速性、挠性附件振动抑制性能、金字塔构型CMG群奇异性及满足各种约束等的优化控制问题；基于Legendre伪谱法，实现对优化控制问题的离散化求解，获得卫星姿态机动的最优状态轨迹及CMG群的最优框架角速度轨迹；

步骤三、将步骤二获得的CMG群的最优框架角速度轨迹带入步骤一建立的离散非线性状态空间方程，并根据当前测量的卫星姿态信息，建立卫星姿态的预测输出方程，实现对预测时域内的卫星姿态信息的预测；

步骤四、根据步骤三建立的预测输出方程，以步骤二获得的卫星姿态机动的最优状态轨迹为跟踪目标，建立兼顾卫星姿态跟踪快速性及鲁棒性的优化控制问题，经过对该优化控制问题的求解，获得反作用飞轮的非线性模型预测机动控制力矩；将规划得到的CMG群的最优框架角速度及飞轮机动控制力矩作用于卫星，驱动卫星姿态运动；

步骤五、采用离散控制技术，在每个采样时刻重复步骤三和步骤四，通过逐步更新卫星姿态信息，实现挠性卫星姿态机动的复合滚动控制。

本发明的有益效果：本发明通过将非线性模型预测控制方法及Legendre伪谱法相互融合，基于轨迹规划及滚动跟踪控制思想，提出一种联合CMG群及反作用飞轮实现挠性卫星姿态大角度快速机动及高精度跟踪的复合控制策略。

一、联合CMG群及反作用飞轮，建立了包含卫星姿态动力学、运动学及挠性附件振动的非线性状态空间方程，用于进行卫星姿态最优规划及对未来姿态信息的准确预测。

二、本发明为兼顾卫星姿态机动性能、附件振动的抑制性能、CMG群奇异性及成像约束等因素，建立了挠性卫星敏捷姿态机动的最优控制问题；以Legendre伪谱法手段，实现挠性卫星姿态机动最优轨迹及CMG群最优框架角速度的求取。

三、根据当前测量的卫星姿态信息，以建立的非线性状态空间方程及规划的CMG群最优框架角速度为依据，通过迭代计算确立卫星姿态的预测输出方程，实现对预测时域内卫星姿态信息的准确预测。

四、基于非线性模型预测控制技术，以规划得到的卫星姿态最优轨迹快速跟踪为目标，建立了兼顾卫星姿态跟踪快速性及鲁棒性的优化控制问题，经过求解获得反作用飞轮的姿态机动控制力矩。综合CMG群最优框架角速度及反作用飞轮控制力矩，实现挠性卫星的姿态控制。

附图说明

图1为本发明所述的挠性敏捷卫星的联合执行机构控制方法的原理框图；

图2为本发明所述的挠性敏捷卫星的联合执行机构控制方法的金字塔构型CMG群坐标系示意图；其中，分别为CMG群的四个框架轴，β为安装倾角，为沿卫星本体坐标系三轴的基矢量，为四个CMG的输出角动量；

图3为本发明所述的挠性敏捷卫星的联合执行机构控制方法的流程图；

图4中(a)(b)(c)和(d)分别为卫星姿态四元数分别为q0、q1、q2和q3时变化曲线示意图；

图5中(a)(b)(c)分别为姿态角速度的变化曲线示意图；

图6中(a)(b)(c)分别为合控制力矩的变化曲线示意图。

具体实施方式

具体实施方式一、结合图1至图4说明本实施方式，挠性敏捷卫星的联合执行机构控制方法，该方法由以下步骤实现：

步骤A:选择惯性坐标系为参考坐标系，建立以金字塔构型CMG群和反作用飞轮为执行机构的挠性卫星姿态动力学及运动学方程；以挠性附件的模态坐标系为基准，建立挠性附件振动的动力学方程；

步骤B:定义新的状态变量，基于建立的三个方程获得用于CMG群框架角速度最优轨迹规划及预测卫星姿态未来信息的非线性连续状态空间方程，通过离散化方法，获得该方程的离散化形式；

步骤C:对步骤B建立的非线性连续状态空间方程，在反作用飞轮控制力矩较小情况下，建立兼顾卫星姿态机动快速性、挠性附件振动抑制性能、金字塔构型CMG群奇异性及满足状态约束等的优化控制问题；

步骤D:基于Legendre伪谱法，采用拉格朗日全局插值多项式在一系列离散的时间节点上近似状态变量和控制变量，并通过插值多项式的导数来逼近动力学方程中状态变量对时间的导数，实现对优化控制问题的离散化求解，进而获得卫星姿态机动的最优状态轨迹及CMG群的最优框架角速度轨迹；

步骤E:将步骤D获得的CMG群最优框架角速度轨迹带入到步骤B建立的离散非线性状态空间方程，并根据当前测量的卫星姿态信息，建立卫星姿态的预测输出方程，实现对预测时域内的卫星姿态信息的预测；

步骤F:根据步骤E建立的预测输出方程，以步骤D规划得到的卫星姿态最优状态轨迹为跟踪目标，建立兼顾卫星姿态跟踪快速性及鲁棒性的优化控制问题，经过对该优化控制问题的求解，获得反作用飞轮的非线性模型预测机动控制力矩；

步骤G:将当前时刻规划得到的CMG群最优框架角速度及飞轮控制力矩作用于卫星，驱动卫星姿态运动，并更新卫星姿态下一时刻测量信息；

步骤H:采用离散控制技术，在每个采样时刻重复步骤E至步骤G，通过滚动优化，实现挠性卫星姿态机动的高精度复合控制。

具体实施方式二、结合图1至图6说明本实施方式，本实施方式为具体实施方式一所述的挠性敏捷卫星的联合执行机构控制方法的实施例，其具体过程为：

一、以惯性坐标系为参考坐标系，建立的具有联合执行机构的挠性卫星姿态动力学为：

$I\stackrel{\·}{w}+{\mathrm{\σ}}^T\stackrel{\·\·}{\mathrm{\η}}+\[w\×\]Iw+\[w\×\]{\mathrm{\σ}}^T\stackrel{\·}{\mathrm{\η}}=-{\stackrel{\·}{H}}_{CMG}-\[w\×\]H_{CMG}-{\stackrel{\·}{H}}_F-\[w\×\]H_F+T_d$

式中，I为卫星的转动惯量矩阵，w为星体三轴姿态角速度，σ为挠性附件与星体的刚柔耦合矩阵，η为挠性附件在模态坐标系下的位移，H_CMG为金字塔构型CMG群的三轴角动量，H_F为飞轮系统的三轴角动量，T_d为空间干扰力矩。定义CMG群的控制力矩为飞轮系统的控制力矩为其中[w×]定义为：

$\[w\×\]=\left[\begin{array}{ccc}0& -w_z& w_y\\ w_z& 0& -w_x\\ -w_y& w_x& 0\end{array}\right]$

在模态坐标系下，挠性附件振动的动力学方程为：

$\stackrel{\·\·}{\mathrm{\η}}+2{\mathrm{\ζ}}_f{w}_f\stackrel{\·}{\mathrm{\η}}+w_f^2\mathrm{\η}+\mathrm{\σ}\stackrel{\·}{w}=0$

式中，ζ_f为挠性附件模态的阻尼比矩阵，w_f为挠性附件模态的振动频率矩阵。

考虑到以欧拉角描述卫星姿态存在奇异性的问题，采用四元数表示的运动学模型来表示卫星的姿态变化。选取惯性系为参考坐标系，则有如下的基于四元数的卫星姿态运动学描述：

$\stackrel{\·}{q}=\frac{1}{2}\left[\begin{array}{cccc}0& -w_x& -w_y& -w_z\\ w_x& 0& w_z& -w_y\\ w_y& -w_z& 0& w_x\\ w_z& w_y& -w_x& 0\end{array}\right]\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\end{array}\right]$

其中，q＝[q₀ q₁ q₂ q₃]^T为姿态四元数，整星的绝对角速度在三轴分量为w_x，w_y和w_z分别为滚动轴姿态角速度、俯仰轴姿态角速度和偏航轴姿态角速度。

考虑到通常卫星惯量主要在其主轴上，即取星体的转动惯量矩阵为I＝diag(I_x,I_y,I_z)。且假设三轴飞轮沿惯量主轴正交安装，取为挠性附件在模态坐标系下的角速度。定义状态变量为忽略空间干扰力矩的影响，此时挠性卫星的动力学及运动学如下：w_η为挠性附件在模态坐标系下的位移；

$\left[\begin{array}{ccccccccc}{I}_x& 0& 0& {\mathrm{\σ}}_1& 0& 0& 0& 0& 0\\ 0& I_y& 0& {\mathrm{\σ}}_2& 0& 0& 0& 0& 0\\ 0& 0& I_z& {\mathrm{\σ}}_3& 0& 0& 0& 0& 0\\ {\mathrm{\σ}}_1& {\mathrm{\σ}}_2& {\mathrm{\σ}}_3& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\stackrel{\·}{w}}_x\\ {\stackrel{\·}{w}}_y\\ {\stackrel{\·}{w}}_z\\ {\stackrel{\·}{w}}_{\mathrm{\η}}\\ {\stackrel{\·}{q}}_0\\ {\stackrel{\·}{q}}_1\\ {\stackrel{\·}{q}}_2\\ {\stackrel{\·}{q}}_3\\ \stackrel{\·}{\mathrm{\η}}\end{array}\right]=\left[\begin{array}{c}(I_y-I_z)w_y{w}_z+({\mathrm{\σ}}_2{w}_z-{\mathrm{\σ}}_3{w}_y)w_{\mathrm{\η}}\\ (I_z-I_x)w_x{w}_z+({\mathrm{\σ}}_3{w}_x-{\mathrm{\σ}}_1{w}_z)w_{\mathrm{\η}}\\ (I_x-I_y)w_x{w}_y+({\mathrm{\σ}}_1{w}_y-{\mathrm{\σ}}_2{w}_x)w_{\mathrm{\η}}\\ -2{\mathrm{\ξ}}_f{w}_f{w}_{\mathrm{\η}}-w_f^2\mathrm{\η}\\ 0.5(-w_x{q}_1-w_y{q}_2-w_z{q}_2)\\ 0.5(w_x{q}_0+w_z{q}_2-w_y{q}_3)\\ 0.5(w_y{q}_0-w_z{q}_1+w_y{q}_3)\\ 0.5(w_z{q}_0+w_y{q}_1-w_x{q}_2)\\ w_{\mathrm{\η}}\end{array}\right]+\left[\begin{array}{c}{T}_{CMGx}\\ T_{CMGy}\\ T_{CMGz}\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]+\left[\begin{array}{c}{T}_{Fx}\\ T_{Fy}\\ T_{Fz}\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right]$

定义矩阵：

$B={\left[\begin{array}{ccccccccc}{I}_x& 0& 0& {\mathrm{\σ}}_1& 0& 0& 0& 0& 0\\ 0& I_y& 0& {\mathrm{\σ}}_2& 0& 0& 0& 0& 0\\ 0& 0& I_z& {\mathrm{\σ}}_3& 0& 0& 0& 0& 0\\ {\mathrm{\σ}}_1& {\mathrm{\σ}}_2& {\mathrm{\σ}}_3& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]}^{-1},f\left(x\right)=B\left[\begin{array}{c}(I_y-I_z)w_y{w}_z+({\mathrm{\σ}}_2{w}_z-{\mathrm{\σ}}_3{w}_y)w_{\mathrm{\η}}\\ (I_z-I_x)w_x{w}_z+({\mathrm{\σ}}_3{w}_x-{\mathrm{\σ}}_1{w}_z)w_{\mathrm{\η}}\\ (I_x-I_y)w_x{w}_y+({\mathrm{\σ}}_1{w}_y-{\mathrm{\σ}}_2{w}_x)w_{\mathrm{\η}}\\ -2{\mathrm{\ξ}}_f{w}_f{w}_{\mathrm{\η}}-w_f^2\mathrm{\η}\\ 0.5(-w_x{q}_1-w_y{q}_2-w_z{q}_2)\\ 0.5(w_x{q}_0+w_z{q}_2-w_y{q}_3)\\ 0.5(w_y{q}_0-w_z{q}_1+w_y{q}_3)\\ 0.5(w_z{q}_0+w_y{q}_1-w_x{q}_2)\\ w_{\mathrm{\η}}\end{array}\right]$

综上，以金字塔构型CMG群及反作用飞轮为执行机构的挠性卫星姿态动力学及运动学可以归结为如下的非线性方程：

$\stackrel{\·}{x}=f\left(x\right)+{\mathrm{BT}}_{CMG}+{\mathrm{BT}}_F$

非线性函数f(x)描述了系统状态向量之间的相互影响及耦合关系。在本实施方式中，假设系统的状态均是可观测的，定义系统的输出为y_c＝x。进而有：

$\left\{\begin{array}{c}\stackrel{\·}{x}=f\left(x\right)+B{T}_{CMG}+B{T}_F\\ y_c=x\end{array}\right.$

通过对以上连续时间系统进行离散化，可获得卫星姿态动力学及运动学的离散时间方程：

$\left\{\begin{array}{c}x(k+1)=f_d\left(x\right(k\left)\right)+B_d{T}_{CMG}\left(k\right)+B_d{T}_F\left(k\right)\\ y_c\left(k\right)=x\left(k\right)\end{array}\right.$

金字塔构型SGCMG系统安装构型如图2所示，四个SGCMG单元对称分布，各框架轴在卫星本体坐标系中可表示为：

${\stackrel{\→}{g}}_1=\stackrel{\→}{i}sin\mathrm{\β}+\stackrel{\→}{k}cos\mathrm{\β},{\stackrel{\→}{g}}_2=\stackrel{\→}{j}sin\mathrm{\β}+\stackrel{\→}{k}cos\mathrm{\β}$

${\stackrel{\→}{g}}_3=-\stackrel{\→}{i}sin\mathrm{\β}+\stackrel{\→}{k}cos\mathrm{\β},{\stackrel{\→}{g}}_4=-\stackrel{\→}{j}sin\mathrm{\β}+\stackrel{\→}{k}cos\mathrm{\β}$

式中，分别为CMG群的四个框架轴，β为安装倾角，为沿卫星本体坐标系三轴的基矢量。

SGCMG系统的角动量H可表示为：

$H=\stackrel{4}{\Σ}{h}_i=h_1+h_2+h_3+h_4$

式中，h₁,h₂,h₃,h₄为四个CMG的输出角动量。

在卫星本体坐标系中有：

$\begin{array}{c}H=H_x\stackrel{\→}{i}+H_y\stackrel{\→}{j}+H_z\stackrel{\→}{k}=\[\begin{array}{ccc}\stackrel{\→}{i}& \stackrel{\→}{j}& \stackrel{\→}{k}\end{array}\]\left[\begin{array}{c}{H}_x\\ H_y\\ H_z\end{array}\right]\\ =h_1\left({\mathrm{\δ}}_1\right)+h_2\left({\mathrm{\δ}}_2\right)+h_3\left({\mathrm{\δ}}_3\right)+h_4\left({\mathrm{\δ}}_4\right)\\ =\left[\begin{array}{c}-{\mathrm{cos\βsin\δ}}_1\\ {\mathrm{cos\δ}}_1\\ {\mathrm{sin\βsin\δ}}_1\end{array}\right]+\left[\begin{array}{c}-{\mathrm{cos\δ}}_2\\ -{\mathrm{cos\βsin\δ}}_2\\ {\mathrm{sin\βsin\δ}}_2\end{array}\right]+\left[\begin{array}{c}{\mathrm{cos\βsin\δ}}_3\\ -{\mathrm{cos\δ}}_3\\ {\mathrm{sin\βsin\δ}}_3\end{array}\right]+\left[\begin{array}{c}{\mathrm{cos\δ}}_4\\ {\mathrm{cos\βsin\δ}}_4\\ {\mathrm{sin\βsin\δ}}_4\end{array}\right]\end{array}$

式中，δ₁,δ₂,δ₃,δ₄分别为四个框架角，H_x,H_y,H_z为CMG群系统输出角动量在卫星本体三轴方向上分量，四个CMG的转子角动量均相同，不失一般性，假设为1。

对CMG系统输出角动量H微分得到：

dH＝J₁dδ₁+J₂dδ₂+J₃dδ₃+J₄dδ₄＝Jdδ

式中，J为雅克比矩阵，形式为：

$\begin{array}{c}J=\[\begin{array}{cccc}{J}_1& J_2& J_3& J_4\end{array}\]\\ =\left[\begin{array}{cccc}-{\mathrm{cos\βcos\δ}}_1& {\mathrm{sin\δ}}_2& {\mathrm{cos\βcos\δ}}_3& -{\mathrm{sin\δ}}_4\\ -{\mathrm{sin\δ}}_1& -{\mathrm{cos\βcos\δ}}_2& {\mathrm{sin\δ}}_3& {\mathrm{cos\βcos\δ}}_4\\ {\mathrm{sin\βcos\δ}}_1& {\mathrm{sin\βcos\δ}}_2& {\mathrm{sin\βcos\δ}}_3& {\mathrm{sin\βcos\δ}}_4\end{array}\right]\end{array}$

其中，J₁,J₂,J₃,J₄分别为四个CMG的输出力矩矢量。

CMG系统的输出力矩为：

$T_{CMG}=\stackrel{\·}{H}\left(\mathrm{\δ}\right)=J\stackrel{\·}{\mathrm{\δ}}$

式中，为SGCMG系统的框架运动角速度。

二、考虑挠性卫星姿态机动控制的一般性描述，即：

$J=\mathrm{\Φ}\left(x\right(t_f),t_f)+\underset{t_0}{\overset{t_f}{\∫}}g\left(x\right(t),u(t),t)dt$

满足非线性模型、初始、终端及过程约束条件：

$\stackrel{\·}{x}=f\left(x\right)+{\mathrm{BT}}_{CMG}$

φ(x(t₀),t₀,x(t_f),t_f)＝0

C(x(t),u(t),t；t₀,t_f)≤0

其中：函数g(·)为综合了卫星姿态机动快速性、挠性附件振动抑制性、CMG群奇异性等的加权优化目标函数；φ(·)为初始和终端时刻的状态约束函数；C(·)为卫星姿态机动过程中的约束函数。

伪谱法的主要思想是采用拉格朗日全局插值多项式在一系列离散的时间节点上近似状态变量和控制变量，并通过插值多项式的导数来逼近动力学方程中状态变量对时间的导数，且在选取的配点上严格满足动力学方程的右函数约束。由以下步骤实现：

(1)状态量和控制量的逼近

Legendre伪谱法采用全局插值多项式在区间τ∈[-1,1]内逼近所需的状态变量，以Lagrange插值多项式作为基函数，利用时间区间τ∈[-1,1]的M个点τ₁,τ₂,…τ_M定义这些多项式。因此，最优控制问题中的状态量，控制量以及协态量均可以被近似为：

$y\left(\mathrm{\τ}\right)\≈Y\left(\mathrm{\τ}\right)=\underset{i=1}{\overset{M}{\Σ}}{L}_i\left(\mathrm{\τ}\right)Y\left({\mathrm{\τ}}_i\right)$

$L_i\left(\mathrm{\τ}\right)=\underset{j=0,j\≠i}{\overset{M}{\Π}}\frac{\mathrm{\τ}-{\mathrm{\τ}}_j}{{\mathrm{\τ}}_i-{\mathrm{\τ}}_j}=\frac{g\left(\mathrm{\τ}\right)}{(\mathrm{\τ}-{\mathrm{\τ}}_i)\stackrel{\·}{g}\left(\mathrm{\τ}\right)}$

式中：

L_i(τ)，i＝1,2,…M——Largrange插值多项式；

g(τ)——时间区间τ∈[-1,1]的M个点位置的基函数；

——g(τ)的时间导数；

Y(τ)——真实值y(τ)的(M-1)阶近似，由于其在插值点上等于真实值，因此有：

$L_i\left({\mathrm{\τ}}_j\right)=\left\{\begin{array}{c}1,i=j\\ 0,i\≠j\end{array}\right.$

(2)全局插值点的选取

在上述全局插值逼近状态变量的过程中，需要选取区间τ∈[-1,1]内的M个点τ₁,τ₂,…τ_M以进行离散化。采用非等间距的LGL(Legendre-Gaus-Lobatto)点作为离散化点。在这种方法中，积分节点包含了两个边界点，减少了两个自由度，可以精确逼近小于等于2M-3次的多项式。M个LGL点可以由的零点来确定，其中为(M-1)阶Legendre多项式的导数。其中，P_M(t)定义为：

$P_M\left(\mathrm{\τ}\right)=\frac{1}{2^{M}M!}\frac{d^M}{{\mathrm{d\τ}}^M}\[{({\mathrm{\τ}}^2-1)}^M\]$

(3)时间变量的映射

由于τ₁,τ₂,…τ_M在区间τ∈[-1,1]内，但通常任务规划时间[t₀,t_f]并不只限于该区间，所以需将时间量映射到该区间，映射关系如下：

$t=\frac{t_f-t_0}{2}\mathrm{\τ}+\frac{t_f+t_0}{2}$

式中：

t₀,t_f——表示起始和终端时间。

通过该变换，可以将卫星姿态机动时间区间[t₀,t_f]映射到区间[-1,1]，进而可以在[-1,1]内考虑优化问题的求解。

(4)数值积分的逼近

在选择用于逼近状态、控制和协态量的插值点之后，Legendre伪光谱算法还需要使用另一组插值点精确的逼近最优问题的动力学部分、最优性能指标以及各项包含积分微分的非线性约束。积分逼近的一般形式为：

$\underset{a}{\overset{b}{\∫}}f\left(\mathrm{\τ}\right)d\mathrm{\τ}\≈\underset{l=1}{\overset{K}{\Σ}}{w}_{l}f\left({\mathrm{\τ}}_l\right)$

式中：

τ₁,...τ_K——时间区间τ∈[-1,1]的积分离散点；

w_l——积分权重，

(5)状态微分的逼近

伪光谱算法对最优约束条件离散化时，需将动力学方程等微分条件等转换为代数条件。伪光谱算法应用正交分配方法，在正交多项式(例如Legendre多项式)的零点将状态逼近值的导数进行分配，也即令离散后系统在这些点上的导数逼近实际的导数。

Legendre伪光谱采用LGL点，插值逼近所用的离散点和正交分配所用的点是相同的。其微分方程离散化的表达式为：

$\stackrel{\·}{x}\left({\mathrm{\τ}}_l\right)\≈\stackrel{\·}{X}\left({\mathrm{\τ}}_l\right)=\underset{s=1}{\overset{K}{\Σ}}{\stackrel{\·}{L}}_s\left({\mathrm{\τ}}_l\right)X\left({\mathrm{\τ}}_s\right)=\underset{s=1}{\overset{K}{\Σ}}{D}_{ls}X\left({\mathrm{\τ}}_s\right),(l=1,...K)$

微分矩阵D∈R^K×K定义为：

$D_{ls}=\frac{{\mathrm{dL}}_s\left({\mathrm{\τ}}_l\right)}{d\mathrm{\τ}}=\left\{\begin{array}{c}\frac{P_{N-1}\left({\mathrm{\τ}}_l\right)}{P_{N-1}\left({\mathrm{\τ}}_s\right)}\frac{1}{{\mathrm{\τ}}_l-{\mathrm{\τ}}_i},l\≠s\\ -\frac{(N-1)N}{4},l=s=1\\ \frac{(N-1)N}{4},l=s=N\\ 0,others\end{array}\right.$

这里可以得到状态近似值导数的代数表达式，通过该表达式就可以用一组在配置点上的代数条件在离散最优化的求解中来代替实际系统中的各个微分方程。连续动力学方程通过正交配置可以转换成如下K个代数方程：

$R_l=\underset{s=1}{\overset{K}{\Σ}}{D}_{ls}X\left({\mathrm{\τ}}_s\right)-\frac{t_f-t_0}{2}f\left(X\right({\mathrm{\τ}}_l),u({\mathrm{\τ}}_l),{\mathrm{\τ}}_l;t_0,t_f)=0,(l=1,...K)$

初末状态的积分约束可直接由边界约束条件中近似多项式的边界点X₁,X_N决定，即：

φ(X₁,τ₀,X_N,τ_f)＝0

(6)终端状态约束的逼近

一般情况下，终端状态x(t_f)也是有约束的，需要对其离散化，这里可以根据系统的动力学，利用数值积分逼近得到终端状态约束的离散化：

$X\left({\mathrm{\τ}}_f\right)=X\left({\mathrm{\τ}}_0\right)+\frac{t_f-t_0}{2}\underset{l=1}{\overset{K}{\Σ}}{w}_{l}f\left(X\right({\mathrm{\τ}}_l),U({\mathrm{\τ}}_l),{\mathrm{\τ}}_l;t_0,t_f)$

通过以上设计，可将原来的连续非线性系统最优问题转化为离散非线性系统动态规划问题(NLP)，进而可以利用序列二次规划算法进行求解。

三、以规划得到的最优卫星姿态机动轨迹为期望跟踪目标，建立滚动优化跟踪控制问题，通过对其求解实现反作用飞轮控制力矩的计算。分解为以下步骤：

(1)、卫星姿态机动优化指标的建立；

对于规划得到的k时刻的期望姿态四元数及期望角速度，表示为：

$q^h\left(k\right)={\left[\begin{array}{cccc}{q}_0^h\left(k\right)& q_1^h\left(k\right)& q_2^h\left(k\right)& q_3^h\left(k\right)\end{array}\right]}^T$

$w^h\left(k\right)={\left[\begin{array}{ccc}{w}_x^h\left(k\right)& w_y^h\left(k\right)& w_z^h\left(k\right)\end{array}\right]}^T$

其中，q^h(k)表示k时刻的期望姿态四元数，w^h(k)表示k时刻的期望角速度；

反作用飞轮参与控制的目的是实现对规划轨迹的高精度跟踪控制，即对卫星姿态及角速度等最优规划的跟踪。在本发明中，将卫星姿态及角速度等的跟踪误差、飞轮控制力矩和终端偏差加权组合作为优化指标，定义如下：

$J=\underset{n=1}{\overset{N_p-1}{\Σ}}e{(k+n)}^{T}Qe(k+n)+\underset{m=0}{\overset{N_p-1}{\Σ}}u{(k+m)}^{T}Ru(k+m)+e{\left(N_p\right)}^{T}Pe\left(N_p\right)$

式中：e(k+n)表示对最优轨迹等的跟踪误差，u(k+m)为反作用飞轮的待设计控制力矩，矩阵Q,R,P为相应的跟踪误差、控制量和终端加权矩阵，N_p为预测时域，N_u为控制时域；

(2)、卫星姿态的预测；

根据k时刻的卫星姿态信息x(k)，基于规划获得的最优CMG群框架轴角速度轨迹，采用离散化的卫星姿态动力学及运动学方程，通过迭代计算，求取未来N_p预测步内的卫星姿态信息：

$\begin{array}{c}\stackrel{\‾}{y}(k+1)=f_d\left(x\left(k\right)\right)+B_d{T}_{CMG}\left(k\right)+B_d{T}_F\left(k\right)\\ \stackrel{\‾}{y}\left(k+2\right)=f_d\left(\stackrel{\‾}{x}\left(k+1\right)\right)+B_d{T}_{CMG}\left(k+1\right)+B_d{T}_F\left(k+1\right)\\ \mathrm{...}\\ \stackrel{\‾}{y}\left(k+N_p\right)=f_d\left(\stackrel{\‾}{x}\left(k+N_p-1\right)\right)+B_d{T}_{CMG}\left(k+N_p-1\right)+B_d{T}_F\left(k+N_p-1\right)\end{array}$

其中，预测时域内的规划CMG群框架轴角速度控制量为将相应时刻规划值折算成CMG群的控制力矩；

从以上N_p步预测时域内的系统状态可知，预测后的待设计自由变量为N_p时域步内的控制量{T_F(k),T_F(k+1),…T_F(k+N_p-1)}，即反作用飞轮的待设计控制力矩u(k+m)；

当预测时域N_p超出控制时域N_u时，设定控制力矩输入在区间[N_u,N_p]保持不变，即：

T_F(k+N_u-1)＝T_F(k+N_u)＝…＝T_F(k+N_p-1)

此时，卫星姿态的预测输出方程为：

$\begin{array}{c}\stackrel{\‾}{y}(k+1)=f_d\left(x\left(k\right)\right)+B_d{T}_{CMG}\left(k\right)+B_d{T}_F\left(k\right)\\ \mathrm{...}\\ \stackrel{\‾}{y}\left(k+N_u\right)=f_d\left(\stackrel{\‾}{x}\left(k+N_u-1\right)\right)+B_d{T}_{CMG}\left(k+N_u-1\right)+B_d{T}_F\left(k+N_u-1\right)\\ \mathrm{...}\\ \stackrel{\‾}{y}\left(k+N_p\right)=f_d\left(\stackrel{\‾}{x}\left(k+N_p-1\right)\right)+B_d{T}_{CMG}\left(k+N_p-1\right)+B_d{T}_F\left(k+N_p-1\right)\end{array};$

(3)、卫星姿态机动律的实现；

挠性敏捷卫星姿态机动控制的优化问题为：

且满足由执行机构的能力带来的时域约束条件：

$T_{F\min }\≤\stackrel{\‾}{T_F}(k+m)\≤T_{Fmax},0\≤m\≤N_u-1$

${\mathrm{\ΔT}}_{Fmin}\≤\mathrm{\Δ}\stackrel{\‾}{T_F}(k+m)\≤{\mathrm{\ΔT}}_{Fmax}$

$\mathrm{\Δ}\stackrel{\‾}{T_F}(k+m)=\stackrel{\‾}{T_F}(k+m)-\stackrel{\‾}{T_F}(k+m-1)$

式中，T_Fmin表示飞轮力矩的下界，T_Fmax表示飞轮力矩的上界，ΔT_Fmin表示飞轮力矩增量的下界，ΔT_Fmax表示飞轮力矩增量的上界；

目标函数J(x(k),T_k)整理为：

$J\left(x\right(k),T_k)=\underset{n=1}{\overset{N_p-1}{\Σ}}\left|\right|{\stackrel{\‾}{y}}_c(k+n)-r(k+n)|{|}_Q^2+\underset{m=0}{\overset{N_p-1}{\Σ}}|\left|\stackrel{\‾}{T_F}(k+m)\right|{|}_R^2+|\[{\stackrel{\‾}{y}}_c(k+N_p)-r(k+N_p)|{|}_P^2$

式中：||·||²表示2-范数，表示预测的卫星姿态输出，r(k+n)为期望的卫星姿态输出，为系统的预测控制输入，即待优化的飞轮控制力矩，定义为：

$\stackrel{\‾}{T_F}(k+m)={\stackrel{\‾}{T_F}}_m$

式中：是互相独立的优化变量，记为：

$T_k={\left[\begin{array}{cccc}{\stackrel{\‾}{T_F}}_0& {\stackrel{\‾}{T_F}}_1& \mathrm{...}& {\stackrel{\‾}{T_F}}_{(N_u-1)}\end{array}\right]}^T$

设定卫星姿态控制的约束优化问题的最优解为：

$T_k^{*}={\left[\begin{array}{cccc}{\stackrel{\‾}{T_F}}_0^{*}& {\stackrel{\‾}{T_F}}_1^{*}& \mathrm{...}& {\stackrel{\‾}{T_F}}_{(N_u-1)}^{*}\end{array}\right]}^T$

根据预测控制原理，最优解的第一个元素，即当前时刻的飞轮最优控制力矩定义为：

四、在当前控制时刻，将规划得到的最优CMG群框架角速度及计算的飞轮最优控制力矩作用于卫星，驱动卫星姿态运动；在下一采样时刻，根据卫星姿态更新的状态信息，重复以上优化过程，实现挠性卫星姿态机动滚动时域控制，结合图3。

在本实施中以某型小卫星为例，其转动惯量矩阵如下所示：

$I=\left(\begin{array}{ccc}104& -2& -0.3\\ -2& 105& 1.8\\ -0.3& 1.8& 147\end{array}\right)(Kg\·m2)$

这里，假设卫星的初始姿态四元数为[1；0；0；0]，目标姿态四元数为[0.9274；0.2226；0.0702；0.2924]，初始姿态角速度均为[0°/s,0°/s,0°/s]，目标姿态角速度为[-0.8°/s,0°/s,0°/s]。考虑的挠性帆板一阶频率为2.23Hz，阻尼为0.032，刚柔耦合系数矩阵为[0.00041,3.833,0]。SGCMG最大框架角速度为3rad/s，转子额定角动量为5Nms。考虑的金字塔构型CMG群的控制力矩约束为[-5Nm，5Nm],飞轮力矩约束为[-0.04Nm,0.04Nm]。机动时间设置为60s。空间扰动采用典型表达形式：

$\left\{\begin{array}{c}{M}_{dx}={10}^{-5}(3{\mathrm{cos\ω}}_{0}t+1)\\ M_{dy}={10}^{-5}(1.5{\mathrm{sin\ω}}_{0}t+3{\mathrm{cos\ω}}_{0}t)\\ M_{dz}={10}^{-5}(3{\mathrm{sin\ω}}_{0}t+1)\end{array}\right.$

为说明本发明方法的有效性，仅考虑卫星三轴的主惯量矩阵进行设计，且在控制方法设计时不考虑空间干扰力矩的影响。图4、图5和图6给出了卫星姿态四元数、姿态角速度的变化曲线及合控制力矩的变化曲线。从仿真结果可以看出，在初始及终端时刻，卫星姿态均满足设计值要求，且在机动过程中的姿态数值也在状态约束要求范围内；CMG群与反作用飞轮的合控制力矩均满足设计要求。

Claims (2)

1.挠性敏捷卫星的联合执行机构控制方法，其特征是，该方法由以下步骤实现：

步骤一、选择惯性坐标系为参考坐标系，建立以金字塔构型CMG群和反作用飞轮为执行机构的挠性卫星姿态动力学及运动学方程；以挠性附件的模态坐标系为基准，建立挠性附件振动的动力学方程，通过定义新的状态变量，基于建立的三个方程获得用于CMG群框架角速度最优轨迹规划及预测卫星姿态未来信息的非线性连续状态空间方程，通过离散化方法，获得离散非线性状态空间方程；

步骤二、对步骤一建立的非线性连续状态空间方程，建立兼顾卫星姿态机动快速性、挠性附件振动抑制性能、金字塔构型CMG群奇异性及满足各种约束等的优化控制问题；基于Legendre伪谱法，实现对优化控制问题的离散化求解，获得卫星姿态机动的最优状态轨迹及CMG群的最优框架角速度轨迹；

步骤三、将步骤二获得的CMG群的最优框架角速度轨迹带入步骤一建立的离散非线性状态空间方程，并根据当前测量的卫星姿态信息，建立卫星姿态的预测输出方程，实现对预测时域内的卫星姿态信息的预测；

步骤四、根据步骤三建立的预测输出方程，以步骤二获得的卫星姿态机动的最优状态轨迹为跟踪目标，建立兼顾卫星姿态跟踪快速性及鲁棒性的优化控制问题，经过对该优化控制问题的求解，获得反作用飞轮的非线性模型预测机动控制力矩；将规划得到的CMG群的最优框架角速度及飞轮机动控制力矩作用于卫星，驱动卫星姿态运动；

步骤五、采用离散控制技术，在每个采样时刻重复步骤三和步骤四，通过逐步更新卫星姿态信息，实现挠性卫星姿态机动的复合滚动控制。

2.根据权利要求1所述的挠性敏捷卫星的联合执行机构控制方法，其特征在于，步骤四中，获得反作用飞轮的非线性模型预测机动控制力矩的具体过程为：

步骤四一、卫星姿态机动优化指标的建立；

对于规划得到的k时刻的期望姿态四元数及期望角速度，表示为：

$q^h\left(k\right)={\left[\begin{array}{cccc}{q}_0^h\left(k\right)& q_1^h\left(k\right)& q_2^h\left(k\right)& q_3^h\left(k\right)\end{array}\right]}^T$

$w^h\left(k\right)={\left[\begin{array}{ccc}{w}_x^h\left(k\right)& w_y^h\left(k\right)& w_z^h\left(k\right)\end{array}\right]}^T$

其中，q^h(k)表示k时刻的期望姿态四元数，w^h(k)表示k时刻的期望角速度；

反作用飞轮参与控制的目的是实现对规划轨迹的高精度跟踪控制，即对卫星姿态及角速度等最优规划的跟踪。在本发明中，将卫星姿态及角速度等的跟踪误差、飞轮控制力矩和终端偏差加权组合作为优化指标，定义如下：

$J=\underset{n=1}{\overset{N_p-1}{\Σ}}e{(k+n)}^{T}Qe(k+n)+\underset{m=0}{\overset{N_p-1}{\Σ}}u{(k+m)}^{T}Ru(k+m)+e{\left(N_p\right)}^{T}Pe\left(N_p\right)$

式中：e(k+n)表示对最优轨迹等的跟踪误差，u(k+m)为反作用飞轮的待设计控制力矩，矩阵Q,R,P为相应的跟踪误差、控制量和终端加权矩阵，N_p为预测时域，N_u为控制时域；

步骤四二、卫星姿态的预测；

根据k时刻的卫星姿态信息x(k)，基于规划获得的最优CMG群框架轴角速度轨迹，采用离散化的卫星姿态动力学及运动学方程，通过迭代计算，求取未来N_p预测步内的卫星姿态信息：

$\begin{array}{c}\stackrel{\‾}{y}(k+1)=f_d\left(x\right(k\left)\right)+B_d{T}_{CMG}\left(k\right)+B_d{T}_F\left(k\right)\\ \stackrel{\‾}{y}(k+2)=f_d\left(\stackrel{\‾}{x}\right(k+1\left)\right)+B_d{T}_{CMG}(k+1)+B_d{T}_F(k+1)\\ \mathrm{...}\\ \stackrel{\‾}{y}(k+N_p)=f_d\left(\stackrel{\‾}{x}\right(k+N_p-1\left)\right)+B_d{T}_{CMG}(k+N_p-1)+B_d{T}_F(k+N_p-1)\end{array}$

其中，预测时域内的规划CMG群框架轴角速度控制量为将相应时刻规划值折算成CMG群的控制力矩；

从以上N_p步预测时域内的系统状态可知，预测后的待设计自由变量为N_p时域步内的控制量{T_F(k),T_F(k+1),…T_F(k+N_p-1)}，即反作用飞轮的待设计控制力矩u(k+m)；

当预测时域N_p超出控制时域N_u时，设定控制力矩输入在区间[N_u,N_p]保持不变，即：

T_F(k+N_u-1)＝T_F(k+N_u)＝…＝T_F(k+N_p-1)

此时，卫星姿态的预测输出方程为：

$\begin{array}{c}\stackrel{\‾}{y}(k+1)=f_d\left(x\right(k\left)\right)+B_d{T}_{CMG}\left(k\right)+B_d{T}_F\left(k\right)\\ \mathrm{...}\\ \stackrel{\‾}{y}(k+N_u)=f_d\left(\stackrel{\‾}{x}\right(k+N_u-1\left)\right)+B_d{T}_{CMG}(k+N_u-1)+B_d{T}_F(k+N_u-1)\\ \mathrm{...}\\ \stackrel{\‾}{y}(k+N_p)=f_d\left(\stackrel{\‾}{x}\right(k+N_p-1\left)\right)+B_d{T}_{CMG}(k+N_p-1)+B_d{T}_F(k+N_p-1)\end{array};$

步骤四三、卫星姿态机动律的实现；

挠性敏捷卫星姿态机动控制的优化问题为：

$\underset{U_k}{min}J\left(x\right(k),T_k)$

且满足由执行机构的能力带来的时域约束条件：

$T_{F\min }\≤\stackrel{\‾}{T_F}(k+m)\≤T_{Fmax},0\≤m\≤N_u-1$

${\mathrm{\ΔT}}_{Fmin}\≤\mathrm{\Δ}\stackrel{\‾}{T_F}(k+m)\≤{\mathrm{\ΔT}}_{Fmax}$

$\mathrm{\Δ}\stackrel{\‾}{T_F}(k+m)=\stackrel{\‾}{T_F}(k+m)-\stackrel{\‾}{T_F}(k+m-1)$

式中，T_Fmin表示飞轮力矩的下界，T_Fmax表示飞轮力矩的上界，ΔT_Fmin表示飞轮力矩增量的下界，ΔT_Fmax表示飞轮力矩增量的上界；

目标函数J(x(k),T_k)整理为：

$J\left(x\right(k),T_k)=\underset{n=1}{\overset{N_p-1}{\Σ}}\left|\right|{\stackrel{\‾}{y}}_c(k+n)-r(k+n)|{|}_Q^2+\underset{m=0}{\overset{N_u-1}{\Σ}}|\left|\stackrel{\‾}{T_F}(k+m)\right|{|}_R^2+\left|\right|{\stackrel{\‾}{y}}_c(k+N_p)-r(k+N_p)|{|}_P^2$

式中：||·||²表示2-范数，表示预测的卫星姿态输出，r(k+n)为期望的卫星姿态输出，为系统的预测控制输入，即待优化的飞轮控制力矩，定义为：

$\stackrel{\‾}{T_F}(k+m)={\stackrel{\‾}{T_F}}_m$

式中：是互相独立的优化变量，记为：

$T_k={\left[\begin{array}{cccc}{\stackrel{\‾}{T_F}}_0& {\stackrel{\‾}{T_F}}_1& \mathrm{...}& {\stackrel{\‾}{T_F}}_{(N_u-1)}\end{array}\right]}^T$

设定卫星姿态控制的约束优化问题的最优解为：

$T_k^{*}={\left[\begin{array}{cccc}{\stackrel{\‾}{T_F}}_0^{*}& {\stackrel{\‾}{T_F}}_1^{*}& \mathrm{...}& {\stackrel{\‾}{T_F}}_{(N_u-1)}^{*}\end{array}\right]}^T$

根据预测控制原理，最优解的第一个元素，即当前时刻的飞轮最优控制力矩定义为：