本发明属于轨道设计领域，涉及一种基于bezier曲线的连续推力机动轨道设计方法。
背景技术：
：小推力轨道的表达是解决小推力轨道优化问题的一个关键技术，也是对其优化方法分类的准则。目前较为通用的优化方法包括间接法、直接法和标称轨道法。间接法视原问题为最优控制问题，基于庞特里亚金极小值原理，用变分法推导最优控制的一阶必要条件和终端横截条件，最终将问题归结为用数值方法求解的两点边值问题。这种方法通过对协状态变量的直接数值积分来表达小推力轨道。其主要技术瓶颈在于协状态变量初始猜测不易确定，而且结果对初始猜测的依赖性高，收敛域小。直接法是通过离散和参数化，将原问题转换成非线性规划问题。目前较为常见的直接法是配点法和伪谱法，它们都是同时离散控制变量和状态变量，用多项式来拟合小推力轨道。它们能避免轨道数值积分，但离散的状态变量的初值不容易简便地给出，对结果有效性和鲁棒性有影响，因此需要设计出一种具有普适性的连续推力机动轨道设计方法。技术实现要素：本发明的目的在于克服上述现有技术的缺点，提供了一种基于bezier曲线的连续推力机动轨道设计方法，该方法得到的轨道具有普适性。为达到上述目的，本发明所述的基于bezier曲线的连续推力机动轨道设计方法，其特征在于，包括以下步骤：1)计算轨道机动问题的起止点约束及机动轨道推力表达式；2)选用n点bezier曲线进行轨道设计，其中，设Pj为bezier曲线的第j个控制点，t为bezier曲线的比例因子，1≤j≤n；3)设定复合函数P＝B1H1+B2H2，其中，H表示开普勒轨道方程，B1为第一个bezier曲线，B2为第一个bezier曲线，得复合函数的一阶导数及二阶导数分别为：P′＝B′1.H1+B1.H′1+B′2.H2+B2.H′2P″＝B1″.H1+2B1′.H1′+B1.H1″+B2″.H2+2B2′.H2′+B2.H2″；4)已知满足轨道机动任务的复合函数P约束为：(P|t＝0)＝H1，(θ|t＝0)＝θ1(P′|t＝0)＝H′1，(P″|t＝0)＝H1′(P|t＝1)＝H2，(θ|t＝1)＝θ2(P′|t＝1)＝H′2，(P″|t＝1)＝H″2由步骤3)得bezier曲线的两点边值问题约束为：(B1|t＝0)＝1，(B1|t＝1)＝0(θ|t＝0)＝θ1，(θ|t＝1)＝θ2(B1′|t＝0)＝0，(B1′|t＝1)＝0(B1″|t＝0)＝0，(B1″|t＝1)＝0(B2|t＝0)＝0，(B2|t＝1)＝1(θ|t＝0)＝θ1，(θ|t＝1)＝θ2(B2′|t＝0)＝0，(B2′|t＝1)＝0(B2″|t＝0)＝0，(B2″|t＝1)＝0；5)对于bezier曲线B1，当n＝7时，设(ri，θi)表示第i个控制点的坐标参数，在bezier曲线B1的起点处，各阶导与控制点坐标的关系为；在bezier曲线B1的终点处，各阶导与控制点坐标的关系为：引入系数ki，i＝2，3，5，6，通过第一个控制点及第七个控制点的参数对第2、3、5、6个控制点的参数进行表述，得P2＝(θ2，r2)＝(k2+θ1，k2·(r′|t＝0)+r1)P6＝(θ6，r6)＝(-k6+θ7，-k6·(r′|t＝1)+r7)将两点边值问题的边界条件代入上述公式中，得P2＝(θ2，r2)＝(k2+θ1，1)P3＝(θ3，r3)＝(k3+(2θ2-θ1)，1)P5＝(θ5，r5)＝(k5+(2θ6-θ7)，0)P6＝(θ6，r6)＝(-k6+θ7，0)同理，得bezier曲线B2的四个控制点为：P9＝(θ9，r9)＝(k9+θ8，0)P10＝(θ10，r10)＝(k10+(2θ9-θ8)，0)P12＝(θ12，r12)＝(k12+(2θ13-θ14)，1)P13＝(θ13，r13)＝(-k13+θ14，1)再添加bezier曲线B1中第四个控制点坐标P4＝(θ4，r4)及bezier曲线B2中第四个控制点坐标P11＝(θ11，r11)作为自由优化变量，然后加上引入的八个系数ki得12个优化变量，最后根据所述12个优化变量的优化结果得bezier曲线B1及bezier曲线B2的方程；6)将bezier曲线B1和bezier曲线B2进行复合，得最终的机动轨道，完成基于bezier曲线的连续推力机动轨道设计。步骤6)的具体操作为：bezier曲线B1和bezier曲线B2的比例因子t从0至1等步长变化过程中，对应的和不同，因此bezier曲线B1和bezier曲线B2不能直接进行叠加，通过等步长获取对应的及进而通过的关系，得再通过反解多项式反解出对应的从而得到然后再将bezier曲线B1与bezier曲线B2进行复合，得最终的机动轨道。步骤2)中bezier曲线的方程及一阶导数、二阶导数及三阶导数分别为：反解多项式为：a1＝1·θ8-6·θ9+15·θ10-20·θ11+15·θ12-6·θ13+1·θ14a2＝-6·θ8+30·θ9-60·θ10+60·θ11-30·θ12+6·θ13a3＝15·θ8-60·θ9+90·θ10-60·θ11+15·θ12a4＝-20·θ8+60·θ9-60·θ10+20·θ11a5＝15·θ8-30·θ9+15·θ10a6＝-6·θ8+6·θ9a7＝1·θ8。步骤1)的具体操作为：由轨道方程可知，初始轨道及目标轨道的椭圆轨道极坐标下的轨道矢径r及其一阶导r′及二阶导r″分别为：其中，θ为轨道的相位角，a为轨道的半长轴，e为轨道的偏心率；通过航迹角γ给出轨道速度方向，其中，Vr及Vθ分别代表径向速度及切向速度，其中，则每一时刻机动轨道推力Ta(θ)为：本发明具有以下有益效果：本发明所述的基于bezier曲线的连续推力机动轨道设计方法在具体操作时，选用7点bezier曲线进行轨道设计，另外，通过两个bezier曲线构建复合函数，通过该复合函数在极坐标系下对机动轨道进行有效的描述，并对两个bezier曲线中的各控制点坐标进行优化，从而得到具有普适性的最终机动轨道，即既能表述轨道高度的渐变特性，也能较好的描述椭圆轨道的周期震荡特性。与现有的设计方法相比，本发明能够灵活的解决动轨道设计问题，设计结果推力消耗远小于现有方法。附图说明图1为四点bezier曲线示意图；图2为极坐标系轨道图；图3为直角坐标系轨道图；图4为推力曲线图；图5为优化结果的两条bezier曲线示意图；图6表示总推力消耗随优化过程变化的趋势图；图7为机动轨道各阶导数与初始轨道和目标轨道匹配情况图。具体实施方式下面结合附图对本发明做进一步详细描述：本发明所述的基于bezier曲线的连续推力机动轨道设计方法包括以下步骤：1)计算轨道机动问题的起止点约束及机动轨道推力表达式，其中，由轨道方程可知，初始轨道及目标轨道的椭圆轨道极坐标下的轨道矢径r及其一阶导r′及二阶导r″″分别为：其中，θ为轨道的相位角，a为轨道的半长轴，e为轨道的偏心率，需要说明，本发明中所有的和A′分别代表A对比例因子t和轨道相位角θ求导；通过航迹角γ给出轨道速度方向，其中，Vr及Vθ分别代表径向速度及切向速度，其中，则每一时刻机动轨道推力Ta(θ)为：由步骤1)可知，对于给定θ和r，若r′相同，则速度方向相同；对于给定θ和r，若r′相同，r″″相同，则速度大小方向均相同；因此，现在问题就转化成对表述函数的设计和对设计的函数进行优化，使函数满足两点边值问题的边界约束条件，即优化变量使函数在起点及终点的函数值、一阶导、二阶导与初始轨道和目标轨道分别匹配，即可解决连续推力机动轨道设计问题。2)根据起止点约束可知，基本的6点bezier曲线由于所有控制点都被边界条件进行了约束，因此形式相对较为固定，自由度较差，而适当增加控制点个数则会使bezier曲线更为灵活，可参数化设计性更强。本发明用7点bezier曲线进行轨道设计，其中，设Pj为bezier曲线的第j个控制点，t为bezier曲线的比例因子，1≤j≤7，则bezier曲线的方程及一阶导数、二阶导数及三阶导数分别为：3)由于设计的函数既要受bezier曲线的调整灵活可变，又要包含椭圆轨道极坐标系下曲线的周期振荡趋势，因此单纯的bezier曲线不足以对轨道特性进行合理的描述，设定复合函数P＝B1H1+B2H2，其中，H表示开普勒轨道方程，B1为第一个bezier曲线，B2为第一个bezier曲线，得复合函数的一阶导数及二阶导数分别为：P′＝B1′.H1+B1.H1′+B2′.H2+B2.H2′P″＝B1″.H1+2B1′.H1′+B1.H1″+B2″.H2+2B2＇.H2′+B2.H2″4)已知轨道机动任务的复合函数P约束为：(P|t＝0)＝H1，(θ|t＝0)＝θ1(P′|t＝0)＝H′1，(P″|t＝0)＝H1′(P|t＝1)＝H2，(θ|t＝1)＝θ2(P′|t＝1)＝H′2，(P″|t＝1)＝H2′由步骤3)可知，bezier曲线的两点边值问题约束为：(B1|t＝0)＝1，(B1|t＝1)＝0(θ|t＝0)＝θ1，(θ|t＝1)＝θ2(B1′|t＝0)＝0，(B1′|t＝1)＝0(B1″|t＝0)＝0，(B1″|t＝1)＝0(B2|t＝0)＝0，(B2|t＝1)＝1(θ|t＝0)＝θ1，(θ|t＝1)＝θ2(B2′|t＝0)＝0，(B2′|t＝1)＝0(B2″|t＝0)＝0，(B2″|t＝1)＝05)对于bezier曲线B1，设(ri，θi)表示第i个控制点的坐标参数，在bezier曲线B1的起点处，各阶导与控制点坐标的关系为；在bezier曲线B1的终点处，各阶导与控制点坐标的关系为：引入系数ki，i＝2，3，5，6，通过第一个控制点及第七个控制点的参数对第2、3、5及6个控制点的参数进行表述，得P2＝(θ2，r2)＝(k2+θ1，k2·(r′|t＝0)+r1)P6＝(θ6，r6)＝(-k6+θ7，-k6·(r′|t＝1)+r7)然后将两点边值问题的边界条件代入上述公式中，得P2＝(θ2，r2)＝(k2+θ1，1)P3＝(θ3，r3)＝(k3+(2θ2-θ1)，1)P5＝(θ5，r5)＝(k5+(2θ6-θ7)，0)P6＝(θ6，r6)＝(-k6+θ7，0)同理，得bezier曲线B2的四个控制点为：P9＝(θ9，r9)＝(k9+θ8，0)P10＝(θ10，r10)＝(k10+(2θ9-θ8)，0)P12＝(θ12，r12)＝(k12+(2θ13-θ14)，1)P13＝(θ13，r13)＝(-k13+θ14，1)再添加bezier曲线B1中第四个控制点坐标P4＝(θ4，r4)及bezier曲线B2中第四个控制点坐标P11＝(θ11，r11)作为自由优化变量，然后加上引入的八个系数ki得12个优化变量，最后根据所述12个优化变量的优化结果得bezier曲线B1及bezier曲线B2的方程；6)bezier曲线B1和bezier曲线B2的比例因子t从0至1等步长变化过程中，对应的和不同，因此bezier曲线B1和bezier曲线B2不能直接进行叠加，通过等步长获取对应的及进而通过的关系，得再通过反解多项式反解出对应的从而得到然后再将bezier曲线B1与bezier曲线B2进行复合，得最终的机动轨道，其中，反解多项式为：a1＝1·θ8-6·θ9+15·θ10-20·θ11+15·θ12-6·θ13+1·θ14a2＝-6·θ8+30·θ9-60·θ10+60·θ11-30·θ12+6·θ13a3＝15·θ8-60·θ9+90·θ10-60·θ11+15·θ12a4＝-20·θ8+60·θ9-60·θ10+20·θ11a5＝15·θ8-30·θ9+15·θ10a6＝-6·θ8+6·θ9a7＝1·θ8。实施例一本发明对某一轨道机动问题实例实施的结果为：优化结果累计推力为0.2044VU，与现有方法的设计结果(0.3231VU)对比如下：VU为无量纲推力单位，单位对照表如表1所示，仿真案例的具体参数如表2所示，设计结果的beizer曲线各控制点坐标和比例系数如表3所示。图2及图3分别表示在极坐标系及直角坐标系下的初始轨道、目标轨道及机动轨道，图4为推力曲线图，图5为优化结果的两条bezier曲线示意图，图6表示总推力消耗随优化过程变化的趋势图，图7为机动轨道各阶导数与初始轨道和目标轨道匹配情况图。表1单位符号数值距离单位DU6378km时间单位TU806s速度单位VU7.9km/s表2初始轨道目标轨道半长轴2DU6DU偏心率0.30.6轨道倾角00升交点赤经00近地点幅角10°30°真近点角-30°120°表3bezier曲线控制点参数表当前第1页1 2 3