-
深空探测器在轨运行时需要进行大量的姿态机动来完成各种跟踪和探测等任务,姿态系统需要时刻规划出有效合理的姿态机动路径。所以,对深空探测器姿态机动过程中约束的建模和分析是设计深空探测器姿态机动路径规划算法的前提。
本文将探测器视为刚体。为避免奇异性,采用四元数形式来描述姿态运动学和动力学约束探测器姿态机动
$$ {\dot{ q}} = \frac{1}{2}{{Q\omega }} = \frac{1}{2}{{\varOmega }}q $$ (1) $$ {{J\omega }} = {{T}} - {{{\omega }}^ \times }{{J\omega }} $$ (2) 其中:
${{q}} = {\left[ {{q_0},{q_1},{q_2},{q_3}} \right]^{\rm{T}}}$ ,满足归一化约束即${\left\| {{q}} \right\|_2} = 1$ ;${{J}} = {\rm{diag}}\left( {{J_1},{J_2},{J_3}} \right)$ 表示探测器相对本体系的惯性矩阵,${{u}} = {\left[ {{T_1},{T_2},{T_3}} \right]^{\rm{T}}}$ 为控制力矩在本体系下的分量;${{\omega}} = {\left[ {{\omega _1},{\omega _2},{\omega _3}} \right]^{\rm{T}}}$ 是探测器相对惯性系的角速度在本体系下的表示,${{{\omega }}^ \times }$ 为${{\omega}} $ 的叉乘矩阵,具体形式为$$ {{{\omega }}^ \times } = \left[ {\begin{array}{*{20}{c}} 0&{ - {\omega _3}}&{{\omega _2}} \\ {{\omega _3}}&0&{ - {\omega _1}} \\ { - {\omega _2}}&{{\omega _{\rm{1}}}}&0 \end{array}} \right] $$ (3) -
在实际深空探测任务中,深空探测器的执行机构提供的力矩幅值是有限的,控制输入有界约束为
$$ \left| {{u_i}} \right| \leqslant {\gamma _1} \;\;\; ,i = 1,2,3 $$ (4) 其中:
${\gamma _1} $ 表示系统控制力矩的幅值。系统敏感器的量程有限,要求深空探测器的角速度必须保持在某个范围内,角速度有界约束为
$$ \left| {{\omega _i}} \right| \leqslant {\gamma _2} \;\;\; ,i = 1,2,3 $$ (5) 其中:
${\gamma _2} $ 表示系统角速度幅值。 -
深空探测器在执行空间任务时,会面临较为复杂的指向约束,这些指向约束极大程度地缩小了姿态机动的可行空间。一旦违反指向约束,将会对深空探测器携带的光学载荷造成严重的损坏,进而影响对应空间执行的任务,所以要分析和处理深空探测器的指向约束。
深空探测器在姿态机动过程中必须避免强光天体进入敏感光学元件的视场角中,这种约束称为“禁止指向约束”。图1为深空探测器指向约束关系图,
${{ r}_{g1}}$ 表示该光敏元件在本体下的方向矢量,${{ r}_I}$ 表示某一强光天体在惯性系下的方向矢量。姿态机动过程中要求强光天体方向矢量和光敏元件的视线方向夹角不能低于视场角,约束表示形式为
$$ {{r}}_{{{gi}}}^{\bf T}\left( {{{{C}}_{{{gI}}}}{{{r}}_{{i}}}} \right) \leqslant \cos {\theta _i} $$ (6) 将指向约束的姿态余弦矩阵展开,将指向约束表示为更加简洁的二次型形式
$$ {{{q}}^{\rm{T}}}{{{K}}_{\rm f}}{{q}} \leqslant 0 $$ (7) 其中
$$ {{{K}}_{\rm f}} = \left[ {\begin{array}{*{20}{c}} {{{r}}_i^{\rm{T}}{{{r}}_{g1}} - \cos \theta }&{{{\left( {{{{r}}_{g1}}^ \times {{r}}_i^{\rm{T}}} \right)}^{\rm{T}}}} \\ {{{{r}}_{g1}}^ \times {{r}}_i^{\rm{T}}}&{{{{r}}_i}{{{r}}_{g1}}^{\rm{T}} + {{{r}}_{g1}}{{r}}_i^{\rm{T}} - \left( {{{r}}_i^{\rm{T}}{{{r}}_{g1}} + \cos \theta } \right){{{I}}_3}} \end{array}} \right] $$ (8) -
姿态机动多目标优化不仅要满足上述复杂约束条件,还需要重点考虑时间和能量参数作为优化性能指标。这里将多约束下深空探测器姿态机动多目标规划问题归纳为以下形式,具体比例参数由第2节给出。
$$ \begin{array}{l} \min:J = p(x) \\ {\rm{s}}{\rm{.t}}{\rm{.}}\quad \\ \left\{ \begin{array}{l} {{J\dot \omega }} = {{u}} - {{{\omega }}^ \times }{{J\omega }} \\ {\dot{ q}} = \dfrac{1}{2}{{Q\omega }} = \dfrac{1}{2}{{\varOmega }}{{q}} \\ - {\gamma _1} \leqslant {u_i} \leqslant {\gamma _1} , i = 1,2,3 \\ - {\gamma _2} \leqslant {\omega _i} \leqslant {\gamma _2} , i = 1,2,3 \\ {{{q}}^{\rm{T}}}{{{K}}_{\rm f}}{{q}} \leqslant 0 \\ {{{q}}^{\rm{T}}}{{q}} = 1 \\ \end{array} \right. \end{array} $$ (9) 其中:
$p(x)$ 表示系统的多目标优化函数。 -
这里采用物理规划方法,将探测器姿态机动的多目标优化问题转化为单目标优化问题。通过分析多目标之间的关系,得到一个在Pareto最优解前沿附近的解[19]。
通过物理规划将不同性质的多个指标融合成一个分段指标,在分段区域根据指标优先级设置不同的权重,既能够根据指标的特性进行针对性的优化,又能够对多目标优化进行解耦,简化路径规划过程和提高规划效率。
$$ \min p(x) = \frac{1}{{{n_{\rm sc}}}}\sum\limits_{i = 1}^{{n_{\rm sc}}} {{p_i}({J_i}(x))} $$ (10) 其中:
${p_i}$ 是分解函数;${n_{\rm sc}}$ 是目标函数的个数。 分解函数被划分成多个区间,在各个区间之中,通过分段曲线拟合能够得到各个性能指标函数的分解函数,用分段样条函数表示。第一段和剩下的多段表示如下$$ {p_{i1}} = {x_{i1}}\exp \Bigg(\frac{{{s_i}}}{{{x_{i1}}}}({x_i} - {x_{i1}})\Bigg),{x_i} < {x_{i1}} $$ (11) $$ \begin{aligned}[b] {p_{ik}} = & {c_1}(r_i^k){x_{i(k - 1)}} + {c_2}(r_i^k){x_{ik}}+ \\ & {c_{11}}(r_i^k,\lambda _i^k){s_{i(k - 1)}} + {c_{12}}(r_i^k,\lambda _i^k){s_{ik}} \end{aligned} $$ (12) 各分解函数参数的选取主要根据各个性能指标的区间满意度[20],再配合多性能指标适应度函数,就得到目标优化函数
$p(x)$ 。 -
使用目标优化函数代替差分进化算法的总体适应度函数,用于改进差分进化算法的选择标准。目标优化函数由分解函数
${p_i}$ 和性能指标适应度函数${J_i}(x)$ 组成。差分进化算法的性能指标适应度函数设置为越小越好型,通过寻找目标优化函数的最小值来变异、交叉和选择得到最优姿态机动路径。性能指标适应度函数除了探测器的性能指标机动时间和能量之外,也要考虑探测器遇到的路径指向约束、有界约束和终端约束。
首先考虑深空探测器姿态机动过程中的时间和能量优化指标参数,转换为多目标优化的时间和能量适应度函数为
$$ {J_t}{\rm{ = }}{t_{\rm f}} $$ (13) $$ {J_{\rm e}}{\rm{ = }}{u^2}\Delta t $$ (14) 姿态路径指向约束的适应度函数针对禁忌约束
$${J_{\rm g}} = \left\{ {\begin{array}{*{20}{l}} {\exp \left( {{{q}}^{\rm T}(k){{K}}{{q}}(k)} \right) - 1}, &{ \text{如果} {{q}}^{\rm T}(k){{K}}{{q}}(k) > 0} \\ 0&{\text{否则}} \end{array}} \right.$$ (15) 探测器有界约束的适应度函数取为
$${J_{\rm u}} = \left\{ {\begin{array}{*{20}{l}} {\exp \left( {\left| {{u_i}\left( k \right)} \right| - {\gamma _1}} \right) - 1}, &{{\text{如果}} \; \left| {{u_i}\left( k \right)} \right| > {\gamma _1}} \\ 0&{\text{否则}} \end{array}} \right.$$ (16) 为了确保算法的准确性与收敛性,将终端约束扩展到所有节点的目标导向约束,这时终端约束的适应度函数为
$$ \begin{aligned}[b] {J_{qi}} =& - {\Bigg(\dfrac{{{k_i}}}{{{k_{\max }}}}\Bigg)^p} \lg \left( {\left| {\dfrac{{{q_{ei}} + 1}}{2}} \right|} \right) \\ {J_q} =& \displaystyle\sum\limits_{k = 1}^{{k_{\max }}} {{J_{qi}}} \end{aligned}$$ (17) 其中:
${{{q}}_{\rm e}}{\rm{ = }}{\left[ {\begin{array}{*{20}{c}} {{q_{{\rm e}0}}}&{{q_{{\rm e}1}}}&{{q_{{\rm e}2}}}&{{q_{{\rm e}3}}} \end{array}} \right]^{\rm{T}}}$ ;${{k}_i}$ 是当前节点数;${k_{\max }}$ 是总节点数,这样即能够保证探测器姿态机动到设定的目标点,同时由于终端约束适应度函数的导向作用,能够大幅度提高差分进化算法后期的收敛速率。得到所有性能指标适应度函数,与分解函数组合得到目标优化函数。
-
为了加快算法收敛,本文将首先生成初始机动路径,并与改进差分进化算法组合,共同生成多目标优化姿态机动路径。
本文通过快速欧拉旋转方法得到从初始到目标的机动路径,并不考虑指向和有界约束,得出转动四元数
$${q_i}{\rm{ = }}\overline {{q_0}} \otimes {q_{\rm f}}$$ (18) 根据欧拉定理,得出欧拉转角
$$\theta {\rm{ = 2}}\operatorname{arcos} {q_i}(1)$$ (19) 再根据初始姿态参数,进行逐步姿态路径节点计算,得到整个姿态路径上的姿态四元数,即得到初始的姿态机动路径,用以引导优化过程的中间姿态机动路径。
-
1)初始化
利用n个向量
${{P}}$ 作为每一代的种群,每个个体向量${{{P}}_{i,G}}$ 表示为$${{{P}}_{i,G}} = (\begin{array}{*{20}{l}} {\omega _1^{\rm T}}&{{t_1}}&{\omega _2^{\rm T}}&{{t_2}}&{ \cdot \cdot \cdot }&{\omega _n^{\rm T}}&{{t_n}} \end{array})$$ (20) 其中:
${{\omega }}$ 为当前节点的姿态角速度;$t$ 为节点的机动时间;i表示个体在种群中的序号;G表示进化代数。初始化的过程主要是生成算法优化的初始种群。通常为了避免算法后期的收敛效率过低,常用的生成初始种群的方法是基于给定约束值进行边界内部随机选择。
在整个优化算法中,假定所有随机初始种群都符合均匀概率分布。设参数变量的界限为
${{P}}_j^{(L)}\! < \!{{{P}}_j}\! < \! {{P}}_j^{(U)}$ ,${{{P}}_j}$ 表示${{P}}$ 的第j维向量,j=1,3,···,D,因此初始化的种群为$$ {{{P}}_{i,0}} = {\rm rand(0,1)} \cdot (P_j^{(U)} - P_j^{(L)}) + P_j^{(L)} $$ (21) 2)组合变异
对于单个个体向量
${{{P}}_{i,G}}$ ,采用组合变异的方式。由欧拉旋转得到的初始机动路径为${{{P}}_{\rm c}}$ ,以设定的概率选择是否将${{{P}}_{\rm c}}$ 放入差分进化的变异过程中,变异个体的产生如式(22)所示$$ {{{v}}_{i,G + 1}} = \left\{ {\begin{array}{*{20}{l}} {{{{P}}_{r1,G}} + F \cdot ({{{P}}_{r2,G}} - {{{P}}_{r3,G}}),} \\ {{{{P}}_c} + F \cdot ({{{P}}_{r2,G}} - {{{P}}_{r3,G}}),} \end{array}} \right. \begin{array}{*{20}{l}} { \text{如果} (ra < f)} \\ {\text{如果}(ra \geqslant f)} \end{array}\!\!\! $$ (22) 其中:
$ra$ 为0~1之间的随机数;$f$ 为相应范围内的判断概率,判断是否将初始路径${{{P}}_{\rm c}}$ 加入变异过程中;${{{P}}_c}$ 是欧拉旋转得到的最短距离姿态机动路径,引导当前的中间姿态机动路径向最短距离姿态机动路径变异,加快算法收敛速度;突变算子$F$ 的取值显著影响算法的搜索范围,$F$ 值越大,算法搜索的范围越广,全局搜索能力强,也能够避免过早收敛。相反,一个较小的值$F$ 可以提高算法的局部搜索能力,加快局部搜索速度和算法的收敛速度。一般的$F$ 值固定选取为1,这里将$F$ 值设置为根据算法进程可变的形式,搜索前期较大,搜索后期较小,使算法同时具有较好的前期全局搜索能力和后期收敛速度$$ F = c{5^{ - {{(\frac{{{k_i}}}{{{k_{\max }}}})}^p}}} $$ (23) 其中:
${k_i}$ 是当前节点数;${k_{\max }}$ 是总节点数;c5和p为控制参数。3)交叉
作为变异策略的补充,利用交叉策略在变异向量的基础上保留原始目标个体的信息,并生成一个测试向量与原始目标个体竞争。交叉向量可以表示为
$$ {{{Y}}_{ij,G + 1}} = \left\{ {\begin{array}{*{20}{l}} {{{{v}}_{ij,G + 1}} , \; {\rm{if}} \; rb \leqslant CR || j = rn} \\ {{X_{ij,G + 1}},\;{\rm{else}}} \end{array}} \right. $$ (24) 其中:CR为交叉概率因子,是0~1的正实数;rb是0~1的随机数;rn为[1,2,···,D]的随机整数。
4)选择
按照目标优化函数
$p(x)$ 最小准则将当前最优个体与当前种群中的个体向量进行比较,决定试验向量是否会成为下一代中的成员。具有目标优化函数的向量将保存在下一代种群中。逐代更新,每一代都优于或等于前一代。最后,将基于多约束下多目标姿态机动问题转化成寻找种群中的最优向量、使得目标优化函数
$p(x)$ 最小的优化问题,算法逐步迭代直到收敛到最优解,进而得到优化后的深空探测器姿态机动路径。 -
在仿真验证中,深空探测器在其z轴方向安装了一个光学相机,方向矢量由
${{ r}_I}$ 表示,4个需要探测器姿态机动过程中躲避的发光天体在惯性系下的指向分别为用${{{r}}_{I1}}$ 、${{{r}}_{I2}}$ 、${{{r}}_{I3}}$ 和${{{r}}_{I4}}$ 。最小视场角分别为${\theta _1}$ 、${\theta _2}$ 、${\theta _3}$ 和${\theta _4}$ ,探测器的初始姿态和初始角速度分别为${{{q}}_0}$ 和${{{\omega}} _0}$ ,目标姿态和目标角速度分别为${{{q}}_{\rm{f}}}$ 和${{{\omega }}_{\rm{f}}}$ 。如表1所示,探测器的姿态机动角速度幅值为${{{\gamma}} _\omega }$ ,姿态控制力矩幅值为${{{\gamma}} _T}$ ,探测器的转动惯量为${{J}}$ ,具体数值见表1。表 1仿真条件
Table 1.Simulation conditions
变量 值 J/(kg·m2) diag(100,100,100) q0 [0.646,0.034,0.722,0.241]T ω0/(rad•s–1) [0,0,0] qf [0.733,0.362,–0.544,0.181]T ωf/(rad•s–1) [0,0,0] ${ {{\gamma} } _\omega }$ 0.05 $ { {{\gamma} } _u} / {\rm Nm}$ 0.1 θ1/(°) 40 θ2/(°) 30 θ3/(°) 30° θ4/(°) 25 rgI [0,0,1]T rI1 [–0.76,0,0.64]T rI2 [0.49,0.85,0.17]T rI3 [0.64,–0.76,0]T rI4 [0.76,0,0.64]T 算法的仿真环境为Matlab 2019b,计算机主频为3.2 Ghz,内存4 G,采用本文提出的多目标组合规划算法进行姿态机动路径求解。
深空探测器姿态机动过程中的姿态四元数、控制力矩和角速度分别如图2~4所示。可以看出,探测器的姿态机动路径满足有界约束,没有超出探测器的角速度和控制力矩上限,同时角速度和控制力矩路径的前后转变平稳,比较方便实际工程实现。
经过组合规划算法得到探测器在天球坐标系下的姿态机动路径,如图5所示,环形区域表示姿态指向约束,红色连线表示姿态机动路径。可以看出,深空探测器运动从起始点成功姿态机动到目标点,并且没有进入约束区域,说明该算法能够得到安全的姿态机动路径。
图 5探测器在天球坐标系下姿态机动路径
Figure 5.Attitude maneuver path of deep space probe in celestial coordinate system
同时,探测器在天球坐标系下可见光相机的空间姿态机动路径如图6所示,其中,圆形区域表示的是姿态指向约束,红线表示空间姿态机动路径。由此可以看出,在探测器机动过程中光学相机成功地规避了4个发光天体,满足了姿态指向约束,得到了安全的姿态机动路径。从上述所有结果可以看出,本文提出的组合规划算法可以满足多种姿态约束,证明了算法的有效性。
多目标组合规划算法的平均运行时长为36 s,姿态机动时长为200 s,能量消耗的相对单位为0.8。基于相同的仿真环境,对比基于时间优化的差分进化算法,时间优化后的标准差分进化算法的平均运行时长为60 s,姿态机动时长为220 s,能量消耗的相对单位为1.5。可以看出,多目标组合规划算法同时提高了时间和能量的优化效率以及程序运行速度,具有较好的时间和能量优化效果。
Combination Planning for Attitude Maneuver of Deep Space Probes Based on Multi-Objective Optimization
-
摘要:深空探测器在执行姿态机动任务中,通常需要优化多个性能指标参数。此外,姿态指向约束和有界约束的存在显著减小了姿态机动路径的可行空间。在复杂的多约束条件下,多目标姿态机动优化对深空探测器系统是一个极大的挑战。提出多目标组合的概念,通过物理规划方法将多目标优化转换为单目标规划问题。考虑到欧拉路径的低能量特性,提出一种多目标组合规划方法,通过欧拉路径改进差分进化算法的变异过程,实现多约束多目标姿态机动路径组合规划。最后,对深空探测器大角度姿态机动进行多目标组合规划数值仿真,验证了该方法的可行性和有效性。Abstract:For deep space detectors attitude maneuver tasks, it is necessary to optimize several performance index parameters. In this paper, the multi-objective planning problem of attitude maneuver for deep space detectors is studied, and a multi- target combination planning method is proposed to calculate the attitude maneuver path. For multi-objective programming, the attitude objective function is designed through physical planning method, which transforms multi-objective optimization into a single objective planning problem, and at the same time replaces fitness function as the population selection criterion. The combined differential evolution algorithm is designed to improve the mutation process. The initial attitude path obtained by fast Euler maneuver is added to the mutation formula with a certain probability, so as to improve the convergence rate and program operation efficiency, and realize the multi constraint and multi-objective attitude maneuver path combination planning. Finally, the feasibility and effectiveness of this method are verified by numerical simulation of large angle attitude maneuver of deep space probes.Highlights
● Multi-objective optimization of time and energy for attitude maneuver of deep space probe is studied. ● The multi-objective optimization is transformed into single objective programming to replace fitness function as the population selection criterion of differential evolution. ● The initial attitude path obtained by the fast Euler maneuver is added to the variation formula with the set probability to improve the variation process of differential evolution. ● The combined differential evolution algorithm is designed to improve the convergence rate and the efficiency of the program. -
表 1仿真条件
Table 1Simulation conditions
变量 值 J/(kg·m2) diag(100,100,100) q0 [0.646,0.034,0.722,0.241]T ω0/(rad•s–1) [0,0,0] qf [0.733,0.362,–0.544,0.181]T ωf/(rad•s–1) [0,0,0] ${ {{\gamma} } _\omega }$ 0.05 $ { {{\gamma} } _u} / {\rm Nm}$ 0.1 θ1/(°) 40 θ2/(°) 30 θ3/(°) 30° θ4/(°) 25 rgI [0,0,1]T rI1 [–0.76,0,0.64]T rI2 [0.49,0.85,0.17]T rI3 [0.64,–0.76,0]T rI4 [0.76,0,0.64]T -
[1] ZHOU H,WANG D W,WU B L,et al. Time-optimal reorientation for rigid satellite with reaction wheels[J]. International Journal of Control,2012,85(10):1452-1463.doi:10.1080/00207179.2012.688873 [2] WANG S,ZHAO L,CHENG J,et al. Task scheduling and attitude planning for agile earth observation satellite with intensive tasks[J]. Aerospace Science and Technology,2019,90(6):23-33. [3] SHE Y,LI S,ZHAO Y. Onboard mission planning for agile satellite using modified mixed-integer linear programming[J]. Aerospace Science and Technology,2018,72(11):204-216. [4] 程小军,崔祜涛,徐瑞,等. 几何约束下的航天器姿态机动控制[J]. 控制与决策,2012,27(5):724-730.CHENG X J,CUI H T,XU R,et al. Attitude maneuver control of spacecraft under geometric constraints[J]. Control and Decision,2012,27(5):724-730. [5] MCINNES C R. Large-angle slew maneuvers with autonomous Sun vector avoidance[J]. Journal of Guidance Control and Dynamics,1994,17(4):875-877.doi:10.2514/3.21283 [6] SINGH G, MACALA G, WONG C E, et al. A constraint monitor algorithm for the Cassini spacecraft[C]//In Proceedings of the AIAA Guidance, Navigation, and Control Conference. [S.l.]: AIAA, 1997. [7] 仲维国,崔平远,崔祜涛. 航天器复杂约束姿态机动的自主规划[J]. 航空学报,2007,25(5):1091-1097.doi:10.3321/j.issn:1000-6893.2007.05.011ZHONG W G,CUI P Y,CUI H T. Autonomous attitude maneuver planning for spacecraft under complex constraints[J]. Acta Aeronautica et Astronautica Sinica,2007,25(5):1091-1097.doi:10.3321/j.issn:1000-6893.2007.05.011 [8] KJELLBERG H C,LIGHTSEY E G. Discretized constrained attitude pathfinding and control for satellites[J]. Journal of Guidance,Control,and Dynamics,2013,36(5):1301-1309. [9] LAI L C,YANG C C,WU C J. Time-optimal maneuvering control of a rigid spacecraft[J]. Acta Astronautica,2007,60(10):791-800. [10] BOYARKO G A,ROMANO M,YAKIMENKO O A. Time-optimal reorientation of a spacecraft using a direct optimization method based on inverse dynamics[J]. Journal of Guidance,Control,and Dynamics,2011,34(4):1197-1208. [11] PANAGOU D. A distributed feedback motion planning protocol for multiple unicycle agents of different classes[J]. IEEE Transactions on Automatic Control,2017,62(3):1178-1193.doi:10.1109/TAC.2016.2576020 [12] SEKHAVAT P,YAN H,FLEMING A,et al. Closed-loop time-optimal attitude maneuvering of magnetically actuated spacecraft[J]. The Journal of the Astronautical Sciences,2011,58(1):81-97.doi:10.1007/BF03321160 [13] KIM Y,MESBAHI M,SINGH G,et al. On the convex parameterization of constrained spacecraft reorientation[J]. Aerospace & Electronic Systems IEEE Transactions,2010,46(3):1097-1109. [14] SPILLER D,MELTON R G,CURTI F. Inverse dynamics particle swarm optimization applied to constrained minimum-time maneuvers using reaction wheels[J]. Aerospace Science & Technology,2018,75(4):1-12. [15] WU C,XU R,ZHU S,et al. Time-optimal spacecraft attitude maneuver path planning under boundary and pointing constraints[J]. Acta Astronautica,2017,137(8):128-127. [16] FAHROO F,ROSS I M. Direct trajectory optimization by a chebyshev pseudospectral method[J]. Journal of Guidance,Control,and Dynamics,2002,25(1):160-166. [17] XIAO G,ZHU M. Direct trajectory optimization based on a mapped Chebyshev pseudospectral method[J]. Chinese Journal of Aeronautics,2013(2):401-412. [18] CRAIN A,ULRICH S. Experimental validation of pseudo spectral based optimal trajectory planning for free-floating robots[J]. Journal of Guidance Control and Dynamics,2019,42(8):1-17. [19] MAGINOT J, GUENOV M, FANTINI P. A method for assisting the study of Pareto solutions in multi-objective optimization[C]//7th AIAA ATIO Conference. Belfast, Northern Ireland: AIAA, 2007. [20] ZHAO Q,HUANG H. Multi-objective optimization of zero propellant maneuver using hybrid programming[J]. Acta Astronautica,2015,116(12):154-160.