-
摘要: 为了满足小推力航天器交会轨迹的快速性设计需求，基于形状逼近理论，设计了一种三维轨迹模型.将轨迹设计问题转换为求解傅里叶级数的系数问题，避免了轨迹运动方程非线性强、难以求解的难题，极大地提高了计算效率.考虑到推力加速度的限制，建立了加速度约束方程，并结合轨迹的运动方程，给出了傅里叶级数的求解过程.同时根据边界条件和最大推力加速度值，定性地分析了傅里叶系数的存在条件.仿真验证了该方法的正确性和可行性，并从计算效率上与高斯伪谱法进行了对比，结果表明本文的方法计算耗时仅为高斯伪谱法的0.67%.Abstract: Trajectory design with low-thrust propulsion needs a method for quickly approximating the spacecraft's trajectory during the preliminary stage phase. To address this requirement, a 3D rendezvous trajectory design model is proposed based on the shape-based theory, by which the trajectory design problem is simplified into solving the coefficients of Fourier series. In the process of solving coefficients, thrust constraint is considered, and the existence conditions of the coefficients are qualitatively analyzed. Finally, the correctness and feasibility of this approach are verified through a numerical example. Results show that the proposed method can not only satisfy the boundary conditions and thrust constraint, but also has a good computational efficiency. Its computation time is only 0.67% of that of the Gauss pseudospectral method.
-
Key words:
- Quickly approximating /
- low-thrust /
- trajectory design /
- thrust constraint /
- Fourier serious /
- spacecraft
-
航天器任务设计的初步阶段,需要快速地设计出一条合理的初始轨迹并评估燃料的消耗量,然而由于小推力轨迹模型的非线性较强,直接以解析解的形式求出运动轨迹非常困难,而且实际工程中推力具有最大限制,极大地增加了任务挑战性 [1].传统的方法中一般将小推力轨迹设计问题转化为最优控制问题来求解,求解方法主要分为两类 [2-4]: 直接法和间接法. 但是这两类方法对初始猜测值比较敏感,而且对于不同的初始值需要对全过程重新进行数值积分,不能满足任务初步阶段的快速性需求. 因此,寻找新的快速性轨迹设计方法非常重要.
为了解决任务设计初步阶段的快速性需求,1999年Petropoulos等 [5]提出了一种形状逼近法,为小推力轨迹设计问题打开了一扇新的大门. 形状逼近法是一种基于逆向设计的思想,它首先假设运动轨迹呈某一特定的曲线,然后用形状曲线拟合的方法设计出满足要求的转移轨迹. 由于该方法简单、快速,迅速引起了学者的广泛关注. Petropoulos等 [6] 用正弦指数函数模型设计了小推力拦截轨迹. Izzo [7] 研究了基于正弦指数函数模型的多圈Lambert 问题. 然而,对于交会问题正弦指数函数模型不能满足终端速度约束条件.随后,Wall等 [8-9] 提出了一种6 阶逆多项式的方法,该方法虽然能满足边界条件约束,但是对于交会问题,只能通过增大转移时间或增加轨迹的圈数来减小最大推力加速度值. 为了解决推力约束的问题,Taheri等 [10-11] 引入了傅里叶函数,设计了一种近似轨迹模型,但是该模型适合于转移圈数给定的情况,而且文中没有给出系数的存在条件,对于快速性轨迹设计需求存在一定的局限性.
针对空间小推力轨道交会问题,本文在前人的基础上引入傅里叶级数模型,设计了一种基于傅里叶级数的空间小推力轨迹设计方法. 通过轨迹模型,将小推力轨迹设计问题转化成了求解傅里叶系数的问题.根据加速度约束方程和轨迹运动约束方程,结合序列二次规划算法给出了系数的求解过程,并定性分析了傅里叶系数解的存在条件. 最后仿真验证了该方法的正确性和可行性. 本文所设计的方法可为工程设计初步阶段提供一定的技术参考,或为更进一步的精确计算提供初值依据.
1. 基于傅里叶级数的轨迹模型
1.1 问题的提出
在轨道交会过程中,航天器在空间的转移轨迹是一条满足约束条件的空间曲线. 通常情况下,对于任意给定的一条空间曲线都可以近似用一个傅里叶级数的形式来表示,所以航天器的轨迹设计问题可以转化为傅里叶级数的设计问题,而傅里叶级数取决于多项式的系数,因此,只要确定了多项式的系数,就可以确定航天器在空间的运动轨迹.
为了便于分析和计算,本文以初始轨道平面为参考面,建立柱坐标系,将航天器的轨迹表示为极角${\theta}$、 半径${r}$和轨道面法向位置${z}$的函数,图中${x}$轴和${y}$轴位于初始轨道面内,与${z}$轴构成右手坐标系,如图 1所示.
根据牛顿引力定律,航天器的运动方程在柱坐标系下可以表示为 [9, 12]
$\left\{ \begin{array}{l} \ddot r - r{{\dot \theta }^2} = - \dfrac{\mu }{{{s^3}}}r + {T_{ain}}\sin \alpha \\dfrac{1}{r}\dfrac{\rm d}{{\rm d}t}\left( {{r^2}\dot \theta } \right) = {T_{ain}}\cos \alpha \\ddot z = - \dfrac{\mu }{{{s^3}}}z + {T_{az}} \end{array} \right. $
(1) 式中,${z}$表示$oz$轴方向的位置,${T_{ain}}$表示航天器在$oxy$平面内的推力加速度大小,${T_{az}}$ 表示$oz$轴方向的推力加速度大小,$\alpha $ 为$oxy$ 平面内的推力方向角,$\mu $ 为地心引力常数,s表示地心距,其中$s = \sqrt {{r^2} + {z^2}} $.
在柱坐标系下,只要确定了航天器每个时刻的状态$\left( {r,{\kern 1pt} {\kern 1pt} \theta ,{\kern 1pt} {\kern 1pt} z} \right)$,航天器的空间运动轨迹便可以唯一确定. 下面采用傅里叶级数的设计思想,分别对$oxy$平面内的运动和z向的运动进行建模分析.
1.2 模型的建立
1.2.1 oxy平面内的运动
在平面内航天器的位置是轨道半径r和极角$\theta$ 的函数,r和$\theta$ 可以表示成两种函数形式 [6, 8],一种是以$\theta $为自变量,用$\theta $ 来表示函数r; 另一种是将r 和$\theta $都表示成关于时间t 的函数. 这里选取第二种形式,航天器在$oxy$ 平面内的运动轨迹模型如下:
$\left\{ \begin{array}{l} r\left( t \right) = \dfrac{{{a_0}}}{2} + \displaystyle\sum\limits_{n = 1}^{{n_r}} {\left\{ {{a_n}\cos \left( {\dfrac{{n\pi }}{T}t} \right) + {b_n}\sin \left( {\dfrac{{n\pi }}{T}t} \right)} \right\}} \\theta \left( t \right) = \dfrac{{{c_0}}}{2} + \displaystyle\sum\limits_{n = 1}^{{n_\theta }} {\left\{ {{c_n}\cos \left( {\dfrac{{n\pi }}{T}t} \right) + {d_n}\sin \left( {\dfrac{{n\pi }}{T}t} \right)} \right\}} \end{array} \right. $
(2) 式中,T表示从离轨点到入轨点之间的飞行时间,${a_0}$,${a_n}$,${b_n}$和${c_0}$,${c_n}$,${d_n}$为待定系数,${n_r},{n_\theta } \in {\text{N}}$且${n_r} \geqslant 2,{n_\theta } \geqslant 2$.
假设平面内推力方向角等于飞行路径角,即
$\alpha = \gamma + b\pi $
(3) 式中,$\gamma $表示飞行路径角,当推力方向与速度方向相同时b等于1,相反时b 等于0.
对式(1)中的第二项进行整理可得
${T_{ain}} = \dfrac{{2\dot r\dot \theta + r\ddot \theta }}{{\cos \alpha }} $
(4) 将上式代入式(1)中的第一项,整理可得
$\begin{align} & \ddot{r}-r{{{\dot{\theta }}}^{2}}+\frac{\mu }{{{s}^{3}}}r=\frac{2\dot{r}\dot{\theta }+r\ddot{\theta }}{\cos \alpha }\sin \alpha = \\ & \left( 2\dot{r}\dot{\theta }+r\ddot{\theta } \right)\tan \alpha \\ \end{align}$
(5) 将式(3)代入式(5)可得
$\ddot{r}-r{{{\dot{\theta }}}^{2}}+\frac{\mu }{{{s}^{3}}}r=\left( 2\dot{r}\dot{\theta }+r\ddot{\theta } \right)\tan \gamma $
(6) 根据切向速度与径向速度之间的关系,可知推力方向角的表达式为
$\tan \gamma = \dfrac{{\dot r}}{{r\dot \theta }} $
(7) 将式(7)代入式(6),可得$oxy$平面内轨迹约束方程为
$\begin{align} & f\left( r,\dot{r},\ddot{r},\dot{\theta },\ddot{\theta },s \right)={{r}^{2}}\left( \frac{\mu }{{{s}^{3}}}\dot{\theta }-{{{\dot{\theta }}}^{3}} \right)+ \\ & r\left( \dot{\theta }\ddot{r}-\dot{r}\ddot{\theta } \right)-2{{{\dot{r}}}^{2}}\dot{\theta }=0 \\ \end{align}$
(8) 1.2.2 平面外的运动
平面外的运动即$oz$轴方向的运动,考虑到转移轨迹具有一定的周期性,为了尽量减少优化参数,同时提高轨迹模型的逼近能力,这里选取极角$\theta $ 作为自变量,将平面外的运动表示成关于余弦函数和高阶多项式的复合函数的形式
$z\left( \theta \right) = {a_z}\cos \theta + {b_z}\theta + {c_z}{\theta ^{q - 1}} + {d_z}{\theta ^q} $
(9) 式中,${a_z}$,${b_z}$,${c_z}$,${d_z}$为多项式的系数,q为大于等于3的正整数. 由于平面外的运动边界条件已知,即初始和末端时刻航天器z向运动的位置和速度参数(${z_0},{z_f},{\dot z_0},{\dot z_f}$)是确定的,所以可以直接求出式(9)中的多项式系数.
将式(9)代入式(1)中的第三项得z轴方向推力加速度值为
$\begin{align} & {{T}_{az}}=\frac{\mu }{{{s}^{3}}}z+\ddot{\theta }\left[ \left( q-1 \right){{c}_{z}}{{\theta }^{q-2}}+q{{d}_{z}}{{\theta }^{q-1}}+ \right. \\ & \left. {{b}_{z}}-{{a}_{z}}\sin \theta \right]+{{{\dot{\theta }}}^{2}}\left[ q\left( q-1 \right){{d}_{z}}{{\theta }^{q-2}} \right.+ \\ & \left. \left( q-1 \right)\left( q-2 \right){{c}_{z}}{{\theta }^{q-3}}-{{a}_{z}}\cos \theta \right] \\ \end{align}$
(10) 由式(4)和式(10)可得总的推力加速度大小为
${T_a} = \sqrt {T_{ain}^2 + T_{az}^2} $
(11) 由于实际工程中,发动机的推力是有限制的,所以在轨迹设计中还需要考虑如下约束方程
${f_a}\left( {r,\dot r,\dot \theta ,\ddot \theta ,\alpha } \right) = \left( {\dfrac{{{T_a}}}{{{T_{a,\max }}}}} \right) \le 1 $
(12) 式中,${T_{a,\max }}$表示轨道机动过程中允许的最大推力加速度值.
2. 傅里叶级数系数的求解
由于傅里叶级数各项系数的好坏直接影响着轨迹模型的逼近能力,因此求解满足边界条件和约束方程的最优系数是本文设计方法的关键.
将式(2)和(9)代入式(8)和(12),可以得到一组仅与未知系数和时间变量相关的非线性代数方程
$f\left( {{a}_{rj}},{{a}_{\theta k}},{{a}_{z}},\cdots ,{{d}_{z}},t \right)=0$
(13) ${{f}_{a}}\left( {{a}_{rj}},{{a}_{\theta k}},{{a}_{z}},\cdots ,{{d}_{z}},t \right)\le 1$
(141) 式中,$j = 0,\cdots ,m$,$k = 0,\cdots ,n$. 式(13)和(14)即为各个时刻航天器的运动轨迹需要满足的所有约束函数. 求解未知系数至少需要$2\left( {{n_r} + {n_\theta } + 1} \right)$ 个方程. 为了得出约束方程,需要先对轨道转移时间进行离散化,然后针对每一个离散点建立满足约束函数的方程.
通常情况下,${n_r}$和${n_\theta }$取10以内的值即可达到较高的精度. 对于带推力约束的轨迹设计问题,${n_r}$和${n_\theta }$的取值大于等于3才能保证推力在约束范围以内,为了保证边界点处的位置和速度精度,未知系数分两部分求取.
2.1 由边界条件直接确定的系数
令初始时刻${t_0} = 0$,将初始时刻${t_0}$和末端时刻${t_f}$的状态参数$\left( {{r_0},{r_f},{{\dot r}_0},{{\dot r}_f}} \right)$分别代入式(13),可得式(2)中函数r 的可确定系数为
$\left\{ \begin{matrix} {{a}_{1}}=\frac{{{r}_{0}}-{{r}_{f}}}{z}-\sum\limits_{{{\lambda }_{1}}=3}^{{{n}_{r}}}{{{a}_{{{\lambda }_{1}}}}} \\ {{a}_{2}}=\frac{{{r}_{0}}+{{r}_{f}}-{{a}_{0}}}{2}-\sum\limits_{{{\lambda }_{2}}=4}^{{{n}_{r}}}{{{a}_{{{\lambda }_{2}}}}} \\ {{b}_{1}}=\frac{T}{2\pi }\left( {{{\dot{r}}}_{0}}-{{{\dot{r}}}_{f}} \right)-\sum\limits_{{{\lambda }_{1}}=3}^{{{n}_{r}}}{\left( {{\lambda }_{1}}{{b}_{{{\lambda }_{1}}}} \right)} \\ {{b}_{2}}=\frac{T}{4\pi }\left( {{{\dot{r}}}_{0}}+{{{\dot{r}}}_{f}} \right)-\frac{1}{2}\sum\limits_{{{\lambda }_{2}}=4}^{{{n}_{r}}}{\left( {{\lambda }_{2}}{{b}_{{{\lambda }_{2}}}} \right)} \\ \end{matrix} \right.$
(15) 式中,${\lambda _1} \ge 3$,${\lambda _1}$ 为奇数,${\lambda _2} \ge 4$,${\lambda _2}$ 为偶数. $r_0$和${r_f}$分别为初始和末端时刻轨道半径,${\dot r_0}$和${\dot r_f}$ 分别为初始和末端时刻径向速度.
同理,将初始时刻${t_0}$和末端时刻${t_f}$的状态参数$\left( {{\theta _0},{\theta _f},{{\dot \theta }_0},{{\dot \theta }_f}} \right)$分别代入式(13),可得式(2)中函数$\theta $的可确定系数为
$\left\{ \begin{matrix} {{c}_{1}}=\frac{{{\theta }_{0}}-{{\theta }_{f}}}{2}-\sum\limits_{{{\tau }_{1}}=3}^{{{n}_{\theta }}}{{{c}_{{{\tau }_{1}}}}} \\ {{c}_{2}}=\frac{{{\theta }_{0}}+{{\theta }_{f}}-{{c}_{0}}}{2}-\sum\limits_{{{\tau }_{2}}=4}^{{{n}_{\theta }}}{{{c}_{{{\tau }_{2}}}}} \\ {{d}_{1}}=\frac{T}{2\pi }\left( {{{\dot{\theta }}}_{0}}-{{{\dot{\theta }}}_{f}} \right)-\sum\limits_{{{\tau }_{1}}=3}^{{{n}_{\theta }}}{\left( {{\tau }_{1}}{{d}_{{{\tau }_{1}}}} \right)}\text{ } \\ {{d}_{2}}=\frac{T}{4\pi }\left( {{{\dot{\theta }}}_{0}}+{{{\dot{\theta }}}_{f}} \right)-\frac{1}{2}\sum\limits_{{{\tau }_{2}}=4}^{{{n}_{\theta }}}{\left( {{\tau }_{2}}{{d}_{{{\tau }_{2}}}} \right)} \\ \end{matrix} \right.$
(16) 式中,${\tau _1} \ge 3$,${\tau _1}$ 为奇数,${\tau _2} \ge 4$,${\tau _2}$ 为偶数. ${\theta _0}$和${\theta _f}$分别为初始和末端时刻的极角,${\dot \theta _0}$ 和${\dot \theta _f}$分别为初始和末端时刻极角的变化率.
同理,将${t_0}$和${t_f}$时刻的状态参数$\left( {{z_0},{z_f},{{\dot z}_0},{{\dot z}_f}} \right)$ 分别代入式(2)和式(9),整理得
$\left[ \begin{matrix} {{a}_{z}} \\ {{b}_{z}} \\ {{c}_{z}} \\ {{d}_{z}} \\ \end{matrix} \right]={{\left[ \begin{matrix} A & B \\ \end{matrix} \right]}^{-1}}\cdot \left[ \begin{matrix} {{z}_{0}} \\ {{{\dot{z}}}_{0}} \\ {{z}_{f}} \\ {{{\dot{z}}}_{f}} \\ \end{matrix} \right]$
(17) 其中
$\begin{align} & A=\left[ \begin{matrix} \cos {{\theta }_{0}} & {{\theta }_{0}} \\ -{{{\dot{\theta }}}_{0}}\sin {{\theta }_{0}} & {{{\dot{\theta }}}_{0}} \\ \cos {{\theta }_{f}} & {{\theta }_{f}} \\ -{{{\dot{\theta }}}_{f}}\sin {{\theta }_{f}} & {{{\dot{\theta }}}_{f}} \\ \end{matrix} \right] \\ & B=\left[ \begin{matrix} \theta _{0}^{q-1} & \theta _{0}^{q} \\ \left( q-1 \right)\theta _{0}^{q-2}{{{\dot{\theta }}}_{0}} & q\theta _{0}^{q-1}{{{\dot{\theta }}}_{0}} \\ \theta _{f}^{q-1} & \theta _{f}^{q} \\ \left( q-1 \right)\theta _{f}^{q-2}{{{\dot{\theta }}}_{f}} & q\theta _{f}^{q-1}{{{\dot{\theta }}}_{f}} \\ \end{matrix} \right] \\ \end{align}$
${z_0}$和${z_f}$分别为初始时刻$oz$轴方向的位置,${\dot z_0}$和${\dot z_f}$ 分别为末端时刻$oz$轴方向的速度值.
需要指出的是,式中的${\theta _f} = \tilde \theta + {N_{rev}} \times 2\pi $,其中$\tilde \theta $ 表示${r_0}$和${r_f}$两矢量之间的夹角,${N_{rev}}$ 为交会过程中航天器转过的圈数,具体如图 2 所示. 对于极角${\theta _f}$没有约束的情况,${\theta _f}$ 可以通过确定${N_{rev}}$的值来获得一个最优值,${N_{rev}}$的求解过程将会在后文中具体给出.
2.2 待优化系数的求解
2.2.1 优化算法和求解步骤
为了使计算过程尽量简单,本文选用序列二次规划(Sequential quadratic programming,SQP)[10]的方法对待求系数进行优化求解. 在计算的过程中为了避免有多个极值点和陷入局部极小值的情况,采用了Lagrange乘子法,对Lagrange函数取二次逼近,将带有非线性约束的优化问题转换为二次规划问题,通过求解指标函数的最小值来获得待求系数的最优值. 本文以燃料消耗作为优化指标函数
$\min J = \int_{{t_0}}^{{t_f}} {{T_a}} {\rm d}t $
(18) 对于给定的轨道转移时间T和极角${\theta _f}$、傅里叶级数${n_r}$和${n_\theta }$,以及离散点$DP$的个数,待优化多项式系数的计算步骤如下: 1)根据初始时刻和终端时刻航天器的状态参数$( {{r_0},{\kern 1pt} {\kern 1pt} {\theta _0},{\kern 1pt} {\kern 1pt} {r_f},{\kern 1pt} {\kern 1pt} {\theta _f},{{\dot r}_0},{{\dot \theta }_0},{{\dot r}_f},{{\dot \theta }_f}} )$,在$oxy$平面内,利用第2.2.2 节中的方法初步拟合出一条满足式(13)的参考轨迹曲线; 2)结合离散点$DP$的个数,将轨道交会时间区间$0 \sim {t_f}$,离散成($DP-1$) 个区间,利用步骤1)中给出参考轨迹模型,计算出每个节点i处的状态参数$( {{r_i},{\kern 1pt} {\kern 1pt} {\theta _i},{\kern 1pt} {\kern 1pt} {{\dot r}_i},{{\dot \theta }_i},{{\ddot \theta }_i},{z_i},{{\dot z}_i}})$,将这些状态参数代入式(13) 中,通过等式约束方程初步估算出多项式系数的猜测值; 3)在猜测值的基础上,根据每个节点处的约束方程,采用序列二次规划的方法,优化求解出满足所有约束方程(13)和(14)的一组最优解,所得的最优解即为所求的待优化多项式系数.
2.2.2 待优化系数的初值猜测
在用SQP方法求解待优化系数时,首先要给出一组初始猜测值. 为了避免自主猜测的盲目性,本文采用初始值自主给定的方式,首先假设存在一条满足等式(13)约束的参考轨迹,然后根据参考轨迹上的离散点,拟合出初始猜测值. 文献[8, 10]给出了逆多项式法和 指数正弦函数两种参考轨迹模型,但是文献中的模型都是以极角$\theta $作为自变量,而本文的模型是以时间为自变量,为了能够简单、快速地得出初始参考轨迹,这里选用了立方多项式的方法.
r和$\theta $的参考轨迹模型如下
$\left\{ \begin{array}{l} r\left( t \right) = a{t^3} + b{t^2} + ct + d\\theta \left( t \right) = e{t^3} + f{t^2} + gt + h \end{array} \right. $
(19) 式中,a,b,c,d和e,f,g,h均为常值,其值可由初始和末端边界状态参数直接得出.
2.3 待优化系数的存在条件
待优化系数的存在性,决定着是否存在满足设计要求的轨迹曲线,在有推力约束轨迹的设计任务中,轨迹设计的难点是寻找满足推力加速限制的曲线. 对于本文的方法,影响轨迹曲线性能的因素,除了轨迹模型的逼近能力外,主要的影响因素就是转移轨迹的圈数${{N_{rev}}}$. 本节将定性给出${{N_{rev}}}$ 的确定方法.
对于给定的轨道交会任务,交会时间是确定的,轨道交会时间${T}$满足以下约束
${P_{t0}} \cdot \dfrac{{{\theta _f}}}{{2\pi }}＜T＜{P_{tf}} \cdot \dfrac{{{\theta _f}}}{{2\pi }} $
(20) 式中${P_{t0}}$和${P_{tf}}$分别表示初始和目标轨道的周期(这里假设${P_{t0}}＜{P_{tf}}$).
为了保证算法的收敛、避免轨迹发生奇异现象,由式(20)可得转移轨迹的圈数${{N_{rev}}}$取值范围为
$\dfrac{T}{{{P_{tf}}}} - \dfrac{{\tilde \theta }}{{2\pi }}＜{N_{rev}} ＜\frac{T}{{{P_{t0}}}} - \frac{{\tilde \theta }}{{2\pi }} $
(21) 此外,为了保证在推力约束范围以内存在满足边界条件的最优解,除了${\theta _f} \ge \tilde \theta $以外,还应满足以下条件
${\theta _f} > \dfrac{{\left| {{v_0} - {v_f}} \right|}}{{{T_{a,\max }}}} \cdot \sqrt {\dfrac{{8\mu }}{{{{\left( {{a_0} + {a_f}} \right)}^3}}}} $
(22) 即
${{N_{rev}} > \dfrac{1}{{2\pi }}\left( {\dfrac{{\left| {{v_0} - {v_f}} \right|}}{{{T_{a,\max }}}} \cdot \sqrt {\dfrac{{8\mu }}{{{{\left( {{a_0} + {a_f}} \right)}^3}}}} - \tilde \theta } \right)} $
(23) 式中${v_0}$和${v_f}$表示离轨点和入轨点的速度大小,${a_0}$和${a_f}$表示初始轨道和目标轨道的半长轴.
对于不同的轨道交会任务,根据给定的设计参数${N_{rev}}$可能存在多个可行值,由于${N_{rev}}$ 的值越大,所需离散点的个数就越多,直接会导致计算量增大,为了满足快速性设计需求,在计算过程中选取${N_{rev}}$的最小值作为最优值.
3. 数值分析
为了论证本文所设计的方法的可行性和有效性,选取如下算例进行仿真验证. 仿真电脑为Windows 7 系统,内存(RAM) 4 GB,仿真平台为Matlab 2014a.
假设航天器初始轨道参数和目标轨道参数如表 1所示,表中的a、 e、 i、 $\omega$、 $\Omega$和f分别为轨道半长轴、 偏心率、轨道倾角、 近地点幅角、升交点赤经和真近点角. 输入参数${n_r} = 4$,${n_\theta } = 5$,q取9,离散点数$DP$取22,整个轨道转移时间为17 449 s,初始质量为4 000 kg,发动机的比冲为3 000 m/s. 根据式(21)和(23)计算可得$1.9306＜{{N}_{rev}}＜2.8830$,由于${N_{rev}}$ 为正整数,所以${N_{rev}}$取2. 为了便于计算,仿真中将参数进行无量纲化,选取参考长度为$DU$,1 $DU$ = 6 378.1 km,参考时间为$TU,1TU\text{=}\sqrt{D{{U}^{3}}/D{{U}^{3}}\mu }$,发动机的推力限制设为${T_{a,\max }} = 0.014 DU/T{U^2}$.
表 1 输入状态参数Table 1 Input boundary parameters参数 a e i $\omega $ $\Omega $ f 初始状态7 178.1km 0 0± 70± 20± 0± 末端状态9 378.1km 0.01 2± 90± 0± 180± 由前文可知,本文方法中设计参数有q、离散点个数$DP$、${n_r}$和${n_\theta }$,为了分析各设计参数对计算效率和精度的影响,这里结合仿真算例,在给定的设计参数基础上,分别以某一设计参数作为变量进行仿真. 从仿真结果来看,初始和末端的速度、位置计算精度均在${10^{ - 16}}$量级,由于篇幅有限,表 2～4 给出了部分仿真结果,其中$\Delta v$ 表示燃料消耗,${T_{\max }}$表示最大推力加速度值,$\Delta t$表示仿真时间.
表 2 q取不同值时的仿真结果Table 2 Results with different qq ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 5 0.2453 0.0141 52.5338 7 0.2221 0.0140 4.8387 9 0.2016 0.0140 0.4418 13 0.1969 0.0138 0.4245 17 0.1862 0.0139 0.3718 21 0.2101 0.0140 13.2845 25 0.2269 0.0140 193.4338 27 0.2972 0.0141 239.1416 表 3 $DP$取不同值时的仿真结果Table 3 Results with different $DP$$DP\cdot 2\pi /{{\theta }_{f}}$ ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 6 0.2088 0.0140 19.6138 10 0.2011 0.0139 0.5181 20 0.2010 0.0140 3.1638 30 0.2012 0.0138 3.6007 40 0.2005 0.0140 10.7402 60 0.2003 0.0140 15.2663 80 0.2002 0.0140 36.1296 90 0.2002 0.0140 47.4839 表 4 ${n_r},{n_\theta }$取不同值时的仿真结果Table 4 Results with different ${n_r}$ and ${n_\theta }$nr=${{n}_{\theta }}$ ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 3 0.2773 0.0140 39.1008 4 1.9671 0.0140 3.6494 5 0.2117 0.0139 6.0813 6 0.2013 0.0136 3.8057 8 0.1974 0.0132 25.8237 10 0.1902 0.0132 45.5559 从表 2～4中的数据分析可知,单独改变q和$DP$的值对最大推力加速度值${T_{\rm max}}$的影响不是很明显. 当q较小或过大时,轨迹模型的逼近能力较差,计算耗时长,而且精度下降. q直接影响的是轨迹模型的逼近能力,根据多项式的特性可知,q取7附近的值时,式(9)的逼近能力最强,所以q 一般取7附近的值. 对于转移轨迹,当每圈离散点的个数在10 以上时,可以达到较高的计算精度,对于本文的算例当每圈离散点的个数超过60时,计算精度基本保持不变,但是由于计算量增加,计算效率明显下降. ${n_r}$ 和${n_\theta }$ 的取值对${T_{\rm max}}$ 影响比较明显,随着${n_r} $、$ {n_\theta }$的增大,${T_{\rm max}}$ 和$\Delta v$ 逐渐下降,并趋于稳定,当${n_r} $、$ {n_\theta }$ 增大到一定值时,计算精度趋于饱和,由于模型中多项式的系数增加,计算时间快速上升. 同时,由表 2和表 4对比分析可知,在轨迹设计过程中选定q,将${n_r} $和$ {n_\theta }$ 作为优化变量,更容易快速得出精确的计算结果.
令柱坐标系下初始轨道面内的x轴、y轴与地心惯性坐标系在赤道面内的两个坐标轴重合,算例的仿真结果如图 3～5.考虑到逆多项式设计法(Inverse polynomial,IP method)是形状逼近法中最经典的方法,为了直观地分析本文方法(Proposed method)的合理性,仿真中同时给出了逆多项式设计法的仿真结果.
图 3为轨道交会过程中航天器的轨迹变化,从图中可以看出采用本文的设计方法,在离轨第一圈时间内,轨迹变化较平缓,基本沿着初始轨道附近运动,随着时间的推进,z方向的运动逐渐变快,随后平滑地切入到目标轨道交会点,进入目标轨道. 仿真得到的离轨点状态值$( {{r_0},{\theta _0},{z_0},{{\dot r}_0},{{\dot \theta }_0},{{\dot z}_0}}$) 为(1.1254,1.2034E-18,2.0311E-18,-3.0251E-19,0.8376,4.3610E-18),入轨点的状态值(${{r_f},{\theta _f},{z_f},{{\dot r}_f},{{\dot \theta }_f},{{\dot z}_f}})$ 为(1.4842,3.1415,-0.0518,-1.5261E-18,0.5501,1.3892E-19),与初始给定的状态值误差均小于10$^{-17}$量级,表明所得轨迹曲线能很好地满足边界约束条件.
为了更直观地分析轨迹变化情况,图 4给出了轨道交会过程中航天器地心距的变化曲线,从曲线的变化趋势可以看出在轨迹末端实线和虚线完全重合,两种方式都能保证航天器精确地到达交会点,但是相比于逆多项式设计法,本文方法在离轨和入轨区间段轨迹变化更为平缓.
图 5给出了推力加速的变化曲线,由于本文方法在算法中考虑了推力加速度的限制,推力加速度的最大值保持在给定的范围以内. 而逆多项式方法假定运动轨迹为某一特定的曲线,限制了模型的逼近能力,最大推力加速度值超出了允许的范围.
考虑到高斯伪谱法是最优控制方法中精度非常高的一种算法,最后采用高斯伪谱法对算例进行了仿真. 由表 5中的数据可知本文的方法燃料消耗非常接近高斯伪谱法,说明本文方法燃料预估精度较高,而且明显高于逆多项式方法. 同时本文方法计算得出的${T_{\max }}$ 值小于0.014,进一步表明本文的方法具有满足推力加速度约束的能力. 从仿真时间来看,由于逆多项式是以解析解的形式给出的结果,计算时间最短,但是相比于逆多项式算法仅多了0.2789 s; 相比于高斯伪谱法,本文方法避免了反复迭代寻找初始协态变量的过程,计算时间减少了61.488 s,仅为高斯伪谱法的0.67 %,极大提高了计算效率.
表 5 三种方法计算结果比较Table 5 Comparison of the results using three methods方法 ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 本文的方法 0.1894 0.0139 0.4130 逆多项式法 0.2093 0.0186 0.1341 高斯伪谱法 0.1849 0.0140 61.901 4. 结论
本文针对小推力航天器空间交会轨迹设计问题,提出了一种基于傅里叶级数的快速轨迹设计方法,并通过仿真验证了该方法的正确性和可行性. 本文的方法主要有以下3个优点: 1) 模型简单、计算快速: 通过形状曲线拟合运动轨迹,克服了非线性运动方程难以求解的困难,直接将轨迹设计问题简化成了求解傅里叶级数系数的问题,提高了计算效率; 2) 求解精度较高: 避免了建模时假设运动轨迹成某一特定形状曲线的限制,增强了形状曲线的逼近能力,改善了求解精度; 3) 可用性好: 设计过程中考虑了发动机的推力限制,所得结果不仅能满足边界条件,而且还能满足推力加速度约束条件. 本文方法可为小推力航天器任务设计初步阶段的交会轨迹设计和燃料预估提供一定的参考,或为更进一步的精确计算提供可靠的初始值.
本文方法旨在为任务设计初步阶段小推力轨迹设计问题提供一种新思路,对于如何定量地选择设计参数(q、${n_r} $、$ {n_\theta }$和$DP$)、考虑摄动干扰和常值小推力等轨迹设计问题,将在后续的工作中做进一步的深入研究.
-
表 1 输入状态参数
Table 1 Input boundary parameters
参数 a e i $\omega $ $\Omega $ f 初始状态7 178.1km 0 0± 70± 20± 0± 末端状态9 378.1km 0.01 2± 90± 0± 180± 
下载: 导出CSV
表 2 q取不同值时的仿真结果
Table 2 Results with different q
q ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 5 0.2453 0.0141 52.5338 7 0.2221 0.0140 4.8387 9 0.2016 0.0140 0.4418 13 0.1969 0.0138 0.4245 17 0.1862 0.0139 0.3718 21 0.2101 0.0140 13.2845 25 0.2269 0.0140 193.4338 27 0.2972 0.0141 239.1416
下载: 导出CSV
表 3 $DP$取不同值时的仿真结果
Table 3 Results with different $DP$
$DP\cdot 2\pi /{{\theta }_{f}}$ ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 6 0.2088 0.0140 19.6138 10 0.2011 0.0139 0.5181 20 0.2010 0.0140 3.1638 30 0.2012 0.0138 3.6007 40 0.2005 0.0140 10.7402 60 0.2003 0.0140 15.2663 80 0.2002 0.0140 36.1296 90 0.2002 0.0140 47.4839
下载: 导出CSV
表 4 ${n_r},{n_\theta }$取不同值时的仿真结果
Table 4 Results with different ${n_r}$ and ${n_\theta }$
nr=${{n}_{\theta }}$ ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 3 0.2773 0.0140 39.1008 4 1.9671 0.0140 3.6494 5 0.2117 0.0139 6.0813 6 0.2013 0.0136 3.8057 8 0.1974 0.0132 25.8237 10 0.1902 0.0132 45.5559
下载: 导出CSV
表 5 三种方法计算结果比较
Table 5 Comparison of the results using three methods
方法 ${{\vartriangle }_{v}}$(DU=TU) Tmax(DU=TU2) ${{\vartriangle }_{t}}$(s) 本文的方法 0.1894 0.0139 0.4130 逆多项式法 0.2093 0.0186 0.1341 高斯伪谱法 0.1849 0.0140 61.901
下载: 导出CSV
-
[1] Laipert F E,Longuski J M.Low-thrust trajectories for human missions to Ceres.Acta Astronautica,2014,95:124-132 doi: 10.1016/j.actaastro.2013.11.003 [2] Graham K F,Rao A V.Minimum-time trajectory optimization of multiple revolution low-thrust earth-orbit transfers.Journal of Spacecraft and Rockets,2015,52(3):711-727 doi: 10.2514/1.A33187 [3] Sun C,Yuan J P,Fang Q.Continuous low thrust trajectory optimization for preliminary design.Proceedings of the Institution of Mechanical Engineers.Part G:Journal of Aerospace Engineering,2016,230(5):921-933 [4] Topputo F,Zhang C.Survey of direct transcription for low-thrust space trajectory optimization with applications.Abstract and Applied Analysis,2014,2014:851720,DOI: 10.1155/2014/851720 [5] Petropoulos A E,Longuski J M,Vinh N X.Shape-based analytic representations of low-thrust trajectories for gravity-assist applications.In:Proceedings of the 2000 Advances in Astronautical Sciences.Girdwood,USA:AAS,2000.563-581 [6] Petropoulos A E,Longuski J M.Shape-based algorithm for the automated design of low-thrust,gravity assist trajectories.Journal of Spacecraft and Rockets,2004,41(5):787-796 doi: 10.2514/1.13095 [7] Izzo D.Lambert's problem for exponential sinusoids.Journal of Guidance,Control,and Dynamics,2006,29(5):1242-1245 doi: 10.2514/1.21796 [8] Wall B J,Conway B A.Shape-based approach to low-thrust rendezvous trajectory design.Journal of Guidance,Control,and Dynamics,2009,32(1):95-101 doi: 10.2514/1.36848 [9] Wall B J.Shape-based approximation method for low-thrust trajectory optimization.In:Proceedings of the 2008 AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Honolulu,Hawaii:AIAA,2008. [10] Taheri E,Abdelkhalik O.Shape based approximation of constrained low-thrust space trajectories using Fourier series.Journal of Spacecraft and Rockets,2012,49(3):535-546 http://cn.bing.com/academic/profile?id=2154522750&encoded=0&v=paper_preview&mkt=zh-cn [11] Taheri E,Abdelkhalik O.Initial three-dimensional low-thrust trajectory design.Advances in Space Research,2016,57(3):889-903 doi: 10.1016/j.asr.2015.11.034 [12] Gondelach D J,Noomen R.Hodographic-shaping method for low-thrust interplanetary trajectory design.Journal of Spacecraft and Rockets,2015,52(3):728-738 doi: 10.2514/1.A32991 期刊类型引用(2)
1. 李俊阳，董立财，林肖. 基于改进傅里叶级数轨迹的运动学参数辨识. 湖南大学学报(自然科学版). 2025(02): 55-63 . 
百度学术2. 李思远，叶东，孙兆伟. 考虑通信延迟的卫星集群改进蜂拥控制. 宇航学报. 2022(08): 1097-1108 . 百度学术其他类型引用(3)
-
下载:
计量
- 文章访问数: 1516
- HTML全文浏览量: 196
- PDF下载量: 861
- 被引次数: 5
