中文核心期刊

中国科学引文数据库(CSCD)来源期刊

中国高校优秀科技期刊

中国宇航学会深空探测技术专业委员会会刊

高级检索

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

行星着陆动力下降轨迹优化的中心差分凸化方法

李明翔 泮斌峰

李明翔, 泮斌峰. 行星着陆动力下降轨迹优化的中心差分凸化方法[J]. 深空探测学报(中英文), 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
引用本文: 李明翔, 泮斌峰. 行星着陆动力下降轨迹优化的中心差分凸化方法[J]. 深空探测学报(中英文), 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
LI Mingxiang, PAN Binfeng. Central Difference Convexification Method for Soft-Landing Trajectory Optimization in Planetary Powered Descending Phase[J]. Journal of Deep Space Exploration, 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
Citation: LI Mingxiang, PAN Binfeng. Central Difference Convexification Method for Soft-Landing Trajectory Optimization in Planetary Powered Descending Phase[J]. Journal of Deep Space Exploration, 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079

行星着陆动力下降轨迹优化的中心差分凸化方法

doi: 10.15982/j.issn.2096-9287.2021.20200079
基金项目: 国防科技重点实验室基金资助项目(6142210200312)
详细信息
    作者简介:

    李明翔(1996-),男,硕士研究生,主要研究方向:飞行器轨迹优化与制导。通讯地址:陕西省西安市友谊西路127号西北工业大学航天学院(710072)电话:18829571869 E-mail:limingxiang@mail.nwpu.edu.cn

    通讯作者:

    泮斌峰(1981– ),男,教授,博士生导师,主要研究方向:航天器计算制导控制、机器学习、轨迹优化、深空探测技术。本文通讯作者。通讯地址:陕西省西安市友谊西路127号西北工业大学航天学院(710072)电话:18691959798 Email:panbinfeng@nwpu.edu.cn

  • ● The workload of code programming and formula derivation can be reduced and the generality of algorithm can be improved by using pre-labeled central difference method instead of analytically computing the Jacobi matrix. ● mf-tf curves of optimal trajectories can be fitted by a cluster of functions which can be used to approximately estimate the optimal terminal time. ● Simulation results show that direct linearization of dynamic constraints is better than customized convexification method with variable substitution for planet powered landing.
  • 中图分类号: V467.3

Central Difference Convexification Method for Soft-Landing Trajectory Optimization in Planetary Powered Descending Phase

  • 摘要: 针对行星定点软着陆实时在线制导的任务需求,设计了基于序列凸优化的动力下降燃料最优轨迹求解算法。算法采用预标记的中心差分法线性化动力学方程,并提出将性能指标相对偏差作为收敛终止条件,能够快速生成燃料最优轨迹。此外,在分析最优轨迹簇剩余燃料和终端时间关系的基础上给出其拟合函数,作为最优终端时间的近似估算,以减少算法求解未知变量的维数。数值仿真结果表明,与对动力学方程先进行变量代换再线性化的传统凸化方法相比,该算法对初始猜想不敏感,收敛性好且终端误差较小。
    Highlights
    ● The workload of code programming and formula derivation can be reduced and the generality of algorithm can be improved by using pre-labeled central difference method instead of analytically computing the Jacobi matrix. ● mf-tf curves of optimal trajectories can be fitted by a cluster of functions which can be used to approximately estimate the optimal terminal time. ● Simulation results show that direct linearization of dynamic constraints is better than customized convexification method with variable substitution for planet powered landing.
  • 图  1  性能指标及其绝对偏差和相对偏差随迭代次数变化图

    Fig.  1  Performance index and its absolute deviation,relative deviation versus iterations

    图  2  求解表1算例500次耗时统计

    Fig.  2  Computing time statistic of solving the case in Table 1 500 times

    图  3  发动机3个方向推力曲线

    Fig.  3  Thrust component in three directions

    图  4  飞行器速度曲线

    Fig.  4  Velocity profiles of optimal trajectory

    图  5  飞行器质量曲线

    Fig.  5  Mass profile of optimal trajectory.

    图  6  飞行器位置相对于着陆点变化曲线

    Fig.  6  Position profiles of optimal trajectory

    图  7  tf = 60 s轨迹随迭代次数变化

    Fig.  7  Trajectory versus iterations for case of tf=60 s

    图  8  tf = 60 s优化轨迹和推力矢量

    Fig.  8  Optimal trajectory and thrust for case of t f = 60 s

    图  9  tf∈[50,190]s最优轨迹簇

    Fig.  9  Optimal trajectories at tf∈[50,190]s

    图  10  不同tf下最优轨迹对应的剩余燃料质量mf

    Fig.  10  mf of optimal trajectories at different tf

    图  11  不同离散点个数对应的mf-tf曲线

    Fig.  11  mf-tf curves under different number of discrete points

    图  12  不同离散点个数与动力学约束下的计算耗时

    Fig.  12  Computing time under different dynamics constraint versus discrete points

    图  13  随机初始状态下最优轨迹簇的mf-tf曲线

    Fig.  13  mf-tf curves of optimal trajectories under random initial state

    图  14  3种凸化方法最优轨迹对比

    Fig.  14  Comparison of optimal trajectories obtained by three convexification methods

    图  15  3种凸化方法求解得到的推力幅值曲线对比

    Fig.  15  Comparison of thrust magnitude obtained by three convexification methods

    图  16  P0、P1、P2求得规划轨迹与积分轨迹的误差随时间变化曲线

    Fig.  16  Error versus time between the integral trajectories and the optimal trajectories obtained by methods P0,P1,P2.

    图  17  3种方法积分所得终端位置误差散布图

    Fig.  17  Terminal position error of integral trajectories

    图  18  P0和P2积分所得终端位置误差散布图

    Fig.  18  Terminal position error of integral trajectories obtained by convexification methods P0 and P2

    图  19  3种方法积分所得终端速度误差散布图

    Fig.  19  Terminal velocity error of integral trajectories obtained by three methods

    图  20  P0和P2积分所得终端位置误差散布图

    Fig.  20  Terminal velocity error of integral trajectories obtained by convexification methods P0 and P2

    图  21  积分得到终端位置误差随离散点个数变化图

    Fig.  21  Terminal position error of integral trajectory versus the number of discrete points

    图  22  积分得到终端速度误差随离散点个数变化图

    Fig.  22  Terminal velocity error of integral trajectory versus the number of discrete points

    图  23  不同固定终端时间下积分所得终端位移误差三维散布图及xy平面图

    Fig.  23  3D and projection on xy plane of terminal position error dispersion of integral trajectories at different fixed terminal times

    图  24  不同固定终端时间下积分所得终端速度误差三维散布图及xy平面图

    Fig.  24  3D and projection on XY plane of terminal velocity error dispersion of integral trajectories at different fixed terminal times

    图  25  P2解得任务181的规划轨迹

    Fig.  25  The optimal trajectory of case 181 by using P2 convexification method

    图  26  P2解得任务181的速度轨线

    Fig.  26  Velocity profiles of case 181 by using convexification method P2

    图  27  200组随机任务打靶耗时统计图

    Fig.  27  Computing time statistics of 200 random cases

    图  28  200组随机任务中指标J的统计图

    Fig.  28  The index J statistics of 200 random cases

    图  29  不同信赖域约束下优化得到的轨线

    Fig.  29  Optimal trajectories under different trust regions

    图  30  不同信赖域下求解耗时统计

    Fig.  30  Computing time statistics under different trust regions

    表  1  仿真设定参数

    Table  1  Simulation parameters

    参数数值
    初始坐标r0 /m(0,0,10 000)
    初始速度v0 /(m·s–1(400,–200,–200)
    初始质量m0 /kg50
    推力幅值[FminFmax]/N[100,1 000]
    终端坐标rf /m(0,0,0)
    终端速度vf /(m·s–1(0,0,0)
    发动机比冲g0Isp /(m·s–12 000
    终端时间tf /s60
    离散点个数N100
    下载: 导出CSV

    表  2  随机仿真参数设置

    Table  2  Simulation parameters

    参数数值
    离散点个数N30
    初始坐标r0 /mrx0∈[–300,300]
    ry0∈[–300,300]
    rz0∈[1500,10000]
    初始速度v0 /(m·s–1vx0∈[–300,300]
    vy0∈[–300,300]
    vz0∈[–300,0]
    终端位置rf /mrf =(0,0,0)
    终端速度vf /(m·s–1vf =(0,0,0)
    tf 步长/s5
    下载: 导出CSV

    表  3  3种凸化方法的计算结果

    Table  3  Results of three convexification methods

    收敛次数剩余质量mf/kg(剩余质量分数/%)最优性
    P0331.338 6(62.68)最优
    P1631.330 4(62.66)最优
    P2230.446 2(60.89)次优
    下载: 导出CSV

    表  4  终端误差散布仿真参数

    Table  4  Simulation parameters

    参数数值
    初始位置r0 /mrx0∈[–3 000,3 000]
    ry0∈[–3 000,3 000]
    rz0∈[1 500,10 000]
    初始速度v0 /(m·s–1vx0∈[–300,300]
    vy0∈[–300,300]
    vz0∈[–300,0]
    终端位置rf /mrf =(0,0,0)
    终端速度vf /(m·s–1vf =(0,0,0)
    终端时间tf /s[60,200]
    离散点N100
    下载: 导出CSV

    表  5  终端误差散布仿真参数

    Table  5  Simulation parameters

    参数数值
    离散点个数N100
    初始坐标r0/mrx0∈[–3 000,3 000]
    ry0∈[–3 000,3 000]
    rz0∈[1 500,10 000]
    初始速度v0/(m·s–1vx0∈[–300,300]
    vy0∈[–300,300]
    vz0∈[–300,0]
    终端位置rf/mrf=(0,0,0)
    终端速度vf/(m·s–1vf=(0,0,0)
    下载: 导出CSV

    表  6  不同信赖域下的求解结果

    Table  6  Results under different trust regions.

    编号信赖域约束${{\delta }}{\rm{ }}({\rm{ }}\left| {{{{x}}^{(k + 1)}} - {{{x}}^{(k)}}} \right| \leqslant {{\delta }}{\rm{ }})$收敛次数平均耗时/s
    1$\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{300}&{300}&{300}&{50} \end{array}} \right]$30.142 4
    2$\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{200}&{200}&{200}&{20} \end{array}} \right]$40.192 5
    3$\left[ {\begin{array}{*{20}{c}} {5000}&{3000}&{200}&{200}&{200}&{200}&{20} \end{array}} \right]$50.255 4
    4$\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{100}&{200}&{200}&{200}&{20} \end{array}} \right]$90.470 2
    530.129 5
    下载: 导出CSV
  • [1] 崔平远,赵泽端,朱圣英. 火星大气进入段轨迹优化与制导技术研究进展[J]. 宇航学报,2019,40(6):611-620.

    CUI P Y,ZHAO Z R,ZHU S Y. Research progress of trajectory optimization and guidance techniques for Mars atmospheric entry[J]. Journal of Astronautics,2019,40(6):611-620.
    [2] 于正湜, 崔平远. 行星着陆自主导航与制导控制研究现状与趋势[J]. 深空探测学报(中英文), 201, 3(4): : 345-355.

    YU Z S, CUI P Y. Research status and developing trend of the autonomous navigation, guidance, and control for planetary landing[J]. Journal of Deep Space Exploration, 2016, 3(4): 345-355.
    [3] MICHAEL S, ACIKMESE B. Successive convexification for 6-DoF Mars rocket powered landing with free-final-time[C]//2018 AIAA Guidance, Navigation, and Control Conference. [S. l.]: AIAA, 2018.
    [4] BLACKMORE L,SCHARF D P. Minimum-landing-error powered-descent guidance for Mars landing using convex optimization[J]. Journal of Guidance Control and Dynamics,2010,33(4):1161-1171. doi:  10.2514/1.47202
    [5] LIU X. Fuel-optimal rocket landing with aerodynamic controls[J]. Journal of Guidance Control & Dynamics,2018,42(46):1-13.
    [6] WOLF A A,GRAVES C,POWELL R,et al. Systems for pinpoint landing at Mars[J]. Advances in the Astronautical Sciences,2005,119:2677-2677.
    [7] 李鑫,欧阳高翔,孙成明,等. 基于凸优化的有限推力远程转移轨迹优化[J]. 航天控制,2016,34(3):19-25. doi:  10.3969/j.issn.1006-3242.2016.03.004

    LI X,OUYANG G X,SUN C M,et al. Optimal remote transfer with finite thrust based on convex optimization[J]. Aerospace Control,2016,34(3):19-25. doi:  10.3969/j.issn.1006-3242.2016.03.004
    [8] LIU X,SHEN Z,LU P. Entry trajectory optimization by second-order cone programming[J]. Journal of Guidance Control and Dynamics,2015,39(2):1-15.
    [9] DOMAHIDI A, CHU E, BOYD S. ECOS: an SOCP solver for embedded systems[C]//Proceedings European Control Conference. Zurich: [s. n.], 2013.
    [10] LU P,LIU X. Autonomous trajectory planning for rendezvous and proximity operations by conic optimization[J]. Journal of Guidance Control & Dynamics,2013,36(2):375-389.
    [11] 林晓辉,于文进. 基于凸优化理论的含约束月球定点着陆轨道优化[J]. 宇航学报,2013,34(7):901-908. doi:  10.3873/j.issn.1000-1328.2013.07.003

    LIN X H,YU W J. Constrained trajectory optimization for lunar pin-point landing based on convex optimization theory[J]. Journal of Astronautics,2013,34(7):901-908. doi:  10.3873/j.issn.1000-1328.2013.07.003
  • [1] 杨霞, 李少腾, 贾贺, 刘乃彬, 王垒.  深空大载荷航天器软着陆气囊制造技术 . 深空探测学报(中英文), 2021, 8(3): 324-332. doi: 10.15982/j.issn.2096-9287.2021.20210010
    [2] 武迪, 程林, 王伟, 李俊峰.  基于切换系统的小推力轨迹优化协态初始化方法 . 深空探测学报(中英文), 2021, 8(4): 1-6. doi: 10.15982/j.issn.2096-9287.2021.20200090
    [3] 陈上上, 关轶峰, 于萍, 李骥, 张晓文.  基于粒子群优化的月球陨石坑探测轨迹规划 . 深空探测学报(中英文), 2020, 7(3): 271-277. doi: 10.15982/j.issn.2095-7777.2020.20191031007
    [4] 张熇, 杜宇, 李飞, 张弘, 马继楠, 盛丽艳, 吴克.  月球南极探测着陆工程选址建议 . 深空探测学报(中英文), 2020, 7(3): 232-240. doi: 10.15982/j.issn.2095-7777.2020.20191003002
    [5] 张磊.  月面上升的轨迹优化与参数分析 . 深空探测学报(中英文), 2019, 6(4): 391-397. doi: 10.15982/j.issn.2095-7777.2019.04.012
    [6] 吴伟仁, 王琼, 唐玉华, 于国斌, 刘继忠, 张玮, 宁远明, 卢亮亮.  “嫦娥4号”月球背面软着陆任务设计 . 深空探测学报(中英文), 2017, 4(2): 111-117. doi: 10.15982/j.issn.2095-7777.2017.02.002
    [7] 李飞, 张熇, 吴学英, 马继楠, 卢亮亮.  月球背面地形对软着陆探测的影响分析 . 深空探测学报(中英文), 2017, 4(2): 143-149. doi: 10.15982/j.issn.2095-7777.2017.02.007
    [8] 闫晓鹏, 孙海滨, 郭雷.  火星着陆器的大气进入段有限时间抗干扰制导律设计 . 深空探测学报(中英文), 2016, 3(1): 61-67. doi: 10.15982/j.issn.2095-7777.2016.01.010
    [9] 袁旭, 朱圣英.  基于伪谱法的小天体最优下降轨迹优化方法 . 深空探测学报(中英文), 2016, 3(1): 51-55. doi: 10.15982/j.issn.2095-7777.2016.01.008
    [10] 胡海静, 朱圣英, 崔平远.  基于Lyapunov函数的小天体软着陆障碍规避控制方法 . 深空探测学报(中英文), 2015, 2(2): 149-154. doi: 10.15982/j.issn.2095-7777.2015.02.008
    [11] 李爽, 陶婷, 江秀强, 张树瑜, 周杰.  月球软着陆动力下降制导控制技术综述与展望 . 深空探测学报(中英文), 2015, 2(2): 111-119. doi: 10.15982/j.issn.2095-7777.2015.02.002
    [12] 郭延宁, 马广富, 曾添一, 崔祜涛.  基于燃料最优解的火星精确着陆制导策略研究 . 深空探测学报(中英文), 2015, 2(1): 61-68. doi: 10.15982/j.issn.2095-7777.2015.01.009
    [13] 秦同, 朱圣英, 崔平远.  火星软着陆能量最优制导律转移能力分析 . 深空探测学报(中英文), 2015, 2(3): 218-223. doi: 10.15982/j.issn.2095-7777.2015.03.005
    [14] 韩松涛, 唐歌实, 曹建峰, 陈略, 任天鹏, 王美.  面向"CE-3号"软着陆过程的深空网干涉测量技术 . 深空探测学报(中英文), 2015, 2(2): 120-124. doi: 10.15982/j.issn.2095-7777.2015.02.003
    [15] 王茜茜, 谢慕君, 李元春.  基于模糊参数优化的小行星软着陆控制方法研究 . 深空探测学报(中英文), 2015, 2(2): 162-167. doi: 10.15982/j.issn.2095-7777.2015.02.010
    [16] 王大铁, 李骥, 黄翔宇, 张洪华.  月球软着陆过程高精度自主导航避障方法 . 深空探测学报(中英文), 2014, 1(1): 44-51.
    [17] 吴伟仁, 于登云.  “嫦娥3号”月球软着陆工程中的关键技术 . 深空探测学报(中英文), 2014, 1(2): 105-109.
    [18] 黄翔宇, 张洪华, 王大铁, 李骥, 关软峰, 王鹏基.  “嫦娥三号”探测器软着陆自主导航与制导技术 . 深空探测学报(中英文), 2014, 1(1): 52-59.
    [19] 魏若岩, 阮晓钢, 庞涛, OuattaraSIE, 武旋, 肖尧.  小天体软着陆中的地面特征区域提取与跟踪算法 . 深空探测学报(中英文), 2014, 1(4): 308-314. doi: 10.15982/j.issn.2095-7777.2014.04.011
    [20] 贺晶, 龚胜平, 李俊峰.  利用逃逸能量的太阳帆最快交会轨迹优化 . 深空探测学报(中英文), 2014, 1(1): 60-66.
  • 加载中
图(30) / 表 (6)
计量
  • 文章访问数:  206
  • HTML全文浏览量:  83
  • PDF下载量:  18
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-12-15
  • 修回日期:  2021-03-08
  • 网络出版日期:  2021-05-12
  • 刊出日期:  2021-04-28

行星着陆动力下降轨迹优化的中心差分凸化方法

doi: 10.15982/j.issn.2096-9287.2021.20200079
    基金项目:  国防科技重点实验室基金资助项目(6142210200312)
    作者简介:

    李明翔(1996-),男,硕士研究生,主要研究方向:飞行器轨迹优化与制导。通讯地址:陕西省西安市友谊西路127号西北工业大学航天学院(710072)电话:18829571869 E-mail:limingxiang@mail.nwpu.edu.cn

    通讯作者: 泮斌峰(1981– ),男,教授,博士生导师,主要研究方向:航天器计算制导控制、机器学习、轨迹优化、深空探测技术。本文通讯作者。通讯地址:陕西省西安市友谊西路127号西北工业大学航天学院(710072)电话:18691959798 Email:panbinfeng@nwpu.edu.cn
  • ● The workload of code programming and formula derivation can be reduced and the generality of algorithm can be improved by using pre-labeled central difference method instead of analytically computing the Jacobi matrix. ● mf-tf curves of optimal trajectories can be fitted by a cluster of functions which can be used to approximately estimate the optimal terminal time. ● Simulation results show that direct linearization of dynamic constraints is better than customized convexification method with variable substitution for planet powered landing.
  • 中图分类号: V467.3

摘要: 针对行星定点软着陆实时在线制导的任务需求,设计了基于序列凸优化的动力下降燃料最优轨迹求解算法。算法采用预标记的中心差分法线性化动力学方程,并提出将性能指标相对偏差作为收敛终止条件,能够快速生成燃料最优轨迹。此外,在分析最优轨迹簇剩余燃料和终端时间关系的基础上给出其拟合函数,作为最优终端时间的近似估算,以减少算法求解未知变量的维数。数值仿真结果表明,与对动力学方程先进行变量代换再线性化的传统凸化方法相比,该算法对初始猜想不敏感,收敛性好且终端误差较小。

注释:
1)  ● The workload of code programming and formula derivation can be reduced and the generality of algorithm can be improved by using pre-labeled central difference method instead of analytically computing the Jacobi matrix. ● mf-tf curves of optimal trajectories can be fitted by a cluster of functions which can be used to approximately estimate the optimal terminal time. ● Simulation results show that direct linearization of dynamic constraints is better than customized convexification method with variable substitution for planet powered landing.

English Abstract

李明翔, 泮斌峰. 行星着陆动力下降轨迹优化的中心差分凸化方法[J]. 深空探测学报(中英文), 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
引用本文: 李明翔, 泮斌峰. 行星着陆动力下降轨迹优化的中心差分凸化方法[J]. 深空探测学报(中英文), 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
LI Mingxiang, PAN Binfeng. Central Difference Convexification Method for Soft-Landing Trajectory Optimization in Planetary Powered Descending Phase[J]. Journal of Deep Space Exploration, 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
Citation: LI Mingxiang, PAN Binfeng. Central Difference Convexification Method for Soft-Landing Trajectory Optimization in Planetary Powered Descending Phase[J]. Journal of Deep Space Exploration, 2021, 8(2): 171-181. doi: 10.15982/j.issn.2096-9287.2021.20200079
    • 行星探测是人类探索太空最重要的途径之一,月球和火星探测是当前行星探测的热点。当目标行星大气比较稀薄或者没有大气时,探测器将无法使用降落伞降至安全速度,基于反推火箭的动力下降成为实现行星表面软着陆的一种重要手段。由于行星动力下降过程存在诸多不确定性因素,其在线轨迹优化是实现探测器安全、精确着陆所需要解决的重要问题。

      近年来,序列凸优化方法成为求解上述在线轨迹优化的有效手段[1-2],其本质是将无限维空间中的最优控制问题通过离散转化为有限维空间中的非线性规划问题,并利用序列凸优化理论通过凸化处理后迭代求解。Szmuk等[3]通过引入虚拟控制量并且在动力学方程中定义了新的时间状态量,提出了一种终端时间自由的序列凸优化方法计算框架,能够较好地得到最优燃料轨迹;其不足之处是对初始猜想敏感且迭代次数多。Blackmore等[4]通过将燃料最优问题凸化为二阶锥规划(Second-Order-Cone-Programming,SOCP)问题,并考虑了坡度约束,能够较快求解出燃料最优的下降轨迹。上述方法对于终端时间固定的情况能够得到较好的解,但由于其终端时间采用线性搜索方法计算,效率较低。此外,Liu等[5]考虑了气动力作用下的一级火箭回收问题,通过无损凸化技术和动力学方程归一化提高了算法收敛性和快速性。该方法在凸化过程中重定义了控制量,公式推导复杂且引入了额外的约束和计算量,只能处理固定终端时间问题。

      有别于上述研究方法,本文探索使用预标记中心差分替代解析表达式计算动力学矩阵,并对动力学方程采用直接线性化替代变量代换后再线性化的凸化方法,最后研究行星动力下降燃料最优控制问题中近似最优终端时间的确定方法,以提高序列凸优化在线轨迹优化方法的工程实用性。

    • 以着陆点为原点建立空间直角坐标系,着陆点一般不直接选在地表,而是距离地面稍有一定距离,方便提供安全缓冲余量。x轴方向预先约定,z轴垂直于地表向上,y轴满足右手螺旋关系。行星动力下降可以表述为经典的Bolza型最优控制问题

      $$ \left\{ \begin{array}{l} \min {{\theta}} ({{{x}}_{{t_{\rm{f}}}}}) + \int_{{t_0}}^{{{\rm{t}}_f}} {\dfrac{{\left\| {{F}} \right\|}}{{{g_0}{I_{{\rm{sp}}}}}}{\rm{d}}t} \\ {\rm{s}}{\rm{.t}}{\rm{. }}\;\;\;\;{{{x}}_0} = {{x}}({t_0}) \\ \;\;\;\;\;\;\;\;\; \dot {{x}} = {{f}}({{x}},{{u}},t) \\ \;\;\;\;\;\;\;\;\;{F_{\min }} \text{≤} \left\| {{F}} \right\| \text{≤} {F_{\max }} \\ \;\;\;\;\;\;\;\;\;{{{r}}_z} \text{≥} 0 \\ \end{array} \right. $$ (1)

      其中终端性能指标为

      $$\theta ({{{x}}_{{t_{\rm{f}}}}}) = \left| {{ w^{\rm{T}}}\left( {{ x_{{t_{\rm{f}}}}} - x({t_{\rm{f}}})} \right)} \right|$$ (2)

      式(1)中约束分别为初值约束、动力学约束、控制量幅值约束和防撞约束,其中反推发动机为了避免多次开关机带来的可靠性问题,在实际工作中一般不将其完全关闭,故需限定最小推力。行星动力下降动力学方程为

      $$\left\{ \begin{array}{l} {\dot{{ r}}} = {{ v}} \\ \dot {{v}} = {{g}} + \dfrac{{{F}}}{m} + {{d}} \\ \dot m = - \dfrac{{\left\| {{F}} \right\|}}{{{g_0}{I_{{\rm{sp}}}}}} \end{array} \right.$$ (3)

      其中:g = [0 0 –grz)]T是行星的重力加速度矢量。

      考虑到动力下降段高度在2~5 km之间,高度相比行星半径尺度较小,重力加速度变化也较小,可以直接通过海拔高度计算。本文以火星引力场为例,其表面重力加速度为gM = 0.377、gE = 3.7 m/s2RM = 3 397 km,距地表任意高度的重力加速度为

      $$g({r_z}) = {g_{\rm{M}}}{\left( {\frac{1}{{1 + {r_z}/{R_{\rm{M}}}}}} \right)^2}$$ (4)

      此外,d = [dx dy dz]T是3个方向的气动加速度,分为可建模的气动(如飞行器的气动阻力)和无法建模的随机风场扰动。由于火星大气较为复杂,不适合在在线轨迹优化中应用[1],因此本文采用文献[3~4]的处理方式,在动力下降段中以忽略气动力作用下简化的动力学进行分析,气动扰动可以通过后期设计控制器加以克服[6]

    • 采用直接法求解优化问题式(1)时,需要先将其离散到有限维空间中。当终端时间tf确定时,将飞行时间[0,tf]等距离散为N个时刻,即N–1个子区间,并在这N个时刻对状态量和控制量进行离散,方便施加过程约束。

      对于控制力幅值约束

      $${F_{\min }} \text{≤} \sqrt {F_x^2 + F_y^2 + F_z^2} \text{≤} {F_{\max }}$$ (5)

      因最小推力幅值约束是二阶锥补集,为将非凸约束式(5)凸化,在各个离散点处引入推力松弛因子Fη

      $$ \sqrt{{F}_{xi}^{2}+{F}_{yi}^{2}+{F}_{zi}^{2}}\text{≤} {F}_{\eta i},i=0,1,\cdots,N-1$$ (6)

      此时,控制幅值约束可以表示为

      $$ {F}_{\mathrm{min}}\text{≤} {F}_{\eta i}\text{≤} {F}_{\mathrm{max}},i=0,1,\cdots,N-1$$ (7)

      至此非凸约束式(5)被等价转换为凸约束式(6)和式(7),积分型性能指标也可以表达为关于Fη的线性组合

      $$\int_{{t_0}}^{{t_f}} {\frac{{\left\| {{ F}} \right\|}}{{{g_0}{I_{{\rm{sp}}}}}}{\rm{d}}t} \approx \frac{{\Delta t}}{{{g_0}{I_{{\rm{sp}}}}}}\sum\limits_{i = 0}^{N - 1} {{F_{\eta i}}} $$ (8)

      其中:Δt =(tft0)/(N–1),即使用矩形公式将积分化为求和。

    • 为将非线性动力学微分方程约束处理为凸约束,需进行凸化处理。文献[5]和文献[7]中处理上述动力学方程通常使用2种凸化技巧:将推力矢量分解为推力大小和方向余弦、对质量取对数重定义新的状态量。前者存在的问题是,方向余弦约束为非凸等式约束,虽然文献[3]中通过极小值原理证明其可无损转化为凸锥约束,但实际数值求解中经常出现不稳定因素,造成收敛速度较慢,且控制变量从原来的3维上升到4维,问题复杂度上升。后者对质量取对数后重定义新的状态量,再对推力约束进行一阶泰勒展开线性化和直接对动力学泰勒展开线性化,本质上都是线性化,并未降低原问题的非线性,仅将非线性从动力学方程转移到过程约束中。取对数后推力不等式约束中存在指数函数,相比原来仅含加、乘法的公式增加了计算量,且人为增加了推导公式和编程的工作量。

      本文中不再对动力学使用多余的凸化技巧,直接在第k+1次迭代时于各个离散时间点处,使用第k次优化得到的状态量和控制量的解对非线性动力学方程进行线性化,使其在每个子区间[nΔt,(n+1)Δt](n = 0,1,···,N–2)内变为线性动力学,即$\dot {{x}} = {{A}}_n^{(k)}{{x}} + $${{B}}_n^{(k)}{{u}} + {{C}}_n^{(k)}$

      对于动力学矩阵ABC的分量使用差分直接计算,第i个微分方程,对第n个离散点处的状态量(控制量)的第j个分量的偏导的差分计算公式为

      $$\left\{ \begin{array}{*{20}{l}} {{{A}}_{ij}}({{x}}_n^{(k)},{{u}}_n^{(k)}) = \dfrac{{\partial {f_i}}}{{\partial {{{x}}_j}}} =\\ \dfrac{{{f_i}({{x}}_n^{(k)} + {{{e}}_j} \times {{\delta}},{{u}}_n^{(k)}) - {f_i}({{x}}_n^{(k)} - {{{e}}_j} \times {{\delta}},{{u}}_n^{(k)})}}{{2{{\delta}} }} \\ {{{B}}_{ij}}({{x}}_n^{(k)},{{u}}_n^{(k)}) = \dfrac{{\partial {f_i}}}{{\partial {{{u}}_j}}} =\\ \dfrac{{{f_i}({{x}}_n^{(k)},{{u}}_n^{(k)} + {{{e}}_j} \times {{\delta}} ) - {f_i}({{x}}_n^{(k)},{{u}}_n^{(k)} - {{{e}}_j} \times {{\delta}} )}}{{2{{\delta}} }} \\ {{C}}({{x}}_n^{(k)},{{u}}_n^{(k)}) = {{f}} - {{{{f}}'}_x}{{x}} - {{{{f}}'}_u}{{u}} = {{f}}({{x}}_n^{(k)},\\ {{u}}_n^{(k)}) - {{A}}({{x}}_n^{(k)},{{u}}_n^{(k)}){{x}}_n^{(k)} - {{B}}({{x}}_n^{(k)},{{u}}_n^{(k)}){{u}}_n^{(k)} \end{array} \right.$$ (9)

      为了消除不必要的计算量,在设计程序时,提前对每个动力学微分方程显含的状态量和控制量索引进行标注,仅计算这些位置的偏导数,其余元素直接赋值为0。将式(1)中的微分方程约束转换为线性代数方程约束,为表示方便,简记为

      $$\left\{ \begin{array}{l}{{{A}}}_{n}^{(k)}={{A}}({{x}}_{n}^{(k)},{{{u}}}_{n}^{(k)}),\\ {{{B}}}_{n}^{(k)}={{B}}({{x}}_{n}^{(k)},{{{u}}}_{n}^{(k)}),\\ {{{C}}}_{n}^{(k)}={{C}}({{x}}_{n}^{(k)},{{{u}}}_{n}^{(k)})\end{array} \right.$$ (10)

      为兼顾计算精度和计算效率,使用梯形法施加动力学约束

      $$\left\{ \begin{array}{l} {{{k}}_1} = {{A}}_n^{(k)}{{{x}}_n} + {{B}}_n^{(k)}{{{u}}_n} + {{C}}_n^{(k)} \\ {{{k}}_2} = {{A}}_{n + 1}^{(k)}{{{x}}_{n + 1}} + {{B}}_{n + 1}^{(k)}{{{u}}_{n + 1}} + {{C}}_{n + 1}^{(k)} \\ {{{x}}_{n + 1}} = {{{x}}_n} + \dfrac{{\Delta t}}{2}({{{k}}_1} + {{{k}}_2}) \\ \end{array}\!\!\!, \right.\;\;n = 0,2, \cdots,N - 2$$ (11)

      启动第一次迭代时,动力学矩阵ABC的计算需要提供k = 0时每个离散点n处的状态量x(0)和控制量u(0),即初始猜想。对于状态量,如果同时存在初值约束和终端约束,则状态量的初始猜想直接设计为两个状态的线性插值

      $${{x}}_n^{(0)} = \left( {1 - \frac{n}{N}} \right){{{x}}_0} + \frac{n}{N}{{{x}}_{{t_{\rm f}}}},n = 0,1, \cdots,N - 1$$ (12)

      如果仅存在初值约束(如质量m),则状态量初始猜想设计为常值向量

      $${{x}}_n^{(0)} = {{{x}}_0},n = 0,1, \cdots,N - 1$$ (13)

      对于控制量,初始猜测全部初始化为0,即

      $${{u}}_n^{(0)} = 0,n = 0,1, \cdots,N - 1$$ (14)

      式(12)~(14)简洁的设计过程直接消除了利用数值积分设计的轨迹引入的人为因素和计算量,虽然初始猜想并不符合动力学方程约束,但在每次迭代求解后得到的轨迹将完全满足动力学约束,并在不断迭代中使轨迹收敛到最优解。

    • 为验证前文所设计算法的有效性,给定飞行状态和航天器参数如表1所示,求解行星动力下降段燃料最优的轨迹优化问题。

      表 1  仿真设定参数

      Table 1.  Simulation parameters

      参数数值
      初始坐标r0 /m(0,0,10 000)
      初始速度v0 /(m·s–1(400,–200,–200)
      初始质量m0 /kg50
      推力幅值[FminFmax]/N[100,1 000]
      终端坐标rf /m(0,0,0)
      终端速度vf /(m·s–1(0,0,0)
      发动机比冲g0Isp /(m·s–12 000
      终端时间tf /s60
      离散点个数N100

      设置迭代次数为5次,每步迭代的性能指标J,性能指标绝对偏差e = |Jk+1Jk|和相对偏差er = e/Jk+1,仿真结果如图1所示。一般地,序列凸优化常以$ \Vert {{{x}}}_{n}^{(k+1)}-{{{x}}}_{n}^{(k)}\Vert \text{≤} \varepsilon,$$n=0,1,\cdots,N-1 $作为迭代收敛的条件[8],但在轨迹震荡时易导致无法快速收敛且收敛指标$\varepsilon $需根据不同物理量进行独立设计,不具有普适性。分析图1可知,er作为判断收敛指标效果最好,也在一定程度上克服了上述缺点。

      图  1  性能指标及其绝对偏差和相对偏差随迭代次数变化图

      Figure 1.  Performance index and its absolute deviation,relative deviation versus iterations

      本文采用C++调用SOCP求解器ECOS[9]进行求解。计算环境建立在CPU为i5-8300H(基频2.2 GHz,睿频3.8 GHz)Windows系统的PC上,耗时统计从进入main函数开始到满足序列凸优化迭代条件er < 1%退出为止(不含优化结果输出时间),求解500次,均为迭代3步收敛,耗时统计如图2所示,平均耗时0.130 1 s,优化所得控制量与状态量如图3~6所示。

      优化得到的推力曲线如图3所示,可以看到推力幅值被精确地约束在[100,1 000]N的区间内,表现为bang-bang控制。速度曲线如图4所示,可看到终端速度被精确约束在0 m/s,图5中质量曲线mt)的导数是控制力幅值函数||Ft)||2,可以看到mt)随||Ft)||2也相应分为3段。

      图  2  求解表1算例500次耗时统计

      Figure 2.  Computing time statistic of solving the case in Table 1 500 times

      图  3  发动机3个方向推力曲线

      Figure 3.  Thrust component in three directions

      图  4  飞行器速度曲线

      Figure 4.  Velocity profiles of optimal trajectory

      图1证明了用中心差分法建立序列凸优化模型的有效性,良好收敛性和对初始猜测弱敏感性,事实上由图1还可以看出如果以相对偏差作为性能指标,可认为第3步已经收敛,轨迹随迭代次数变化如图7所示,3~5次迭代产生的轨迹基本重合,算例对应的三维轨迹和各离散点处推力矢量如图8所示。

      图  5  飞行器质量曲线

      Figure 5.  Mass profile of optimal trajectory.

      图  6  飞行器位置相对于着陆点变化曲线

      Figure 6.  Position profiles of optimal trajectory

      图  7  tf = 60 s轨迹随迭代次数变化

      Figure 7.  Trajectory versus iterations for case of tf=60 s

      图  8  tf = 60 s优化轨迹和推力矢量

      Figure 8.  Optimal trajectory and thrust for case of t f = 60 s

    • 如果选取不同的终端时间tf,可以求解出对应的最优控制轨迹如图9所示,并同时可以得到飞行器剩余质量如图10所示,可见存在最优终端时间$t_{\rm{f}}^*$使得在所有最优轨迹中消耗燃料最少,在mf -tf散点图中,可发现右端具有明显的线性递减的特征,左端递增且增速快于右端,极值点必然存在。经过大量拟合测试,发现采用次数分别为–3,0,1的组合多项式${m_{\rm{f}}}({t_{\rm{f}}}) = a{t_{\rm{f}}} + b/t_{\rm{f}}^3 + c$进行拟合后,均方误差和确定系数均较理想,且兼顾了极值点计算的简易性。上述问题中拟合参数a = –0.045 04,b = –6.72 × 105c = 37.1,拟合结果的标准差为0.042 2,确定系数为0.999 1。通过对该函数求极值,可以直接得到其近似最优终端时间

      图  9  tf∈[50,190]s最优轨迹簇

      Figure 9.  Optimal trajectories at tf∈[50,190]s

      图  10  不同tf下最优轨迹对应的剩余燃料质量mf

      Figure 10.  mf of optimal trajectories at different tf

      $$t_{\rm{f}}^{\rm{*}} = \sqrt[\root{10}{{4}}] {{\frac{{3b}}{a}}} = 81.794\rm \;s$$ (15)

      图11可见,降低离散点个数对最优终端时间的确定几乎不产生影响,但可以显著降低计算量,如图12所示。在实际工程应用中,可以通过降低离散点个数的方式或者使用低精度的动力学约束进行最优控制模型的求解,以快速获取(tfmf)样本点。样本点个数不能少于3个,使用最小二乘逼近算法得到abc 3个参数,利用式(15)计算出最优终端时间后,再增大离散点个数进行精规划。

      图  11  不同离散点个数对应的mf-tf曲线

      Figure 11.  mf-tf curves under different number of discrete points

      图  12  不同离散点个数与动力学约束下的计算耗时

      Figure 12.  Computing time under different dynamics constraint versus discrete points

      为验证上述方法的有效性,用计算机随机产生100组初始状态和终端状态,计算出每组条件下的mf -tf曲线,仿真参数如表2所示。仿真结果如图13所示,在100个飞行条件下,除去2个苛刻的飞行情况无解外,其余98个条件下都找到了最优燃料控制的规划轨迹,且有较高的规划成功率。

      表 2  随机仿真参数设置

      Table 2.  Simulation parameters

      参数数值
      离散点个数N30
      初始坐标r0 /mrx0∈[–300,300]
      ry0∈[–300,300]
      rz0∈[1500,10000]
      初始速度v0 /(m·s–1vx0∈[–300,300]
      vy0∈[–300,300]
      vz0∈[–300,0]
      终端位置rf /mrf =(0,0,0)
      终端速度vf /(m·s–1vf =(0,0,0)
      tf 步长/s5

      图  13  随机初始状态下最优轨迹簇的mf-tf曲线

      Figure 13.  mf-tf curves of optimal trajectories under random initial state

    • 记本文中不对动力学改造直接采用中心差分进行线性化建立的序列凸优化算法为P0,另外两种凸化方法分别为P1[10]和P2[11]

      $$\left\{ \begin{array}{l} {\rm{P1}}:\;\;\;\min \dfrac{{\Delta t{T_{\max }}}}{{{g_0}{I_{{\rm{sp}}}}}} \displaystyle\sum\limits_{i = 0}^{N - 1} u + \left| {{{{w}}^{\rm{T}}}({{x}}{\rm{(}}{t_f}{\rm{)}} - {{{x}}_{_f}})} \right|\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;\;\;\dfrac{{{T_{\min }}}}{{{T_{\max }}}} \text{≤} u \text{≤} 1,\;{{{r}}_z} \text{≥} 0\\ \;\;\;\;\;\;\;\;\;\;\;\;\;n_x^2 + n_y^2 + n_z^2 \text{≤} 1\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\dot{ r}} = {{v}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\dot{ v}} = {{g}} + \dfrac{{u{T_{\max }}}}{m}{{n}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;m = - \dfrac{{{T_{\max }}}}{{{g_0}{I_{{\rm{sp}}}}}}u \end{array} \right.$$ (16)
      $$\left\{ \!\!\!\! \begin{array}{l} {\rm{P2}}:\min \;\; - {z_{\rm{f}}} + \left| {{{{w}}^{\rm{T}}}({{x}}({t_{\rm{f}}}) - {{{x}}_{_{\rm{f}}}})} \right|\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\sqrt {a_x^2 + a_y^2 + a_z^2} \text{≤} u,\;{{{r}}_z} \text{≥} 0\\ \;\;\;\;\;\;{\dot{ r}} = {{v}}\\ \;\;\;\;\;\;{\dot{ v}} = {{g}} + {{a}}\\ \;\;\;\;\;\;z = - \dfrac{u}{{{g_0}{I_{{\rm{sp}}}}}}\\ \;\;\;\;\;\;{T_{\min }}{{\rm{e}}^{ - {z^{(k)}}}}\left( {1 - (z - {z^{(k)}})} \right) \text{≤} u \text{≤} {T_{\max }}{{\rm{e}}^{ - {z^{(k)}}}}\left( {1 - (z - {z^{(k)}})} \right) \end{array} (17) \right. $$

      P1中将推力分解为推力大小u和推力方向n共计4个控制量进行优化,并使用无损凸化理论将方向余弦等式约束$n_x^2 + n_y^2 + n_z^2 = 1$松弛为$n_x^2 + n_y^2 + n_z^2 \text{≤} 1$不等式约束,P2定义了新的状态量为z = ln m,这将原问题动力学的一部分非线性转移到不等式约束中。设置3个模型状态量的初始猜想一致(P2中新状态量z全部初始化为ln m0),控制量都全部初始化为0,使用表1中的任务参数求解P0、P1和P2。3种方法对比结果如图14~16所示。结果显示P0和P1求解结果基本一致,仅推力分量曲线有较小差别,P2收敛得更快,但是并未收敛到全局最优解,P2的最优性依赖于更好的初始猜想。

      图  14  3种凸化方法最优轨迹对比

      Figure 14.  Comparison of optimal trajectories obtained by three convexification methods

      图  15  3种凸化方法求解得到的推力幅值曲线对比

      Figure 15.  Comparison of thrust magnitude obtained by three convexification methods

      表 3  3种凸化方法的计算结果

      Table 3.  Results of three convexification methods

      收敛次数剩余质量mf/kg(剩余质量分数/%)最优性
      P0331.338 6(62.68)最优
      P1631.330 4(62.66)最优
      P2230.446 2(60.89)次优

      为将非线性动力学模型转化为凸问题进行求解,三者都使用了等距离散的分段线性动力学进行近似逼近,线性化误差和离散误差不可避免。使用控制量积分得到的轨迹和规划所得轨迹进行对比,积分步长为0.01 s。发现3种凸化方式规划解和积分解误差随时间变化都在同一数量级内。

      图  16  P0、P1、P2求得规划轨迹与积分轨迹的误差随时间变化曲线

      Figure 16.  Error versus time between the integral trajectories and the optimal trajectories obtained by methods P0,P1,P2.

      表4设置的参数下,产生500组随机动力下降任务,分别使用3种方法进行求解,观察终端误差散布。仿真结果如图17~20所示。

      表 4  终端误差散布仿真参数

      Table 4.  Simulation parameters

      参数数值
      初始位置r0 /mrx0∈[–3 000,3 000]
      ry0∈[–3 000,3 000]
      rz0∈[1 500,10 000]
      初始速度v0 /(m·s–1vx0∈[–300,300]
      vy0∈[–300,300]
      vz0∈[–300,0]
      终端位置rf /mrf =(0,0,0)
      终端速度vf /(m·s–1vf =(0,0,0)
      终端时间tf /s[60,200]
      离散点N100

      图  17  3种方法积分所得终端位置误差散布图

      Figure 17.  Terminal position error of integral trajectories

      图  18  P0和P2积分所得终端位置误差散布图

      Figure 18.  Terminal position error of integral trajectories obtained by convexification methods P0 and P2

      图  19  3种方法积分所得终端速度误差散布图

      Figure 19.  Terminal velocity error of integral trajectories obtained by three methods

      图  20  P0和P2积分所得终端位置误差散布图

      Figure 20.  Terminal velocity error of integral trajectories obtained by convexification methods P0 and P2

      根据对比,可以发现P1方法的数值稳定性较差,部分算例达到最大迭代次数20次后仍然没有收敛,导致终端误差较大,P0和P2方法终端误差数量级接近。P0相比于P2的位置终端误差散布更加集中,且速度终端误差在z方向上更小。

      另外,离散步长对终端误差也都有一定影响。以P0为例,计算表1任务在不同离散点个数下得到的规划轨迹,并通过积分得到终端误差如图21图22所示。可以看到,终端误差随离散点个数的增加而不断减少。

      图  21  积分得到终端位置误差随离散点个数变化图

      Figure 21.  Terminal position error of integral trajectory versus the number of discrete points

      图  22  积分得到终端速度误差随离散点个数变化图

      Figure 22.  Terminal velocity error of integral trajectory versus the number of discrete points

      固定离散点个数,针对表5设置的随机参数,采用P0在不同终端时间下各打靶200次,如图23图24所示,积分所得终端误差散布也处于可接受范围内。

      表 5  终端误差散布仿真参数

      Table 5.  Simulation parameters

      参数数值
      离散点个数N100
      初始坐标r0/mrx0∈[–3 000,3 000]
      ry0∈[–3 000,3 000]
      rz0∈[1 500,10 000]
      初始速度v0/(m·s–1vx0∈[–300,300]
      vy0∈[–300,300]
      vz0∈[–300,0]
      终端位置rf/mrf=(0,0,0)
      终端速度vf/(m·s–1vf=(0,0,0)

      图  23  不同固定终端时间下积分所得终端位移误差三维散布图及xy平面图

      Figure 23.  3D and projection on xy plane of terminal position error dispersion of integral trajectories at different fixed terminal times

      图  24  不同固定终端时间下积分所得终端速度误差三维散布图及xy平面图

      Figure 24.  3D and projection on XY plane of terminal velocity error dispersion of integral trajectories at different fixed terminal times

      对于收敛性,200个随机初始状态中,P0、P1都求出198个解(第43和181个失败),P2求出199个解(第43个失败),三者收敛性无显著差别。两任务的初始速度大,高度低,规划也失败属于合理现象。虽然第181个任务的解被P2找到,但这条轨迹十分危险,在18.18 s处离地面较近且仍然有较大的水平速度,考虑存在离散误差和环境扰动的情况下可能会导致任务失败,两任务参数如下。

      $$ \begin{array}{l} 43:\;{{{r}}_0} = {\left[ {\begin{array}{*{20}{l}} {317.972}&{117.466}&{2\;182.24} \end{array}} \right]^{\rm{T}}},\\ \;\;\;\;\;{{{v}}_0} = {\left[ {\begin{array}{*{20}{c}} {70.836\;5}&{ - 274.877}&{ - 279.446} \end{array}} \right]^{\rm{T}}} \end{array} $$
      $$\begin{array}{l} 181:\;{{{r}}_0} = {\left[ {\begin{array}{*{20}{l}} { - 2\;725.7}&{ - 790.948}&{2\;673.82} \end{array}} \right]^{\rm{T}}},\\ \;\;\;\;\;\;\;{{{v}}_0} = {\left[ {\begin{array}{*{20}{c}} {228.953}&{ - 182.388}&{ - 298.709} \end{array}} \right]^{\rm{T}}} \end{array}$$

      图  25  P2解得任务181的规划轨迹

      Figure 25.  The optimal trajectory of case 181 by using P2 convexification method

      使用100个离散点求解,所得的仿真结果如图26~28所示。算法P0平均耗时0.129 2 s,平均3.04步收敛。P2平均耗时0.115 7 s,2步收敛。单步迭代P0耗时更少,总耗时P2略少但优化结果P0更优。若定义P0,P2收敛后剩余燃料质量分别为${m_{P0{\rm{f}}}}$${m_{P2{\rm{f}}}}$,定义指标$J = ({m_{P0{\rm{f}}}} - {m_{P2{\rm{f}}}})/{m_0}$,统计打靶结果显示J总为正,即P0总是收敛到相对于P2更好的燃料最优解。

      图  26  P2解得任务181的速度轨线

      Figure 26.  Velocity profiles of case 181 by using convexification method P2

      图  27  200组随机任务打靶耗时统计图

      Figure 27.  Computing time statistics of 200 random cases

      图  28  200组随机任务中指标J的统计图

      Figure 28.  The index J statistics of 200 random cases

      以上求解过程都未引入信赖域约束,若增加状态量信赖域约束,在表1仿真参数下采用P0进行求解,每种信赖域各计算100次,图29显示对解最优性并未有明显改善。图30为不同信赖域下求解耗时统计。前3条轨线和无信赖域约束求解结果基本重合,第4条轨迹消耗燃料相比前3条节省了0.032 7 kg,仅占总质量的0.065 4%,但耗时相对较长。在本问题中基本可视信赖域为冗余约束。

      图  29  不同信赖域约束下优化得到的轨线

      Figure 29.  Optimal trajectories under different trust regions

      图  30  不同信赖域下求解耗时统计

      Figure 30.  Computing time statistics under different trust regions

      表 6  不同信赖域下的求解结果

      Table 6.  Results under different trust regions.

      编号信赖域约束${{\delta }}{\rm{ }}({\rm{ }}\left| {{{{x}}^{(k + 1)}} - {{{x}}^{(k)}}} \right| \leqslant {{\delta }}{\rm{ }})$收敛次数平均耗时/s
      1$\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{300}&{300}&{300}&{50} \end{array}} \right]$30.142 4
      2$\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{5000}&{200}&{200}&{200}&{20} \end{array}} \right]$40.192 5
      3$\left[ {\begin{array}{*{20}{c}} {5000}&{3000}&{200}&{200}&{200}&{200}&{20} \end{array}} \right]$50.255 4
      4$\left[ {\begin{array}{*{20}{c}} {5000}&{5000}&{100}&{200}&{200}&{200}&{20} \end{array}} \right]$90.470 2
      530.129 5
    • 本文使用中心差分逐次线性化动力学,提高了序列凸优化算法的规范性和工程性,避免因使用推力方向余弦无损凸化和质量方程取对数凸化而引入额外的优化变量,降低了公式推导工作量和编程复杂性。通过对比证明,相对于传统动力学凸化技巧,直接线性化的求解结果并未对算法收敛性、解的最优性和终端误差造成显著影响。同时,采用性能指标相对偏差作为迭代收敛条件,使收敛条件和具体模型解耦,具有更好的普适性。

      针对动力下降段最优控制问题,给出最优轨迹簇剩余燃料和终端时间的拟合函数,可解析求出近似最优终端时间,能够有效提高在线轨迹方法的效率。

      本文使用的动力学模型较为简单,未考虑更多的过程约束和环境条件,如发动机推力方向约束、落地倾角约束和气动力作用,但可为复杂约束下的轨迹规划算法提供参考初始轨迹,也可为轨迹跟踪控制算法提供标称轨迹。

      由于存在线性化误差和离散误差二者的复合作用,主要采用大批次数值仿真计算来评估算法收敛能力,对序列凸优化算法收敛性的严格数学证明尚比较困难,需要进一步进行深入研究。

参考文献 (11)

目录

    /

    返回文章
    返回