CN109459929A - 火星大气进入段纵向可达区生成的解析同伦法 - Google Patents
火星大气进入段纵向可达区生成的解析同伦法 Download PDFInfo
- Publication number
- CN109459929A CN109459929A CN201811470234.7A CN201811470234A CN109459929A CN 109459929 A CN109459929 A CN 109459929A CN 201811470234 A CN201811470234 A CN 201811470234A CN 109459929 A CN109459929 A CN 109459929A
- Authority
- CN
- China
- Prior art keywords
- formula
- parsing
- state
- homotopy
- condition
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
- G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
- G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
- G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
本发明公开的火星大气进入段纵向可达区生成的解析同伦法,属于深空探测领域。本发明实现方法包括如下步骤:建立火星大气进入段探测器纵向动力学模型;基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓出最大航程;基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最小航程;基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最大开伞高度;基于解析同伦法,根据延拓得到的最大航程、最小航程和最大开伞高度,通过构造同伦参数,延拓出纵向可达区。本发明能够避免协态初值猜测和内部点可达性验证,进而提高可达区生成效率。
Description
技术领域
本发明涉及一种火星大气进入段纵向可达区生成的解析同伦法,属于深空探测领域。
背景技术
火星大气进入段指的是探测器从进入火星大气直至降落伞展开的飞行阶段,火星大气进入段可达区分析对于任务设计、着陆点选取及风险评估等均具有重要意义。可达区求解涉及到进入段轨迹优化问题,可利用直接法或间接法求解。对于直接法,已有学者通过网格划分,结合数值优化方法求解大气进入段可达区。该方法需要对每一个网格点单独求取最优化问题,计算量大,求解依赖于初值猜测,求解过程的收敛性及可靠性难以保证。间接法的理论基础是庞特里亚金极小值原理,通过将原优化问题转换为一个等价的两点边值问题求解,由于极小值理论中的协状态缺乏物理意义,且优化问题对协状态初值比较敏感,故间接法存在协态初值猜测困难的问题。特别是对于可以横纵向解耦的小升阻比火星进入探测器,基于其纵向动力学方程的最优控制问题一般均属于bang-bang控制,初值收敛域小且对初值敏感。目前,同伦法在协状态初值猜测方面已有一定的理论基础,在求解小推力星际转移轨道最优问题方面已经有了较为深入的研究。在火星大气进入段轨迹优化方面,已有学者利用同伦法求解了最大开伞高度问题。
发明内容
针对火星大气进入段飞行器的飞行能力分析问题,本发明公开的火星大气进入段纵向可达区生成的解析同伦法要解决的技术问题是:基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓出纵向可达区(开伞点高度-航程剖面),所述同伦过程能够避免协态初值猜测和内部点可达性验证,进而提高可达区生成效率。
本发明目的是通过下述技术方案实现的。
本发明公开的火星大气进入段纵向可达区生成的解析同伦法,
本发明公开的火星大气进入段纵向可达区生成的解析同伦法,包括如下步骤:
步骤一、建立火星大气进入段探测器纵向动力学模型;
在火星惯性坐标系下,忽略火星自转,取探测器的纵向平面内运动状态为x=[r,V,γ,s]T,其中,r探测器质心到火星质心的距离,V为探测器速度大小,γ为飞行路径角,s为航程,则大气进入段无量纲的纵向动力学模型为:
式(1)中,τ为无量纲时间,σ为倾侧角。在无量纲化过程中,长度的量纲单位为火星半径R0,速度的无量纲单位为其中为火表引力加速度,μ为火星引力常数。时间的无量纲单位为角度的单位为弧度,不需要无量纲化处理。式(1)中,L和D分别为探测器受到的无量纲升力和阻力加速度,分别具有如下形式:
L=D·L/D (3)
式(1)中,B为探测器的弹道系数,L/D为探测器的升阻比,ρ为火星大气密度,采用如下指数模型:
式(4)中,ρ0为参考密度,r0为参考半径,h为探测器的飞行高度,hs为标高。
步骤二、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓出最大航程;
取最大航程的辅助优化问题如式(5)、式(6)所示。
问题寻找最优控制使得
满足约束:
其中,σmax和σmin分别为倾侧角可以取的最大值及最小值,x0为探测器在大气进入点的初始状态,Vf为探测器顺利开伞所需满足的末端速度约束。
根据庞特里亚金极小值原理,引入协状态λ=[λr λV λγ λs]T,问题对应的哈密顿函数为
问题的协状态微分方程为
问题的横截条件如下
由于中不显含时间项,故
由于0≤umax-u≤umax-umin,且终端条件只有速度约束,故将动力学方程按照常值倾侧角σmin正向积分直到满足终端速度约束为止,即得到问题对应的最优状态和飞行时间。
由式(10)知,将式(9)代入式(10),得终端速度协状态值
λV(τf)=0 (11)
根据式(11)和式(9)的协状态终端值、式(8)的协状态微分方程以及正向积分得到的最优状态和飞行时间逆向积分即可得到辅助优化问题对应的最优轨迹协状态值,即辅助优化问题的最优解是解析已知的。
由辅助优化问题延拓得到原最优问题解的关键在于构建同伦参数。将同伦参数ε置于性能指标中,构建如下的优化问题
问题寻找最优控制使得
满足约束:
问题的哈密顿函数为
由于中不显含时间,故成立,则
控制量u以一阶形式出现,故问题的最优解是bang-bang形式,最优解具有如下的形式
问题的协状态微分方程与问题Ps,max相同。横截条件如下
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(13)、式(11)、式(8)和(17)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件。
由式(14)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最大航程问题。构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1。从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最大航程问题的最优解,即延拓得到最大航程。
步骤三、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最小航程。
构造辅助约束优化问题如式(18)和式(19)所示。
问题寻找最优控制使得
满足约束:
问题的哈密顿函数如下
由于中不显含时间项,故
由于0≤u-umin≤umax-umin,且终端条件只有速度约束,故将动力学方程按照常值倾侧角σmax正向积分直到满足终端速度约束为止,即得到问题对应的最优状态和飞行时间,即辅助优化问题的最优解是解析已知的。
问题为拉格朗日问题,性能指标中只包含与状态变量无关的积分项,且中不显含时间项。故问题的横截条件与问题的横截条件相同,即
从式(22)所示的末端条件出发,根据式(8)的微分方程及积分得到的最优状态和飞行时间,逆向积分即得到该辅助优化问题所对应的协状态值。
构建如公式(24)所示的优化问题
问题寻找最优控制使得
满足约束:
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(8)、式(11)、式(22)、和式(23)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件。
由式(23)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最小航程问题。构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1。从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最小航程问题的最优解,即延拓得到最小航程。
步骤四、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最大开伞高度;
取性能指标为构建辅助优化问题如公式(25)、(26)所示。
问题寻找最优控制使
满足约束:
由于0≤u2,故问题的最优解为常值σ=π/2。获取最优状态的方法为从初始状态开始,以常值倾侧角正向积分飞行轨迹。直到满足末端速度天剑为止,得到辅助最优问题的最优状态和飞行时间,即辅助优化问题的最优解是解析已知的。
问题的协态微分方程与问题相同,横截条件与问题相同。根据协状态微分方程逆向积分即得到最优轨迹对应的协状态值。
为从辅助优化问题得到最大开伞高度问题的最优解,取性能指标为构建如公式(27)、(28)所示的优化问题。
问题寻找最优控制使得
满足约束:
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(8)、式(11)、式(22)、和式(28)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件。
由式(27)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最大开伞高度问题。构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1。从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最大开伞高度问题的最优解,即延拓得到最小航程。
步骤五、基于解析同伦法,根据步骤二延拓得到的最大航程、步骤三延拓得到的最小航程和步骤四延拓得到的最大开伞高度,通过构造同伦参数,延拓出纵向可达区,所述纵向可达区为开伞点高度-航程剖面。
取性能指标为Jsk=-r(τf),构建如公式(29)、(30)所示的优化模型。
问题RA:寻找使得
满足约束:
问题RA与问题Ph,max哈密顿函数相同,故问题RA的横截条件为
构造两个单调航程序列:{s0,...si,si+1,...sI}与{s0,...sj,sj+1,...sJ},满足s0为问题Ph,max对应的航程,sI为问题Ps,min对应的最小航程,sJ为问题Ps,max对应的最大航程,分别从问题Ph,max的最优解出发,按照序列依次求解问题RA,其中第k次的最优解zk作为第k+1次问题求解初值,延拓出纵向可达区,所述纵向可达区为开伞点高度-航程剖面。
有益效果:
1、本发明公开的火星大气进入段纵向可达区生成的解析同伦法,在构造最优解已知的辅助优化问题基础上,通过构造合适的同伦参数,能够快速求解最大/小航程问题及最大开伞高度问题,有效避免最优问题协态初值猜测过程,同伦求解过程稳定。
2、本发明公开的火星大气进入段纵向可达区生成的解析同伦法,根据步骤二延拓得到的最大航程、步骤三延拓得到的最小航程和步骤四延拓得到的最大开伞高度,通过构造合适的同伦参数,延拓出纵向可达区,能够避免对可达区内部点可达性的验证,进而提高可达区生成效率。
附图说明
图1为火星大气进入段纵向可达区生成的解析同伦法流程图。
图2为问题中初始协状态值随参数ε的变化。
图3为问题中航程随参数ε的变化关系。
图4为问题中初始协状态值随参数ε的变化。
图5为问题中航程随参数ε的变化关系。
图6为问题中初始协状态值随参数ε的变化。
图7为问题中航程随参数ε的变化关系。
图8为协状态初值随航程的变化关系。
图9为不同初始航迹角下的纵向可达区。
具体实施方式
为了更好的说明本发明的目的和优点,下面结合附图和实施实例对发明内容做进一步说明。
本实例为针对火星大气进入段纵向可达区求解问题,利用解析同伦法,在求解最大/小航程和最大开伞高度问题的基础上,进一步利用延拓的方法得到纵向可达区。
火星大气进入段纵向可达区生成的解析同伦法,如图1所示,具体步骤如下:
步骤一、建立火星大气进入段探测器纵向动力学模型;
在火星惯性坐标系下,忽略火星自转,取探测器的纵向平面内运动状态为x=[r,V,γ,s]T,其中,r探测器质心到火星质心的距离,V为探测器速度大小,γ为飞行路径角,s为航程,则大气进入段无量纲的纵向动力学模型为:
式(1)中,τ为无量纲时间,σ为倾侧角。在无量纲化过程中,长度的量纲单位为火星半径R0,速度的无量纲单位为其中为火表引力加速度,μ为火星引力常数。时间的无量纲单位为角度的单位为弧度,不需要无量纲化处理。式(1)中,L和D分别为探测器受到的无量纲升力和阻力加速度,分别具有如下形式:
L=D·L/D (3)
式(1)中,B为探测器的弹道系数,L/D为探测器的升阻比,ρ为火星大气密度,采用如下指数模型:
式(4)中,ρ0为参考密度,r0为参考半径,h为探测器的飞行高度,hs为标高。
步骤二、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓出最大航程;
取最大航程的辅助优化问题如式(5)、式(6)所示。
问题寻找最优控制使得
满足约束:其中,σmax和σmin分别为倾侧角可以取的最大值及最小值,x0为探测器在大气进入点的初始状态,Vf为探测器顺利开伞所需满足的末端速度约束。
根据庞特里亚金极小值原理,引入协状态λ=[λr λV λγ λs]T,问题对应的哈密顿函数为
问题的协状态微分方程为
问题的横截条件如下
由于中不显含时间项,故
由于0≤umax-u≤umax-umin,且终端条件只有速度约束,故将动力学方程按照常值倾侧角σmin正向积分直到满足终端速度约束为止,即得到问题对应的最优状态和飞行时间。
由式(10)知,将式(9)代入式(10),得终端速度协状态值
λV(τf)=0 (11)
根据式(11)和式(9)的协状态终端值、式(8)的协状态微分方程以及正向积分得到的最优状态和飞行时间逆向积分即可得到辅助优化问题对应的最优轨迹协状态值,即辅助优化问题的最优解是解析已知的。
由辅助优化问题延拓得到原最优问题解的关键在于构建同伦参数。将同伦参数ε置于性能指标中,构建如下的优化问题
问题寻找最优控制使得
满足约束:
问题的哈密顿函数为
由于中不显含时间,故成立,则
控制量u以一阶形式出现,故问题的最优解是bang-bang形式,最优解具有如下的形式
问题的协状态微分方程与问题Ps,max相同。横截条件如下
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(13)、式(11)、式(8)和(17)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件。
由式(14)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最大航程问题。构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1。从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最大航程问题的最优解,即延拓得到最大航程。
步骤三、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最小航程。
构造辅助约束优化问题如式(18)和式(19)所示。
问题寻找最优控制使得
满足约束:
问题的哈密顿函数如下
由于中不显含时间项,故
由于0≤u-umin≤umax-umin,且终端条件只有速度约束,故将动力学方程按照常值倾侧角σmax正向积分直到满足终端速度约束为止,即得到问题对应的最优状态和飞行时间,即辅助优化问题的最优解是解析已知的。
问题为拉格朗日问题,性能指标中只包含与状态变量无关的积分项,且中不显含时间项。故问题的横截条件与问题的横截条件相同,即
从式(22)所示的末端条件出发,根据式(8)的微分方程及积分得到的最优状态和飞行时间,逆向积分即得到该辅助优化问题所对应的协状态值。
构建如公式(24)所示的优化问题
问题寻找最优控制使得
满足约束:
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(8)、式(11)、式(22)、和式(23)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件。
由式(23)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最小航程问题。构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1。从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最小航程问题的最优解,即延拓得到最小航程。
步骤四、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最大开伞高度;
取性能指标为构建辅助优化问题如公式(25)、(26)所示。
问题寻找最优控制使
满足约束:
由于0≤u2,故问题的最优解为常值σ=π/2。获取最优状态的方法为从初始状态开始,以常值倾侧角正向积分飞行轨迹。直到满足末端速度天剑为止,得到辅助最优问题的最优状态和飞行时间,即辅助优化问题的最优解是解析已知的。
问题的协态微分方程与问题相同,横截条件与问题相同。根据协状态微分方程逆向积分即得到最优轨迹对应的协状态值。
为从辅助优化问题得到最大开伞高度问题的最优解,取性能指标为构建如公式(27)、(28)所示的优化问题。
问题寻找最优控制使得
满足约束:
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(8)、式(11)、式(22)、和式(28)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件。
由式(27)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最大开伞高度问题。构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1。从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最大开伞高度问题的最优解,即延拓得到最小航程。
步骤五、基于解析同伦法,根据步骤二延拓得到的最大航程、步骤三延拓得到得最小航程和步骤四延拓得到的最大开伞高度,通过构造同伦参数,延拓出纵向可达区,所述纵向可达区为开伞点高度-航程剖面。
取性能指标为构建如公式(29)、(30)所示的优化模型。
问题RA:寻找使得
满足约束:
问题RA与问题Ph,max哈密顿函数相同,故问题RA的横截条件为
构造两个单调航程序列:{s0,...si,si+1,...sI}与{s0,...sj,sj+1,...sJ},满足s0为问题Ph,max对应的航程,sI为问题Ps,min对应的最小航程,sJ为问题Ps,max对应的最大航程,分别从问题Ph,max的最优解出发,按照序列依次求解问题RA,其中第k次的最优解zk作为第k+1次问题求解初值,延拓出纵向可达区,所述纵向可达区为开伞点高度-航程剖面。为了保证探测器的安全,约束可达区的最低高度为6km。
步骤六、在MATLAB环境下对上述算法进行仿真分析;
探测器的物理参数如下:质量m为2800kg,参考面积Sref为15.9m2,升阻比L/D为0.24,阻力系数CD为1.45,初始条件及开伞条件具体如表1所示。为了给制导控制系统留有控制余量,取最大倾侧角为σmax为120°,最小倾侧角σmin为30°。在同伦延拓过程中,以下三个优化子问题的ε均从0开始,其中ε=0代表辅助优化问题的解。
表1探测器初始条件及开伞条件
参数 | r<sub>0</sub> | V<sub>0</sub> | γ<sub>0</sub> | s<sub>0</sub> | V<sub>f</sub> |
数值 | 3.5222×10<sup>6</sup> | 5800 | -0.2182 | 0 | 450 |
单位 | m | m/s | rad | m | m/s |
具体仿真结果如下。
对于问题协状态初值随参数ε:0→1的变化如图2所示。由图2可知,协状态随参数ε的增加连续变化,一定程度上代表了求解过程的稳定性,证明了该同伦求解过程的可行性。相比于σ≡σmin的进入轨迹,航程在ε由0变化到1时增加了约150m,出现这种现象的主要原因在于飞行时间的增加,相比于ε由0.1变化到1的绝大部分过程,最大航程问题对应的末端航程仅增加了3m,如图3所示。故在任务的快速分析设计中,可用σ≡σmin的进入轨迹近似代替最大航程轨迹来分析最大航程。
对于问题协状态初值随参数ε:0→1的变化分别如图4所示。由图4可知,协状态随参数ε的增加连续变化,一定程度上代表了求解过程的稳定性,证明了该同伦求解过程的可行性。
最小航程轨迹比倾侧角保持为σ≡σmax的进入轨迹末端航程高135m,如图5所示,这种现象应该是辅助问题求解时采取0.5s步长,积分精度较低引起的。在任务的快速分析设计中,可用σ≡σmax的进入轨迹近似代替最小航程轨迹。
对于问题协状态初值随参数ε:0→1的变化分别如图6所示。由图6可知,协状态随参数ε的增加连续变化,一定程度上代表了求解过程的稳定性,证明了该同伦求解过程的可行性。由图7可知,在ε由0变化到1的过程中,相比于辅助优化问题,开伞高度由2.65km增至10.23km。
另外选取两个不同的初始航迹角:γ0=12°和γ0=13°。在分别求取三个子优化问题的基础上分别求解问题RA。其中γ0=12°和γ0=13°的子问题既可以分别从辅助优化问题出发求解,也可以从γ0=12.5°的三个子问题最优解出发,构造{γk1,γk2,...,γkN}的序列,采用同伦法求解,其中γk1=12.5°,γkN为12°或13°,也可以是其他任务设计时感兴趣的初始航迹角。其中,γ0=12.5°时,求解问题RA得到的初始协状态值随航程的变化关系如图8所示,均为光滑曲线,一定程度上代表了求解过程的稳定性和方法的可行性。三个不同初始航迹角下,纵向可达区分别如图9所示。可以看出,进入角更陡时,航程范围更窄,末端可达最大开伞高度更高。
综上,本发明所提出的火星大气进入段纵向可达区生成的解析同伦法,可从最优解已知的辅助优化问题出发,最终得到火星大气进入段纵向可达区。
以上所述的具体描述,对发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (6)
1.火星大气进入段纵向可达区生成的解析同伦法,其特征在于:包括如下步骤,
步骤一、建立火星大气进入段探测器纵向动力学模型;
步骤二、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓出最大航程;
步骤三、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最小航程;
步骤四、基于解析同伦法,从最优解解析已知的辅助问题出发,通过构造同伦参数,延拓得到最大开伞高度;
步骤五、基于解析同伦法,根据步骤二延拓得到的最大航程、步骤三延拓得到的最小航程和步骤四延拓得到的最大开伞高度,通过构造同伦参数,延拓出纵向可达区,所述纵向可达区为开伞点高度-航程剖面。
2.如权利要求1所述的火星大气进入段纵向可达区生成的解析同伦法,其特征在于:步骤一具体实现方法为,
在火星惯性坐标系下,忽略火星自转,取探测器的纵向平面内运动状态为x=[r,V,γ,s]T,其中,r探测器质心到火星质心的距离,V为探测器速度大小,γ为飞行路径角,s为航程,则大气进入段无量纲的纵向动力学模型为:
式(1)中,τ为无量纲时间,σ为倾侧角;在无量纲化过程中,长度的量纲单位为火星半径R0,速度的无量纲单位为其中为火表引力加速度,μ为火星引力常数;时间的无量纲单位为角度的单位为弧度,不需要无量纲化处理;式(1)中,L和D分别为探测器受到的无量纲升力和阻力加速度,分别具有如下形式:
L=D·L/D (3)
式(1)中,B为探测器的弹道系数,L/D为探测器的升阻比,ρ为火星大气密度,采用如下指数模型:
式(4)中,ρ0为参考密度,r0为参考半径,h为探测器的飞行高度,hs为标高。
3.如权利要求2所述的火星大气进入段纵向可达区生成的解析同伦法,其特征在于:步骤二具体实现方法为,
取最大航程的辅助优化问题如式(5)、式(6)所示;
问题寻找最优控制使得
满足约束:
其中,σmax和σmin分别为倾侧角可以取的最大值及最小值,x0为探测器在大气进入点的初始状态,Vf为探测器顺利开伞所需满足的末端速度约束;
根据庞特里亚金极小值原理,引入协状态λ=[λr λV λγ λs]T,问题对应的哈密顿函数为
问题的协状态微分方程为
问题的横截条件如下
由于中不显含时间项,故
由于0≤umax-u≤umax-umin,且终端条件只有速度约束,故将动力学方程按照常值倾侧角σmin正向积分直到满足终端速度约束为止,即得到问题对应的最优状态和飞行时间;
由式(10)知,将式(9)代入式(10),得终端速度协状态值
λV(τf)=0 (11)
根据式(11)和式(9)的协状态终端值、式(8)的协状态微分方程以及正向积分得到的最优状态和飞行时间逆向积分即可得到辅助优化问题对应的最优轨迹协状态值,即辅助优化问题的最优解是解析已知的;
由辅助优化问题延拓得到原最优问题解的关键在于构建同伦参数;将同伦参数ε置于性能指标中,构建如下的优化问题
问题寻找最优控制使得
满足约束:
问题的哈密顿函数为
由于中不显含时间,故成立,则
控制量u以一阶形式出现,故问题的最优解是bang-bang形式,最优解具有如下的形式
问题的协状态微分方程与问题Ps,max相同;横截条件如下
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(13)、式(11)、式(8)和(17)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件;
由式(14)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最大航程问题;构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1;从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最大航程问题的最优解,即延拓得到最大航程。
4.如权利要求3所述的火星大气进入段纵向可达区生成的解析同伦法,其特征在于:步骤三具体实现方法为,
构造辅助约束优化问题如式(18)和式(19)所示;
问题寻找最优控制使得
满足约束:
问题的哈密顿函数如下
由于中不显含时间项,故
由于0≤u-umin≤umax-umin,且终端条件只有速度约束,故将动力学方程按照常值倾侧角σmax正向积分直到满足终端速度约束为止,即得到问题对应的最优状态和飞行时间,即辅助优化问题的最优解是解析已知的;
问题为拉格朗日问题,性能指标中只包含与状态变量无关的积分项,且中不显含时间项;故问题的横截条件与问题的横截条件相同,即
从式(22)所示的末端条件出发,根据式(8)的微分方程及积分得到的最优状态和飞行时间,逆向积分即得到该辅助优化问题所对应的协状态值;
构建如公式(24)所示的优化问题;
问题寻找最优控制使得
满足约束:
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(8)、式(11)、式(22)、和式(23)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件;
由式(23)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最小航程问题;构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1;从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最小航程问题的最优解,即延拓得到最小航程。
5.如权利要求4所述的火星大气进入段纵向可达区生成的解析同伦法,其特征在于:步骤四具体实现方法为,
取性能指标为构建辅助优化问题如公式(25)、(26)所示;
问题寻找最优控制使
满足约束:
由于0≤u2,故问题的最优解为常值σ=π/2;获取最优状态的方法为从初始状态开始,以常值倾侧角正向积分飞行轨迹;直到满足末端速度天剑为止,得到辅助最优问题的最优状态和飞行时间,即辅助优化问题的最优解是解析已知的;
问题的协态微分方程与问题相同,横截条件与问题相同;根据协状态微分方程逆向积分即得到最优轨迹对应的协状态值;
为从辅助优化问题得到最大开伞高度问题的最优解,取性能指标为构建如公式(27)、(28)所示的优化问题;问题寻找最优控制使得
满足约束:
问题的求解转换为求解两点边值问题寻找由协状态初值和终端飞行时间组成的变量z,满足式(8)、式(11)、式(22)、和式(28)描述的边值条件和微分方程条件,式(15)中的哈密顿函数条件及式(16)中的控制切换条件;
由式(27)知,当ε=0时,问题等价于问题当ε=1时,问题等价于最大开伞高度问题;构造递增序列{ε1,ε2,...εk,εk+1,...εN},其中ε1=0,εN=1;从ε1开始,以第k步的解zk作为第k+1步的问题的初值,依次求解两点边值问题,即能够最终得到原最大开伞高度问题的最优解,即延拓得到最小航程。
6.如权利要求5所述的火星大气进入段纵向可达区生成的解析同伦法,其特征在于:步骤五具体实现方法为,
取性能指标为构建如公式(29)、(30)所示的优化模型;
问题RA:寻找使得
满足约束:
问题RA与问题Ph,max哈密顿函数相同,故问题RA的横截条件为
构造两个单调航程序列:{s0,...si,si+1,...sI}与{s0,...sj,sj+1,...sJ},满足s0为问题Ph,max对应的航程,sI为问题Ps,min对应的最小航程,sJ为问题Ps,max对应的最大航程,分别从问题Ph,max的最优解出发,按照序列依次求解问题RA,其中第k次的最优解zk作为第k+1次问题求解初值,延拓出纵向可达区,所述纵向可达区为开伞点高度-航程剖面。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811470234.7A CN109459929B (zh) | 2018-12-04 | 2018-12-04 | 火星大气进入段纵向可达区生成的解析同伦法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201811470234.7A CN109459929B (zh) | 2018-12-04 | 2018-12-04 | 火星大气进入段纵向可达区生成的解析同伦法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109459929A true CN109459929A (zh) | 2019-03-12 |
CN109459929B CN109459929B (zh) | 2020-06-16 |
Family
ID=65612330
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201811470234.7A Active CN109459929B (zh) | 2018-12-04 | 2018-12-04 | 火星大气进入段纵向可达区生成的解析同伦法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109459929B (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113377128A (zh) * | 2021-06-10 | 2021-09-10 | 北京空天技术研究所 | 一种飞行器可达区域估算方法 |
CN114280934A (zh) * | 2021-12-15 | 2022-04-05 | 北京航天自动控制研究所 | 一种可重复使用运载火箭全程轨迹规划方法 |
CN115618171A (zh) * | 2022-06-06 | 2023-01-17 | 北京理工大学 | 一种基于同伦算法的推进剂燃烧平衡产物求解方法 |
CN116880527A (zh) * | 2023-07-20 | 2023-10-13 | 中国空气动力研究与发展中心空天技术研究所 | 高超声速飞行器跳跃滑翔飞行射程最大的控制方法及系统 |
CN118092502A (zh) * | 2024-04-28 | 2024-05-28 | 北京航空航天大学 | 一种基于摄动、同伦理论的解析轨迹预测方法及系统 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20100006704A1 (en) * | 2008-07-08 | 2010-01-14 | Thales | Method for Lightening the Weight of Fuel Stowed Onboard During an Interplanetary Mission |
CN106295000A (zh) * | 2016-08-10 | 2017-01-04 | 北京理工大学 | 一种考虑不确定性影响的火星大气进入段轨迹优化方法 |
CN106742069A (zh) * | 2016-12-29 | 2017-05-31 | 北京理工大学 | 一种火星大气进入段最优预测制导方法 |
CN107323691A (zh) * | 2017-07-04 | 2017-11-07 | 北京理工大学 | 一种多约束火星大气进入预测制导方法 |
CN107942673A (zh) * | 2017-12-11 | 2018-04-20 | 北京理工大学 | 一种开伞点高度跟踪的火星大气进入段解析制导方法 |
-
2018
- 2018-12-04 CN CN201811470234.7A patent/CN109459929B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20100006704A1 (en) * | 2008-07-08 | 2010-01-14 | Thales | Method for Lightening the Weight of Fuel Stowed Onboard During an Interplanetary Mission |
CN106295000A (zh) * | 2016-08-10 | 2017-01-04 | 北京理工大学 | 一种考虑不确定性影响的火星大气进入段轨迹优化方法 |
CN106742069A (zh) * | 2016-12-29 | 2017-05-31 | 北京理工大学 | 一种火星大气进入段最优预测制导方法 |
CN107323691A (zh) * | 2017-07-04 | 2017-11-07 | 北京理工大学 | 一种多约束火星大气进入预测制导方法 |
CN107942673A (zh) * | 2017-12-11 | 2018-04-20 | 北京理工大学 | 一种开伞点高度跟踪的火星大气进入段解析制导方法 |
Non-Patent Citations (4)
Title |
---|
AMITABH SARAF,ET AL.: "landing footprint computation for entry vehicles", 《AIAA GUIDANCE,NAVIGATION,AND CONTROL CONFERENCE AND EXHIBIT,PROVIDENCE,RHODE ISLAND,USA》 * |
JIANG FANGHUA,ET AL.: "Practical techniques for low-thrust trajectory optimization with homotopic approach", 《JOURNAL OF GUIDANCE,CONTROL,AND DYNAMICS》 * |
JOEL BENITO,ET AL.: "Reachable and Controllable Sets for Planetary Entry and Landing", 《JOURNAL OF GUIDANCE,CONTROL,AND DYNAMICS》 * |
ZHENG YIYU,ET AL.: "Indirect trajectory optimization for Mars entry with maximum terminal altitude", 《JORRNAL OF SPACECRAFT AND ROCKETS》 * |
Cited By (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113377128A (zh) * | 2021-06-10 | 2021-09-10 | 北京空天技术研究所 | 一种飞行器可达区域估算方法 |
CN113377128B (zh) * | 2021-06-10 | 2022-12-09 | 北京空天技术研究所 | 一种飞行器可达区域估算方法 |
CN114280934A (zh) * | 2021-12-15 | 2022-04-05 | 北京航天自动控制研究所 | 一种可重复使用运载火箭全程轨迹规划方法 |
CN114280934B (zh) * | 2021-12-15 | 2023-08-15 | 北京航天自动控制研究所 | 一种可重复使用运载火箭全程轨迹规划方法 |
CN115618171A (zh) * | 2022-06-06 | 2023-01-17 | 北京理工大学 | 一种基于同伦算法的推进剂燃烧平衡产物求解方法 |
CN115618171B (zh) * | 2022-06-06 | 2023-10-24 | 北京理工大学 | 一种基于同伦算法的推进剂燃烧平衡产物求解方法 |
CN116880527A (zh) * | 2023-07-20 | 2023-10-13 | 中国空气动力研究与发展中心空天技术研究所 | 高超声速飞行器跳跃滑翔飞行射程最大的控制方法及系统 |
CN116880527B (zh) * | 2023-07-20 | 2024-02-23 | 中国空气动力研究与发展中心空天技术研究所 | 高超声速飞行器跳跃滑翔飞行射程最大的控制方法及系统 |
CN118092502A (zh) * | 2024-04-28 | 2024-05-28 | 北京航空航天大学 | 一种基于摄动、同伦理论的解析轨迹预测方法及系统 |
Also Published As
Publication number | Publication date |
---|---|
CN109459929B (zh) | 2020-06-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109459929A (zh) | 火星大气进入段纵向可达区生成的解析同伦法 | |
CN110930770B (zh) | 一种基于管制意图和飞机性能模型的四维航迹预测方法 | |
CN110908396B (zh) | 可重复使用运载器全阶段再入返回制导方法 | |
CN108958285B (zh) | 一种基于分解思想的高效多无人机协同航迹规划方法 | |
CN108614420B (zh) | 基于非线性规划的星簇级卫星容错控制方法 | |
CN105953800A (zh) | 一种无人飞行器航迹规划栅格空间划分方法 | |
CN116257082B (zh) | 多无人机分布式主动协同探测方法 | |
CN106768123A (zh) | 一种无人直升机燃油预估方法 | |
CN104648695A (zh) | 一种基于倾侧角可用性的再入走廊最优规划方法 | |
CN108491650B (zh) | 一种基于智能学习的火星进入终端状态高效评估方法 | |
CN103198187A (zh) | 基于微分修正的深空探测器的轨道设计方法 | |
Malaek et al. | Dynamic based cost functions for TF/TA flights | |
CN103324993A (zh) | 一种基于多飞行器协同作战的轨迹优化方法 | |
CN116187199A (zh) | 融合智能技术的非定常气动力建模方法 | |
CN106873615A (zh) | 应急返场着陆速度指令集设计方法 | |
CN113093789B (zh) | 一种基于路径点优选的飞行器禁飞区规避轨迹规划方法 | |
CN108168558A (zh) | 应用于河流目标搜索任务的无人机航迹规划算法 | |
Tian et al. | Wake encounter simulation and flight validation with UAV close formation flight | |
CN115542939A (zh) | 空空导弹分布式协同中制导律分析方法及其导引系统 | |
CN115454145A (zh) | 飞行速度不可控条件下高速飞行器协同制导方法及系统 | |
Pokhrel et al. | Extremum Seeking in Nature: Examination of Soaring Birds Flights Maneuver | |
Cheng et al. | Cross-cycle iterative unmanned aerial vehicle reentry guidance based on reinforcement learning | |
Li et al. | Modeling and control for morphing unmanned aircraft vehicle | |
Huang et al. | On the 3D Track Planning for Electric Power Inspection Based on the Improved Ant Colony Optimization and A∗ Algorithm | |
CN112364433B (zh) | 一种高效的固定翼飞行器飞行动力学模型配平方法 |
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 |