¹中国人民解放军32021部队，北京

²南京信息工程大学，江苏 南京

收稿日期：2022年2月16日；录用日期：2022年3月10日；发布日期：2022年3月17日

摘要

以采用有限推力的航天器轨道机动问题为研究对象，根据航天器的运动模型，构建了轨道转移的两点边值问题，给出了发动机开关函数，设计了开关逻辑，应用非线性规划求解两点边值问题，得到了精确的发动机开关时刻和推力方向。最后，用仿真算例验证了本文方法的有效性。

关键词

有限推力，非线性规划，轨道机动，两点边值问题

A Orbit Maneuver Method of the Spacecraft with Finite Thrust

Yafeng Wang¹, Tailong Chen², Kaiguo Fan¹

¹32021 Troops, PLA, Beijing

²Nanjing University of Information Engineering, Nanjing Jiangsu

Received: Feb. 16^th, 2022; accepted: Mar. 10^th, 2022; published: Mar. 17^th, 2022

ABSTRACT

Taking the orbit maneuver problem of spacecraft with finite thrust as the research object, the two-point boundary value problem of orbit transfer is constructed according to the motion model of spacecraft. The engine switching function is given and the switching logic is designed. Then, the two-point boundary value problem is solved by nonlinear programming, and the accurate engine switching time and thrust direction are obtained. Finally, a simulation example is given to verify the effectiveness of the proposed method.

Keywords:Finite Thrust, Nonlinear Programming, Orbit Maneuver, Two-Point Boundary Value Problem

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

卫星轨道机动可采用脉冲轨道机动 [1] [2] 和有限推力轨道机动 [3] - [8] 两种方式。有限推力轨道机动问题可归结为求解非线性两点边值问题，只要解出两点边值问题，即可得到发动机开关逻辑和推力方向的变化规律。为了避开困难的两点边值问题，一些学者用直接法或混合法对最优轨道转移问题进行求解，文献 [3] 用分段多项式来表示整个轨道的状态和控制向量，将最优控制问题转化为非线性规划问题，文献 [4] 采用改进的遗传算法得到了远程导引变轨策略，文献 [5] 在选定变轨策略基础上，以燃料消耗最小为目标，优化给出了GEO卫星快速入轨定点标称轨迹，文献 [6] 采用脉冲式匀速和有限推力式变速直线制导律，获得了燃料最省的制导条件。文献 [7] 将两点边值问题转化为针对协状态初值等的参数优化问题，把开关机时刻作为未知参数去优化，最后应用非线性规划方法求解形成的参数优化问题。文献 [8] 基于任务规划分析，在变轨时间固定、连续变推力的情况下，以总冲最小、满足终端位置和速度约束为指标，对推力大小和方向进行优化。本文在文献 [7] 的基础上，应用非线性规划求解异面最优轨道转移中的两点边值问题，特点在于能得到精确的开关函数，进而确定最优的开关机逻辑，对于有限推力下最优轨道转移问题提供了一种较为有效的方法。

2. 运动模型

航天器的运动，可以通过位置矢量(三个分量)、速度矢量(三个分量)共六个参数描述，也可用六个轨道根素来描述。为避免经典轨道根素的更新在轨道倾角为零或偏心率为零时存在奇异的问题，本文采用改进春分点根素。定义经典轨道根素：a为半长轴；e为偏心率；i为轨道倾角； $\omega$ 为近地点幅角； $\Omega$ 为升交点赤经；v为真近点角。经典轨道根素向改进春分点根素的转化如下：

$p=a\left(1-e^2\right)$ ; $f=e\cos \left(\omega +\Omega \right)$ ; $g=e\sin \left(\omega +\Omega \right)$ ;

$h=\tan \left(i/2\right)\cos \Omega$ ; $k=\tan \left(i/2\right)\sin \Omega$ ; $L=\Omega +\omega +v$.

不考虑地球非球形引力、空气阻力和太阳光压等摄动力，以改进春分点根素表示的运动模型可描述为：

$\stackrel{\dot{}}{X}=Bd\frac{T}{m}+D$ (1)

其中， $X={\left[\begin{array}{cccccc}p& f& g& h& k& L\end{array}\right]}^{\text{T}}$ 为改进春分点根素状态向量；T为发动机推力；m为飞行器质量； $d={\left[\begin{array}{ccc}{d}_r& d_n& d_t\end{array}\right]}^{\text{T}}$ 为推力方向单位矢量，满足 $\Vert d\Vert =1$，$d_r$ 从地心指向飞行器， $d_n$ 沿密切轨道角动量方向， $d_t$ 在密切轨道面内与 $d_r$ 、 $d_n$ 构成右手定则； $D={\left[\begin{array}{cccccc}0& 0& 0& 0& 0& D\end{array}\right]}^{\text{T}}$，$D=\sqrt{p\mu }{\left(w/p\right)}^2$ ； $B$ 的形式如下：

$B=\left[\begin{array}{ccc}0& B_{12}& 0\\ B_{21}& B_{22}& B_{23}\\ B_{31}& B_{32}& B_{33}\\ 0& 0& B_{43}\\ 0& 0& B_{53}\\ 0& 0& B_{63}\end{array}\right]$

各元素的表达式为：

$B_{12}=\sqrt{\frac{p}{\mu }}\frac{2p}{w}$ ; $B_{21}=\sqrt{\frac{p}{\mu }}\sin L$ ; $B_{22}=\sqrt{\frac{p}{\mu }}\frac{\left(1+w\right)\cos L+f}{w}$ ;

$B_{23}=-\sqrt{\frac{p}{\mu }}\frac{\left(h\sin L-k\cos L\right)g}{w}$ ; $B_{31}=-\sqrt{\frac{p}{\mu }}\cos L$ ;

$B_{32}=\sqrt{\frac{p}{\mu }}\frac{\left(1+w\right)\sin L+g}{w}$ ; $B_{33}=\sqrt{\frac{p}{\mu }}\frac{\left(h\sin L-k\cos L\right)f}{w}$ ;

$B_{43}=\sqrt{\frac{p}{\mu }}\frac{s^2\cos L}{2w}$ ; $B_{53}=\sqrt{\frac{p}{\mu }}\frac{s^2\sin L}{2w}$ ; $B_{63}=\sqrt{\frac{p}{\mu }}\frac{\left(h\sin L-k\cos L\right)}{w}$.

其中， $\mu =3.986012\times {10}^{14}{\text{m}}^{\text{3}}/{\text{s}}^{\text{2}}$ 为地球引力系数； $w=1+f\cos L+g\sin L$，$s^2=1+h^2+k^2$ 为中间变量。

有限推力交会变轨过程中，质量的变化是必须考虑的，为此，引入质量变化方程：

$\stackrel{\dot{}}{m}=-\frac{T}{I_{sp}{g}_0}$ (2)

其中： $I_{sp}$ 为发动机比冲， $g_0$ 为海平面处地球重力加速度。

3. 两点边值问题

对于有限推力作用下的轨道优化问题，要解决的问题就是给出某一性能指标，确定喷气发动机何时进行开机关机操作及发动机开机时的喷气方向(最优推力方向)。

3.1. 性能指标

定义交会初始时刻为 $t_0$，终端时刻为 $t_f$。在轨道优化问题中，性能指标常见为燃料消耗最少或交会变轨时间最短，本文选取燃料消耗最少为性能指标：

$J=-m\left(t_f\right)$ (3)

进行轨道优化的目标就是使J取得极小值，即明确初始时刻与终端时刻的约束条件，以发动机的开关时刻和喷气方向为待优化参数，求解优化问题 $\min \left(J\right)$。由于开关逻辑的不确定， $\min \left(J\right)$ 的直接求解是非常困难的，可以引入协状态概念，构建哈密顿函数，设计开关逻辑，在改进春分点根素的运动方程上增加协状态的运动方程构成正则方程，形成两点边值问题，通过以协状态为待优化参数求解两点边值问题，得到协状态初始值，进而得到喷气发动机的开关时刻和喷气方向。

3.2. 协状态的运动方程

定义 $\lambda ={\left[\begin{array}{cccccc}{\lambda }_p& {\lambda }_f& {\lambda }_g& {\lambda }_h& {\lambda }_k& {\lambda }_L\end{array}\right]}^{\text{T}}$ 为改进春分点根素的协状态， ${\lambda }_m$ 为对应于质量的协状态。构建

系统哈密顿函数为：

$H={\lambda }^{\text{T}}Bd\frac{T}{m}+{\lambda }_{L}D-{\lambda }_m\frac{T}{I_{sp}{g}_0}$ (4)

轨道优化的目标要使H取得极小值，最优推力方向矢量应取如下形式：

$d^{*}=-\frac{{\left({\lambda }^{\text{T}}B\right)}^{\text{T}}}{\Vert {\lambda }^{\text{T}}B\Vert }$ (5)

对于有限推力交会变轨，推力T应使H取极小值，由(4)可知，H关于T是线性的：

$\frac{\partial H}{\partial T}=-\frac{\Vert {\lambda }^{\text{T}}B\Vert }{m}-{\lambda }_m\frac{1}{I_{sp}{g}_0}$ (6)

取发动机开关函数S为：

$S=\frac{\Vert {\lambda }^{\text{T}}B\Vert }{m}+{\lambda }_m\frac{1}{I_{sp}{g}_0}$ (7)

发动机的开关逻辑可设计为：

1) 若： $S<0$，则： $T=0$ ；

2) 若： $S>0$，则： $T=T_{\max }$ ；

3) 若： $S=0$，则： $T=0$。

其中： $T_{\max }$ 为推力最大值。

此开关逻辑的物理意思在于：当 $S<0$ 时，令 $T=0$，使得H保持当前值不增大；当 $S>0$ 时，令 $T=T_{\max }$，使得H往小的方向变化；当 $S=0$ 时，T的变化对H没有影响，为节省能量，可取 $T=0$。

协状态的变化方程为：

$\stackrel{\dot{}}{\lambda }=-{\lambda }^{\text{T}}\frac{\partial B}{\partial X}{d}^{*}\frac{T}{m}-{\lambda }_L\frac{\partial D}{\partial X}$ (8)

${\stackrel{\dot{}}{\lambda }}_m=-\Vert {\lambda }^{\text{T}}B\Vert \frac{T}{m^2}$ (9)

其中：

$\frac{\partial B}{\partial p}=\left[\begin{array}{ccc}0& \frac{3B_{12}}{2p}& 0\\ \frac{B_{21}}{2p}& \frac{B_{22}}{2p}& \frac{B_{23}}{2p}\\ \frac{B_{31}}{2p}& \frac{B_{32}}{2p}& \frac{B_{33}}{2p}\\ 0& 0& \frac{B_{43}}{2p}\\ 0& 0& \frac{B_{53}}{2p}\\ 0& 0& \frac{B_{53}}{2p}\end{array}\right]$ ; $\frac{\partial B}{\partial f}=\left[\begin{array}{ccc}0& -\frac{\cos L}{w}{B}_{12}& 0\\ 0& \frac{B_{12}}{2p}\left({\cos }^{2}L+1\right)-\frac{B_{22}}{w}\cos L& -\frac{\cos L}{w}{B}_{23}\\ 0& \frac{\cos L}{w}\left(B_{21}-B_{32}\right)& B_{33}\left(\frac{1}{f}-\frac{\cos L}{w}\right)\\ 0& 0& -\frac{\cos L}{w}{B}_{43}\\ 0& 0& -\frac{\cos L}{w}{B}_{53}\\ 0& 0& -\frac{\cos L}{w}{B}_{63}\end{array}\right]$

$\frac{\partial B}{\partial g}=\left[\begin{array}{ccc}0& -\frac{\sin L}{w}{B}_{12}& 0\\ 0& \cos L\frac{B_{21}}{w}-\frac{B_{22}}{w}\sin L& B_{23}\left(\frac{1}{g}-\frac{\sin L}{w}\right)\\ 0& \frac{B_{12}}{2p}\left({\sin }^{2}L+1\right)-\frac{B_{32}}{w}\sin L& -\frac{\sin L}{w}{B}_{33}\\ 0& 0& -\frac{\sin L}{w}{B}_{43}\\ 0& 0& -\frac{\sin L}{w}{B}_{53}\\ 0& 0& -\frac{\sin L}{w}{B}_{63}\end{array}\right]$ ; $\frac{\partial B}{\partial h}=\left[\begin{array}{ccc}0& 0& 0\\ 0& \frac{B_{12}}{2p}\left({\cos }^{2}L+1\right)-\frac{B_{22}}{w}\cos L& -B_{21}\frac{g}{w}\\ 0& 0& B_{21}\frac{f}{w}\\ 0& 0& -B_{31}\frac{h}{w}\\ 0& 0& B_{21}\frac{h}{w}\\ 0& 0& B_{21}\frac{1}{w}\end{array}\right]$

$\frac{\partial B}{\partial k}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -B_{31}\frac{g}{w}\\ 0& 0& B_{31}\frac{f}{w}\\ 0& 0& -B_{31}\frac{k}{w}\\ 0& 0& B_{21}\frac{k}{w}\\ 0& 0& B_{31}\frac{1}{w}\end{array}\right]$ ; $\frac{\partial B}{\partial L}=\left[\begin{array}{ccc}0& -\frac{B_{12}}{w}\frac{\partial w}{\partial L}& 0\\ -B_{31}& -\frac{1}{w}\left(B_{31}\frac{\partial w}{\partial L}+B_{21}\left(w+1\right)+B_{22}\frac{\partial w}{\partial L}\right)& -bg-\frac{B_{23}}{w}\frac{\partial w}{\partial L}\\ B_{21}& \frac{1}{w}\left(B_{21}\frac{\partial w}{\partial L}-B_{31}\left(w+1\right)-B_{32}\frac{\partial w}{\partial L}\right)& bf-\frac{B_{33}}{w}\frac{\partial w}{\partial L}\\ 0& 0& -\left(B_{53}+\frac{B_{43}}{w}\frac{\partial w}{\partial L}\right)\\ 0& 0& B_{43}-\frac{B_{53}}{w}\frac{\partial w}{\partial L}\\ 0& 0& b-\frac{B_{63}}{w}\frac{\partial w}{\partial L}\end{array}\right]$ ;

$\frac{\partial D}{\partial X}={\left[\begin{array}{cccccc}-\frac{3w}{2p{B}_{12}}& \frac{4\cos L}{B_{12}}& \frac{4\sin L}{B_{12}}& 0& 0& \frac{4}{B_{12}}\frac{\partial w}{\partial L}\end{array}\right]}^{\text{T}}$ ;

$\frac{\partial w}{\partial L}=-f\sin L+g\cos L$,$b=\sqrt{\frac{p}{\mu }}\left(h\cos L+k\sin L\right)\frac{1}{w}$.

至此，得到了系统的正则方程组(1)、(2)和(8)、(9)。

3.3. 边值约束

优化问题的始端状态约束直接由交会变轨的起始状态决定，终端状态约束为追踪航天器与目标航天器的末端状态一致，即

$X\left(t_f\right)-X_{obj}\left(t_f\right)=0$,(10)

其中： $X\left(t_f\right)$ 为追踪航天器末端状态； $X_{obj}\left(t_f\right)$ 为目标航天器末端状态。

协状态无明确物理意思，初始值为两点边值问题的待优化参数，终端约束可由横截条件计算得出，具体过程此处不再详述，只给出结果：

${\left[\begin{array}{ccc}{\lambda }_L\left(t_f\right)& {\lambda }_m\left(t_f\right)& H\left(t_f\right)\end{array}\right]}^{\text{T}}={\left[\begin{array}{ccc}0& -1& 0\end{array}\right]}^{\text{T}}$. (11)

4. 两点边值问题的求解

上述内容将有限推力作用下最优轨道转移问题转变为两点边值问题，下面进一步将其整理为可以使用非线性规划方法求解的形式。

待优化参数为七个协状态初值和终端时刻：

${\lambda }_p\left(t_0\right),{\lambda }_f\left(t_0\right),{\lambda }_g\left(t_0\right),{\lambda }_h\left(t_0\right),{\lambda }_k\left(t_0\right),{\lambda }_L\left(t_0\right),{\lambda }_m\left(t_0\right),t_f$.

优化目标为最小化性能指标函数： $J=-m\left(t_f\right)$。

参数的约束条件为：满足 $t_0<t_f<t_{\max }$，$t_{\max }$ 为完成轨道转移的最大时限；满足正则方程(1)、(2)和(8)、(9)。

边界条件为：满足式(10)和(11)。

两点边值问题转化为非线性规划问题之后，就可以选择合适的非线性规划算法对其进行求解。非线性规划算法有很多，本文采用内点算法，求解得到七个协状态初值和终端时刻 $t_f$ 等参数后，根据追踪航天器的初始条件，就可对正则方程组(1)、(2)和(8)、(9)进行数值积分，由开关函数(7)控制发动机的开关机操作，由推力矢量控制逻辑控制其推力方向，从而得到各个状态和协状态的时间历程。

至此，我们就将有限推力交会变轨的整个优化控制过程给出，下面将结合某一算例对此优化控制作仿真验证。

5. 仿真算例

假定追踪航天器和被动航天器的初始经典轨道根素如表1所示，追踪航天器的初始推重比为0.05，发动机比冲为2000s，初始时刻为 $t_0=0$，性能指标为燃料最省。各单位量均用长度单位1 DU = 6578.145 km 与时间单位1 TU = 806.8119 s进行无量纲化处理。仿真平台为配置2.90 GHz CPU和8.00 GB RAM的DESKTOP-MP3FGKJ PC，仿真软件为matlab2016b。

表1. 追踪航天器和被动航天器的初始经典轨道根素

整个交会过程耗时27.3492 TU，得到的仿真结果如下：

图1. 哈密顿函数变化曲线

图2. 开关函数变化曲线

图3. 推力切换曲线

图4. 推力方向变化曲线

图5. 两航天器交会曲线

由图1可以看出，哈密顿函数在轨道转移末端时刻收敛到零，因为末端时刻不再喷气且协状态 ${\lambda }_L\left(t_f\right)=0$。从图2和图3可以看出，当开关函数值大于零时，才有推力作用，推力方向由图4确定。图5显示了交会变轨过程中，追踪航天器和目标航天器的运行轨迹。

表2. 交会结束时两航天器相对位置误差

表2给出了交会结束时追踪航天器和目标航天器的相对位置误差，可以看出，追踪航天器逼近被动航天器的精度较高。

6. 结论

研究了一种应用非线性规化求解有限推力作用下异面最优轨道转移问题的方法。采用改进春分点根素形式的运动模型，将有限推力作用下异面最优轨道转移问题转化为两点边值问题。在考虑边界条件、横截条件及开关函数的前提下，将两点边值问题转化为针对协状态初值等的参数优化问题，最后应用非线性规划方法求解形成的参数优化问题。本方法特点是能得到开关函数从而得到最优发动机开关机逻辑，无需对发动机开关机做事先假设。

文章引用

王亚锋,陈泰龙,范开国. 一种采用有限推力的航天器轨道机动方法
A Orbit Maneuver Method of the Spacecraft with Finite Thrust[J]. 应用数学进展, 2022, 11(03): 1044-1052. https://doi.org/10.12677/AAM.2022.113113

参考文献