CN108846162B - 一种基于bezier曲线的连续推力机动轨道设计方法 - Google Patents
一种基于bezier曲线的连续推力机动轨道设计方法 Download PDFInfo
- Publication number
- CN108846162B CN108846162B CN201810433526.7A CN201810433526A CN108846162B CN 108846162 B CN108846162 B CN 108846162B CN 201810433526 A CN201810433526 A CN 201810433526A CN 108846162 B CN108846162 B CN 108846162B
- Authority
- CN
- China
- Prior art keywords
- bezier curve
- track
- orbit
- point
- thrust
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Image Generation (AREA)
Abstract
本发明公开了一种基于bezier曲线的连续推力机动轨道设计方法,包括以下步骤:1)计算轨道机动问题的起止点约束及推力表达式;2)选用7点bezier曲线进行轨道设计;3)设定复合函数P=B1H1+B2H2,得复合函数的一阶导数及二阶导数;4)由步骤3)得满足轨道机动任务的复合函数P约束;5)得bezier曲线B1及bezier曲线B2的方程;6)将bezier曲线B1和bezier曲线B2进行复合,得最终的机动轨道,该方法得到的轨道具有普适性。
Description
技术领域
本发明属于轨道设计领域,涉及一种基于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′2
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)=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·θ14
a2=-6·θ8+30·θ9-60·θ10+60·θ11-30·θ12+6·θ13
a3=15·θ8-60·θ9+90·θ10-60·θ11+15·θ12
a4=-20·θ8+60·θ9-60·θ10+20·θ11
a5=15·θ8-30·θ9+15·θ10
a6=-6·θ8+6·θ9
a7=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″″分别为:
通过航迹角γ给出轨道速度方向,其中,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)=0
5)对于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和bezi er曲线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·θ14
a2=-6·θ8+30·θ9-60·θ10+60·θ11-30·θ12+6·θ13
a3=15·θ8-60·θ9+90·θ10-60·θ11+15·θ12
a4=-20·θ8+60·θ9-60·θ10+20·θ11
a5=15·θ8-30·θ9+15·θ10
a6=-6·θ8+6·θ9
a7=1·θ8。
实施例一
本发明对某一轨道机动问题实例实施的结果为:优化结果累计推力为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 (3)
1.一种基于bezier曲线的连续推力机动轨道设计方法,其特征在于,包括以下步骤:
1)计算轨道机动问题的起止点约束及机动轨道推力表达式;
2)选用n点bezier曲线进行轨道设计,其中,设Pj为bezier曲线的第j个控制点,t为bezier曲线的比例因子,1≤j≤n;
3)设定复合函数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)=H″1
(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个控制点的坐标参数,r为初始轨道的椭圆轨道极坐标下的轨道矢径,θ为轨道的相位角,在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曲线的连续推力机动轨道设计。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810433526.7A CN108846162B (zh) | 2018-05-08 | 2018-05-08 | 一种基于bezier曲线的连续推力机动轨道设计方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810433526.7A CN108846162B (zh) | 2018-05-08 | 2018-05-08 | 一种基于bezier曲线的连续推力机动轨道设计方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108846162A CN108846162A (zh) | 2018-11-20 |
CN108846162B true CN108846162B (zh) | 2022-11-22 |
Family
ID=64212829
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810433526.7A Active CN108846162B (zh) | 2018-05-08 | 2018-05-08 | 一种基于bezier曲线的连续推力机动轨道设计方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108846162B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111639403A (zh) * | 2020-05-28 | 2020-09-08 | 上海理工大学 | 一种基于Bézier曲线的燃料棒形状优化方法 |
CN112269965B (zh) * | 2020-08-10 | 2024-04-05 | 中国北方车辆研究所 | 一种非完整约束条件下的连续曲率路径优化方法 |
CN114742975B (zh) * | 2022-06-10 | 2022-09-02 | 西南交通大学 | 一种车载图像铁轨曲线建模方法 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105676640A (zh) * | 2016-01-13 | 2016-06-15 | 南京航空航天大学 | 基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
FR3009104B1 (fr) * | 2013-07-29 | 2015-08-21 | Snecma | Procede de modelisation d'une pale d'une helice non-carenee |
CN104965519B (zh) * | 2015-06-10 | 2018-08-10 | 北京理工大学 | 一种基于贝塞尔曲线的带落角约束的末制导方法 |
CN107490965B (zh) * | 2017-08-21 | 2020-02-07 | 西北工业大学 | 一种空间自由漂浮机械臂的多约束轨迹规划方法 |
-
2018
- 2018-05-08 CN CN201810433526.7A patent/CN108846162B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105676640A (zh) * | 2016-01-13 | 2016-06-15 | 南京航空航天大学 | 基于贝塞尔曲线的涡扇发动机加速过程控制律设计方法 |
Non-Patent Citations (1)
Title |
---|
飞行器轨迹优化的造型及多分辨率方法研究;张松;《中国博士学位论文全文数据库电子期刊 工程科技II辑》;20170115;第2017年卷(第1期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN108846162A (zh) | 2018-11-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108846162B (zh) | 一种基于bezier曲线的连续推力机动轨道设计方法 | |
Zhang et al. | Design of ship course-keeping autopilot using a sine function-based nonlinear feedback technique | |
CN110018687B (zh) | 基于强化学习方法的无人水面船最优轨迹跟踪控制方法 | |
CN109283842B (zh) | 一种无人艇航迹跟踪智能学习控制方法 | |
Shao et al. | Active disturbance rejection based trajectory linearization control for hypersonic reentry vehicle with bounded uncertainties | |
CN102929283B (zh) | 基于sdre的再入飞行器自适应最优滑模姿态控制方法 | |
CN110221541B (zh) | 伺服***中前馈控制器的前馈系数获取装置、方法 | |
CN108717265B (zh) | 一种基于控制变量参数化的无人飞行器巡航跟踪控制***及控制方法 | |
CN106788036A (zh) | 一种直流电机的改进型自抗扰位置控制器及其设计方法 | |
CN111650832B (zh) | 一种水下多足步行机器人机械足姿态跟踪控制方法 | |
CN105867139B (zh) | 一种模型参考自适应的飞行器鲁棒控制方法 | |
Mu et al. | Tracking control of podded propulsion unmanned surface vehicle with unknown dynamics and disturbance under input saturation | |
CN110162073B (zh) | 一种适应边界约束的分段正弦姿态机动轨迹规划方法 | |
CN109782759B (zh) | 一种轮式移动机器人的近似解耦、快速轨迹跟踪控制方法 | |
Yuan et al. | Terminal sliding mode fuzzy control based on multiple sliding surfaces for nonlinear ship autopilot systems | |
Wu et al. | Adaptive robust course-tracking control of time-varying uncertain ships with disturbances | |
Liu et al. | High AOA decoupling control for aircraft based on ADRC | |
CN107728596A (zh) | 一种内燃机车电控***的模糊控制方法 | |
CN113110527B (zh) | 一种自主水下航行器有限时间路径跟踪的级联控制方法 | |
Shojaei | An adaptive output feedback proportional-integral-derivative controller for n-link type (m, s) electrically driven mobile manipulators | |
CN111007867B (zh) | 一种可预设调整时间的高超声速飞行器姿态控制设计方法 | |
CN104122793A (zh) | 一种满足预设性能的导弹过载控制方法 | |
CN108536164B (zh) | 挠性航天器无角速度测量的姿态控制方法 | |
CN112363538B (zh) | 一种速度信息不完备下的auv区域跟踪控制方法 | |
Kamath et al. | Vision Augmented 3 DoF Quadrotor Control using a Non-singular Fast-terminal Sliding Mode Modified Super-twisting Controller |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |