CN111580535B - 一种基于凸优化的再入轨迹三维剖面规划方法及系统 - Google Patents

一种基于凸优化的再入轨迹三维剖面规划方法及系统 Download PDF

Info

Publication number
CN111580535B
CN111580535B CN202010419723.0A CN202010419723A CN111580535B CN 111580535 B CN111580535 B CN 111580535B CN 202010419723 A CN202010419723 A CN 202010419723A CN 111580535 B CN111580535 B CN 111580535B
Authority
CN
China
Prior art keywords
track
current
max
angle
convex optimization
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
Application number
CN202010419723.0A
Other languages
English (en)
Other versions
CN111580535A (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202010419723.0A priority Critical patent/CN111580535B/zh
Publication of CN111580535A publication Critical patent/CN111580535A/zh
Application granted granted Critical
Publication of CN111580535B publication Critical patent/CN111580535B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
    • G05D1/08—Control of attitude, i.e. control of roll, pitch, or yaw
    • G05D1/0808—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
    • G05D1/10—Simultaneous control of position or course in three dimensions
    • G05D1/101—Simultaneous control of position or course in three dimensions specially adapted for aircraft
    • 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
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Feedback Control In General (AREA)

Abstract

本发明公开了一种基于凸优化的再入三维剖面规划方法及系统，针对再入飞行器在滑翔段的轨迹规划问题，提出了一种基于凸优化的三维剖面规划方法。首先，建立飞行器的三自由度运动模型，并选择合适的中间变量表征三维剖面，以此作为轨迹规划的控制量。其次，通过合适的凸化与离散方法，将三维剖面规划问题建模为凸优化问题，利用凸优化求解器进行求解。最后，利用凸优化问题的最优解反算得到攻角和倾侧角指令。本发明能有效解决滑翔轨迹的三维剖面规划求解效率低的难题。

Description

一种基于凸优化的再入轨迹三维剖面规划方法及系统
技术领域
本发明涉及飞行器制导控制技术领域，可用于再入飞行器滑翔段轨迹规划与制导，尤其涉及一种基于凸优化的再入轨迹三维剖面规划方法及系统。
背景技术
高超声速滑翔飞行器是一种主要在临近空间飞行，且飞行速度大于5马赫的飞行器。相比于传统弹道式飞行器，此类飞行器具有机动能力强、飞行速域宽、弹道全程可控等优点，在军事领域具有广阔应用价值。该类飞行器的再入轨迹可分为初始下降段和滑翔段，滑翔段由于飞行时间长、航程远、过程约束复杂等，是轨迹规划问题研究的重点和难点。传统轨迹规划方法的不足主要包括：1)基于极大值原理的间接法需要复杂公式推导，难以处理复杂约束；2)以伪谱法为代表的直接法，计算效率低，对初值比较敏感；3)传统滑翔段轨迹规划方法多是基于二维剖面，即固定攻角剖面，将倾侧角作为轨迹规划的唯一控制量，由于参数优化范围小，所以限制了飞行器的可用机动能力。在基于三维剖面的轨迹规划中，将攻角和倾侧角同时作为控制量，增大了参数优化范围，从而充分释放了飞行器的可用机动能力。但传统基于三维剖面的规划方法为了降低计算量采用了大量模型简化方法，降低了解的精度，未能很好解决精度与计算效率之间的矛盾。为了克服该不足，研究能够充分释放飞行器机动能力，且满足在线快速计算要求的滑翔段轨迹规划方法尤为重要。
凸优化是一类具备良好理论性质和高效求解效率的数学规划问题。由于该问题能够在多项式时间内收敛到全局最优解，目前已成为在线求解轨迹规划等实时最优控制问题的常用方法之一。基于凸优化求解轨迹规划问题主要包括三个步骤：1)将原始轨迹规划问题转化为一个连续的凸优化问题，包括对目标函数、非线性动力学、控制量约束和状态量约束等原问题中非凸项的凸化处理；2)将上述连续凸优化问题进行离散化，得到一个以有限个参数为优化变量的凸优化问题；3)利用凸优化求解器进行问题求解，得到控制量。凸优化方法的难点主要在于第一步，如何对原问题进行凸化处理，目前常用的方法包括变量替换、松弛等，针对实际问题需要灵活采用多种凸化方法进行处理。特别是对于一些具有强非线性动力学或状态约束的轨迹规划问题，就目前所知，连续线性化是处理这类约束的仅有方法，这种情况下由于需要求解一个凸优化子问题序列，以这个序列的收敛解作为原问题的最优解，因此该方法又称为序列凸优化。由于再入轨迹规划问题具有强非线性动力学与非线性状态约束，因此序列凸优化是求解这类问题的主要手段。原始对偶内点法是目前求解凸优化问题最有效的算法之一，很多凸优化求解器均以该算法为基础进行二次开发，该算法对凸优化问题的求解效率明显优于其它轨迹规划方法，为轨迹规划问题的在线求解提供了技术保证。
参考文献
[1]Liu X,Shen Z,Lu P.Entry Trajectory optimization by second-ordercone programming.Journal of guidance control and dynamics,2015.DOI:10.2514/1.G001210.
[2]Phillips T H.A common aero vehicle model,description,andemployment guide.http://www.dtic.Mil/matris/sbir041/srch/af03 1a.doc,2003.
发明内容
本发明所要解决的技术问题是，针对现有技术不足，提供一种基于凸优化的再入轨迹三维剖面规划方法及系统，克服传统方法在求解三维剖面规划问题时计算效率低、精度低的问题，实现飞行器在大机动任务要求下的轨迹在线规划。
为解决上述技术问题，本发明所采用的技术方案是：一种基于凸优化的再入轨迹三维剖面规划方法，包括以下步骤：
1)给定飞行器的当前归一化机械能E0、终端归一化机械能Ef、当前运动状态x0＝[r0λ0φ0θ0σ0]和预期终端状态xf＝[rfλfφfθfσf]，其中，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；
2)对滑翔段参考轨迹进行初始化，具体方法为：将地心距从当前值r0到终端值rf的变化规律假设为线性的，将从当前值r0到终端值rf的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x(0)，并令参考轨迹x(k)＝x(0)，上标k表示当前迭代次数；
3)利用所述参考轨迹x(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x(k+1)(i)；
4)以x(k+1)(i)作为输出的最优轨迹，以第k+1次迭代的攻角指令和倾侧角指令(α(k+1)(i),υ(k+1)(i))作为输出的最优三维剖面。
借由上述过程，将三维剖面规划问题转化为一个凸优化问题，克服传统方法在求解三维剖面规划问题时计算效率低、精度低的问题。
上述步骤3)的具体实现过程包括：
3a)给定最小攻角αmin、最大攻角αmax、最小倾侧角υmin和最大倾侧角υmax；根据参考轨迹x(k)计算升阻比下界(u3)min和升阻比上界(u3)max；基于参考轨迹x(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数
Figure GDA0004071086540000031
rq、rn；
3b)将上述基于当前参考轨迹x(k)所得的连续变量(u3)min、(u3)max、(A,B,C)、
Figure GDA0004071086540000032
rq、rn进行等间距离散，得到每个离散点处的取值，表示为(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000033
rq(i)、rn(i)，其中，i表示离散点序号，变化范围从1到N；
3c)将离散化后的参数：(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000034
rq(i)、rn(i)以及预期终端状态xf作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x(k+1)(i),u(k+1)(i))；
3d)根据(x(k+1)(i),u(k+1)(i))计算当前攻角和倾侧角指令，其中，
Figure GDA0004071086540000035
Figure GDA0004071086540000036
分别表示在离散点i处的三个控制量分量；
3e)令x(k)＝x(k+1)(i)，重复上述步骤3a)～3d)，直到下列收敛准则成立：
Figure GDA0004071086540000041
其中，ε为一个常值矢量，维数为5×1。
上述求解过程计算量相对较小，求解效率高。
本发明所述升阻比下界(u3)min和升阻比上界(u3)max分别表示的计算公式表示如下：
Figure GDA0004071086540000042
Figure GDA0004071086540000043
其中，
Figure GDA0004071086540000044
表示对括号中函数在攻角α的取值范围[αmin,αmax]内求极小，
Figure GDA0004071086540000045
表示对括号中函数在攻角α的取值范围[αmin,αmax]内求极大；G表示升阻比计算函数；Mar是由参考轨迹x(k)根据大气模型得到的参考马赫数。
本发明矩阵A、B和矩阵C的计算公式表示如下：
Figure GDA0004071086540000046
Figure GDA0004071086540000047
C＝f(x(k),u(k),E)-Ax(k)+h(x(k))
其中，E表示当前离散点处的归一化机械能，所述当前离散点处的归一化机械能通过对归一化机械能变化域[E0,Ef]进行等间距离散获得；
Figure GDA0004071086540000051
Figure GDA0004071086540000052
上标k表示当前迭代次数，r(k)表示当前参考轨迹x(k)对应的地心距，φ(k)表示当前参考轨迹x(k)对应的纬度，V(k)表示当前参考轨迹x(k)对应的速度，θ(k)表示当前参考轨迹x(k)对应的速度倾角，σ(k)表示当前参考轨迹x(k)对应的航迹偏航角，D(k)表示当前参考轨迹x(k)对应的阻力加速度，
Figure GDA0004071086540000053
表示当前参考轨迹x(k)对应的阻力加速度关于地心距的偏导数，且
Figure GDA0004071086540000054
R0＝6371004m，u(k)表示当前参考轨迹x(k)对应的中间控制量，f(x(k),u(k),E)和h(x(k))表示如下：
Figure GDA0004071086540000055
Figure GDA0004071086540000061
ωe为地球自转角速度，ωe＝7.292115×10-5rad/s。
本发明中，热流、过载、动压三个约束线性化后的系数
Figure GDA0004071086540000062
rq、rn计算公式如下：
Figure GDA0004071086540000063
Figure GDA0004071086540000064
Figure GDA0004071086540000065
其中，
Figure GDA0004071086540000066
qmax和nmax分别为最大热流、最大动压、最大过载；
Figure GDA0004071086540000067
qk和nk分别为当前参考轨迹x(k)确定的热流、动压和过载；rk为当前参考轨迹x(k)确定的地心距；β为大气密度指数模型中的常数。
所述攻角指令α(k+1)(i)由
Figure GDA0004071086540000068
进行升阻比反插值得到，倾侧角指令由下式得到：
Figure GDA0004071086540000069
可见，攻角指令和倾侧角指令的计算简单，容易实现。
本发明中，为了进一步简化计算，常值矢量ε中的元素对应小于[5000.20.20.20.5]T中的元素，其中[5000.20.20.20.5]T第一项单位为米，后面四项单位为度。
相应地，本发明还提供了一种基于凸优化的再入轨迹三维剖面规划系统，包括：
初始化单元，对滑翔段参考轨迹进行初始化，具体执行如下操作：将地心距从当前值r0到终端值rf的变化规律假设为线性的，将从当前值r0到终端值rf的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x(0)，并令参考轨迹x(k)＝x(0)，上标k表示当前迭代次数；其中，所述归一化机械能包括飞行器的当前归一化机械能E0、终端归一化机械能Ef；飞行器的当前运动状态x0＝[r0λ0φ0θ0σ0]和预期终端状态xf＝[rfλfφfθfσf]，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；
计算单元，用于根据所述参考轨迹x(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x(k+1)(i)；
输出单元，用于输出最优轨迹x(k+1)(i)和最优三维剖面；所述最优三维剖面为第k+1次迭代的攻角指令和倾侧角指令(α(k+1)(i),υ(k+1)(i))。
计算单元包括：
第一计算模块，用于根据参考轨迹x(k)计算升阻比下界(u3)min和升阻比上界(u3)max；基于参考轨迹x(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数
Figure GDA0004071086540000071
rq、rn；
离散模块，用于将基于当前参考轨迹x(k)所得的连续变量(u3)min、(u3)max、(A,B,C)、
Figure GDA0004071086540000072
、rq、rn进行等间距离散，得到每个离散点处的取值，表示为(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000073
rq(i)、rn(i)，其中，i表示离散点序号，变化范围从1到N；
第二计算模块，用于将离散化后的参数：(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000081
rq(i)、rn(i)以及预期终端状态xf作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x(k+1)(i),u(k+1)(i))；
第三计算模块，用于根据(x(k+1)(i),u(k+1)(i))计算当前攻角和倾侧角指令，其中，
Figure GDA0004071086540000082
Figure GDA0004071086540000083
分别表示在离散点i处的三个控制量分量；
迭代模块，用于判断
Figure GDA0004071086540000084
是否成立，若成立，则输出x(k+1)(i)；否则，重复第一计算模块、离散模块、第二计算模块和第三计算模块的操作。
本发明相对现有技术的有益效果为：本发明将三维剖面规划问题转化为一个凸优化问题，利用这类问题的良好理论性质，实现对原问题高效精确的求解，既充分释放了飞行器的机动能力，又克服了传统方法求解三维剖面效率较低的难题，实现了飞行器在大机动任务要求下的轨迹在线规划。
附图说明
图1为本发明地面轨迹对比结果；
图2为本发明空间三维轨迹曲线；
图3为本发明速度倾角和航迹偏航角随时间的变化曲线；
图4为本发明攻角和倾侧角随时间的变化曲线；
图5为本发明热流、过载和动压随时间的变化曲线；
图6为本发明高度、经纬度、速度倾角和航迹偏航角随迭代次数的收敛曲线；
图7为本发明方法流程图；
图8为本发明系统结构框图；
图9为本发明计算单元结构框图。
具体实施方式
下面将结合附图和具体实例对本发明做进一步详细说明。
该实例中以美国CAV-H为仿真对象[2]，质量m＝907kg，参考面积Sm＝0.4839m2，仿真初始条件设置为：初始高度h0＝80km，速度V0＝7000m/s，速度倾角θ0＝-2°，初始经纬度均为0°，初始航迹偏航角为90°。预期终端条件设置为：高度hf＝25km，速度Vf＝2000m/s，速度倾角θf＝-4°，终端经度为100°，终端纬度为50°，终端航迹偏航角为30°。攻角取值范围为10到20度，倾侧角取值范围为-85到85度。最大热流
Figure GDA0004071086540000091
qmax＝90kPa，nmax＝4.5。计算机硬件条件为Inter Core i5-3470 CPU@3.2GHz，内存RAM为8G。
如图7所示，具体包括如下步骤：
1)步骤1：由公式E＝V2/2-μ/r计算飞行器当前归一化机械能
E0＝-0.5955和终端归一化机械能Ef＝-0.9641，其中，V表示飞行器速度，μ表示引力常数，μ＝3.986005×1014m3/s2。根据所给参数确定当前运动状态
x0＝[1.012600-0.03491.5708]和预期终端状态
xf＝[1.00391.74530.8727-0.06980.5236]；确定离散点数目N＝300；
2)步骤2：对滑翔段参考轨迹进行初始化，具体方法为：将地心距
r0＝1.0126到rf＝1.0039的变化范围均匀离散为300个区间，则初始化的地心距可以表示为
Figure GDA0004071086540000092
采用同样的方法得到初始化经度
Figure GDA0004071086540000093
初始化纬度
Figure GDA0004071086540000094
初始化速度倾角
Figure GDA0004071086540000095
初始化航迹偏航角
Figure GDA0004071086540000096
和离散后的归一化机械能{E0,E1,…,EN}，其中EN＝Ef；其中，上标0表示初始化参考轨迹的状态值，下标表示对应的离散点顺序。将这些初始化状态量联合在一起记为初始轨迹x(0)，并令
x(k)＝x(0)，上标k表示当前迭代次数；
步骤3：迭代执行以下步骤：
3a)根据参考轨迹x(k)计算升阻比下界(u3)min和升阻比上界(u3)max；
3b)基于参考轨迹x(k)，计算动力学方程的线性化矩阵(A,B,C)，矩阵A、B和C的计算公式表示如下：
Figure GDA0004071086540000101
Figure GDA0004071086540000102
C＝f(x(k),u(k),E)-Ax(k)+h(x(k))
其中，E表示当前离散点对应的归一化机械能，当前离散点处的归一化机械能通过对归一化机械能变化域[E0,Ef]进行等间距离散获得，对归一化机械能变化域[E0,Ef]进行等间距离散得到的能量依次取值为{E0,E1,…,EN}，
Figure GDA0004071086540000103
Figure GDA0004071086540000104
上标k表示当前迭代次数，r(k)表示当前参考轨迹x(k)对应的地心距，φ(k)表示当前参考轨迹x(k)对应的纬度，V(k)表示当前参考轨迹x(k)对应的速度，θ(k)表示当前参考轨迹x(k)对应的速度倾角，σ(k)表示当前参考轨迹x(k)对应的航迹偏航角，D(k)表示当前参考轨迹x(k)对应的阻力加速度，
Figure GDA0004071086540000117
表示当前参考轨迹x(k)对应的阻力加速度关于地心距的偏导数，且
Figure GDA0004071086540000118
R0＝6371004m，u(k)表示当前参考轨迹x(k)对应的中间控制量，f(x(k),u(k),E)和h(x(k))的表达式如下所示
Figure GDA0004071086540000111
Figure GDA0004071086540000112
ωe为地球自转角速度，ωe＝7.292115×10-5rad/s。
3c)基于参考轨迹x(k)，计算热流、过载、动压三个约束线性化后的系数
Figure GDA0004071086540000113
rq、rn；
3d)将上述基于当前参考轨迹x(k)所得的连续变量(u3)min、(u3)max、(A,B,C)、
Figure GDA0004071086540000119
rq、rn进行均匀离散，分别得到(u3)min、(u3)max、(A,B,C)、
Figure GDA0004071086540000114
rq、rn在每个离散点处的取值，表示为(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000115
rq(i)、rn(i)，其中，i表示离散点序号，变化范围从1到N。
3e)将离散化后的参数：(u3)min、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000116
rq(i)、rn(i)以及预期终端状态xf作为凸优化求解器的输入并进行求解，得到当前迭代的优化解(x(k+1)(i),u(k+1)(i))，其中，凸优化求解器可以通过开源方式获得或购买商用软件；
3f)根据(x(k+1)(i),u(k+1)(i))计算当前攻角和倾侧角指令，其中，
Figure GDA0004071086540000121
攻角指令α(k+1)(i)由u3 (k+1)(i)进行升阻比反插值得到，倾侧角指令
Figure GDA0004071086540000122
3g)令x(k)＝x(k+1)(i)，重复上述步骤3a)～3f)，直到下式收敛准则成立
Figure GDA0004071086540000123
其中，ε为一个常值矢量，本实施例中取值为ε＝[10001115]T。
步骤4：以x(k+1)(i)作为输出的最优轨迹，以(α(k+1)(i),υ(k+1)(i))作为输出的最优三维剖面。
根据上述步骤，本发明所提方法规划一条满足各项约束的可行轨迹计算耗时为12.6秒，而参考方法(即伪谱法)在同样参数设置和硬件条件下的计算耗时为235.00秒，由此证明了本发明中所提方法在同等计算条件下具有求解效率高的有益效果。
如图8，本发明的规划系统包括：
初始化单元，具体执行如下操作：将地心距从当前值r0到终端值rf的变化规律假设为线性的，将从当前值r0到终端值rf的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x(0)，并令参考轨迹x(k)＝x(0)，上标k表示当前迭代次数；其中，所述归一化机械能包括飞行器的当前归一化机械能E0、终端归一化机械能Ef；飞行器的当前运动状态x0＝[r0λ0φ0θ0σ0]和预期终端状态xf＝[rfλfφfθfσf]，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；
计算单元，用于根据所述参考轨迹x(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x(k+1)(i)；
输出单元，用于输出最优轨迹x(k+1)(i)和最优三维剖面；所述最优三维剖面为第k+1次迭代的攻角指令和倾侧角指令(α(k+1)(i),υ(k+1)(i))。
如图9，本发明的计算单元包括三个计算模块和一个离散模块，分别执行如下操作：
第一计算模块，用于根据参考轨迹x(k)计算升阻比下界(u3)min和升阻比上界(u3)max；基于参考轨迹x(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数
Figure GDA0004071086540000131
rq、rn；
离散模块，用于将基于当前参考轨迹x(k)所得的连续变量(u3)min、(u3)max、(A,B,C)、
Figure GDA0004071086540000132
rq、rn进行等间距离散，分别得到(u3)min、(u3)max、(A,B,C)、
Figure GDA0004071086540000138
rq、rn在每个离散点处的取值，表示为(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000133
rq(i)、rn(i)，其中，i表示离散点序号，变化范围从1到N；
第二计算模块，用于将离散化后的参数：(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure GDA0004071086540000134
rq(i)、rn(i)以及预期终端状态xf作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x(k+1)(i),u(k+1)(i))；
第三计算模块，用于根据(x(k+1)(i),u(k+1)(i))计算当前攻角和倾侧角指令，其中，
Figure GDA0004071086540000135
Figure GDA0004071086540000136
分别表示在离散点i处的三个控制量分量；
迭代模块，用于判断
Figure GDA0004071086540000137
是否成立，若成立，则输出x(k+1)(i)；否则，重复第一计算模块、离散模块、第二计算模块和第三计算模块的操作。
为了突出本发明中所提方法可以充分释放飞行器机动能力的有益效果，将所提方法得到的轨迹与基于二维剖面方法得到的轨迹进行对比。图1给出了地面轨迹的对比结果，可以看出，在给定再入点和预期目标点后，基于三维剖面的方法可以得到一条到达该目标的可行轨迹，但基于二维剖面的方法无法得到一条到达该目标的可行轨迹，因此基于三维剖面的规划方法可以使得飞行器的可达目标“更远”。由于任务参数设置是相同的，只有控制量选择的差异，说明本发明所提的三维剖面规划方法可以充分释放飞行器的机动能力，证明了本方法的在释放飞行器机动能力方面的有益效果。图2至图4给出了由本发明所提方法得到轨迹对应的空间飞行轨迹及其它参数与控制量的变化规律，均满足给定约束。图5说明了所得轨迹满足热流、动压和过载三项过程约束，最大值均未超过约束要求。图6说明了本发明所提方法是有效收敛的。

Claims (8)

1.一种基于凸优化的再入轨迹三维剖面规划方法，其特征在于，包括以下步骤：
1)给定飞行器的当前归一化机械能E0、终端归一化机械能Ef、当前运动状态x0＝[r0λ0φ0θ0σ0]和预期终端状态xf＝[rfλfφfθfσf]，其中，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；
2)对滑翔段参考轨迹进行初始化，具体方法为：将地心距从当前值r0到终端值rf的变化规律假设为线性的，将从当前值r0到终端值rf的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x(0)，并令参考轨迹x(k)＝x(0)，上标k表示当前迭代次数；
3)利用所述参考轨迹x(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x(k+1)(i)；
4)以x(k+1)(i)作为输出的最优轨迹，以第k+1次迭代的攻角指令和倾侧角指令(α(k+1)(i),υ(k+1)(i))作为输出的最优三维剖面；
步骤3)的具体实现过程包括：
3a)给定最小攻角αmin、最大攻角αmax、最小倾侧角υmin和最大倾侧角υmax；根据参考轨迹x(k)计算升阻比下界(u3)min和升阻比上界(u3)max；基于参考轨迹x(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数
Figure QLYQS_1
rq、rn；
3b)将上述基于当前参考轨迹x(k)所得的连续变量(u3)min、(u3)max、(A,B,C)、
Figure QLYQS_2
rq、rn进行等间距离散，得到(u3)min、(u3)max、(A,B,C)、
Figure QLYQS_3
rq、rn在每个离散点处的取值，分别表示为(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure QLYQS_4
rq(i)、rn(i)，其中，i表示离散点序号，变化范围从1到N；
3c)将离散化后的参数：(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure QLYQS_5
rq(i)、rn(i)以及预期终端状态xf作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x(k +1)(i),u(k+1)(i))；
3d)根据(x(k+1)(i),u(k+1)(i))计算当前攻角和倾侧角指令，其中，
Figure QLYQS_6
分别表示在离散点i处的三个控制量分量；
3e)令x(k)＝x(k+1)(i)，重复上述步骤3a)～3d)，直到下列收敛准则成立：
Figure QLYQS_7
其中，ε为一个常值矢量，维数为5×1。
2.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，所述升阻比下界(u3)min和升阻比上界(u3)max分别表示的计算公式表示如下：
Figure QLYQS_8
Figure QLYQS_9
其中，
Figure QLYQS_10
表示对括号中函数在攻角α的取值范围[αmin,αmax]内求极小，
Figure QLYQS_11
表示对括号中函数在攻角α的取值范围[αmin,αmax]内求极大；G表示升阻比计算函数；Mar是由参考轨迹x(k)根据大气模型得到的参考马赫数。
3.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，矩阵A、B、C的计算公式表示如下：
Figure QLYQS_12
Figure QLYQS_13
C＝f(x(k),u(k),E)-Ax(k)+h(x(k))；
其中，E表示当前离散点处的归一化机械能，所述当前离散点处的归一化机械能通过对归一化机械能变化域[E0,Ef]进行等间距离散获得；
Figure QLYQS_14
上标k表示当前迭代次数，r(k)表示当前参考轨迹x(k)对应的地心距，φ(k)表示当前参考轨迹x(k)对应的纬度，V(k)表示当前参考轨迹x(k)对应的速度，θ(k)表示当前参考轨迹x(k)对应的速度倾角，σ(k)表示当前参考轨迹x(k)对应的航迹偏航角，D(k)表示当前参考轨迹x(k)对应的阻力加速度，
Figure QLYQS_15
表示当前参考轨迹x(k)对应的阻力加速度关于地心距的偏导数，且
Figure QLYQS_16
R0＝6371004m，β为大气密度指数模型中的常数；u(k)表示当前参考轨迹x(k)对应的中间控制量，f(x(k),u(k),E)和h(x(k))，表示如下：
Figure QLYQS_17
Figure QLYQS_18
ωe为地球自转角速度，ωe＝7.292115×10-5rad/s。
4.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，热流、过载、动压三个约束线性化后的系数
Figure QLYQS_19
rq、rn计算公式如下：
Figure QLYQS_20
Figure QLYQS_21
Figure QLYQS_22
其中，
Figure QLYQS_23
qmax和nmax分别为最大热流、最大动压、最大过载；
Figure QLYQS_24
qk和nk分别为当前参考轨迹x(k)确定的热流、动压和过载；rk为当前参考轨迹x(k)确定的地心距；β为大气密度指数模型中的常数。
5.根据权利要求4所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，β＝1/7110。
6.根据权利要求1～5之一所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，所述攻角指令α(k+1)(i)由
Figure QLYQS_25
进行升阻比反插值得到，倾侧角指令由下式得到：
Figure QLYQS_26
7.根据权利要求1所述的基于凸优化的再入轨迹三维剖面规划方法，其特征在于，常值矢量ε中的元素对应小于[5000.20.20.20.5]T中的元素，其中[5000.20.20.20.5]T第一项单位为米，后面四项单位为度。
8.一种基于凸优化的再入轨迹三维剖面规划系统，其特征在于，包括：
初始化单元，对滑翔段参考轨迹进行初始化，具体执行如下操作：将地心距从当前值r0到终端值rf的变化规律假设为线性的，将从当前值r0到终端值rf的变化范围均分为N个区间，得到每个离散点上的初始地心距，采用同样方式得到经度、纬度、速度倾角、航迹偏航角、归一化机械能在每个离散点上的初始值，将这些初始值联合在一起记为初始轨迹x(0)，并令参考轨迹x(k)＝x(0)，上标k表示当前迭代次数；其中，所述归一化机械能包括飞行器的当前归一化机械能E0、终端归一化机械能Ef；飞行器的当前运动状态x0＝[r0λ0φ0θ0σ0]和预期终端状态xf＝[rfλfφfθfσf]，r为地心距，λ为经度，φ为纬度，θ为速度倾角，σ为航迹偏航角，下标为0表示当前状态的取值，下标为f表示终端状态的取值；给定离散点数目N，N为正整数；
计算单元，用于根据所述参考轨迹x(k)，结合凸优化求解器，得到第k+1次迭代的参考轨迹x(k+1)(i)；
输出单元，用于输出最优轨迹x(k+1)(i)和最优三维剖面；所述最优三维剖面为第k+1次迭代的攻角指令和倾侧角指令(α(k+1)(i),υ(k+1)(i))；
所述计算单元包括：
第一计算模块，用于根据参考轨迹x(k)计算升阻比下界(u3)min和升阻比上界(u3)max；基于参考轨迹x(k)，计算动力学方程的线性化矩阵(A,B,C)、热流、过载、动压三个约束线性化后的系数
Figure QLYQS_27
rq、rn；
离散模块，用于将基于当前参考轨迹x(k)所得的连续变量(u3)min、(u3)max、(A,B,C)、
Figure QLYQS_28
rq、rn进行等间距离散，分别得到(u3)min、(u3)max、(A,B,C)、
Figure QLYQS_29
rq、rn在每个离散点处的取值，表示为(u3)min(i)、(u3)max(i)、(A(i),B(i),C(i))、
Figure QLYQS_30
rq(i)、rn(i)，其中，i表示离散点序号，变化范围从1到N；
第二计算模块，用于将离散化后的参数：(u3)min、(u3)max(i)、(A(i),B(i),C(i))、rQ(i)、rq(i)、rn(i)以及预期终端状态xf作为凸优化求解器的输入并进行求解，从而得到当前迭代的优化解(x(k+1)(i),u(k+1)(i))；
第三计算模块，用于根据(x(k+1)(i),u(k+1)(i))计算当前攻角和倾侧角指令，其中，
Figure QLYQS_31
分别表示在离散点i处的三个控制量分量；
迭代模块，用于判断
Figure QLYQS_32
是否成立，若成立，则输出x(k+1)(i)；否则，重复第一计算模块、离散模块、第二计算模块和第三计算模块的操作。
CN202010419723.0A 2020-05-18 2020-05-18 一种基于凸优化的再入轨迹三维剖面规划方法及系统 Active CN111580535B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010419723.0A CN111580535B (zh) 2020-05-18 2020-05-18 一种基于凸优化的再入轨迹三维剖面规划方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010419723.0A CN111580535B (zh) 2020-05-18 2020-05-18 一种基于凸优化的再入轨迹三维剖面规划方法及系统

Publications (2)

Publication Number Publication Date
CN111580535A CN111580535A (zh) 2020-08-25
CN111580535B true CN111580535B (zh) 2023-06-06

Family

ID=72126924

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010419723.0A Active CN111580535B (zh) 2020-05-18 2020-05-18 一种基于凸优化的再入轨迹三维剖面规划方法及系统

Country Status (1)

Country Link
CN (1) CN111580535B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112380692B (zh) * 2020-11-12 2022-10-11 北京航天自动控制研究所 一种运载火箭的大气层内在线轨迹规划方法
CN112693631B (zh) * 2020-11-27 2022-07-29 中国人民解放军国防科技大学 飞行器在线序列凸优化中的初始轨迹生成方法及系统
CN112596537B (zh) * 2020-11-27 2022-03-29 中国人民解放军国防科技大学 用于在线轨迹规划的模型误差补偿方法、系统及存储介质
CN112965532B (zh) * 2021-03-22 2022-03-15 北京航空航天大学 一种基于路径优选的飞行器绕多禁飞区轨迹优化方法
CN113589838B (zh) * 2021-05-31 2023-08-01 南京航空航天大学 一种基于圆柱位置离散化的三维轨迹调度方法
CN113671826B (zh) * 2021-07-18 2023-10-13 北京理工大学 一种跨大气层飞行器气动辅助轨道可达能力快速评估方法
CN113960926B (zh) * 2021-10-18 2024-04-16 北京理工大学 一种气动捕获制导参数边界的自适应调节方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106020231A (zh) * 2016-05-30 2016-10-12 中国人民解放军国防科学技术大学 基于再入点参数的高超声速飞行器再入轨迹优化方法
CN107941087A (zh) * 2017-10-18 2018-04-20 北京航空航天大学 一种基于阻力剖面的高升阻比高超平稳滑翔再入制导方法
CN109740198A (zh) * 2018-12-14 2019-05-10 中国人民解放军国防科技大学 一种基于解析预测的滑翔飞行器三维再入制导方法
CN109976154A (zh) * 2019-03-04 2019-07-05 北京理工大学 一种基于混沌多项式和序列凸优化的飞行器轨迹优化方法
CN110750850A (zh) * 2019-07-12 2020-02-04 中国人民解放军国防科技大学 强约束复杂任务条件下的三维剖面优化设计方法、系统及介质

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8719212B2 (en) * 2011-05-09 2014-05-06 King Fahd University Of Petroleum And Minerals Parallel kinematic machine trajectory planning method
FR3053781B1 (fr) * 2016-07-07 2018-08-17 Thales Procede de calcul par un systeme de gestion de vol d'une trajectoire presentant des transitions ameliorees

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106020231A (zh) * 2016-05-30 2016-10-12 中国人民解放军国防科学技术大学 基于再入点参数的高超声速飞行器再入轨迹优化方法
CN107941087A (zh) * 2017-10-18 2018-04-20 北京航空航天大学 一种基于阻力剖面的高升阻比高超平稳滑翔再入制导方法
CN109740198A (zh) * 2018-12-14 2019-05-10 中国人民解放军国防科技大学 一种基于解析预测的滑翔飞行器三维再入制导方法
CN109976154A (zh) * 2019-03-04 2019-07-05 北京理工大学 一种基于混沌多项式和序列凸优化的飞行器轨迹优化方法
CN110750850A (zh) * 2019-07-12 2020-02-04 中国人民解放军国防科技大学 强约束复杂任务条件下的三维剖面优化设计方法、系统及介质

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张共济 ; 孙春贞 ; .基于拟平衡滑翔的再入轨迹快速规划方法.兵器装备工程学报.2016,(第01期),全文. *
郑总准 ; 吴浩 ; 王永骥 ; .基于序列二次规划算法的再入轨迹优化研究.航天控制.2009,(第06期),全文. *

Also Published As

Publication number Publication date
CN111580535A (zh) 2020-08-25

Similar Documents

Publication Publication Date Title
CN111580535B (zh) 一种基于凸优化的再入轨迹三维剖面规划方法及系统
Li et al. Stochastic gradient particle swarm optimization based entry trajectory rapid planning for hypersonic glide vehicles
Liu et al. Immersion and invariance-based output feedback control of air-breathing hypersonic vehicles
Seigler et al. Modeling and flight control of large-scale morphing aircraft
CN102073755B (zh) 近空间高超声速飞行器运动控制仿真方法
Tang et al. Nonlinear dynamic modeling and hybrid control design with dynamic compensator for a small-scale UAV quadrotor
Luo et al. Accurate flight path tracking control for powered parafoil aerial vehicle using ADRC-based wind feedforward compensation
Chen et al. Preliminary design of multirotor UAVs with tilted-rotors for improved disturbance rejection capability
CN109446582B (zh) 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
Lhachemi et al. A structured H∞-based optimization approach for integrated plant and self-scheduled flight control system design
Luo et al. On decoupling trajectory tracking control of unmanned powered parafoil using ADRC-based coupling analysis and dynamic feedforward compensation
CN105912025B (zh) 一种基于特征模型的高空飞艇水平位置控制方法
Zandavi et al. Multidisciplinary design of a guided flying vehicle using simplex nondominated sorting genetic algorithm II
Tripathi et al. Fast terminal sliding mode super twisting controller for position and altitude tracking of the quadrotor
Zheng et al. Constrained trajectory optimization with flexible final time for autonomous vehicles
Wu et al. Integrated optimization design using improved pigeon-inspired algorithm for a hypersonic vehicle model
Mohamad et al. Dynamic aerodynamic parameter estimation using a dynamic particle swarm optimization algorithm for rolling airframes
Zhao et al. Dynamic modelling of parafoil system based on aerodynamic coefficients identification
CN117850446A (zh) 一种基于改进序列凸规划的变翼飞行器在线轨迹优化方法
Li et al. Trajectory optimization for hypersonic boost‐glide missile considering aeroheating
Chen et al. Switching multi-model predictive control for hypersonic vehicle
Sattar et al. Disturbance rejection enhancement using predictive control for the fixed‐wing UAV with multiple ailerons
Li et al. The fuel optimal control problem of a hypersonic aircraft with periodic cruising mode
Su et al. Modeling and robust decoupling control for hypersonic scramjet vehicle
Dzul et al. Helicopter's nonlinear control via backstepping techniques

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