下载:
-
以电推进为代表的小推力推进系统近年来受到了广泛的关注，相较于传统化学推进，电推进具有更高的推进效率，即更能节省燃料[1]。自“深空1号”（Deep Space 1）探测器成功验证电推进技术后，“黎明号”（Dawn）、“隼鸟号”（Hayabusa）和正在进行的“贝皮科伦布”探测器（BepiColombo）等深空任务均采用了电推进发动机。而在近地任务中，波音公司的Boeing 902SP和中国的“实践十七”“实践二十”将电推进技术应用于地球同步轨道的入轨和轨道保持任务[2]。因此，电推进技术由于其出色的推进效率，在未来航天任务设计中将逐渐获得更普遍的应用。
在轨道设计方面，由于电推进所产生的推力较小，需要长时间的开机工作才能产生足够的速度增量，以电推进驱动的航天器小推力最优轨迹将和脉冲推进的最优轨迹产生显著的差别。圆锥曲线拼接[1]、基于形状的标称轨道方法[3]等相对解析的设计方法难以实现对小推力最优轨迹的准确近似。而对于一般的深空探测轨道设计，在连续推力作用下的轨道无法采用类似解析的方法得到最优解，而是采用较为耗时的数值优化的方法（主要为直接法和间接法）[4]。
对于探测任务尤其是多目标任务的全局规划问题[5-6]，需要快速且准确地确定小推力轨迹的燃料消耗、转移时间等参数。因此，实现对小推力轨迹的快速优化具有重要意义。
目前，主要的数值优化方法分为直接法和间接法[7]。
直接法将小推力优化问题进行离散，转化为非线性规划问题求解，通常需要较多的离散点以提高解的精度和最优性[8]。相较于直接法，间接法利用极大值原理构建两点边值问题求解，能够保证最优性的一阶必要条件[9]。但是这些数值方法均需要提供初值进行优化，直接法需要提供离散点的状态和控制初值，比较容易利用设计经验或标称轨道方法得到，而且直接法的收敛域相对较大，其求解对初值敏感性较弱[8]，因此其计算性能主要取决于离散点数。间接法需要提供协态的初值，相较于具有实际物理意义较容易猜测的状态，协态并无明确的物理意义，其取值范围也难以确定，因此较难提供初值[1]。同伦方法[9]将小推力燃料最优问题通过同伦参数联系转换为求解相对较容易的能量最优问题求解，提高了初值猜测成功的概率。协态归一化方法[10]通过引入归一化参数将协态的猜测范围限制在高维单位球面上，提供了初值猜测的范围，进一步提高了初值猜测的效率。基于标称轨道的协态初值估计方法[1]能够提供能量最优问题的初值猜测，但是其猜测的成功率依赖于标称轨道的设计，对于较复杂的任务，标称轨道与最优轨迹偏离较大，协态估计的效果较差。因此针对小推力间接法轨迹优化，解析的保证收敛的协态初值能够显著降低求解失败的概率，以提高小推力轨迹优化的效率。
本研究提出了一种解析的协态初始化方法，引入切换系统，将小推力燃料最优问题同一个具有解析初值的轨迹优化问题相联系，通过迭代方法求解得到小推力燃料最优轨迹。本文首先将小推力优化模型嵌入到更为广义的切换系统模型中，给出了切换系统的构建形式，并设计了相应的切换函数，从而将小推力优化问题转化为更简单优化模型的优化问题；其次，给出了具有解析协态初值的简单系统优化模型，能够通过切换系统与小推力优化模型相联系，进一步给出了计算迭代算法，并通过数值仿真验证了该方法相较于传统同伦方法，具有更高的求解效率。
-
针对常见的小推力轨迹优化问题，本节首先给出了一般指标形式下的小推力轨迹优化模型，其中动力学模型也采用了一般的矩阵形式描述。在此优化模型的基础上，进一步介绍了切换系统的概念，从而将小推力优化模型嵌入进更广义的切换系统优化模型中，同时给出了切换函数设计。
-
一般而言，考虑如下小推力轨迹优化模型，优化指标为
$$ {J_1} = \int_{{t_0}}^{{t_f}} {{L_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right)\;{\rm{d}}t} $$ (1) 其中：指标J1通常取为总的燃料消耗[9]、能量消耗[1]、飞行时间[8]或混合指标；
${\boldsymbol{\xi}} $ 表示状态矢量，通常由航天器的位置、速度（或轨道要素）和质量组成，即${\boldsymbol{\xi}} = {\left[ {x,m} \right]^{\rm{T}}}$ ；${{\boldsymbol{u}}_1}$ 表示控制量，通常包括控制量幅值${\tau _1} \in \left[ {0,1} \right]$ 及其方向${\hat {\boldsymbol{\alpha}} _1}$ 。动力学方程为$$ \dot {\boldsymbol{x}} = {\boldsymbol{A}}\left( {{\boldsymbol{x}},t} \right) + \frac{{{T_{\max }}}}{m}{\boldsymbol{B}}\left( {{\boldsymbol{x}},t} \right){\tau _1}{\hat {\boldsymbol{\alpha}} _1},\;\;\;\;{\boldsymbol{x}}\left( {{t_0}} \right) = {{\boldsymbol{x}}_0} $$ (2) $$ \dot m =- \frac{{{T_{\max }}}}{{{I_{{\rm{sp}}}}{g_0}}}{\tau _1},\;\;\;\;m\left( {{t_0}} \right) = {m_0} $$ (3) 其中：假设小推力推进系统具有恒定的最大推力
${T_{\max }}$ 和恒定的比冲${I_{{\rm{sp}}}}$ ；矩阵${\boldsymbol{A}}$ 和${\boldsymbol{B}}$ 分别为$6 \times 1$ 和$6 \times 3$ 的矩阵，其具体形式依赖于坐标系的选取，在本文中采用春分点轨道根数形式描述动力学方程；${g_0}$ 表示海平面重力加速度常量。此优化模型可以转化为两点边值问题并采用间接法求解。相比于直接法，间接法能够保证最优的一阶必要条件，因此一旦求解成功即可得到局部最优的轨道。但另一方面，在间接法的求解过程中需要人为提供无物理意义的协态初值猜测作为求解的初始化，而且间接法的收敛域也较直接法更小，难以对估计不准的协态初值实现成功求解。因此，优化时的初值猜测情况会显著影响间接法的求解成功率和效率，当前主要采用随机猜测的方式多次猜测求解，猜测成功率较低且求解效率较慢。
基于协态初值估计的方法，本文将小推力优化模型嵌入进切换系统中，同时设计切换系统形式使其具有解析的协态初值，避免了初值猜测困难。以此解析协态初值作为初始化，通过切换函数迭代实现了小推力最优轨迹的快速求解。
-
基于小推力轨迹优化模型，引入如下一般形式的切换系统[11]优化模型。优化指标为
$$ {J_1} = \int\limits_{{t_0}}^{{t_f}} {\left\{ {\left[ {1 - \upsilon \left( {t,\varepsilon } \right)} \right]{L_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right) + \upsilon \left( {t,\varepsilon } \right){L_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right)} \right\}{\rm{d}}t} $$ (4) 动力学方程为
$$ \dot {\boldsymbol{\xi}} = \left[ {1 - \upsilon \left( {t,\varepsilon } \right)} \right]{{\boldsymbol{f}}_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right) + \upsilon \left( {t,\varepsilon } \right){{\boldsymbol{f}}_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right) $$ (5) 其中：
$\upsilon \left( {t,\varepsilon } \right) \in \left[ {0,1} \right]$ 为切换函数，通过此函数即可实现系统优化模型的切换，$\varepsilon \in \left[ {0,1} \right]$ 为同伦参数以控制切换函数的具体形式。一般而言，传统的从能量最优问题同伦求解燃料最优问题的算法，将同伦参数添加到优化指标中，而动力学方程保持不变。此时，首先求解能量最优问题的协态初值，然后逐渐改变同伦参数值，以得到的协态初值作为猜测值，求解对应同伦参数下的优化问题，最终得到燃料最优解。
与传统方法不同，本文的优化模型将同伦参数引入到两个系统之间的切换函数中，系统1的优化指标为
${L_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right)$ ，动力学方程的右函数为${{\boldsymbol{f}}_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right)$ 。系统2的优化指标为${L_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right)$ ，动力学方程的右函数为${{\boldsymbol{f}}_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right)$ ，采用同样的同伦过程即可得到燃料最优解。在本文中，系统1设计为小推力燃料最优轨迹优化模型，系统2及其解析协态初值的设计将在下一节中详细给出，切换函数设计为$$ \upsilon \left( {t,\varepsilon } \right) \triangleq \left\{ \begin{array}{l} 0,\;\;t \leqslant {t_f} - \varepsilon \left( {{t_f} - {t_0}} \right) \\ {\rm{1}},\;\;t > {t_f} - \varepsilon \left( {{t_f} - {t_0}} \right) \end{array} \right. $$ (6) 接下来，首先推导此切换系统的最优控制律。引入协态变量
${{\boldsymbol{\lambda}} _{{\xi}} }$ ，系统的哈密顿函数为$$ \begin{aligned}[b] H = & \left[ {1 - \upsilon \left( {t,\varepsilon } \right)} \right]\left[ {{{\boldsymbol{\lambda}} _{{\xi}} } \cdot {{\boldsymbol{f}}_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right) + {L_1}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_1}} \right)} \right] + \\ & \upsilon \left( {t,\varepsilon } \right)\left[ {{{\boldsymbol{\lambda}} _{{\xi}} } \cdot {{\boldsymbol{f}}_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right) + {L_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right)} \right] \end{aligned} $$ (7) 根据极大值原理，系统的最优控制使哈密顿函数取极值。对于确定的同伦参数
$\varepsilon $ ，系数$1 - \upsilon \left( {t,\varepsilon } \right)$ 和$\upsilon \left( {t,\varepsilon } \right)$ 均为大于等于0的常数，因此哈密顿函数取极值即为函数中的第一项和第二项分别取极值。这意味着控制${{\boldsymbol{u}}_1}$ 和${{\boldsymbol{u}}_2}$ 是解耦的，相互无关。确定系统1和系统2的形式后，便可得到最优控制${{\boldsymbol{u}}_1}^ * $ 和${{\boldsymbol{u}}_2}^ * $ 的形式。 -
基于标称轨道线性化的方法，本节首先给出了系统2的模型，并进一步给出了确定小推力燃料最优轨迹的求解算法。
-
系统2考虑如下线性优化模型。优化指标为
$$ {J_2} = \int_{{t_0}}^{{t_f}} {{L_2}\left( {t,{\boldsymbol{\xi}} ,{{\boldsymbol{u}}_2}} \right)\;{\rm{d}}t} = \beta \int_{{t_0}}^{{t_f}} {\tau _2^2\;{\rm{d}}t} $$ (8) 其中：
$\beta \in {{\mathbb{R}}^ + }$ 为协态变量归一化参数，或简单取为和燃料最优问题的比例系数一致即可。动力学方程为$$ \dot {\boldsymbol{x}} = {\boldsymbol{C}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right){\boldsymbol{x}} + \frac{{{T_{\max }}}}{{{m_0}}}{\boldsymbol{B}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right){\tau _{\rm{2}}}{\hat {\boldsymbol{\alpha}} _{\rm{2}}}{\rm{ + }}{\boldsymbol{D}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right),\;\;{\boldsymbol{x}}\left( {{t_0}} \right) = {{\boldsymbol{x}}_0} $$ (9) 其中：
${{\boldsymbol{x}}_{\rm{n}}}$ 表示选取的标称轨道；航天器质量为常数${m_0}$ ，即质量变化率恒为0。并且在此模型中，控制${{\boldsymbol{u}}_2}$ 的幅值${\tau _2}$ 无约束。矩阵${\boldsymbol{C}}$ 和${\boldsymbol{D}}$ 可以通过矩阵${\boldsymbol{A}}$ 和标称轨道给出$$ {\boldsymbol{C}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right) \triangleq {\left. {\frac{{\partial {\boldsymbol{A}}\left( {{\boldsymbol{x}},t} \right)}}{{\partial {\boldsymbol{x}}}}} \right|_{{\boldsymbol{x}} = {{\boldsymbol{x}}_{\rm{n}}}}} $$ (10) $$ {\boldsymbol{D}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right) \triangleq {\boldsymbol{A}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right) - {\boldsymbol{C}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right){{\boldsymbol{x}}_{\rm{n}}} $$ (11) 另外，对于标称轨道的选取，在本文中仅考虑的简单线性的标称轨道，具体形式为
$$ {{\boldsymbol{x}}_{\rm{n}}}\left( t \right) = {{\boldsymbol{x}}_0} + \frac{{{{\boldsymbol{x}}_f} - {{\boldsymbol{x}}_0}}}{{{t_f} - {t_0}}}\left( {t - {t_0}} \right) $$ (12) 由于系统2设计即为简单线性的优化模型，此种标称轨道能够进一步简化解析解求解的难度，快速给出解析协态初值。本文的数值仿真证明此标称轨道能够处理较为复杂的燃料最优问题。
对于系统2，同理可以引入协态变量
${{\boldsymbol{\lambda}} _x} \in {{\mathbb{R}}^6}$ 。则系统2的哈密顿函数为$$ {H_2} = {{\boldsymbol{\lambda}} _x} \cdot \left[ {{\boldsymbol{C}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right){\boldsymbol{x}} + \frac{{{T_{\max }}}}{{{m_0}}}{\boldsymbol{B}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right){\tau _{\rm{2}}}{{\hat {\boldsymbol{\alpha}} }_{\rm{2}}}{\rm{ + }}{\boldsymbol{D}}\left( {{{\boldsymbol{x}}_{\rm{n}}},t} \right)} \right] + \beta \tau _2^2 $$ (13) 根据极大值原理可得最优控制律为
$$ {\left( {{\tau _2}{{\hat {\boldsymbol{\alpha}} }_2}} \right)^ * } = - \frac{{{T_{\max }}}}{{2{m_0}\beta }}{{\boldsymbol{B}}^{\rm{T}}}{{\boldsymbol{\lambda}} _x} $$ (14) 此控制律即为最小化
${H_2}$ 所得。利用式（12）对状态变量${\boldsymbol{x}}$ 求偏导，并将式（13）代入式（9）可得状态微分方程和协态微分方程如下$$ \left[ \begin{array}{l} {\dot {\boldsymbol{x}}} \\ {{\dot {\boldsymbol{\lambda}} }_x} \end{array} \right] = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{C}}&{\dfrac{{T_{\max }^{\rm{2}}}}{{{\rm{2}}m_0^2\beta }}{\boldsymbol{B}}{{\boldsymbol{B}}^{\rm{T}}}} \\ 0&{ - {{\boldsymbol{C}}^{\rm{T}}}} \end{array}} \right]\left[ \begin{array}{l} {\boldsymbol{x}} \\ {{\boldsymbol{\lambda}} _x} \end{array} \right]{\rm{ + }}\left[ \begin{array}{l} {\boldsymbol{D}} \\ 0 \end{array} \right] $$ (15) 因此，得到了一个线性时变的动力学微分方程，此方程中系数仅和标称轨道以及当前时刻有关，而和状态、协态无关。此方程的解可以解析得到
$$ \left[ \begin{array}{l} {\boldsymbol{x}}\left( {{t_f}} \right) \\ {{\boldsymbol{\lambda}} _x}\left( {{t_f}} \right) \end{array} \right] = {\boldsymbol{\varPhi}} \left( {{t_f},{t_0}} \right)\left[ \begin{array}{l} {\boldsymbol{x}}\left( {{t_0}} \right) \\ {{\boldsymbol{\lambda}} _x}\left( {{t_0}} \right) \end{array} \right] + \left[ \begin{array}{l} {\boldsymbol{z}}\left( {{t_f},{t_0}} \right) \\ 0 \end{array} \right] $$ (16) 其中：状态转移矩阵
${\boldsymbol{\varPhi}} \left( {{t_f},{t_0}} \right) = \left[ {{{\boldsymbol{\varPhi}} _{11}},{{\boldsymbol{\varPhi}} _{12}};0,{{\boldsymbol{\varPhi}} _{22}}} \right]$ 是由4个$6 \times 6$ 分块矩阵组成的$12 \times 12$ 的矩阵，状态量${\boldsymbol{z}}\left( {{t_f},{t_0}} \right)$ 由矩阵${\boldsymbol{D}}$ 的存在产生的积分量。这两项均可以通过积分求解得到。因此，式（15）表示了初始状态、初始协态、末端状态、末端协态之间的仿射关系，通过此式并联立边界条件即可解析求解协态初值。对于深空探测中的交会问题，通常考虑如下末端条件
$$ {\boldsymbol{x}}\left( {{t_f}} \right) = {{\boldsymbol{x}}_f},\;\;\;\;m\left( {{t_f}} \right) = {\rm{free}} $$ (17) 由于此系统中质量为常数，质量协态的初值并不影响此系统的解，因此在计算中将质量协态的初值设为0，即
${{\boldsymbol{\lambda}} _m}\left( {{t_0}} \right) = 0$ 。将式（16）代入式（15），其它协态初值的解析表达式为$$ {{\boldsymbol{\lambda}} _x}\left( {{t_0}} \right) = {\boldsymbol{\varPhi}} _{12}^{ - 1}\left( {{{\boldsymbol{x}}_f} - {{\boldsymbol{x}}_0} - {\boldsymbol{z}}} \right) $$ (18) -
对于燃料最优问题，系统1的指标取为燃料消耗，系统2的指标中
$\beta $ 取值和系统1保持一致，即$$ {L_1} = \frac{{{T_{\max }}}}{{{I_{{\rm{sp}}}}{g_0}}}{\tau _1},\;\;\;\;\beta = \frac{{{T_{\max }}}}{{{I_{{\rm{sp}}}}{g_0}}} $$ (19) 将式（18）代入切换系统的表达式，可得切换系统的优化指标为
$$ J = \frac{{{T_{\max }}}}{{{I_{{\rm{sp}}}}{g_0}}}\left[ {\int_{{t_0}}^{{t_s}} {{\tau _1}\;{\rm{d}}t} + \int_{{t_s}}^{{t_f}} {\tau _2^2\;{\rm{d}}t} } \right] $$ (20) 其中：
${t_s}$ 即为切换时刻，可以看到此指标和传统方法燃料和能量最优控制问题的同伦方法基本一致，不同的是此处控制量${\tau _2}$ 为可以解析求解的系统2的控制幅值。同理可以得到切换系统的动力学方程和协态微分方程，即可构建切换系统的两点边值问题。由于基本方程和2.1节内容基本一致，在此不再赘述。此时，当
$\varepsilon = 1$ 时，系统可以完全解析求解，协态初值为式（17）以及${{{\lambda}} _m} = 0$ 。随后，逐渐减小$\varepsilon $ 并将上一步收敛的结果作为初值求解此步的两点边值问题，求解的打靶方程为$$ {\boldsymbol{\phi}} \left( {{\boldsymbol{\lambda}} \left( {{t_0}} \right)} \right) = {\left[ {{{\boldsymbol{x}}}^{\rm{T}}\left( {{t_f}} \right) - {\boldsymbol{x}_f^{\rm{T}}},{{{\lambda}} _m}\left( {{t_f}} \right)} \right]^{\rm{T}}}{\rm{ = 0}} $$ (21) 其中：
${\boldsymbol{\lambda}} \triangleq {\left[ {{\boldsymbol{\lambda}} _x^{\rm{T}},{{{\lambda}} _m}} \right]^{\rm{T}}}$ 为打靶变量。基于解析初始化的迭代过程的详细算法如下1）按照系统2的解析解初始化打靶变量
${\boldsymbol{\lambda}} {\rm{ = }} $ ${\left[ {{{\left[ {{\boldsymbol{\varPhi}} _{12}^{ - 1}\left( {{{\boldsymbol{x}}_f} - {{\boldsymbol{x}}_0} - {\boldsymbol{z}}} \right)} \right]}^{\rm{T}}},{\rm{0}}} \right]^{\rm{T}}}$ ，${\varepsilon _1} = 1$ ；2）以固定的步长
${\rm{d}}\varepsilon = 0.2$ 减小$\varepsilon $ ，即${\varepsilon _{n + 1}} = {\varepsilon _n} - 0.2$ ，若${\varepsilon _n} < 0.2$ ，则${\varepsilon _{n + 1}} = 0$ ；3）以上一步收敛的协态初值作为打靶量的初值进行打靶；
4）若打靶收敛，则转到步骤2继续求解过程，若打靶不收敛，则以更小的步长
${\rm{d}}\varepsilon = 0.5{\rm{d}}\varepsilon $ 重新转到步骤3计算；5）当同伦参数为0时，输出收敛的小推力燃料最优轨迹的结果。
在此切换过程中，若每步均能够收敛，则最终的同伦过程为
${\varepsilon _n} = 1.0,\;\;0.8,\;\;0.6,\;\;0.4,\;\;0.2,\;\;0.0$ ，迭代步数为5，若存在不收敛情况，则迭代步数将增加。因此，在仿真中，将迭代步数和仿真所需时间作为迭代过程计算效率的重要指标。 -
本节仿真试验了自地球出发前往不同目标[1]的深空探测轨迹，以验证方法的有效性。在仿真中，采用了ODE45积分器，其相对误差和绝对误差设置为
${10^{ - 13}}$ ，打靶算法采用Minpack-1，其参数设置参考相关文献[10]。为加快积分并降低数值敏感性，归一化参数设置为长度单位1 AU（天文单位），时间单位1年，质量单位归一化为航天器初始质量。此时，太阳的引力常数为$\mu = 39.476\;926\;{\rm{A}}{{\rm{U}}^{\rm{3}}}{\rm{/}}{{\text{年}}^{\rm{2}}}$ 。所有算例均以C++语言编写（Microsoft Visual Studio 2012编译），并在台式机使用单核运行（Intel Core i7-7700 CPU，主频3.6 GHz，内存8.0 GB）。算例任务参数如表1所示，由于本文采用春分点轨道根数描述，因此各个算例的航天器轨道圈数同时给定。相应时刻的地球和目标的位置速度可根据JPL系统查询得到。
表 1数值仿真任务参数表
Table 1.Parameters for the numerical simulation missions
目标 出发
时间到达
时间运行
圈数比冲/
s推力幅
值/N初始
质量/kg金星 2006.08.03 2008.07.03 2 3 800 0.33 1 500 坦普尔1号 2015.08.09 2016.10.02 1 3 000 0.60 1 000 小行星3671 2012.12.23 2022.08.27 5 3 000 0.32 4 000 -
本文方法的计算效率和传统能量到燃料最优同伦方法的计算效率进行了对比。这两种方法均应用了前文所述的迭代过程。由于能量最优求解需要进行随机猜测初值，因此其每个算例进行了1 000次随机猜测并记录收敛个数。而本文方法采用线性的标称轨道，且由于具有解析初值仅进行单次计算即可。相应的计算收敛率、所需的迭代步数和得到单个解的平均计算时间如表2所示。本文方法对3个算例均能够实现100%收敛，而传统同伦方法的收敛率分别为86.7%、83.6%和39.3%，尤其是对第3个算例，传统方法收敛率较低，计算效率差别较大。
表 2该方法与传统同伦方法计算效率对比
Table 2.Table1 Comparisons of proposed method with the traditional homotopic method
目标 本文方法 传统同伦方法 收敛率 迭代
步数计算时
间/ms收敛率 迭代
步数计算时
间/ms金星 1/1 5 224.0 867/1 000 5 222.5 “坦普尔1号” 1/1 6 291.0 836/1 000 6.3 574.2 小行星3671 1/1 5 360.0 393/1 000 5 1 928.2 -
对每个任务，迭代过程中的控制曲线如图1所示。对于地球至金星的任务，飞行时间为700 d，燃料最优控制包含6个关机段和5个开机段，总的开机时长为275.1 d，消耗燃料210.45 kg。对于地球至“坦普尔1号”（Tempel 1）任务，燃料最优控制包含2个关机段和3个开机段，最短的关机段为初始的弧段，时长0.63 d。任务中总开机时间为197.7 d，消耗燃料348.36 kg。地球至小行星3671（Dionysus）的任务总时长为3 534 d，燃料最优控制包含7个关机段和6个开机段，总计开机时间1 361.5 d，总计燃料消耗1 279.93 kg。3个任务场景的燃料最优轨迹如图2所示，地球至金星任务为2圈近圆轨迹，轨道倾角改变较小，求解相对较为容易，因此本文方法和传统同伦方法均具有较高的求解效率。Tempel 1的轨道为椭圆轨道，因此，此任务的控制曲线较为复杂，而计算所需的迭代步数较多。Dionysus的轨道为倾角较大的椭圆轨道，本文方法求解此问题所需时间仅为传统方法的1/5。
图 1各任务迭代过程的推力曲线图
Figure 1.History of the thrust magnitude of different missions
图 2各任务场景的燃料最优轨迹图
Figure 2.Fuel-optimal trajectories of different missions
-
该研究提出了基于切换系统的解析初始化方法，将小推力燃料最优轨迹优化问题同具有解析协态初值的线性优化问题联系起来，能够解析初始化，并提高间接法小推力燃料最优轨迹的求解效率。仿真表明：本方法和从能量最优同伦至燃料最优轨迹的方法所需的迭代步数基本相同，而其协态初值能够解析得到，因而具有更高的求解效率。
Analytical Initialization for Low-thrust Trajectory Optimization Based on Switching System
-
摘要:传统同伦方法通常将小推力燃料最优问题转化为能量最优问题求解，以增加间接法求解的收敛率，但求解能量最优问题仍需要猜测协态初值以初始化求解算法。该研究针对小推力燃料最优问题，将其优化模型嵌入到切换系统，从而将该问题转化为具有解析协态初值的优化问题进行求解，进一步提高了求解的收敛率，并能够以解析协态初值初始化求解算法。首先，将切换系统优化模型与小推力优化模型相结合，常规的切换系统的切换函数是由最优控制得出的，而本研究则采用了人为设计给定的切换函数，实现了不同系统之间的切换和联系；其次，基于标称轨道线性化的方法，给出了具有解析协态初值的目标系统的设计，无需复杂标称轨道即可实现协态初始化；最后，数值仿真验证了该方法的有效性，相比于传统同伦方法具有更高的求解效率。Abstract:The traditional homotopy method usually transforms the low-thrust fuel-optimal control problem into the energy-optimal problem to increase the convergence rate of the indirect method. However, it is still necessary to guess the initial values of the co-states to initialize the solving algorithm. In this paper, the optimization model of the fuel-optimal problem is embedded in the switching system with the analytical initial co-states, which further improves the convergence rate with the analytical initialization. Firstly, the switching system is introduced with the embedded fuel-optimal problem. The switching function of the conventional switching system is derived from the optimal control, but in this paper, the given switching function is designed artificially to realize the switching and continuation among different systems. Secondly, based on the linearization technique, the target system is designed with analytical initial co-states, initializing the solving algorithm by a simple nominal trajectory. Finally, the numerical simulation verifies the effectiveness of the proposed method, which is more efficient than the traditional homotopy method.Highlights
● The low-thrust fuel-optimal control problem is solved from an analytical solution. ● A switching system is established with a new type of switching function formulation. ● High efficiency solving algorithm is achieved via homotopy to the switching system. -
表 1数值仿真任务参数表
Table 1Parameters for the numerical simulation missions
目标 出发
时间到达
时间运行
圈数比冲/
s推力幅
值/N初始
质量/kg金星 2006.08.03 2008.07.03 2 3 800 0.33 1 500 坦普尔1号 2015.08.09 2016.10.02 1 3 000 0.60 1 000 小行星3671 2012.12.23 2022.08.27 5 3 000 0.32 4 000 
下载:
导出CSV
表 2该方法与传统同伦方法计算效率对比
Table 2Table1 Comparisons of proposed method with the traditional homotopic method
目标 本文方法 传统同伦方法 收敛率 迭代
步数计算时
间/ms收敛率 迭代
步数计算时
间/ms金星 1/1 5 224.0 867/1 000 5 222.5 “坦普尔1号” 1/1 6 291.0 836/1 000 6.3 574.2 小行星3671 1/1 5 360.0 393/1 000 5 1 928.2 下载:
导出CSV
-
[1] JIANG F,TANG G,LI J. Improving low-thrust trajectory optimization by adjoint estimation with shape-based path[J]. Journal of Guidance,Control,and Dynamics,2017,40(12):3282-3289. [2] KITAMURA K,YAMADA K,SHIMA T. Minimum energy coplanar orbit transfer of geostationary spacecraft using time-averaged Hamiltonian[J]. Acta Astronautica,2019,160:270-279.doi:10.1016/j.actaastro.2019.04.033 [3] WALL B J,CONWAY B A. Shape-based approach to low-thrust rendezvous trajectory design[J]. Journal of Guidance,Control,and Dynamics,2009,32(1):95-101. [4] BETTS J T. Very low-thrust trajectory optimization using a direct sqp method[J]. Journal of Computational and Applied Mathematics,2000,120(1-2):27-40.doi:10.1016/S0377-0427(00)00301-0 [5] 刘俊丽,高扬. 十年一剑刃锋利,苦寒方得梅花香——全国空间轨道设计竞赛发展历程回顾[J]. 力学与实践,2019,41(4):488-497.doi:10.6052/1000-0879-19-283LIU J,GAO Y. With ten years’ effort the sword edge is grinded sharp,after the cold winter plum flower blossoms—a review of China Trajectory Optimization Competition[J]. Mechanics in Engineering,2019,41(4):488-497.doi:10.6052/1000-0879-19-283 [6] YANG H,TANG G,JIANG F. Optimization of observing sequence based on nominal trajectories of symmetric observing configuration[J]. Astrodynamics,2018,2(1):25-37.doi:10.1007/s42064-017-0009-2 [7] BETTS J T. Survey of numerical methods for trajectory optimization[J]. Journal of guidance,control,and dynamics,1998,21(2):193-207. [8] GAO Y. Near-optimal very low-thrust earth-orbit transfers and guidance schemes[J]. Journal of guidance,control,and dynamics,2007,30(2):529-539. [9] HABERKORN T,MARTINON P,GERGAUD J. Low thrust minimum-fuel orbital transfer:a homotopic approach[J]. Journal of Guidance,Control,and Dynamics,2004,27(6):1046-1060. [10] JIANG F,BAOYIN H,LI J. Practical techniques for low-thrust trajectory optimization with homotopic approach[J]. Journal of Guidance,Control,and Dynamics,2012,35(1):245-258. [11] BENGEA S C,DECARLO R A. Optimal control of switching systems[J]. Automatica,2005,41(1):11-27. -
点击查看大图
图(2)/
表 (2)
计量
- 文章访问数:422
- HTML全文浏览量:246
- PDF下载量:46
- 被引次数:0
