本发明公开了一种考虑反作用飞轮特性及功率受限的航天器姿态控制方法，包括以下步骤：建立四元数描述方式下的航天器运动学和动力学模型；建立航天器体装式太阳能帆板输出功率与太阳入射角的关系；考虑反作用飞轮特性，建立角动量约束、控制力矩约束以及体装式太阳能帆板故障时反作用飞轮的实际可用功率约束；根据任务需求设计目标性能函数；基于模型预测控制策略完成上述约束条件下的航天器姿态控制任务。本方法能够很好的解决航天器在考虑反作用飞轮特性以及体装式太阳能帆板发生故障时的姿态控制问题，通过对目标函数的设计达到控制精度和能量消耗的综合最优，保证航天器在反作用飞轮特性及功率受限的情况下仍能完成高精度姿态控制任务。
1.一种考虑反作用飞轮特性及功率受限的航天器姿态控制方法，其特征在于，包括以下步骤：S1：针对航天器姿态控制问题，建立四元数描述方式下的航天器姿态运动学和动力学模型，将所述航天器姿态运动学与动力学模型离散化，建立模型预测控制方法的预测模型；
S2：建立航天器体装式太阳能帆板输出功率与太阳入射角的关系，并将其转化为与航天器姿态角有关的表达式；
S3：基于航天器执行机构反作用飞轮特性，建立反作用飞轮的角动量约束和控制力矩约束，同时建立航天器体装式太阳能帆板故障时反作用飞轮的实际可用功率约束；
S5：基于步骤S1中建立的模型预测控制方法的预测模型，完成步骤S3中建立的约束下的航天器姿态控制任务；
步骤S1中的四元数描述方式下的航天器姿态运动学和动力学模型如下：其中，J∈R3×3为航天器惯量矩阵，R3×3表示3×3矩阵集合；ω＝(ω1,ω2,ω3)T为航天器
本体三轴角速度；H∈R为航天器总角动量，R 表示3维实向量空间，其具体形式为H＝Jω+HRWA，HRWA∈R3为反作用飞轮组合角动量在航天器惯量主轴上的投影；uc∈R3为作用于航天器的姿态控制力矩；q＝(q0,q1,q2,q3)T为描述航天器姿态的四元数，且满足Ξ(ω)是斜反对称矩阵，其具体形式如下：对于任意一个向量x∈R3，其对应的斜反对称阵x×形式如下：在忽略航天器轨道角速度对反作用飞轮角动量的影响时，反作用飞轮组合的动力学模型如下：其中，τRWA∈RL为由L个反作用飞轮产生的控制力矩，RL表示L维实向量空间；hRWA∈RL为反作用飞轮组合的角动量， 为反作用飞轮组合角动量关于时间的导数；JRWA∈RL×L为反作用飞轮组合的转动惯量矩阵，其形式为为各个反作用飞轮的转动惯
量，RL×L表示L×L矩阵集合； 为各个反作用飞轮的转速， 为各个反作用飞轮转速的关于时间的导数；
基于反作用飞轮组合的动力学模型及反作用飞轮组合构型，得到反作用飞轮组合作用在航天器上的姿态控制力矩为：其中，C∈R3×L为反作用飞轮组合安装矩阵，R3×L表示3×L矩阵集合；
根据角动量守恒定理，作用于航天器的姿态控制力矩uc为将式(5)代入式(1)，得到由反作用飞轮组合作为执行机构的航天器姿态运动学与动力学模型：将所述航天器姿态运动学与动力学模型离散化，取采样时间间隔为Δt，以将所述航天器姿态运动学与动力学模型作为模型预测控制方法的预测模型，离散化后的航天器姿态运动学与动力学模型为：其中，下标k和k+1分别表示k时刻和k+1时刻；Hk为k时刻的航天器总角动量；
令xk＝(qk,ωk)T为第k时刻航天器姿态控制系统状态量，表征航天器的姿态和角速度，为第k时刻系统控制量，将离散化后的航天器姿态运动学与动力学模型(7)简化为如下形式：
步骤S2中建立的航天器体装式太阳能帆板输出功率与太阳入射角的关系为：其中，Ptotal为航天器上m个体装式太阳能帆板的输出总功率； 为第i个太阳能帆板的输出功率； 为第i个太阳能帆板的有效面积；θi为太阳入射角，表示第i个太阳能帆板法线与太阳光的夹角；PEOL(θi)为太阳入射角θi的函数，表示第i个太阳能帆板的单位面积功率，其具体关系式如下：其中，Ps为入射太阳辐射；η为太阳能帆板光电转换率；Id为太阳能帆板固有退化因子；
太阳入射角随航天器运行轨道、航天器当前时刻的姿态角以及所处的轨道位置变化而变化，为了获得航天器在轨道上任意位置的瞬时可用总功率，建立太阳入射角模型如下：θi＝arccos,i＝1,2,...,m (13)其中，ni∈R3,|ni|＝1为第i个太阳能帆板在航天器本体坐标系中的单位法向量；
<·,·>表示内积运算；rsb∈R3为太阳单位矢量在航天器本体坐标系中的投影，其有如下转换关系：rsb＝RbIrsI (14)其中，rsI＝rs/|rs|∈R3为太阳在地心惯性坐标系下的位置矢量rs∈R3的单位向量；RbI为地心惯性坐标系到航天器本体坐标系的转换矩阵，与欧拉转动顺序以及姿态欧拉角有关；
将航天器体装式太阳能帆板输出功率与太阳入射角的关系转化为与航天器姿态角相关的表达式：其中，Psolar∈Rm为体装式太阳能帆板输出功率向量，Rm表示m维实向量空间。
2.根据权利要求1所述的方法，其特征在于，步骤S3具体过程如下：航天器执行机构反作用飞轮的角动量约束以及控制力矩约束的具体形式如下：
由于航天器姿态控制系统为离散系统形式，故控制力矩约束用如下角加速度 的离散形式来近似表示：
3)体装式太阳能帆板故障时反作用飞轮的实际可用功率约束：在假设航天器上的体装式太阳能帆板总是能受太阳光照射的情况下，体装式太阳能帆板故障时可供反作用飞轮的实际功率约束的表达式如下：其中， 为k时刻各体装式太阳能帆板能提供的瞬时功率向量；κ为可供执行机构使用的功率占比；Ek＝(e1,e2,…,em)为k时刻各太阳能帆板的故障因子序列，其中0≤ei≤1,i＝1,2,…,m； 为k时刻反作用飞轮组合为了提供期望姿态控制力矩而所需的总功率，其具体形式如下：其中， 为k时刻反作用飞轮组合角动量关于时间的导数， 为k时刻反作用飞轮组合的角加速度，上标T表示转置，
将所述反作用飞轮的角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束整理成如下形式：其中， 04L+1表示4L+1维零向量；L为反作用飞轮组合中的反作用飞轮总个数。
3.根据权利要求2所述的方法，其特征在于，在航天器姿态控制过程中在满足高控制精度的同时尽可能减少能量消耗情况下，步骤S4中设计的目标性能函数为：其中，N为预测步长，下标j|k表示第k时刻下的第j预测时刻；x0|k为第k时刻航天器姿态控制系统预测步长内的初始状态量，xj|k＝(qj|k,ωj|k)T为第k时刻预测步长内第j步航天器姿态控制系统状态量，表征航天器的姿态和角速度； 为控制系统的期望状态；
为第k时刻预测步长内第j步系统控制量；Q,P为状态量和控制量的常值正定权重矩阵。
4.根据权利要求3所述的方法，其特征在于，步骤S5具体过程如下：基于步骤S1中建立的模型预测控制方法的预测模型、步骤S4中设计的目标性能函数和步骤S3中建立的反作用飞轮的角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束条件，将航天器姿态控制问题转化为如下优化问题：满足：
xj+1|k＝F(xj|k,uj|k) (27)K(xj|k,uj|k)≤04L+1 (28)其中，xj+1|k＝F(xj|k,uj|k)为航天器姿态控制系统状态方程；K(xj|k,uj|k)为包括反作用飞轮的角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束的系统约束；航天器姿态控制系统第k时刻预测步长内第j步的状态方程和系统约束的具体形式如下：求解所述优化问题得到系统在k时刻下满足约束条件的最优控制序列：其中，上标*表示最优解，
将所述最优控制序列中的第一项作为k时刻航天器姿态控制系统的最优控制输入量，即最后，令 求解下一时刻航天器姿态控制系统的最优控制输入量，完成步骤S3中建立的约束下的航天器姿态控制任务。
考虑反作用飞轮特性及功率受限的航天器姿态控制方法
技术领域
[0001] 本发明属于航天器姿态控制技术领域，具体地涉及一种考虑反作用飞轮特性及功率受限的航天器姿态控制方法，主要应用于考虑太阳入射角变化、飞轮角动量约束、控制力矩约束以及体装式太阳能帆板故障时可用功率受限的航天器姿态控制问题。
背景技术
[0002] 近年来，随着空间技术的飞速发展以及空间飞行任务越发趋于多样化和复杂化，航天器在轨姿态控制是实现各类包括姿态调整操作的航天器工程任务的重要保证，是一项关键技术。反作用飞轮由于其具有控制力矩精细化且不消耗星上工质的优点，已逐渐成为目前航天器的主要执行机构之一，但却不可忽略其在工作过程中容易出现的角动量饱和以及控制力矩饱和等问题对航天器姿态控制性能产生的影响。其次，对于在轨执行空间任务的航天器而言，其完成科学探测及姿态控制等任务所需的能源均来源于航天器上的太阳能帆板，但复杂的太空环境可能导致太阳能帆板发生故障。同时，考虑到太阳能帆板产生的能源需要供有效负载、姿态控制系统及其他分系统共同使用，航天器所能提供给姿态控制系统使用的功率是有限且变化的，主要与航天器在轨道上所处的空间位置、帆板的太阳入射角以及航天器本体姿态角有关。因此，研究航天器在反作用飞轮特性及功率受限等约束下完成姿态高精度控制的同时满足功率消耗要求并使姿态控制系统消耗的能量尽可能小，从而确保航天器完成科学任务，这将具有一定的工程意义。
[0003] 针对目前的考虑反作用飞轮特性的航天器姿态控制方法，专利 201610196190.8提出一种针对飞轮饱和及摩擦特性的挠性航天器姿态控制方法，但该方法只考虑了飞轮的控制力矩饱和特性，没有考虑飞轮功率受限问题；同样地，专利201810298997.1提出一种考虑飞轮动态特性的航天器姿态规避的模型预测控制方法，但该方法也只考虑了飞轮的角动量及控制力矩饱和特性，并没有考虑到飞轮可用功率受限的情况，在实际情况中飞轮的实际功率受限将影响航天器的姿态控制性能。
[0004] 针对目前的考虑反作用飞轮功率的航天器姿态控制问题研究成果，大部分现有的控制策略都是致力于实现反作用飞轮组合在提供姿态控制力矩时的瞬时总功率最优，但忽视了最优瞬时总功率仍大于最小可用功率的情况，在这样的情况下，反作用飞轮组合仍将无法提供所需的控制力矩，进而影响航天器的姿态控制性能。
发明内容
[0005] 本发明的技术解决问题是：在航天器进行姿态机动的过程中，反作用飞轮的角动量饱和、控制力矩饱和不可忽视，同时航天器体装式太阳能帆板也可能发生故障，从而帆板提供的可用功率是有限的。由于该有限功率需要供有效负载、姿态控制系统及其他分系统共同使用，因此可供姿态控制系统使用的功率仅为故障帆板输出功率中的一小部分。为了确保航天器正常运行，使姿态控制系统在如此小且变化的可用功率要求下，能够完成高精度且低能耗的姿态控制任务。本发明基于模型预测控制策略，提出一种考虑反作用飞轮特性及功率受限的航天器姿态控制方法。该方法能够同时处理多种实际工程约束且其鲁棒性较好，能使航天器在给定功率要求下完成姿态控制任务，并能够根据任务需求设计目标性能函数使得航天器在姿态控制过程中实现控制精度和能量消耗的综合优化。
[0006] 本发明提供了一种考虑反作用飞轮特性及功率受限的航天器姿态控制方法，包括如下步骤：
[0007] S1：针对航天器姿态控制问题，建立四元数描述方式下的航天器姿态运动学和动力学模型，将所述航天器姿态运动学与动力学模型离散化，建立模型预测控制方法的预测模型；
[0008] S2：建立航天器体装式太阳能帆板输出功率与太阳入射角的关系，并将其转化为与航天器姿态角有关的表达式；
[0009] S3：基于航天器执行机构反作用飞轮特性，建立反作用飞轮的角动量约束和控制力矩约束，同时建立航天器体装式太阳能帆板故障时反作用飞轮的实际可用功率约束；
[0010] S4：根据任务需求设计目标性能函数；
[0011] S5：基于步骤S1中建立的模型预测控制方法的预测模型，完成步骤S3中建立的约束下的航天器姿态控制任务。
[0012] 进一步，所述步骤S1中的四元数描述方式下的航天器姿态运动学和动力学模型如下：
[0013]
[0014] 其中，J∈R3×3为航天器惯量矩阵；ω＝(ω1，ω2，ω3)T为航天器本体三轴角速度；H∈R3为航天器总角动量，其具体形式为H＝Jω+HRWA，HRWA∈R3为反作用飞轮组合角动量在航3 T
天器惯量主轴上的投影；uc∈R为作用于航天器的姿态控制力矩；q＝(q0，q1，q2，q3)为描述航天器姿态的四元数，且满足 Ξ(ω)是斜反对称矩阵，其具体形式如
下：
[0015]
[0016] 对于任意一个向量x∈R3，其对应的斜反对称阵x×形式如下：
[0017]
[0018] 在忽略航天器轨道角速度对反作用飞轮角动量的影响时，反作用飞轮组合的动力学模型如下：
[0019]
[0020] 其中，τRWA∈RL为由L个反作用飞轮产生的控制力矩；hRWA∈RL为反作用飞轮组合的角动量， 为反作用飞轮组合角动量关于时间的导数； JRWA∈RL×L为反作用飞轮组合的转动惯量矩阵，其形式为 为各个反作用飞轮的转动惯量； 为各个反作用飞轮的转
速， 为各个反作用飞轮转速的关于时间的导数；
[0021] 基于反作用飞轮组合的动力学模型及反作用飞轮组合构型，得到反作用飞轮组合作用在航天器上的姿态控制力矩为：
[0022]
[0023] 其中，C∈R3×L为反作用飞轮组合安装矩阵；
[0024] 根据角动量守恒定理，作用于航天器的姿态控制力矩uc为
[0025]
[0026] 将式(5)代入式(1)，得到由反作用飞轮组合作为执行机构的航天器姿态运动学与动力学模型：
[0027]
[0028] 将所述航天器姿态运动学与动力学模型离散化，取采样时间间隔为Δt，以将所述航天器姿态运动学与动力学模型作为模型预测控制方法的预测模型，离散化后的航天器姿态运动学与动力学模型为：
[0029]
[0030] 其中，下标k和k+1表示k时刻和k+1时刻；Hk为k时刻的航天器总角动量；I4×4为四阶单位矩阵；上标x为斜反对称矩阵；
[0031] 令xk＝(qk，ωk)T为第k时刻航天器姿态控制系统状态量，表征航天器的姿态和角速度， 为第k时刻系统控制量，将离散化后的航天器姿态运动学与动力学模型(7)简化为如下形式：
[0032]
[0033] 定义
[0034]
[0035] 则
[0036] xk+1＝F(xk，uk) (10)
[0037] 进一步，步骤S2中建立的航天器体装式太阳能帆板输出功率与太阳入射角的关系为：
[0038]
[0039] 其中，Ptotal为航天器上m个体装式太阳能帆板的输出总功率； 为第i 个太阳能帆板的输出功率； 为第i个太阳能帆板的有效面积；θi为太阳入射角，表示第i个太阳能帆板法线与太阳光的夹角；PEOL(θi)为太阳入射角θi的函数，表示第i个太阳能帆板的单位面积功率，其具体关系式如下：
[0040]
[0041] 其中，Ps为入射太阳辐射；η为太阳能帆板光电转换率；Id为太阳能帆板固有退化因子；Ld为太阳能帆板寿命退化因子；
[0042] 但考虑到航天器在轨运行时，太阳入射角不是一个固定值，它随航天器运行轨道、航天器当前时刻的姿态角以及所处的轨道位置变化而变化，为了获得航天器在轨道上任意位置的瞬时可用总功率，建立太阳入射角模型如下：
[0043] θi＝arccos，i＝1，2，...，m (13)
[0044] 其中，ni∈R3，|ni|＝1为第i个太阳能帆板在航天器本体坐标系中的单位法向量；<·，·>表示内积运算；rsb∈R3为太阳单位矢量在航天器本体坐标系中的投影，其有如下转换关系：
[0045] rsb＝RbIrsI (14)
[0046] 其中，rsI＝rs/|rs|∈R3为太阳在地心惯性坐标系下的位置矢量rs∈R3的单位向量；RbI为地心惯性坐标系到航天器本体坐标系的转换矩阵，与欧拉转动顺序以及姿态欧拉角有关；
[0047] 将航天器体装式太阳能帆板输出功率与太阳入射角的关系转化为与航天器姿态角相关的表达式：
[0048]
[0049] 其中，Psolar∈Rm为体装式太阳能帆板输出功率向量。
[0050] 进一步，步骤S3具体过程如下：
[0051] 航天器执行机构反作用飞轮的角动量约束以及控制力矩约束的具体形式如下：
[0052] 1)角动量约束用转速约束表示为：
[0053]
[0054] 其中， 为反作用飞轮组合的最大转速向量；
[0055] 将式(16)变形为增量形式：
[0056]
[0057] 2)控制力矩约束表示为：
[0058]
[0059] 其中， 为反作用飞轮组合的最大控制力矩向量；
[0060] 由于航天器姿态控制系统为离散系统形式，故控制力矩约束用如下角加速度的离散形式来近似表示：
[0061]
[0062] 3)体装式太阳能帆板故障时反作用飞轮的实际可用功率约束：
[0063] 首先，假设航天器进入地影区的时间可以忽略不计，即近似认为航天器上的体装式太阳能帆板总是能受太阳光照射；
[0064] 其次，考虑到航天器上体装式太阳能帆板的输出功率需要供有效负载、执行机构及其他分系统共同使用；同时还要考虑当航天器运行到不同的轨道位置时，太阳能帆板能提供的输出功率是随太阳入射角或者说随航天器卫星姿态角变化的，而且太阳能帆板也有可能因复杂的太空环境而发生故障，因此可供执行机构使用的实际功率是有限的；
[0065] 因此，体装式太阳能帆板故障时可供反作用飞轮的实际功率约束的表达式如下：
[0066]
[0067] 其中， 为k时刻各体装式太阳能帆板能提供的瞬时功率向量；κ为可供执行机构使用的功率占比；Ek＝(e1，e2，…，em)为k时刻各太阳能帆板的故障因子序列，其中0≤ei≤1，i＝1，2，…，m； 为k时刻反作用飞轮组合为了提供期望姿态控制力矩而所需的总功率，其具体形式如下：
[0068]
[0069] 其中， 为k时刻反作用飞轮组合角动量关于时间的导数， 为 k时刻反作用飞轮组合的角加速度，上标T表示转置，
[0070] 将式(21)离散化后得到
[0071]
[0072] 将式(22)代入式(20)中得到：
[0073]
[0074] 将所述反作用飞轮的角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束整理成如下形式：
[0075]
[0076] 其中， 04L+1表示4L+1维零向量；L为反作用飞轮组合中的反作用飞轮总个数。
[0077] 进一步，在航天器姿态控制过程中在满足高控制精度的同时尽可能减少能量消耗情况下，步骤S4中设计的目标性能函数为：
[0078]
[0079] 其中，N为预测步长，下标j|k表示第k时刻下的第j预测时刻；x0|k为第k时刻航天器姿态控制系统预测步长内的初始状态量，xj|k＝(qj|k，ωj|k)T为第 k时刻预测步长内第j步航天器姿态控制系统状态量，表征航天器的姿态和角速度； 为控制系统的期望状态； 为第k时刻预测步长内第j步系统控制量；Q，P为状态量和控
制量的常值正定权重矩阵。
[0080] 进一步，步骤S5具体过程如下：
[0081] 基于步骤S1中建立的模型预测控制方法的预测模型、步骤S4中设计的目标性能函数和步骤S3中建立的反作用飞轮的角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束条件，将航天器姿态控制问题转化为如下优化问题：
[0082]
[0083] subject to xj+1|k＝F(xj|k，uj|k) (27)
[0084] K(xj|k，uj|k)≤04L+1 (28)
[0085] 其中，xj+1|k＝F(xj|k，uj|k)为航天器姿态控制系统状态方程；K(xj|k，uj|k)为包括反作用飞轮的角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束的系统约束；航天器姿态控制系统第k时刻预测步长内第j步的状态方程和系统约束的具体形式如下：
[0086]
[0087]
[0088] 求解所述优化问题得到系统在k时刻下满足约束条件的最优控制序列：
[0089]
[0090] 其中，上标*表示最优解，
[0091] 将所述最优控制序列中的第一项作为k时刻航天器姿态控制系统的最优控制输入量，即
[0092]
[0093] 最后，令 求解下一时刻航天器姿态控制系统的最优控制输入量，完成步骤S3中建立的约束下的航天器姿态控制任务。
[0094] 本发明的有益效果：
[0095] (1)与现有的考虑反作用飞轮特性的航天器姿态控制方法相比，本发明在考虑反作用飞轮角动量饱和、力矩饱和的基础上，同时考虑了体装式太阳能帆板故障时引起的反作用飞轮实际可用功率受限问题，更全面地分析反作用飞轮在太空中为航天器提供姿态控制力矩时可能遇到的实际约束；
[0096] (2)与现有的反作用飞轮瞬时功率最优的航天器姿态控制方法相比，现有最优方法求得的瞬时最优功率可能仍大于太阳能帆板故障时所能提供的最大可用功率，本发明能够克服这一缺点，不论太阳能帆板所能提供的可用功率有多大，总能使航天器用于姿态控制的需用功率满足系统要求，更符合实际工程需要。
[0097] (3)本发明中描述的反作用飞轮的实际可用功率约束是一个动态约束，它不仅与太阳能帆板的可能发生的故障相关，同时也与航天器的当前所处的轨道位置、当前姿态角密切相关，求解这样的动态约束下的航天器姿态控制问题，更具有工程意义。
附图说明
[0098] 图1为本发明的考虑反作用飞轮特性及功率受限的航天器姿态控制方法的流程图。
具体实施方式
[0099] 下面结合附图和实施方式对本发明作进一步详细说明。
[0100] 本发明的考虑反作用飞轮特性及功率受限的航天器姿态控制方法基于模型预测控制策略，考虑的是在反作用飞轮角动量、控制力矩以及功率受限情况下的航天器姿态控制问题。如图1所示，本发明首先建立四元数描述方式下的航天器运动学与动力学模型；其次建立航天器体装式太阳能帆板输出功率与太阳入射角的关系；然后考虑反作用飞轮特性建立角动量约束和控制力矩约束，以及根据航天器体装式太阳能帆板输出功率与太阳入射角的关系建立体装式太阳能帆板故障时反作用飞轮的实际可用功率约束；然后设计目标性能函数使得航天器在姿态机动过程中满足高精度控制的同时尽可能减少能量消耗；最后，基于模型预测控制策略求解最优控制力矩，从而完成上述约束条件下的航天器姿态控制任务。
[0101] 其具体操作步骤如下：
[0102] 第一步，首先建立四元数描述方式下的航天器姿态运动学与动力学模型如下：
[0103]
[0104] 其中，J＝diag(1，0.5，0.8)kg·m2∈R3×3为航天器惯量矩阵；ω＝(ω1，ω2，ω3)T为航天器本体三轴角速度；H∈R3为航天器总角动量，其具体形式为 H＝Jω+HRWA，而HRWA∈R3为反作用飞轮组合角动量在航天器惯量主轴上的投影；uc∈R3为作用于航天器的姿态控制力矩，本发明中的姿态控制力矩由反作用飞轮组合提供；q＝(q0，q1，q2，q3)T为描述航天器姿态的四元数，且满足 Ξ(ω)是斜反对称矩阵，其具体形式如下：
[0105]
[0106] 对于任意一个向量x∈R3，其对应的斜反对称阵x×形式如下：
[0107]
[0108] 若忽略航天器轨道角速度对飞轮角动量的影响，给出反作用飞轮组合动力学模型：
[0109]
[0110] 其中，τRWA∈RL为由L个反作用飞轮产生的控制力矩；hRWA∈RL为反作用飞轮组合的角动量， 为反作用飞轮组合角动量关于时间的导数； JRWA∈RL×L为反作用飞轮组合的转动惯量矩阵，其形式为为各个反作用飞轮
的转动惯量； 为各个反作用飞轮的转速。在本
实施例中，飞轮数量L＝4，且反作用飞轮组合构型为三正交一斜装。
[0111] 由反作用飞轮组合的动力学模型及反作用飞轮组合构型可知，反作用飞轮组合作用在航天器上的姿态控制力矩为：
[0112]
[0113] 其中，C∈R3×4为反作用飞轮组合安装矩阵，具体形式为；
[0114]
[0115] 再由角动量守恒定理可知，在不受外力矩作用下，航天器系统总角动量守恒，因此作用于航天器的姿态控制力矩为
[0116]
[0117] 将反作用飞轮提供的姿态控制力矩代入航天器姿态动力学模型可得由反作用飞轮组合作为执行机构的航天器姿态运动学与动力学模型：
[0118]
[0119] 为了将上述动力学模型作为模型预测控制方法的预测模型，需要将其离散化，取采样时间间隔为Δt＝0.4s，则离散化后的动力学模型为：
[0120]
[0121] 其中，下标k和k+1表示k时刻和k+1时刻；Hk为k时刻的航天器总角动量；I4×4为四阶单位矩阵；上标x为斜反对称矩阵。
[0122] 令xk＝(qk，ωk)T为第k时刻航天器姿态控制系统状态量，表征航天器的姿态和角速度， 为第k时刻系统控制量，可将上述离散动力学模型简化为如下形式：
[0123]
[0124] 若定义
[0125]
[0126] 则有
[0127] xk+1＝F(xk，uk) (10)
[0128] 第二步，建立航天器体装式太阳能帆板输出功率与太阳入射角的关系：
[0129]
[0130] 其中，Ptotal为航天器上m个体装式太阳能帆板的输出总功率，为航天器提供能源完成科学探测及姿态控制等任务，这里取m＝4； 为第i个太阳能帆板的输出功率； 为第i个太阳能帆板的有效面积；θi为太阳入射角，表示第i个太阳能帆板法线与太阳光的夹角；PEOL(θi)为太阳入射角θi的函数，表示第i个太阳能帆板的单位面积功率，其具体关系式如下：
[0131]
[0132] 其中，Ps＝1358W/m2为入射太阳辐射；η＝0.14为太阳能帆板光电转换率； Id＝0.77为太阳能帆板固有退化因子；Ld＝0.625为太阳能帆板寿命退化因子。
[0133] 但考虑到航天器在轨运行时，太阳入射角不是一个固定值，它随航天器运行轨道、航天器当前时刻的姿态角以及所处的轨道位置变化而变化，为了获得航天器在轨道上任意位置的瞬时可用总功率，建立太阳入射角模型如下：
[0134] θi＝arccos，i＝1，2，...，m (13)
[0135] 其中，ni∈R3，|ni|＝1为第i个太阳能帆板在航天器本体坐标系中的单位法向量；<·，·>表示内积运算；rsb∈R3为太阳单位矢量在航天器本体坐标系中的投影，有如下转换关系：
[0136] rsb＝RbIrsI (14)
[0137] 其中，rsI＝rs/|rs|∈R3为太阳在地心惯性坐标系下的位置矢量rs∈R3的单位向量；RbI为地心惯性坐标系到航天器本体坐标系的转换矩阵，与欧拉转动顺序以及姿态欧拉角有关；
[0138] 将航天器体装式太阳能帆板输出功率与太阳入射角的关系转化为与航天器姿态角相关的表达式：
[0139]
[0140] 其中，Psolar∈Rm为体装式太阳能帆板输出功率向量。
[0141] 若欧拉转动顺序为“2-3-1”，则地心惯性坐标系到航天器本体坐标系的转换矩阵RbI的形式如下：
[0142]
[0143] 其中，Rx(γ)为绕ox轴转过滚转角γ的基元变换矩阵；Ry(ψ)为绕oy轴转过偏航角ψ的基元变换矩阵； 为绕oz轴转过俯仰角 的基元变换矩阵；上述三个基元变换矩阵形式如下：
[0144]
[0145] 展开可得：
[0146]
[0147] 此时，航天器体装式太阳能帆板输出功率与太阳入射角的关系表达式 (15)可以转化为与航天器姿态角相关的表达式：
[0148]
[0149] 其中，ψ， γ分别表示航天器的偏航、俯仰、滚转姿态欧拉角。
[0150] 第三步，基于第二步航天器体装式太阳能帆板输出功率与航天器姿态角的关系，进一步考虑航天器执行机构反作用飞轮的特性，并建立角动量约束、控制力矩约束以及体装式太阳能帆板故障时反作用飞轮的实际可用功率约束：
[0151] 由于反作用飞轮属机械结构，其存在角动量约束以及控制力矩约束，具体形式如下：
[0152] 1)角动量约束：反作用飞轮的角动量体现在飞轮转速，则角动量约束可以用转速约束来表示，即
[0153]
[0154] 其中， 为反作用飞轮组合的最大转速向量，这里取
[0155] 将上式变形为增量形式：
[0156]
[0157] 2)控制力矩约束：
[0158]
[0159] 其中， 为反作用飞轮组合的最大控制力矩向量，这里取
[0160] 由于航天器姿态控制系统为离散系统形式，故控制力矩约束可用如下角加速度的离散形式来近似表示：
[0161]
[0162] 3)体装式太阳能帆板故障时反作用飞轮的实际可用功率约束：
[0163] 首先，假设航天器进入地影区的时间可以忽略不计，即近似认为航天器上的体装式太阳能帆板总是能受太阳光照射；
[0164] 其次，考虑到航天器上体装式太阳能帆板的输出功率需要供有效负载、执行机构及其他分系统共同使用；同时还要考虑当航天器运行到不同的轨道位置时，太阳能帆板能提供的输出功率是随太阳入射角或者说随航天器姿态角变化的，而且太阳能帆板也有可能因复杂的太空环境而发生故障，因此可供执行机构使用的实际功率是有限的；
[0165] 因此体装式太阳能帆板故障时可供反作用飞轮组合使用的实际功率约束的表达式如下：
[0166]
[0167] 其中， 为k时刻的各体装式太阳能帆板能提供的瞬时功率向量，将 k时刻下的偏航角ψk、俯仰角 和滚转角γk具体代入，可得到当前k时刻、当前姿态下航天器各体装式太阳能帆板能提供的瞬时功率向量 具体形式如式(15-3)所示；κ为可供执行机构使用的功率占比； Ek＝(e1，e2，…，em)为k时刻各太阳能帆板的故障因子序列，其中 0≤ei≤1，i＝1，2，...，m； 为k时刻反作用飞轮组合产生期望姿态控制力矩所需的总功率，其具体形式如下：
[0168]
[0169] 其中， 为k时刻反作用飞轮组合角动量关于时间的导数， 为k 时刻反作用飞轮组合的角加速度，上标T表示转置，
[0170] 将上式离散化后可得
[0171]
[0172] 整理后得体装式太阳能帆板故障时可供反作用飞轮组合使用的实际功率约束为：
[0173]
[0174] 其中，欧拉角ψk， γk与四元数qk之间存在相互转换关系，因此上式可表示为：
[0175]
[0176] 根据上述定义 将上述反作用飞轮角动量约束、控制力矩约束和体装式太阳能帆板故障时可供反作用飞轮组合使用的实际功率约束整理成如下形式：
[0177]
[0178] 其中，04L+1表示4L+1维零向量；L＝4为反作用飞轮组合中的反作用飞轮总个数。
[0179] 第四步，考虑到航天器姿态控制过程中要在满足高控制精度的同时尽可能减少能量消耗，设计步骤S4中的目标性能函数为：
[0180]
[0181] 其中，N＝5为预测步长，下标j|k表示第k时刻下的第j预测时刻；x0|k为第k时刻航天器姿态控制系统预测步长内的初始状态量，xj|k＝(qj|k，ωj|k)T为第k时刻预测步长内第j步航天器姿态控制系统状态量，表征航天器的姿态和角速度； 为控制系统的期望状态，由具体任务决定其形式； 为第k时刻预测步长内第j步系
统控制量； Q＝0.7I7×7，P＝0.3I4×4为状态量和控制量的常值正定权重矩阵。
[0182] 第五步，基于第三步的约束条件以及第四步设计的目标性能函数，基于模型预测控制策略完成有约束的航天器姿态控制任务：
[0183] 由离散化动力学模型、目标性能函数和约束条件，可以将控制问题转化为如下优化问题：
[0184]
[0185] subject to xj+1|k＝F(xj|k，uj|k) (27)
[0186] K(xj|k，uj|k)≤04L+1 (28)
[0187] 其中，xj+1|k＝F(xj|k，uj|k)为航天器姿态控制系统状态方程；K(xj|k，uj|k)为系统约束，包括反作用飞轮角动量约束、控制力矩约束和体装式太阳能帆板故障时反作用飞轮的实际可用功率约束；航天器姿态控制系统第k时刻预测步长内第j步的状态方程和系统约束的具体形式如下：
[0188]
[0189]
[0190] 求解上述优化问题可得系统在k时刻下满足约束条件的最优控制序列：
[0191]
[0192] 将该最优控制序列中的第一项作为k时刻航天器姿态控制系统的最优控制输入量，即
[0193]
[0194] 最后，令 求解下一时刻航天器姿态控制系统的最优控制输入量，即可完成有约束的航天器姿态控制任务。
[0195] 本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。
