CN114117631A - 一种带有最优终端时间估计的火箭回收轨迹优化方法 - Google Patents

一种带有最优终端时间估计的火箭回收轨迹优化方法 Download PDF

Info

Publication number
CN114117631A
CN114117631A CN202111358336.1A CN202111358336A CN114117631A CN 114117631 A CN114117631 A CN 114117631A CN 202111358336 A CN202111358336 A CN 202111358336A CN 114117631 A CN114117631 A CN 114117631A
Authority
CN
China
Prior art keywords
acceleration
rocket
optimal
terminal
terminal time
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
Application number
CN202111358336.1A
Other languages
English (en)
Other versions
CN114117631B (zh
Inventor
郭杰
相岩
王肖
史鹏飞
唐胜景
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202111358336.1A priority Critical patent/CN114117631B/zh
Publication of CN114117631A publication Critical patent/CN114117631A/zh
Application granted granted Critical
Publication of CN114117631B publication Critical patent/CN114117631B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/62—Systems for re-entry into the earth's atmosphere; Retarding or landing devices
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00—Road transport of goods or passengers
    • Y02T10/10—Internal combustion engine [ICE] based vehicles
    • Y02T10/40—Engine management systems

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Testing Of Engines (AREA)

Abstract

本发明公开了一种带有最优终端时间估计的火箭回收轨迹优化方法，首先基于最优控制理论论证了火箭子级加速度变化趋势，进一步利用终端状态约束给出了终端加速度上界的解析形式，同时考虑到火箭子级着陆段飞行状态、发动机推力调节能力以及喷管质量流量等限制条件，垂直着陆过程中加速度变化率相对其速度变化率(加速度)为小量，选择最大加速度作为
Figure DDA0003358112550000011
下界判断指标。最后，根据火箭子级加速度的连续变化特性，引入n阶多项式形式的终端加速度估计公式，通过比较多项式形式的终端加速度估计值与终端加速度上界给出最优终端时间上下界；本发明能够快速估计最优终端时间，收敛性能良好，具有较高的精度和计算效率，具备在线应用的潜力。

Description

一种带有最优终端时间估计的火箭回收轨迹优化方法
技术领域
本发明属于火箭回收的技术领域，具体涉及一种带有最优终端时间估计的火箭回收轨迹优化方法。
背景技术
自人类开展空间探索活动以来，如何低成本快速进入太空一直是航天运载领域中一大重要研究课题。随着Space X公司猎鹰火箭回收和复用的多次成功，垂直起降重复使用运载火箭方案受到了广泛关注，这一方案采用两级入轨、一子级垂直回收的方式实现子级箭体的重复使用，其典型飞行剖面如图1所示。
根据可重复使用火箭子级的飞行环境与动力学特性，可将垂直回收过程划分为调姿段、修航段、高空下降段、动力减速段、气动减速段以及垂直着陆段。其中着陆段约束复杂、外界扰动强、动力学高动态变化，同时需要考虑经济效益、实现燃料消耗最小飞行，对制导算法在线处理多约束优化问题的能力提出了很高的要求。
现有的轨迹优化方法一般可分为直接法和间接法。间接法基于Pontryagin极小值原理，构造Hamiltonian函数，推导一阶必要条件，进而将问题转化为两点-多点边值问题求解，求解精度高，且最优性可以得到理论证明。但边值问题中的协态变量物理意义不明确且对初始猜测敏感，求解困难。直接法将最优控制问题转化为非线性规划问题，进而采用非线性规划算法对目标函数直接寻优，实现简单，数值解法成熟，目前在轨迹优化领域得到了广泛的应用。
在直接法中，凸优化方法求解速度快、具备确定的收敛准则的特点使其在数值计算方面具有明显的优势。凸优化方法通常先将原问题离散后再进行凸化处理。目前常用的凸化方法主要有序列凸化和无损凸化。
现有的无损凸化方法在确定最优终端时间方面研究较少，多将终端时间tf设置为固定值，无法保证开机－终端时刻组合的最优性。现有方法采用变质量动力学方程估计终端时间的上下界，利用黄金分割法对终端时间进行一维线性搜索，完成了最优终端时间估计。现有的无损凸化方法研究中在确定最优终端时间方面研究较少，估计最优终端时间的计算效率较差。
发明内容
有鉴于此，本发明提供了一种带有最优终端时间估计的火箭回收轨迹优化方法，能够快速估计最优终端时间，收敛性能良好，具有较高的精度和计算效率，具备在线应用的潜力。
实现本发明的技术方案如下：
一种带有最优终端时间估计的火箭回收轨迹优化方法，包括以下步骤：
步骤一：估算终端加速度上界||amax||，求解火箭子级加速度的估计值
Figure BDA0003358112530000021
步骤二：求解最优终端时间
Figure BDA0003358112530000022
的下界对应的多项式次数iL，进而利用iL和iL+1计算得到最优终端时间
Figure BDA0003358112530000023
的上下界
Figure BDA0003358112530000024
步骤三：初始化试探点
Figure BDA0003358112530000025
将t0f、
Figure BDA0003358112530000026
分别代入轨迹优化问题P2(如式(31)所示)，采用hp伪谱策略离散，进而基于内点法求解得到对应的性能指标函数J0f、
Figure BDA0003358112530000027
目标函数：
J＝min-z(tf) (23)
过程约束：
Figure BDA0003358112530000031
边界约束：
Figure BDA0003358112530000032
推力方向及变化率约束：
Figure BDA0003358112530000033
推力幅值约束：
Figure BDA0003358112530000034
动力学方程：
Figure BDA0003358112530000035
轨迹优化问题P2：
min(23)
s.t.(9),(19),(22),(25)～(26) (31)
其中，r、v为火箭子级的位置、速度矢量；Isp为火箭子级发动机比冲；g＝[0,-g0，0]T为重力加速度矢量；T为发动机推力矢量，Γ为松弛变量，
Figure BDA0003358112530000036
Tmin、Tmax为推力幅值上下界；m为火箭子级质量，mdry为结构质量，z＝lnm，
Figure BDA0003358112530000041
t为飞行时间，z(tf)为终端时刻的z值；γ为推力摆角；ρ为大气密度；qmax为最大动压；下标0、f分别表示初始和终端值，而下标x、y、z则表示沿着陆点坐标系三轴方向的分量；问题中，状态变量x＝[r,v,z]T，控制变量u＝[uc,σ]T。
步骤四：基于(t0f,J0f)、
Figure BDA0003358112530000042
三点信息采用二次插值法寻优，得到
Figure BDA0003358112530000043
将
Figure BDA0003358112530000044
代入轨迹优化问题P2，采用hp伪谱策略离散，进而基于内点法求解得到启动解
Figure BDA0003358112530000045
步骤五：初始化同伦参数ε＝0，置迭代次数k＝1，读取
Figure BDA0003358112530000046
与启动解
Figure BDA0003358112530000047
采用hp伪谱策略离散，进而基于内点法求解同伦迭代轨迹优化问题P3(如式(32)所示)，得到最优解
Figure BDA0003358112530000048
同伦迭代动力学方程：
Figure BDA0003358112530000049
推力幅值约束：
Figure BDA00033581125300000410
过程约束：
Figure BDA00033581125300000411
轨迹优化问题P3：
Figure BDA0003358112530000051
其中，k为迭代次数，ρk、vk、zk为第k次迭代解，
Figure BDA0003358112530000052
为第k次迭代求得的阻力加速度和升力加速度。
步骤六：令ε＝ε+hε，hε为同伦迭代步长，将
Figure BDA0003358112530000053
作为启动解，采用hp伪谱策略离散，进而基于内点法求解轨迹优化问题P3，令k＝k+1。
步骤七：判断同伦参数ε≥1是否满足，如是则同伦过程完成，令ε＝1并执行步骤八，否则返回执行步骤六。
步骤八：判断最优解收敛准则
Figure BDA0003358112530000054
(xε为收敛常数)是否满足，如是执行步骤九，否则将当前结果作为启动解，继续迭代求解轨迹优化问题P3，直至收敛。
步骤九：输出最优轨迹
Figure BDA0003358112530000055
更新计算时间。
进一步地，步骤一具体为：
根据式(56)估算终端加速度上界||amax||，根据式(59)和式(61)求解火箭子级加速度的估计值
Figure BDA0003358112530000056
Figure BDA0003358112530000057
Figure BDA0003358112530000058
Figure BDA0003358112530000059
其中，Tmax为火箭子级最大推力；g0为重力加速度；mdry为结构质量；i为多项式阶数；tif为采用i阶多项式估计的终端时间；||r0||为垂直着陆段火箭子级初始位置；||rf||为火箭着陆的终端位置；||v0||为垂直着陆段火箭子级初始速度。
进一步地，步骤二具体为：
根据式(62)得到最优终端时间
Figure BDA0003358112530000061
的下界对应的多项式次数iL，进而将iL和iL+1代入式(61)计算得到最优终端时间
Figure BDA0003358112530000062
的上下界
Figure BDA0003358112530000063
Figure BDA0003358112530000064
其中，N*为正整数集；函数argmin{A|B}表示满足条件B的A的最小值。
进一步地，所述hp伪谱策略离散，具体为：
时域变换：
Figure BDA0003358112530000065
动力学离散形式：
Figure BDA0003358112530000066
连接条件：
Figure BDA0003358112530000067
其中，τ为变换后的时间变量；h为子区间个数；p为Lagrange多项式次数；Dlj为FRPM微分矩阵第l行j列的元素；f1[·]为系统动力学方程；上标n表示第n个子区间，下标j表示子区间内的第j个离散点。
有益效果：
1、本发明提出了一套基于解析推导与二次插值的最优终端时间快速估计策略。
针对实际飞行任务中应用无损凸化方法存在最优终端时间上下界估计精度差、线性搜索方法收敛速度慢、需要多次迭代求解等问题，本发明在最优终端时间估计过程中充分分析了垂直着陆段火箭子级的飞行特性，首先基于最优控制理论论证了火箭子级加速度变化趋势，进一步利用终端状态约束给出了终端加速度上界的解析形式，同时考虑到火箭子级着陆段飞行状态、发动机推力调节能力以及喷管质量流量等限制条件，垂直着陆过程中加速度变化率相对其速度变化率(加速度)为小量，选择最大加速度作为
Figure BDA0003358112530000071
下界判断指标。最后，根据火箭子级加速度的连续变化特性，引入n阶多项式形式的终端加速度估计公式，通过比较多项式形式的终端加速度估计值与终端加速度上界给出最优终端时间上下界。
另一方面，本发明针对线性搜索方法收敛速度慢的问题，在得到相对精确的最优终端时间上下界的条件下，结合凸优化问题中凸函数优良的单峰性质，基于二次插值法对最优终端时间进行寻优，二次插值法自身具备的超线性收敛的特性有效提升了最优终端时间的估计速度。
2、本发明基于hp多阶段伪谱离散策略设计了轨迹优化方法，在保证一定的求解精度的条件下有效提升了计算速度，使得该方法具备在线应用的潜力。
本发明针对火箭回收任务着陆段燃料最省轨迹优化问题的在线求解需求，在伪谱凸优化方法的基础上增加hp多阶段策略，使得精度控制参数由原来的全局多项式次数一个参数拓展为子区间个数h和多项式次数p两个参数，增强了轨迹优化方法的精度控制能力，提升了方法的灵活性。
另一方面，相较于全局伪谱离散策略，基于计算流体力学的网格划分思想的hp多阶段伪谱离散策略使得轨迹优化问题的一阶偏微分矩阵(Jacobi矩阵)与二阶偏微分矩阵(Hessian矩阵)的稀疏性得到了极大的提升，更稀疏的问题形式有效缩减了优化方法的计算时间。同时采用伪谱法选取离散点能够有效避免数值近似的龙格现象，抑制了多项式阶数较高时出现的局部饱和现象。此外，伪谱法引入高斯型积分权重可以在相同离散点数量的条件下达到最高的近似精度。综上所述，hp多阶段伪谱离散策略兼具了hp网格细化的稀疏特性与伪谱离散的高精度特性，使得基于hp多阶段伪谱离散策略设计的轨迹优化方法能够在保证一定的求解精度的条件下有效缩减了计算时间，具备在线求解轨迹优化问题的潜力。
附图说明
图1为火箭回收飞行剖面图。
图2为着陆点坐标系示意图。
图3为火箭回收着陆段所受作用力幅值示意图。
图4为最优终端时间上下界解析多项式估计示意图。
图5为带有最优终端时间估计的火箭回收轨迹优化方法流程图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本发明提供了一种带有最优终端时间估计的火箭回收轨迹优化方法，可以快速估计最优终端时间，收敛性能良好，具有较高的精度和计算效率，具备在线应用的潜力。
1问题描述
(1)动力学模型
不考虑地球自转，将地球视为均匀圆球，忽略哥氏惯性力的影响，考虑阻力和升力，定义着陆点坐标系如图2所示。图中，坐标系原点O为着陆点，Oy轴垂直地面向上，Ox轴与初始时刻的速度矢量在地面上的投影的方向一致，Oz轴与Ox、Oy轴呈右手坐标系。
在着陆点坐标系下建立可重复使用运载火箭子级着陆段三自由度动力学方程如下
Figure BDA0003358112530000091
式中：r＝[rx,ry,rz]T和v＝[vx,vy,vz]T分别为着陆点坐标系下的火箭子级质心位置和速度矢量，T＝[Tx,Ty,Tz]T为发动机推力矢量，
Figure BDA0003358112530000092
m为火箭子级质量，Isp为火箭子级发动机比冲，g＝[0,-g0，0]T为重力加速度矢量，其中g0＝9.81m/s2；aL和aD分别为升力加速度矢量和阻力加速度矢量，其计算公式为
Figure BDA0003358112530000093
式中：CL和CD分别为升力系数和阻力系数，Sref为火箭子级参考特征面积，ρ为大气密度，其计算公式为
ρ＝ρ0exp(-ry/Href) (3)
式中：ρ0为海平面密度，Href为参考高度，ry为火箭子级飞行高度。
对于火箭垂直回收过程，令状态变量
Figure BDA0003358112530000094
控制变量
Figure BDA0003358112530000095
则系统状态空间表达式为
Figure BDA0003358112530000096
式中：
Figure BDA0003358112530000097
Figure BDA0003358112530000101
(2)轨迹优化问题
火箭子级垂直回收需满足定点软着陆的任务需求，因此定义轨迹优化问题的初始和终端约束为
Figure BDA0003358112530000102
式中：t0和tf分别为着陆段初始和终端时刻；m0和mdry分别为火箭子级着陆段初始质量和结构干重；对于运载火箭子级一般有vf＝0，rf＝0。进一步考虑控制量约束，其中推力幅值和方向约束为
Figure BDA0003358112530000103
Figure BDA0003358112530000104
式中：Tmax和Tmin为总推力幅值的上下限，γxmax和γzmax为最大推力方向角。
此外，考虑到火箭子级发动机推力调节能力有限，还需要对推力变化率施加约束
Figure BDA0003358112530000105
另一方面，由于着陆段起始高度较低，火箭子级飞行速度相对较小，而空气密度较大，因此仅考虑动压作为过程约束
Figure BDA0003358112530000106
选择燃料消耗最少即终端剩余质量最大为性能指标
J＝-m(tf) (10)
综上所述，可重复使用运载火箭子级着陆段轨迹优化问题可描述为
P0：
min(10)
s.t.(4)～(9) (11)
2带有最优终端时间估计的轨迹优化方法
(1)问题凸化
①无损凸化
针对P0中非凸形式的动力学方程和约束，无法直接应用凸优化方法求解。本节主要考虑通过无损凸化将非凸问题转化为凸问题。
首先针对非凸推力幅值约束，引入变量Γ，则式(6)可以转化为
Tmin≤Γ≤Tmax (12)
||T||＝Γ (13)
基于凸优化理论和问题的具体形式，对式(13)进行松弛变换，变换后的约束为
||T||≤Γ (14)
则动力学方程可变换为
Figure BDA0003358112530000111
基于上述变换，P0中控制变量可增广变换为uaug＝[T,Γ]T，定义变换后的问题为P1，具体形式如下
P1：
min(10)
s.t.(5),(7)～(9),(12),(14)～(15) (16)
进一步对上述动力学方程中的非线性项进行处理，作如下变量替换
Figure BDA0003358112530000121
变换后的动力学方程如下
Figure BDA0003358112530000122
式中：
Figure BDA0003358112530000123
相应地，边界约束、推力约束和目标函数变换为
Figure BDA0003358112530000124
Tmin·e-z≤σ≤Tmax·e-z (20)
||uc||≤σ (21)
Figure BDA0003358112530000125
J＝min-z(tf) (23)
根据凸优化和解析几何理论，易知式(22)中第一式所定义的可行域为二阶锥，是一种典型的凸约束。
②同伦方法
对于式(18)中的非凸气动力项和式(20)所表示的非凸约束，采用同伦方法(Homotopic approach)对其进行凸化。在数学上，两个拓扑空间如果可以通过一系列连续变化从一个变到另一个，那么就称这两个拓扑空间同伦。对于函数而言，若存在连续映射h(x,ε)，使得连续函数f(x)、g(x)满足
Figure BDA0003358112530000126
则称f(x)与g(x)同伦，h(x,ε)为同伦映射，ε为同伦因子，由上式可知，随着同伦因子ε从0变化到1，函数f(x)逐渐延拓为g(x)。
值得注意的是，应用同伦方法的前提条件是同伦过程中不改变系统的性质，下面对本问题中同伦方法的适用性进行简要分析。
不同于以气动力作为主要控制力的滑翔再入飞行器，垂直回收运载火箭的主要控制力为火箭发动机推力。此外，在垂直回收过程中火箭子级飞行速度相对较小，气动力相对发动机推力要小很多，对应的典型作用力曲线如图3所示。
从图3中可以看出，气动力小于推力，且由于推力对飞行速度的主要控制作用，气动力曲线也在一定程度上受到了推力曲线的影响。
通过上述分析可知，气动力为火箭子级的非主要控制力，可以应用同伦方法凸化问题。同伦凸优化算法的具体设计过程如下。
首先不考虑气动力的影响，即在式(18)中删去气动力加速度项，对式(20)进行Taylor展开，将推力幅值约束变换为
Figure BDA0003358112530000131
式中：
Figure BDA0003358112530000132
t为飞行时间。
相应地，动力学方程可变换为
Figure BDA0003358112530000133
由式(25)可得，不等式左侧为二阶锥约束、右侧为线性约束，对于凸化后的问题应用原-对偶内点法(Primal-Dual Interior-Point Method,PDIPM)求解凸优化问题获得初始解。随后，在动力学方程中引入同伦因子ε，将动力学方程、推力约束以及过程约束变换为
Figure BDA0003358112530000141
Figure BDA0003358112530000142
Figure BDA0003358112530000143
式中：k为迭代次数，vk、zk为第k次迭代解，
Figure BDA0003358112530000144
为第k+1次迭代所需的阻力加速度和升力加速度。即利用上一次的迭代结果构造剖面近似动力学方程和约束中的非线性项，使得式(27)～(28)转化为线性约束、式(29)转化为旋转二次锥约束，进而完成动力学方程及约束的凸化。随着同伦因子ε从0变化到1，完成从无气动力轨迹优化问题至考虑气动力轨迹优化问题的过渡。
定义变换后的动力学系统的状态变量x＝[r,v,z]T，控制变量u＝[uc,σ]T，系统状态空间表达式为
Figure BDA0003358112530000145
式中：
Figure BDA0003358112530000146
Figure BDA0003358112530000147
需要说明的是，上述状态空间表达式中的系数矩阵A1与控制矩阵B1在同伦过程中为线性时变剖面形式，具体取值由上一次优化结果决定。
综上所述，凸化后的问题可以描述为
无气动力轨迹优化问题P2：
min(23)
s.t.(9),(19),(22),(25)～(26) (31)
同伦迭代轨迹优化问题P3：
min(23)
s.t.(19),(22),(27)～(29) (32)
(2)问题离散
本方案采用flipped Radau伪谱法(flipped Radau Pseudospectral Method,FRPM)离散连续最优控制问题P2与P3，将离散后的参数优化问题分别记为P4和P5。FRPM是一种非对称的多区间伪谱法，不同于RPM，FRPM将终端时间设置为配点，可以更好地定义终端约束。
在FRPM中，首先定义新的时间变量τ，将时间区间从[0,tf]映射到[-1,1]，定义仿射变换关系如下
Figure BDA0003358112530000151
然后选取N个Flipped Legendre Gauss Radau(FLGR)配点以及τ0＝-1作为离散点，故在τ∈[-1,1]上共有N个配点和N+1个离散点。定义N阶Lagrange多项式Ll(τ)为如下形式
Figure BDA0003358112530000152
相应地，状态量和控制量可以近似为
Figure BDA0003358112530000153
对上式中的状态量求导即可得到状态量微分的近似形式
Figure BDA0003358112530000161
因此动力学方程式(26)与式(27)的离散形式分别为
Figure BDA0003358112530000162
Figure BDA0003358112530000163
式中：Djl为FRPM微分矩阵第j行l列的元素，f1[·]、f2[·]分别为动力学方程式(26)与式(27)。此外，目标函数为Mayer型性能指标，其离散形式为
Figure BDA0003358112530000164
相应地，初始和终端约束变换为
Figure BDA0003358112530000165
推力方向及变化率约束变换为
Figure BDA0003358112530000166
对于无气动力轨迹优化问题P2，推力幅值约束可变换为
Figure BDA0003358112530000167
过程约束可变换为
Figure BDA0003358112530000168
而对于同伦迭代轨迹优化问题P3，推力幅值约束可变换为
Figure BDA0003358112530000171
过程约束变换为
Figure BDA0003358112530000172
进一步地，为加快求解速度、提升求解效率，采用hp伪谱法代替全局伪谱法对问题进行离散化处理。采用hp-FRPM离散轨迹优化问题增加了雅各比矩阵的稀疏性，加快了轨迹优化问题的求解速度。
在本发明方案中用上标n表示第n个子区间，下标j表示子区间内的第j个离散点，则动力学方程式(37)与式(38)的hp离散形式分别为
Figure BDA0003358112530000173
Figure BDA0003358112530000174
此外，考虑各子区间之间的连接条件
Figure BDA0003358112530000175
综上所述，采用hp-FRPM离散后的参数优化问题可描述为
P4：
min(39)
s.t.(40)～(43),(46),(48) (49)
P5：
min(39)
s.t.(40)～(41),(44)～(45),(47)～(48) (50)
(3)最优终端时间快速估计策略
现有的无损凸化方法在确定最优终端时间方面研究较少，多将tf设置为固定值，在一定程度上对求解结果施加了限制条件，不能完全保证求解结果的最优性。
但是，若直接将tf纳入轨迹优化问题框架则会破坏原问题良好的收敛特性。另一方面，若采用传统的线性搜索算法对tf寻优则存在需要多次求解轨迹优化问题，导致搜索算法收敛速度慢，难以在线应用。
基于上述分析，本方案基于最优控制理论和垂直起降可重复使用运载火箭的运动学特性采用解析方法估计包含最优解区间的上下界，进一步设计二次插值型搜索策略完成最优终端时间的快速估计。
首先根据图3中对火箭垂直回收问题中作用力的分析结果，利用气动力为非主要控制力的特性，将无气动力轨迹优化问题对应的最优终端时间
Figure BDA0003358112530000181
作为考虑气动力的轨迹优化问题对应的最优终端时间
Figure BDA0003358112530000182
的有效近似，即认为
Figure BDA0003358112530000183
然后确定最优终端时间
Figure BDA0003358112530000184
所在区间的上下界。由于存在多种严格约束，仅采用质量流量公式估计
Figure BDA0003358112530000185
的上下界容易造成问题不可行，进而导致线性搜索算法初始化困难。此外，仅采用质量流量公式估计
Figure BDA0003358112530000186
的上下界还存在估计区间长度过大，影响线性搜索算法收敛速度的问题。因此本方案将加速度最大值作为判断指标，基于运动学公式和质量流量公式综合估计
Figure BDA0003358112530000187
的上下界。根据最优控制原理，燃料最优控制问题解的形式如下
Figure BDA0003358112530000188
式中：k为切换系数。
进一步地，由火箭子级动力学式(1)中第二式可得，不考虑气动力时火箭子级加速度幅值为
Figure BDA0003358112530000191
虽然火箭子级加速度幅值为时变值，难以分析其具体变化趋势，但结合式(52)中的燃料最优控制形式和火箭子级变质量飞行特性，可以通过定性分析式(53)得出，在垂直回收过程中加速度总体呈增大趋势，终端时刻加速度最大。同时，根据终端时刻火箭子级姿态约束和最优控制理论可知，终端时刻火箭子级推力方向沿着陆点坐标系y轴方向且幅值为Tmax。通过上述分析可得，火箭子级终端加速度具体形式为
Figure BDA0003358112530000192
对于上式中的终端质量mf，结合动力学方程式(1)中的第三式和燃料最优控制问题解式(52)估计其下界
Figure BDA0003358112530000193
则终端加速度的上界为
Figure BDA0003358112530000194
至此，完成终端加速度上界的解析估计。另一方面，由于火箭子级初始位置和速度已知，根据运动学关系，构建多项式函数来估计速度，进而求解最优终端时间上下界，原理如图4所示。
上图中实际速度曲线为基于典型火箭子级参数得到的仿真结果，i次多项式分别对应以(0,||v0||)为顶点的i次多项式函数，其函数表达式为
Figure BDA0003358112530000195
式中：
Figure BDA0003358112530000196
为t时刻多项式函数对速度的估计值；tif为i阶多项式形式的速度对应的终端时间。对上式求导得到加速度的估计
Figure BDA0003358112530000201
由上式可知，终端时刻tif火箭子级加速度的估计值为
Figure BDA0003358112530000202
进一步根据运动学公式，速度曲线对时间的积分为位置，即
Figure BDA0003358112530000203
将式(57)代入上式求解可得
Figure BDA0003358112530000204
另一方面，从图4中可以看出，随着多项式次数的增加，对应的终端时间逐渐缩短。同时由于火箭子级加速度总体呈增大趋势，易知t1f为
Figure BDA0003358112530000205
的确保可行的上界。另一方面，通过比较多项式形式的终端加速度式(59)与式(56)的大小即可得出
Figure BDA0003358112530000206
的下界对应的多项式次数iL，即
Figure BDA0003358112530000207
式中：N*为正整数集，通过上式即可求出最优终端时间
Figure BDA0003358112530000208
的上下界
Figure BDA0003358112530000209
需要说明的是，考虑到火箭子级着陆段飞行状态、发动机推力调节能力以及喷管质量流量等限制条件，垂直着陆过程中加速度变化率相对其速度变化率(加速度)为小量，这也是本方案选择最大加速度作为
Figure BDA00033581125300002010
下界判断指标的主要原因。
进一步利用凸函数的单峰性质，采用二次插值法对
Figure BDA00033581125300002011
寻优。二次插值法将目标函数拟合为二次函数，进而利用二次函数的性质求解最优解。
至此，完成最优终端时间
Figure BDA0003358112530000211
的快速估计。值得注意的是，本方案关于最优终端时间的估计推导为全解析形式，计算量小；同时采用二次插值法可以有效利用凸函数优良的收敛特性，相较于一般的线性搜索方法收敛速度更快。
综上所述，带有最优终端时间估计的火箭回收轨迹优化方法流程如图5所示。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种带有最优终端时间估计的火箭回收轨迹优化方法，其特征在于，包括以下步骤：
步骤一：估算终端加速度上界||amax||，求解火箭子级加速度的估计值
Figure FDA0003358112520000011
步骤二：求解最优终端时间
Figure FDA0003358112520000012
的下界对应的多项式次数iL，进而利用iL和iL+1计算得到最优终端时间
Figure FDA0003358112520000013
的上下界
Figure FDA0003358112520000014
步骤三：初始化试探点
Figure FDA0003358112520000015
将t0f、
Figure FDA0003358112520000016
分别代入轨迹优化问题P2，采用hp伪谱策略离散，进而基于内点法求解得到对应的性能指标函数J0f、
Figure FDA0003358112520000017
目标函数：
J＝min-z(tf) (23)
过程约束：
Figure FDA0003358112520000018
边界约束：
Figure FDA0003358112520000019
推力方向及变化率约束：
Figure FDA00033581125200000110
推力幅值约束：
Figure FDA0003358112520000021
动力学方程：
Figure FDA0003358112520000022
轨迹优化问题P2：
min(23)
s.t.(9),(19),(22),(25)～(26) (31)
其中，r、v为火箭子级的位置、速度矢量；Isp为火箭子级发动机比冲；g＝[0,-g0，0]T为重力加速度矢量；T为发动机推力矢量，Γ为松弛变量，
Figure FDA0003358112520000023
Tmin、Tmax为推力幅值上下界；m为火箭子级质量，mdry为结构质量，z＝lnm，
Figure FDA0003358112520000024
t为飞行时间，z(tf)为终端时刻的z值；γ为推力摆角；ρ为大气密度；qmax为最大动压；下标0、f分别表示初始和终端值，而下标x、y、z则表示沿着陆点坐标系三轴方向的分量；问题中，状态变量x＝[r,v,z]T，控制变量u＝[uc,σ]T；
步骤四：基于(t0f,J0f)、
Figure FDA0003358112520000025
三点信息采用二次插值法寻优，得到
Figure FDA0003358112520000026
将
Figure FDA0003358112520000027
代入轨迹优化问题P2，采用hp伪谱策略离散，进而基于内点法求解得到启动解
Figure FDA0003358112520000028
步骤五：初始化同伦参数ε＝0，置迭代次数k＝1，读取
Figure FDA0003358112520000029
与启动解
Figure FDA0003358112520000031
采用hp伪谱策略离散，进而基于内点法求解同伦迭代轨迹优化问题P3，得到最优解
Figure FDA0003358112520000032
同伦迭代动力学方程：
Figure FDA0003358112520000033
推力幅值约束：
Figure FDA0003358112520000034
过程约束：
Figure FDA0003358112520000035
轨迹优化问题P3：
min(23)
s.t.(19),(22),(27)～(29) (32)
其中，k为迭代次数，ρk、vk、zk为第k次迭代解，
Figure FDA0003358112520000036
为第k次迭代求得的阻力加速度和升力加速度；
步骤六：令ε＝ε+hε，hε为同伦迭代步长，将
Figure FDA0003358112520000037
作为启动解，采用hp伪谱策略离散，进而基于内点法求解轨迹优化问题P3，令k＝k+1；
步骤七：判断同伦参数ε≥1是否满足，如是则同伦过程完成，令ε＝1并执行步骤八，否则返回执行步骤六；
步骤八：判断最优解收敛准则
Figure FDA0003358112520000038
是否满足，xε为收敛常数，如是执行步骤九，否则将当前结果作为启动解，继续迭代求解轨迹优化问题P3，直至收敛；
步骤九：输出最优轨迹
Figure FDA0003358112520000041
更新计算时间。
2.如权利要求1所述的一种带有最优终端时间估计的火箭回收轨迹优化方法，其特征在于，步骤一具体为：
根据式(56)估算终端加速度上界||amax||，根据式(59)和式(61)求解火箭子级加速度的估计值
Figure FDA0003358112520000042
Figure FDA0003358112520000043
Figure FDA0003358112520000044
Figure FDA0003358112520000045
其中，Tmax为火箭子级最大推力；g0为重力加速度；mdry为结构质量；i为多项式阶数；tif为采用i阶多项式估计的终端时间；||r0||为垂直着陆段火箭子级初始位置；||rf||为火箭着陆的终端位置；||v0||为垂直着陆段火箭子级初始速度。
3.如权利要求1所述的一种带有最优终端时间估计的火箭回收轨迹优化方法，其特征在于，步骤二具体为：
根据式(62)得到最优终端时间
Figure FDA0003358112520000046
的下界对应的多项式次数iL，进而将iL和iL+1代入式(61)计算得到最优终端时间
Figure FDA0003358112520000047
的上下界
Figure FDA0003358112520000048
Figure FDA0003358112520000049
其中，N*为正整数集；函数argmin{A|B}表示满足条件B的A的最小值。
4.如权利要求1所述的一种带有最优终端时间估计的火箭回收轨迹优化方法，其特征在于，所述hp伪谱策略离散，具体为：
时域变换：
Figure FDA0003358112520000051
动力学离散形式：
Figure FDA0003358112520000052
连接条件：
Figure FDA0003358112520000053
其中，τ为变换后的时间变量；h为子区间个数；p为Lagrange多项式次数；Dlj为FRPM微分矩阵第l行j列的元素；f1[·]为系统动力学方程；上标n表示第n个子区间，下标j表示子区间内的第j个离散点。
CN202111358336.1A 2021-11-16 2021-11-16 一种带有最优终端时间估计的火箭回收轨迹优化方法 Active CN114117631B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111358336.1A CN114117631B (zh) 2021-11-16 2021-11-16 一种带有最优终端时间估计的火箭回收轨迹优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111358336.1A CN114117631B (zh) 2021-11-16 2021-11-16 一种带有最优终端时间估计的火箭回收轨迹优化方法

Publications (2)

Publication Number Publication Date
CN114117631A true CN114117631A (zh) 2022-03-01
CN114117631B CN114117631B (zh) 2024-06-28

Family

ID=80396833

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111358336.1A Active CN114117631B (zh) 2021-11-16 2021-11-16 一种带有最优终端时间估计的火箭回收轨迹优化方法

Country Status (1)

Country Link
CN (1) CN114117631B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116892866A (zh) * 2023-07-25 2023-10-17 东方空间技术(山东)有限公司 一种火箭子级回收轨迹规划方法、设备及存储介质

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3506041A1 (en) * 2017-12-29 2019-07-03 Deutsches Zentrum für Luft- und Raumfahrt e.V. Method, apparatus and spacecraft for constrained atmospheric entry
CN110466804A (zh) * 2019-08-30 2019-11-19 北京理工大学 火箭动力下降着陆过程快速轨迹优化方法
CN111196382A (zh) * 2019-12-25 2020-05-26 北京理工大学 保证收敛的火箭动力下降段实时轨迹规划方法
CN111783232A (zh) * 2020-07-23 2020-10-16 中国人民解放军国防科技大学 一种基于聚类分析的可回收火箭返回段自适应优化方法
CN112001029A (zh) * 2020-07-28 2020-11-27 清华大学 一种基于凸优化的火箭在线轨迹优化定制化求解器
CN112507461A (zh) * 2020-12-15 2021-03-16 北京航天自动控制研究所 一种运载火箭动力软着陆段发动机开机方法
CN112520071A (zh) * 2020-12-17 2021-03-19 清华大学 一种可回收火箭动力段燃料最优着陆轨迹快速规划方法
CN112550770A (zh) * 2020-12-15 2021-03-26 北京航天自动控制研究所 一种基于凸优化的火箭软着陆轨迹规划方法
CN112629339A (zh) * 2020-12-15 2021-04-09 北京航天自动控制研究所 一种基于直接法的火箭软着陆轨迹规划方法
US20210164783A1 (en) * 2019-08-29 2021-06-03 Dalian University Of Technology Method for directly planning reentry trajectory in height-velocity profile
CN113479347A (zh) * 2021-07-13 2021-10-08 北京理工大学 一种火箭垂直回收着陆段轨迹控制方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3506041A1 (en) * 2017-12-29 2019-07-03 Deutsches Zentrum für Luft- und Raumfahrt e.V. Method, apparatus and spacecraft for constrained atmospheric entry
US20210164783A1 (en) * 2019-08-29 2021-06-03 Dalian University Of Technology Method for directly planning reentry trajectory in height-velocity profile
CN110466804A (zh) * 2019-08-30 2019-11-19 北京理工大学 火箭动力下降着陆过程快速轨迹优化方法
CN111196382A (zh) * 2019-12-25 2020-05-26 北京理工大学 保证收敛的火箭动力下降段实时轨迹规划方法
CN111783232A (zh) * 2020-07-23 2020-10-16 中国人民解放军国防科技大学 一种基于聚类分析的可回收火箭返回段自适应优化方法
CN112001029A (zh) * 2020-07-28 2020-11-27 清华大学 一种基于凸优化的火箭在线轨迹优化定制化求解器
CN112507461A (zh) * 2020-12-15 2021-03-16 北京航天自动控制研究所 一种运载火箭动力软着陆段发动机开机方法
CN112550770A (zh) * 2020-12-15 2021-03-26 北京航天自动控制研究所 一种基于凸优化的火箭软着陆轨迹规划方法
CN112629339A (zh) * 2020-12-15 2021-04-09 北京航天自动控制研究所 一种基于直接法的火箭软着陆轨迹规划方法
CN112520071A (zh) * 2020-12-17 2021-03-19 清华大学 一种可回收火箭动力段燃料最优着陆轨迹快速规划方法
CN113479347A (zh) * 2021-07-13 2021-10-08 北京理工大学 一种火箭垂直回收着陆段轨迹控制方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116892866A (zh) * 2023-07-25 2023-10-17 东方空间技术(山东)有限公司 一种火箭子级回收轨迹规划方法、设备及存储介质
CN116892866B (zh) * 2023-07-25 2024-01-23 东方空间技术(山东)有限公司 一种火箭子级回收轨迹规划方法、设备及存储介质

Also Published As

Publication number Publication date
CN114117631B (zh) 2024-06-28

Similar Documents

Publication Publication Date Title
CN111897214B (zh) 一种基于序列凸优化的高超声速飞行器轨迹规划方法
Hesse et al. Reduced-order aeroelastic models for dynamics of maneuvering flexible aircraft
Mir et al. Optimization of dynamic soaring maneuvers to enhance endurance of a versatile UAV
CN110309571B (zh) 基于径向基函数模型的翼身融合水下滑翔机外形优化方法
Lane et al. Inverse airfoil design utilizing CST parameterization
Fujikawa et al. Multi-objective, multidisciplinary design optimization of TSTO space planes with RBCC engines
HOLST et al. The NASA Computational Aerosciences Program-Toward Teraflops Computing
Peng et al. Analysis of morphing modes of hypersonic morphing aircraft and multiobjective trajectory optimization
Leifsson et al. Variable-fidelity aerodynamic shape optimization
CN114117631A (zh) 一种带有最优终端时间估计的火箭回收轨迹优化方法
Li et al. A segmented and weighted adaptive predictor-corrector guidance method for the ascent phase of hypersonic vehicle
Jung et al. Low Reynolds number airfoil design for a mars exploration airplane using a transition model
Pescetelli et al. Ascent trajectory optimisation for a single-stage-to-orbit vehicle with hybrid propulsion
Cameron et al. Metamodel assisted multi-objective global optimisation of natural laminar flow aerofoils
Leifsson et al. Multi-fidelity design optimization of transonic airfoils using shape-preserving response prediction
Zhang et al. Integral performance optimization for the two-stage-to-orbit RBCC-RKT launch vehicle based on GPM
Yao et al. Trajectory tracking control of a stratospheric airship with fuzzy sliding mode control
CN109117544B (zh) 一种天地往返飞行器全轨迹的优化方法
Klock Efficient Numerical Simulation of Aerothermoelastic Hypersonic Vehicles
Kier et al. Development of aircraft flight loads analysis models with uncertainties for pre-design studies
Du et al. Application of iterative learning tuning to a dragonfly-like flapping wing micro aerial vehicle
Zhu et al. Air-breathing combined power flight vehicle trajectory optimization method based on the aerodynamic surrogate model
Saranathan Algorithmic Advances to Increase the Fidelity of Conceptual Hypersonic Mission Design
Avanzini et al. Evolutionary design of a full-envelope flight control system for an unstable fighter aircraft
Chen et al. Robust neighboring optimal guidance for endoatmospheric powered descent under uncertain wind fields

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