CN108846162A - 一种基于bezier曲线的连续推力机动轨道设计方法 - Google Patents

Abstract

本发明公开了一种基于bezier曲线的连续推力机动轨道设计方法，包括以下步骤：1)计算轨道机动问题的起止点约束及推力表达式；2)选用7点bezier曲线进行轨道设计；3)设定复合函数P＝B₁H₁+B₂H₂，得复合函数的一阶导数及二阶导数；4)由步骤3)得满足轨道机动任务的复合函数P约束；5)得bezier曲线B₁及bezier曲线B₂的方程；6)将bezier曲线B₁和bezier曲线B₂进行复合，得最终的机动轨道，该方法得到的轨道具有普适性。

Description

一种基于bezier曲线的连续推力机动轨道设计方法

技术领域

本发明属于轨道设计领域，涉及一种基于bezier曲线的连续推力机动轨道设计方法。

背景技术

小推力轨道的表达是解决小推力轨道优化问题的一个关键技术，也是对其优化方法分类的准则。目前较为通用的优化方法包括间接法、直接法和标称轨道法。间接法视原问题为最优控制问题，基于庞特里亚金极小值原理，用变分法推导最优控制的一阶必要条件和终端横截条件，最终将问题归结为用数值方法求解的两点边值问题。这种方法通过对协状态变量的直接数值积分来表达小推力轨道。其主要技术瓶颈在于协状态变量初始猜测不易确定，而且结果对初始猜测的依赖性高，收敛域小。直接法是通过离散和参数化，将原问题转换成非线性规划问题。目前较为常见的直接法是配点法和伪谱法，它们都是同时离散控制变量和状态变量，用多项式来拟合小推力轨道。它们能避免轨道数值积分，但离散的状态变量的初值不容易简便地给出，对结果有效性和鲁棒性有影响，因此需要设计出一种具有普适性的连续推力机动轨道设计方法。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供了一种基于 bezier曲线的连续推力机动轨道设计方法，该方法得到的轨道具有普适性。

为达到上述目的，本发明所述的基于bezier曲线的连续推力机动轨道设计方法，其特征在于，包括以下步骤：

1)计算轨道机动问题的起止点约束及机动轨道推力表达式；

2)选用n点bezier曲线进行轨道设计，其中，设P_j为bezier曲线的第j个控制点，t为bezier曲线的比例因子，1≤j≤n；

3)设定复合函数P＝B₁H₁+B₂H₂，其中，H表示开普勒轨道方程，B₁为第一个bezier曲线，B₂为第一个bezier曲线，得复合函数的一阶导数及二阶导数分别为：

P′＝B′₁.H₁+B₁.H′₁+B′₂.H₂+B₂.H′₂

P″＝B₁″.H₁+2B₁′.H₁′+B₁.H₁″+B₂″.H₂+2B₂′.H₂′+B₂.H₂″；

4)已知满足轨道机动任务的复合函数P约束为：

(P|t＝0)＝H₁，(θ|t＝0)＝θ₁

(P′|t＝0)＝H′₁，(P″|t＝0)＝H₁′

(P|t＝1)＝H₂，(θ|t＝1)＝θ₂

(P′|t＝1)＝H′₂，(P″|t＝1)＝H″₂

由步骤3)得bezier曲线的两点边值问题约束为：

(B₁|t＝0)＝1，(B₁|t＝1)＝0

(θ|t＝0)＝θ₁，(θ|t＝1)＝θ₂

(B₁′|t＝0)＝0，(B₁′|t＝1)＝0

(B₁″|t＝0)＝0，(B₁″|t＝1)＝0

(B₂|t＝0)＝0，(B₂|t＝1)＝1

(θ|t＝0)＝θ₁，(θ|t＝1)＝θ₂

(B₂′|t＝0)＝0，(B₂′|t＝1)＝0

(B₂″|t＝0)＝0，(B₂″|t＝1)＝0；

5)对于bezier曲线B₁，当n＝7时，设(r_i，θ_i)表示第i个控制点的坐标参数，在bezier曲线B₁的起点处，各阶导与控制点坐标的关系为；

在bezier曲线B₁的终点处，各阶导与控制点坐标的关系为：

引入系数k_i，i＝2，3，5，6，通过第一个控制点及第七个控制点的参数对第2、3、5、6个控制点的参数进行表述，得

P₂＝(θ₂，r₂)＝(k₂+θ₁，k₂·(r′|t＝0)+r₁)

P₆＝(θ₆，r₆)＝(-k₆+θ₇，-k₆·(r′|t＝1)+r₇)

将两点边值问题的边界条件代入上述公式中，得

P₂＝(θ₂，r₂)＝(k₂+θ₁，1)

P₃＝(θ₃，r₃)＝(k₃+(2θ₂-θ₁)，1)

P₅＝(θ₅，r₅)＝(k₅+(2θ₆-θ₇)，0)

P₆＝(θ₆，r₆)＝(-k₆+θ₇，0)

同理，得bezier曲线B₂的四个控制点为：

P₉＝(θ₉，r₉)＝(k₉+θ₈，0)

P₁₀＝(θ₁₀，r₁₀)＝(k₁₀+(2θ₉-θ₈)，0)

P₁₂＝(θ₁₂，r₁₂)＝(k₁₂+(2θ₁₃-θ₁₄)，1)

P₁₃＝(θ₁₃，r₁₃)＝(-k₁₃+θ₁₄，1)

再添加bezier曲线B₁中第四个控制点坐标P₄＝(θ₄，r₄)及bezier曲线 B₂中第四个控制点坐标P₁₁＝(θ₁₁，r₁₁)作为自由优化变量，然后加上引入的八个系数k_i得12个优化变量，最后根据所述12个优化变量的优化结果得bezier曲线B₁及bezier曲线B₂的方程；

6)将bezier曲线B₁和bezier曲线B₂进行复合，得最终的机动轨道，完成基于bezier曲线的连续推力机动轨道设计。

步骤6)的具体操作为：bezier曲线B₁和bezier曲线B₂的比例因子 t从0至1等步长变化过程中，对应的和不同，因此bezier曲线B₁和bezier曲线B₂不能直接进行叠加，通过等步长获取对应的及进而通过的关系，得再通过反解多项式反解出对应的从而得到然后再将bezier曲线B₁与bezier曲线B₂进行复合，得最终的机动轨道。

步骤2)中bezier曲线的方程及一阶导数、二阶导数及三阶导数分别为：

反解多项式为：

a₁＝1·θ₈-6·θ₉+15·θ₁₀-20·θ₁₁+15·θ₁₂-6·θ₁₃+1·θ₁₄

a₂＝-6·θ₈+30·θ₉-60·θ₁₀+60·θ₁₁-30·θ₁₂+6·θ₁₃

a₃＝15·θ₈-60·θ₉+90·θ₁₀-60·θ₁₁+15·θ₁₂

a₄＝-20·θ₈+60·θ₉-60·θ₁₀+20·θ₁₁

a₅＝15·θ₈-30·θ₉+15·θ₁₀

a₆＝-6·θ₈+6·θ₉

a₇＝1·θ₈。

步骤1)的具体操作为：

由轨道方程可知，初始轨道及目标轨道的椭圆轨道极坐标下的轨道矢径r及其一阶导r′及二阶导r″分别为：

其中，θ为轨道的相位角，a为轨道的半长轴，e为轨道的偏心率；

通过航迹角γ给出轨道速度方向，其中，V_r及V_θ分别代表径向速度及切向速度，其中，

则每一时刻机动轨道推力T_a(θ)为：

本发明具有以下有益效果：

本发明所述的基于bezier曲线的连续推力机动轨道设计方法在具体操作时，选用7点bezier曲线进行轨道设计，另外，通过两个bezier 曲线构建复合函数，通过该复合函数在极坐标系下对机动轨道进行有效的描述，并对两个bezier曲线中的各控制点坐标进行优化，从而得到具有普适性的最终机动轨道，即既能表述轨道高度的渐变特性，也能较好的描述椭圆轨道的周期震荡特性。与现有的设计方法相比，本发明能够灵活的解决动轨道设计问题，设计结果推力消耗远小于现有方法。

附图说明

图1为四点bezier曲线示意图；

图2为极坐标系轨道图；

图3为直角坐标系轨道图；

图4为推力曲线图；

图5为优化结果的两条bezier曲线示意图；

图6表示总推力消耗随优化过程变化的趋势图；

图7为机动轨道各阶导数与初始轨道和目标轨道匹配情况图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

本发明所述的基于bezier曲线的连续推力机动轨道设计方法包括以下步骤：

1)计算轨道机动问题的起止点约束及机动轨道推力表达式，其中，

由轨道方程可知，初始轨道及目标轨道的椭圆轨道极坐标下的轨道矢径r及其一阶导r′及二阶导r″″分别为：

其中，θ为轨道的相位角，a为轨道的半长轴，e为轨道的偏心率，需要说明，本发明中所有的和A′分别代表A对比例因子t和轨道相位角θ求导；

通过航迹角γ给出轨道速度方向，其中，V_r及V_θ分别代表径向速度及切向速度，其中，

则每一时刻机动轨道推力T_a(θ)为：

由步骤1)可知，对于给定θ和r，若r′相同，则速度方向相同；对于给定θ和r，若r′相同，r″″相同，则速度大小方向均相同；

因此，现在问题就转化成对表述函数的设计和对设计的函数进行优化，使函数满足两点边值问题的边界约束条件，即优化变量使函数在起点及终点的函数值、一阶导、二阶导与初始轨道和目标轨道分别匹配，即可解决连续推力机动轨道设计问题。

2)根据起止点约束可知，基本的6点bezier曲线由于所有控制点都被边界条件进行了约束，因此形式相对较为固定，自由度较差，而适当增加控制点个数则会使bezier曲线更为灵活，可参数化设计性更强。本发明用7点bezier曲线进行轨道设计，其中，设P_j为bezier曲线的第j个控制点，t为bezier曲线的比例因子，1≤j≤7，则bezier曲线的方程及一阶导数、二阶导数及三阶导数分别为：

3)由于设计的函数既要受bezier曲线的调整灵活可变，又要包含椭圆轨道极坐标系下曲线的周期振荡趋势，因此单纯的bezier曲线不足以对轨道特性进行合理的描述，设定复合函数P＝B₁H₁+B₂H₂，其中，H 表示开普勒轨道方程，B₁为第一个bezier曲线，B₂为第一个bezier曲线，得复合函数的一阶导数及二阶导数分别为：

P′＝B1′.H1+B1.H₁′+B₂′.H₂+B₂.H₂′

P″＝B₁″.H₁+2B₁′.H₁′+B₁.H₁″+B₂″.H₂+2B₂＇.H₂′+B₂.H₂″

4)已知轨道机动任务的复合函数P约束为：

(P|t＝0)＝H₁，(θ|t＝0)＝θ₁

(P′|t＝0)＝H′₁，(P″|t＝0)＝H₁′

(P|t＝1)＝H₂，(θ|t＝1)＝θ₂

(P′|t＝1)＝H′₂，(P″|t＝1)＝H₂′

由步骤3)可知，bezier曲线的两点边值问题约束为：

(B₁|t＝0)＝1，(B₁|t＝1)＝0

(θ|t＝0)＝θ₁，(θ|t＝1)＝θ₂

(B₁′|t＝0)＝0，(B₁′|t＝1)＝0

(B₁″|t＝0)＝0，(B₁″|t＝1)＝0

(B₂|t＝0)＝0，(B₂|t＝1)＝1

(θ|t＝0)＝θ₁，(θ|t＝1)＝θ₂

(B₂′|t＝0)＝0，(B₂′|t＝1)＝0

(B₂″|t＝0)＝0，(B₂″|t＝1)＝0

5)对于bezier曲线B₁，设(r_i，θ_i)表示第i个控制点的坐标参数，在bezier曲线B₁的起点处，各阶导与控制点坐标的关系为；

在bezier曲线B₁的终点处，各阶导与控制点坐标的关系为：

引入系数k_i，i＝2，3，5，6，通过第一个控制点及第七个控制点的参数对第2、3、5及6个控制点的参数进行表述，得

P₂＝(θ₂，r₂)＝(k₂+θ₁，k₂·(r′|t＝0)+r₁)

P₆＝(θ₆，r₆)＝(-k₆+θ₇，-k₆·(r′|t＝1)+r₇)

然后将两点边值问题的边界条件代入上述公式中，得

P₂＝(θ₂，r₂)＝(k₂+θ₁，1)

P₃＝(θ₃，r₃)＝(k₃+(2θ₂-θ₁)，1)

P₅＝(θ₅，r₅)＝(k₅+(2θ₆-θ₇)，0)

P₆＝(θ₆，r₆)＝(-k₆+θ₇，0)

同理，得bezier曲线B₂的四个控制点为：

P₉＝(θ₉，r₉)＝(k₉+θ₈，0)

P₁₀＝(θ₁₀，r₁₀)＝(k₁₀+(2θ₉-θ₈)，0)

P₁₂＝(θ₁₂，r₁₂)＝(k₁₂+(2θ₁₃-θ₁₄)，1)

P₁₃＝(θ₁₃，r₁₃)＝(-k₁₃+θ₁₄，1)

再添加bezier曲线B₁中第四个控制点坐标P₄＝(θ₄，r₄)及bezier曲线 B₂中第四个控制点坐标P₁₁＝(θ₁₁，r₁₁)作为自由优化变量，然后加上引入的八个系数k_i得12个优化变量，最后根据所述12个优化变量的优化结果得bezier曲线B₁及bezier曲线B₂的方程；

6)bezier曲线B₁和bezi er曲线B₂的比例因子t从0至1等步长变化过程中，对应的和不同，因此bezier曲线B₁和bezier曲线B₂不能直接进行叠加，通过等步长获取对应的及进而通过的关系，得再通过反解多项式反解出对应的从而得到然后再将bezier曲线B₁与bezier曲线B₂进行复合，得最终的机动轨道，其中，反解多项式为：

a₁＝1·θ₈-6·θ₉+15·θ₁₀-20·θ₁₁+15·θ₁₂-6·θ₁₃+1·θ₁₄

a₂＝-6·θ₈+30·θ₉-60·θ₁₀+60·θ₁₁-30·θ₁₂+6·θ₁₃

a₃＝15·θ₈-60·θ₉+90·θ₁₀-60·θ₁₁+15·θ₁₂

a₄＝-20·θ₈+60·θ₉-60·θ₁₀+20·θ₁₁

a₅＝15·θ₈-30·θ₉+15·θ₁₀

a₆＝-6·θ₈+6·θ₉

a₇＝1·θ₈。

实施例一

本发明对某一轨道机动问题实例实施的结果为：优化结果累计推力为0.2044VU，与现有方法的设计结果(0.3231VU)对比如下：VU为无量纲推力单位，单位对照表如表1所示，仿真案例的具体参数如表2所示，设计结果的beizer曲线各控制点坐标和比例系数如表3所示。

图2及图3分别表示在极坐标系及直角坐标系下的初始轨道、目标轨道及机动轨道，图4为推力曲线图，图5为优化结果的两条bezier 曲线示意图，图6表示总推力消耗随优化过程变化的趋势图，图7为机动轨道各阶导数与初始轨道和目标轨道匹配情况图。

表1

单位	符号	数值
			距离单位	DU	6378km
时间单位	TU	806s
			速度单位	VU	7.9km/s

表2

	初始轨道	目标轨道
			半长轴	2DU	6DU
偏心率	0.3	0.6
			轨道倾角	0	0
升交点赤经	0	0
			近地点幅角	10°	30°
真近点角	-30°	120°

表3bezier曲线控制点参数表

Claims (4)

1.一种基于bezier曲线的连续推力机动轨道设计方法，其特征在于，包括以下步骤：

1)计算轨道机动问题的起止点约束及机动轨道推力表达式；

2)选用n点bezier曲线进行轨道设计，其中，设P_j为bezier曲线的第j个控制点，t为bezier曲线的比例因子，1≤j≤n；

3)设定复合函数P＝B₁H₁+B₂H₂，其中，H表示开普勒轨道方程，B₁为第一个bezier曲线，B₂为第一个bezier曲线，得复合函数的一阶导数及二阶导数分别为：

P′＝B₁′.H₁+B₁.H₁′+B₂′.H₂+B₂.H₂′

P″＝B₁″.H₁+2B₁′.H₁′+B₁.H₁″+B₂″.H₂+2B₂＇.H₂′+B₂.H₂″；

4)已知满足轨道机动任务的复合函数P约束为：

(P|t＝0)＝H₁，(θ|t＝0)＝θ₁

(P′|t＝0)＝H′₁，(P″|t＝0)＝H₁″

(P|t＝1)＝H₂，(θ|t＝1)＝θ₂

(P′|t＝1)＝H′₂，(P″|t＝1)＝H₂″

由步骤3)得bezier曲线的两点边值问题约束为：

(B₁|t＝0)＝1，(B₁|t＝1)＝0

(θ|t＝0)＝θ₁，(θ|t＝1)＝θ₂

(B₁′|t＝0)＝0，(B₁′|t＝1)＝0

(B₁″|t＝0)＝0，(B₁″|t＝1)＝0

(B₂|t＝0)＝0，(B₂|t＝1)＝1

(θ|t＝0)＝θ₁，(θ|t＝1)＝θ₂

(B₂′|t＝0)＝0，(B₂′|t＝1)＝0

(B₂″|t＝0)＝0，(B₂″|t＝1)＝0；

5)对于bezier曲线B₁，当n等于7时，设(r_i，θ_i)表示第i个控制点的坐标参数，在bezier曲线B₁的起点处，各阶导与控制点坐标的关系为；

在bezier曲线B₁的终点处，各阶导与控制点坐标的关系为：

引入系数k_i，i＝2，3，5，6，通过第一个控制点及第七个控制点的参数对第2、3、5、6个控制点的参数进行表述，得

P₂＝(θ₂，r₂)＝(k₂+θ₁，k₂·(r′|t＝0)+r₁)

P₆＝(θ₆，r₆)＝(-k₆+θ₇，-k₆·(r′|t＝1)+r₇)

将两点边值问题的边界条件代入上述公式中，得

P₂＝(θ₂，r₂)＝(k₂+θ₁，1)

P₃＝(θ₃，r₃)＝(k₃+(2θ₂-θ₁)，1)

P₅＝(θ₅，r₅)＝(k₅+(2θ₆-θ₇)，0)

P₆＝(θ₆，r₆)＝(-k₆+θ₇，0)

同理，得bezier曲线B₂的四个控制点为：

P₉＝(θ₉，r₉)＝(k₉+θ₈，0)

P₁₀＝(θ₁₀，r₁₀)＝(k₁₀+(2θ₉-θ₈)，0)

P₁₂＝(θ₁₂，r₁₂)＝(k₁₂+(2θ₁₃-θ₁₄)，1)

P₁₃＝(θ₁₃，r₁₃)＝(-k₁₃+θ₁₄，1)

再添加bezier曲线B₁中第四个控制点坐标P₄＝(θ₄，r₄)及bezier曲线B₂中第四个控制点坐标P₁₁＝(θ₁₁，r₁₁)作为自由优化变量，然后加上引入的八个系数k_i得12个优化变量，然后根据所述12个优化变量的优化结果得bezier曲线B₁及bezier曲线B₂的方程；

6)将bezier曲线B₁和bezier曲线B₂进行复合，得最终的机动轨道，并计算相应的推力消耗，完成基于bezier曲线的连续推力机动轨道设计。

2.根据权利要求1所述的基于bezier曲线的连续推力机动轨道设计方法，其特征在于，步骤6)的具体操作为：bezier曲线B₁和bezier曲线B₂的比例因子t从0至1等步长变化过程中，对应的和不同，因此bezier曲线B₁和bezier曲线B₂不能直接进行叠加，通过等步长获取对应的及进而通过的关系，得再通过反解多项式反解出对应的从而得到然后再将bezier曲线B₁与bezier曲线B₂进行复合，得最终的机动轨道，其中，反解多项式为：

a₁＝1·θ₈-6·θ₉+15·θ₁₀-20·θ₁₁+15·θ₁₂-6·θ₁₃+1·θ₁₄

a₂＝-6·θ₈+30·θ₉-60·θ₁₀+60·θ₁₁-30·θ₁₂+6·θ₁₃

a₃＝15·θ₈-60·θ₉+90·θ₁₀-60·θ₁₁+15·θ₁₂

a₄＝一20·θ₈+60·θ₉-60·θ₁₀+20·θ₁₁

a₅＝15·θ₈-30·θ₉+15·θ₁₀

a₆＝-6·θ₈+6·θ₉

a₇＝1·θ₈。

3.根据权利要求1所述的基于bezier曲线的连续推力机动轨道设计方法，其特征在于，步骤2)中bezier曲线的方程及一阶导数、二阶导数及三阶导数分别为：

4.根据权利要求1所述的基于bezier曲线的连续推力机动轨道设计方法，其特征在于，步骤1)的具体操作为：

由轨道方程可知，初始轨道及目标轨道的椭圆轨道极坐标下的轨道矢径r及其一阶导r′及二阶导r″分别为：

其中，θ为轨道的相位角，a为轨道的半长轴，e为轨道的偏心率；

通过航迹角γ给出轨道速度方向，其中，V_r及V_θ分别代表径向速度及切向速度，其中，

则每一时刻机动轨道推力T_a(θ)为：