基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法

797.)，并将该算法应用于高超声速飞行器的最省燃料爬升问题，可以以较少的时间代价获得精度较高的解。
6.现有论文文献，如：颜楚雄,童轶男,宋加洪,卢宝刚,赵亮博.基于多区间radau伪谱法带末端约束的轨迹优化[j].战术导弹技术,2021(02):94-100.doi:10.16358/j.issn.1009-1300.2021.9.093.采用多区间radau伪谱法设计多重约束的导弹末端轨迹优化，设计得到的弹道曲线光滑，较易实现。洪蓓等人将hp自适应伪谱法应用于滑翔弹道快速优化问题和再入轨迹优化问题(洪蓓,辛万青.hp自适应伪谱法在滑翔弹道快速优化中的应用[j].计算机测量与控制,2012,20(05):1283-1286.doi:10.16526/j.cnki.11-4762/tp.2012.05.018.)。苏志娟等人提出了一种基于密度函数的改进自适应伪谱法(苏志娟.基于改进自适应伪谱法的中低速磁浮列车节能驾驶[d].北京交通大学,2021.doi:10.26944/d.cnki.gbfju.2021.001740.)，得到了中低速磁浮列车最优节能驾驶曲线，但是未考虑控制力率影响，会出现控制力曲线脉冲现象。
[0007]
尽管现有hp自适应法的研究取得了很大进展，但在追求高精度求解情况下，降低计算复杂度仍需进一步探索，尤其是在将传统hp自适应伪谱法应用在处理多约束遥感卫星的最优姿态规划问题时，随着精度的不断提升，导致网格细化次数大量增加，不可避免地增加迭代次数，耗费更多的星上计算机资源，降低卫星在轨应用效能。
[0008]
本发明提出的改进hp自适应方法可以在满足精度要求的前提下，通过改进采样点和迭代判据来降低算法计算量，显著提高了求解效率。仿真结果表明，与目前已有的hp自适应伪谱法相比，改进的hp自适应伪谱法能够显著减少网格迭代次数，以较少的时间代价获得较高精度的解，验证了改进算法的优越性。


技术实现要素：

[0009]
本发明为解决现有hp自适应法在处理多约束遥感卫星的最优姿态规划问题中，易导致网格细化次数大量增加，不可避免地增加迭代次数，耗费更多的星上计算机资源，进而降低卫星在轨应用效能等问题，提供一种基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法。
[0010]
基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，该方法由以下步骤实现：
[0011]
步骤一、建立连续最优控制问题模型，将各类约束转化为最优控制问题，并在满足多约束的前提下，使得优化性能指标最小；
[0012]
步骤二、将所述最优控制问题离散化，转化为非线性规划问题进行求解，利用改进hp自适应radau伪谱法法进行配点；具体配点过程为：
[0013]
步骤二一、采用误差评估准则进行判定：
[0014]
采用高斯分布选择采样点，判断误差判定准则εk是否大于微分代数方程的误差门限ε；若是，采用基于改进hp自适应伪谱法重新细化网络，执行步骤二二；并设定误差判定准则作为采样点在动态约束方程中的残差εk，否则执行步骤二五；
[0015]
步骤二二、采用网格分类细化准则进行判定：
[0016]
判断状态曲率标准ρk是否大于比例因子ρ
max
；若是，细化区间提高精度，执行步骤二三；若否，增加配点数提高精度；执行步骤二三；
[0017]
步骤二三、获得新的网格点和配点；
[0018]
步骤二四、将步骤二三获得的新的网格点和配点代入步骤二中，进行新的求解；
[0019]
步骤二五、完成解的迭代，以当前数值解作为最终输出。
[0020]
本发明的有益效果：本发明方法首先结合卫星动力学，运动学和挠性附件的振动方程，建立了面向最优控制方法的刚挠耦合状态空间方程，其次在建立了包含代数-微分约束，控制力矩约束，角动量约束，路径约束以及以时间最短为优化目标函数的基础上，采用改进hp自适应radau伪谱法将最优控制问题转化为一般的非线性规划问题进行求解。
[0021]
本发明方法针对高分辨率挠性遥感卫星多约束姿态最优机动规划问题，为了提高求解最优问题的效率，提出了一种改进的hp自适应radau伪谱法。本方法在网格细化过程中对易发生突变点的区域分配更多采样点，同时修改了曲率迭代判据，以保证突变点出分配更多网格点，该方法大大减少了网格迭代次数，缩短了计算时间。
[0022]
本发明方法以挠性卫星姿态机动任务需求为依据，建立了刚挠耦合模型，在保证星载设备正常运行、不超过卫星执行机构能力的前提下，以某一实际任务性能为优化指标，设计了最优的挠性遥感卫星姿态机动规划方法，相比于传统hp自适应伪谱法，在高精度下该方法生成满足所设需求的最优轨迹的仿真时间更少，效率更高。
[0023]
仿真结果表明，本发明方法可以在满足卫星机动任务需求下规划出最优的姿态机动轨迹，能够实现快速且高精度的机动，本发明方法为高分辨率遥感卫星执行复杂空间任务提供了一定参考。
附图说明
[0024]
图1为本发明所述的基于改进hp自适应伪谱法的挠性遥感卫星姿态机动规划方法的流程图；
[0025]
图2为采用本发明所述的基于改进hp自适应伪谱法的挠性遥感卫星姿态机动规划方法进行仿真的效果图；其中，(a)为三轴欧拉角变化效果图；(b)为三轴角速度变化效果图；(c)为三轴飞轮控制力矩变化效果图；(d)为挠性附件振动位移变化效果图；(e)挠性附件振动速度变化效果图；
[0026]
图3为改进前后卫星三轴姿态角对比效果图；
[0027]
图4为改进前后卫星三轴角速度对比效果图；
[0028]
图5为改进前后卫星三轴力矩对比效果图；
[0029]
图6中(a)和(b)分别为改进前后采样点分布对比效果图。
具体实施方式
[0030]
具体实施方式一、结合图1说明本实施方式，基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，该方法由以下步骤实现：
[0031]
步骤一、建立连续最优控制问题模型，将各类约束转化为最优控制问题，并在满足多约束的前提下，使得优化性能指标最小；
[0032]
本实施方式中，最优控制问题模型具体为：挠性卫星姿态动力学及运动学建模；
[0033]
根据动量和角动量定理，以惯性坐标系为参考坐标系，由卫星本体及其挠性部件组成的系统刚挠耦合运力学方程为：
[0034][0035]
式中，为卫星转动惯量矩阵；为卫星相对惯性系的角速度；为挠性附件与卫星的耦合系数；η为挠性附件在模态坐标系下的位移；为执行机构的角动量；为卫星受到的空间干扰力矩。符号定义为：
[0036][0037]
对于模型中存在的空间干扰力矩td可以通过姿态控制算法进行补偿，在设计卫星姿态轨迹时先不考虑空间干扰的影响。由欧拉角描述的卫星运动学模型，定义θ
x
,θy,θz为3-2-1转序下的滚动角、俯仰角和偏航角。且考虑到卫星惯性定向及初始角度较小，因此卫星角速度可近似为
[0038]
挠性附件的振动方程在模态坐标系下可以表示为:
[0039][0040]
式中，为挠性附件模态的阻尼比矩阵；为挠性附件模态的振动频率矩阵。
[0041]
定义状态变量x＝[θ η ω ω
η
]
t
，θ为卫星角度，ω
η
为挠性附件在模态坐标系下的速度，联立卫星姿态动力学和运动学方程(1)和挠性附件的振动方程(2)，得到状态空间方程如下所示：
[0042][0043]
式中：
[0044][0045]
式中，i为单位矩阵，该方程描述了卫星刚体及其挠性部件的姿态变化，表现了这两者的强非线性耦合关系。
[0046]
本实施方式中，将微分-代数约束、路径约束、执行机构力矩约束及角速度测量能力约束的轨迹优化问题转化为有限时间内的最优控制问题，即在时间区间t∈[t0,tf]内已知状态变量初值x(t0)，通过设置控制力矩t(t)，使得卫星姿态机动至目标姿态x(t0)，在满
足约束的前提下，使得卫星姿态的机动时间最优，即时间段(t
f-t0)最小。t0表示开始机动时刻；tf表示机动完成时刻。
[0047]
卫星在姿态机动过程中，需要考虑的约束主要为来自卫星自身的动力学约束和来自外部环境的几何约束。
[0048]
动力学约束为星上控制力矩t＝[t
x t
y tz]
t
有界和ω＝[ω
x ω
y ωz]
t
有界。
[0049]
由于卫星机动过程中需要规避一些危险的指向，将这类危险的姿态指向称为姿态禁区；同时在进行对地通信中卫星需要保持天线对准地面站，这类区域称为强制区域。将其抽象成几何约束主要包括“禁止区域约束”和“强制区域约束”两种。以卫星光敏器件设备视场不得靠近明亮天体的“禁止区域”为例：
[0050]
设卫星此时机动到姿态[a
i β
i γi]
t
，卫星指向明亮天体的矢量在地心赤道坐标系中为[x
i y
i zi]
t
；视场角为λ，则卫星在机动过程中应满足如下约束：
[0051][0052]
同理，若是“强制区域”则在机动过程应该满足
[0053][0054]
因此，具有bolza型的最优控制问题可以定义如下：
[0055][0056]
式中：表示当前状态值；φ表示边界性能值；g表示系统动态性能值。λ1和λ2分别为禁止区域和强制区域的视场角。
[0057]
约束条件为：
[0058][0059]
其中，f、c以及分别为微分-代数约束、路径约束以及边界约束。
[0060]
步骤二、将所述最优控制问题离散化，转化为非线性规划问题进行求解，利用改进自适应radau伪谱法法进行配点；
[0061]
本实施方式中，伪谱法常用于航天器轨迹优化问题，采用自适应radau伪谱法求解最优轨迹规划问题。其配点采用legendre-gauss-radau(lgr)点。通过将整个时间区间划分为多个子区间，在每个子区间内通过lgr配点的方式实现将连续的最优控制问题离散化。获得状态量和控制量的数值解；将非线性最优控制问题转化为非线性规划问题，最终通过非线性规划求解器进行求解，例如工具包snopt。
[0062]
本实施方式中，采用改进hp自适应多区间配点方法实现配点，对于约束简单且光滑的最优控制问题，自适应p法能以较少的计算代价得到高精度的解，但在解决不存在高阶维数的不连续或者非光滑最优控制问题，自适应h法更能得到高精度的解，而针对一些在某些区段表现为平滑，某些区段表现为不平滑时，两种方式的结合更能得到最佳的近似解。本实施方式中的一种基于radau伪谱法，结合重新细化区间的自适应h法和增加插值多项式位
数的自适应p法的优化采样点选取及迭代判据选取的改进hp自适应伪谱方法，其拥有自适应h法的计算稀疏性和自适应p法的快速收敛性，能更好找到不连续点和奇点的位置，以达到求解精度和计算代价的平衡。具体配点过程为：
[0063]
1、radau伪谱法首先将问题的整体时间区间划分为k个子区间，每个子区间sk＝[t
k-1
,tk],k＝1,2,3,...k通过变换映射到区间τ∈[-1,1]上，t
k-1
和tk为第k个区间的初始时刻和末端时刻；在子区间sk＝[t
k-1
,tk]内选择初始网格点[u0,u1,...,uk]及配点
[0064]
2、采用步骤二一中网格点及配点，利用nlp求解器将连续的最优控制问题离散化，得到状态量和控制量数值解
[0065]
3、采用误差评估准则进行判定：利用高斯分布的采样点，得到动态约束方程的残差εk，若εk≤ε，则进入步骤7，否则进入步骤4，定义ε为微分代数方程的所设误差门限；
[0066]
4、利用网格分类细化准则进行判定：若ρk≤ρ
max
，则通过增加配点数提高精度，否则通过细化区间提高精度；
[0067]
5、得到新的网格点和配点。
[0068]
6、将二五中重新细化的网格点和配点代入到步骤2中进行新一轮求解。
[0069]
7、完成解的迭代，以当前数值解作为最终输出。
[0070]
本实施方式中，hp自适应伪谱法首先将多区间连续的最优控制问题转换为非线性规划问题，参照全局伪谱法解的误差判定准则，以区段[t
k-1
,tk]为例，令nk为第k区段中的lgr配点数，k∈[1,
…
,k]，传统的hp自适应法在进行误差评估时，会选取相邻lgr配点的中心点作为采样点，而在各区间端点处是配点分布最多的地方，这就导致易出现突变点的区间中部分配到的采样点少，为了达到预设的精度网格就需要继续细化，增加网格迭代次数，计算效率大大降低。
[0071]
本实施方式中按照高斯分布选取的采样点，相比于传统采样点，针对lgr配点采样分布更加合理，有效地弥补了易发生突变点的区间中段采样值过少导致求解效率下降的情况。在第k个区间内以高斯分布选取500个点作为采样点。采取此种采样的原因是因为在段内中部产生突变点的概率远远大于段两端。采样点改进前后对比如图6所示。
[0072]
本实施方式的步骤3中，误差判定准则εk为采样点在动态约束方程的残差，可表示为：
[0073][0074]
其中，为第s个采样点处的第v个状态变量的一阶微分残差；为第s个采样点处的第w个路径约束残差。
[0075][0076]
[0077]
式中，为第s个采样点的一阶微分残差，为第s个采样点的路径约束残差，为状态变量在采样点s处的导数，[t
k-1
,tk]为采样点时间区间，f为状态方程微分约束，c为路径约束，为第s个采样点处的状态变量值，为第s个采样点处的控制变量值。
[0078]
本实施方式中的步骤4中，采用基于状态曲率的网格分类细化准则：
[0079]
定义为最大残差所在的状态变量所构成的列向量，即
[0080][0081]
同时，第l个状态变量在第p个采样点的曲率可以表示为：
[0082][0083]
对进行排序：
[0084][0085]
为了进一步提高网格细化效率，结合曲率，将分类细化的标准设置为各采样点处曲率最大值与标准差比值的四分之一：
[0086][0087]
式中，为第k个区间各采样点的曲率最大值，n
k-1为采样点数，为第l个状态在第i个采样点处的曲率，为各采样点处曲率的平均值。
[0088]
传统hp自适应法的分类细化标准为各采样点处最大值与平均值之比，主要反应离散数据偏离平均值的偏离度，本文选择标准差代替平均值，更能体现离散点的偏离度，根据对采样点的改进图可以看出，采样点的区间缩减为原区间的1/4，采样点处的细化标准也相应缩减为1/4倍。
[0089]
本实施方式中，定义ρ
max
为判断细化分区还是提高插值多项式维数的比例因子，为给定常数。
[0090]
1)当ρk中元素都比ρ
max
小时，此时，可认为轨迹相对比较平滑，采用增加区间内配点数的方式来提高求解精度。
[0091]
设和为更新前后每子区间内的配点数，可以由下式求得：
[0092][0093]
其中c1为任意整数常量，控制每一段内配点数的增长数，ceil为向上取整函数。
[0094]
2)当ρk中有比ρ
max
大的元素时，可认为轨迹不平滑，采用细化区间的方式来提高求解精度。分段数seg可由下式求得。
[0095][0096]
其中，c2为任意的整数常数，来控制子区间段的个数。分段至少要分为两段，故seg不能少于2。新增分段点位置为ρk中元素的局部最大值所对应的采样点处。
[0097]
采用改进自适应伪谱法中新增分区和插值多项式维数是变化的，这能大大提高计算效率。
[0098]
具体实施方式二、结合图2至图6说明本实施方式，本实施方式为采用具体实施方式一所述的基于改进hp自适应伪谱法的挠性遥感卫星姿态机动规划方法进行仿真的实施例：
[0099]
以某型具有挠性帆板的光学遥感卫星进行终端角度和角速度不为零的姿态机动场景为例，来说明本实施方式的规划方法的有效性和在计算耗时等方面的性能提升。具体地，考虑的卫星转动惯量如下：
[0100][0101]
边界约束为：
[0102][0103]
路径约束：禁止区域的矢量在地心赤道坐标系中为[15.7
°
,-9.7
°
,14.6
°
]，视场角λ1为10
°
；强制区域的矢量在地心赤道坐标系中为[17.8
°
,0.55
°
,19.26
°
]，视场角λ2为30
°
。
[0104]
为验证本方法在遥感卫星快速机动任务中的有效性，最优问题目标函数为时间最优。
[0105]
minj＝tf[0106]
考虑的挠性帆板一阶频率为2.23hz，阻尼为0.032，刚柔耦合系数矩阵为[0.01,0.0383,0.01]。飞轮最大角速度为5rad/s，三轴控制力矩范围为[-0.25nm0.25nm]。
[0107]
采用改进hp自适应伪谱法对此最优控制问题进行求解，将原问题转化为非线性规划问题，利用matlab外接的软件包gpops2进行求解，采用的nlp求解器为snopt。
[0108]
优化计算过程中软件环境为：windows 10专业版64位造作系统和matlab r2020b；硬件环境为：intel(r)core(tm)i5-8400处理器，12g内存。
[0109]
在相同的误差门限ε下，采用改进hp自适应radau伪谱法解决该问题，并与传统hp自适应法进行比较，最大误差为在误差门限内达到的最大误差，可以在一定程度上衡量结果的精度。仿真相关参数设置：误差门限ε＝10-5
，比例因子ρ
max
＝2，初始网格点选取[-1,1]，初始lgr配点数为3个。利用改进hp自适应伪谱法对所有状态变量及控制变量进行仿真，得到结果如图2所示。
[0110]
在10-5
精度下改进方法仿真结果表明，在仿真时间2.28275s时间时，实现了挠性遥感卫星在25s内机动27
°
，末端角速度不为0，规避姿态禁区，且挠性部件在模态坐标下振动位移小于6*10-3
m，角速度小于8*10-3
m/s的性能。本文提出的方法可以有效解决所提出的刚挠耦合系统存在复杂约束的最优控制问题。
[0111]
在相同的误差门限ε下，采用本文提出改进hp自适应radau伪谱法解决该问题，并与传统hp自适应方法进行比较。得到结果如图3、4、5所示。
[0112]
在不同求解精度下，分别采用hp自适应伪谱法和改进hp自适应伪谱法的优化结果如表1所示，表1为不同精度下伪谱法的优化比较结果。
[0113]
表1
[0114][0115]
从对比仿真结果可以看出本发明方法可以更高效地解决挠性遥感卫星快速机动问题。与传统hp自适应伪谱法相比，如图3-图5所示，两种方法的求解结果基本相同。但本发明方法在网格细化方面具有明显优势，网格迭代次数较少、最大误差较小并且求解效率更高，随着精度的增高，这些优势也在不断的变大，可以在一定程度上消解利用传统伪谱法求解高精度姿态规划问题时求解时间过长的影响。
[0116]
本实施方式所述的方法针对高分辨率挠性遥感卫星在进行对地观测任务时的多约束姿态最优机动规划问题，提出了一种改进的hp自适应radau伪谱法。本文提出改进方法采用双层优化迭代策略，以微分-代数约束在以高斯分布的采样点处的残差为误差评估准则进行优化，以各采样点处曲率最大值与标准差比值的四分之一为细化准则进行二次优化，这两点共同提高了该方法的求解效率。
[0117]
本方法能够在最短时间内生成满足各项约束要求的最优轨迹，且相比于传统hp自适应伪谱法在同精度下耗时更短，具有更高的求解效率。本方法为高分辨率遥感卫星执行复杂空间任务提供了设计依据，具有工程意义。
[0118]
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0119]
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

技术特征：
1.基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，其特征是：该方法由以下步骤实现：步骤一、建立连续最优控制问题模型，将各类约束转化为最优控制问题，并在满足多约束的前提下，使得优化性能指标最小；步骤二、将所述最优控制问题离散化，转化为非线性规划问题进行求解，利用改进自适应radau伪谱法法进行配点；具体配点过程为：步骤二一、采用误差评估准则进行判定：采用高斯分布选择采样点，判断误差判定准则ε
k
是否大于微分代数方程的误差门限ε；若是，采用基于改进hp自适应伪谱法重新细化网络，执行步骤二二；并设定误差判定准则作为采样点在动态约束方程中的残差ε
k
，否则执行步骤二五；步骤二二、采用网格分类细化准则进行判定：判断状态曲率标准ρ
k
是否大于比例因子ρ
max
；若是，细化区间提高精度，执行步骤二三；若否，增加配点数提高精度；执行步骤二三；步骤二三、获得新的网格点和配点；步骤二四、将步骤二三获得的新的网格点和配点代入步骤二中，进行新的求解；步骤二五、完成解的迭代，以当前数值解作为最终输出。2.根据权利要求1所述的基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，其特征在于：步骤二中，采用自适应radau伪谱法求解最优控制问题，具体为：利用伪谱法将状态量和控制量离散，获得状态量和控制量的数值解；将非线性最优控制问题转化为非线性规划问题，通过非线性规划求解器进行求解。3.根据权利要求1所述的基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，其特征在于：步骤二二中，设定所述状态曲率标准ρ
k
为各采样点处曲率最大值与标准差比值的四分之一；用下式表示为：式中，为第k个区间各采样点的曲率最大值，n
k-1为采样点数，为第l个状态在第i个采样点处的曲率，为各采样点处曲率的平均值。4.根据权利要求1所述的基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，其特征在于：步骤二二中，采用增加区间内配点数的方式提高求解精度，用下式表示为：设定和为更新前后每个子区间内的配点数，由下式求得：式中，c1为任意整数常量，控制每一段内配点数的增长数，ceil()为向上取整函数。5.根据权利要求1所述的基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，其特征在于：步骤二二中，采用细化区间的方式来提高求解精度；分段数seg用下式表示为：
式中，c2为任意的整数常数，用于控制子区间段的分段数；新增分段点位置为ρ
k
中元素的局部最大值所对应的采样点处。

技术总结
基于改进自适应伪谱法的挠性遥感卫星姿态机动规划方法，涉及挠性遥感卫星姿态控制技术领域，解决现有方法易导致网格细化次数大量增加，不可避免地增加迭代次数，耗费更多的星上计算机资源，进而降低卫星在轨应用效能等问题，首先结合卫星动力学，运动学和挠性附件的振动方程，建立面向最优控制方法的刚挠耦合状态空间方程，其次在建立包含代数-微分约束，控制力矩约束，角动量约束，路径约束以及以时间最短为优化目标函数的基础上，采用改进hp自适应伪谱法将最优控制问题转化为非线性规划问题进行求解。本方法能够在最短时间内生成满足各项约束要求的最优轨迹，且相比于传统hp自适应伪谱法在同精度下耗时更短，具有更高的求解效率。效率。效率。


技术研发人员：范国伟 高宇 刘泉志 肖阳 吕雪莹 张乐 张刘
受保护的技术使用者：吉林大学
技术研发日：2023.04.24
技术公布日：2023/7/12