多飞行器协同再入轨迹规划方法、系统、电子设备及介质

1.本发明涉及飞行器轨迹规划技术领域，特别是涉及一种多飞行器协同再入轨迹规划方法、系统、电子设备及介质。


背景技术：

2.高超声速滑翔飞行器(hypersonic glidingvehicle，hgv)是一种临近空间飞行器，具有高速、大规模机动和远程滑翔能力，具备作战范围广、突防能力强等特点，有广阔的应用前景。为了对抗导弹防御系统，有必要提高hgv的突防能力。但随着防空武器系统的发展，目前对单一高超声速飞行器的防御和拦截能力越来越强，对其打击效果产生了较大的限制。在此种背景下，发展多枚高超声速飞行器协同作战技术变得更加必要。多飞行器协同作战打击可以大大增加防御系统的拦截难度，实现对目标的饱和攻击。而对于高超声速滑翔飞行器这样飞行环境不确定性大、参数剧烈时变、耦合严重的控制对象，对其多约束协同打击情况下的轨迹规划和跟踪控制更加具有挑战性。
3.高超声速飞行器再入段的轨迹规划研究是一个具有复杂多约束限制的非线性规划问题。在早期的研究中，规划轨迹的状态剖面是常用的方法，并在不断改进中取得了许多成果。值得注意的是，实时性在hgv协同轨迹优化中至关重要。凸优化方法由于具有良好的收敛性能和广泛的适用性，在轨迹优化领域得到了广泛的应用。具有复杂约束和非线性动力学的轨迹优化问题可以通过凸化技术转化为一系列凸优化问题。
4.关于多高超声速飞行器协同再入制导问题的研究大多集中在于实现再入段的时间协同。目前的主要方法是分析和设计再入段飞行时间与各状态量的对应关系，校正飞行轨迹实现期望的时间协同。但其不能直接控制末端时间，浪费了求解次数。而且其目标函数致力于保证终端约束，无法满足如最小热流过载等更好的性能指标。因此，针对多飞行器的协同再入轨迹优化问题需要进一步关注和研究。


技术实现要素：

5.为解决或至少缓解上述问题，本发明提出一种多飞行器协同再入轨迹规划方法、系统、电子设备及介质，以提高多飞行器协同再入轨迹规划效率和精度。
6.为实现上述目的，本发明提供了如下方案：
7.一方面，本发明提供一种多飞行器协同再入轨迹规划方法，包括：
8.建立飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型；
9.根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型；
10.对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型；
11.根据期望的初始和终端条件生成可用的初始参考轨迹；
12.根据所述处理后的轨迹规划问题模型和所述初始参考轨迹，采用序列凸优化方法
确定多飞行器协同再入轨迹。
13.可选地，所述建立高超声速滑翔飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型，具体包括：
14.建立高超声速滑翔飞行器的动力学模型；所述动力学模型采用以时间t为自变量的运动学方程描述飞行器再入段无动力滑翔过程；
15.通过重新定义归一化时间τ∈[0,1]为自变量，对所述动力学模型进行改进，生成改进后的飞行器动力学模型其中x＝[r,θ,φ,v,γ,ψ,σ,t]
t
为飞行状态向量；r为地心到飞行器的距离，θ和φ分别表示经度和纬度，v是飞行器相对于地球的飞行速度，γ和ψ分别是弹道倾角和速度航向角，σ是倾侧角；u
σ
为倾侧角的变化率；时间控制量f(x,u
σ
)和f
ω
(x,u
σ
)分别表示所述改进后的飞行器动力学模型中不带有和带有地球自转角速度的项。
[0016]
可选地，所述根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型，具体包括：
[0017]
根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型p0其中j表示优化函数，表示热流密度；分别是飞行器能承受的热流密度过载n和动压q的无量纲化最大值；g
1-g3表示约束的函数表达式，kq是热流密度常数，ρ是大气密度，l和d分别是升力和阻力加速度；σ
max
是倾侧角的限幅值；u
σmax
是倾侧角变化率的限幅值；是期望的总飞行时间；r0,θ0,φ0,v0,γ0,ψ0分别表示初始状态量；rf,θf,φf,vf,γf,ψf分别表示终端时刻状态量；v
fmin
，v
fmax
分别表示期望终端速度的最小值和最大值。
[0018]
可选地，所述对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型，具体包括：
[0019]
基于小扰动线性化方法对所述轨迹规划问题模型进行线性化处理和凸化处理，将所述轨迹规划问题模型p0转换为连续凸优化问题模型p1；
[0020]
对所述连续凸优化问题模型p1进行离散化处理，得到处理后的轨迹规划问题模型p2。
[0021]
可选地，所述根据期望的初始和终端条件生成可用的初始参考轨迹，具体包括：
[0022]
根据期望的初始和终端条件确定控制量初始序列其中和分别表示i离散点的倾侧角变化率控制量和时间控制量；
[0023]
将所述控制量初始序列代入所述轨迹规划问题模型p0，通过四阶龙格库塔数值积分生成可用的初始参考轨迹。
[0024]
另一方面，本发明还提供一种多飞行器协同再入轨迹规划系统，包括：
[0025]
飞行器建模及模型改进模块，用于建立飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型；
[0026]
轨迹规划问题模型建立模块，用于根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型；
[0027]
轨迹规划问题模型处理模块，用于对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型；
[0028]
初始参考轨迹生成模块，用于根据期望的初始和终端条件生成可用的初始参考轨迹；
[0029]
协同再入轨迹求解模块，用于根据所述处理后的轨迹规划问题模型和所述初始参考轨迹，采用序列凸优化方法确定多飞行器协同再入轨迹。
[0030]
另一方面，本发明还提供一种电子设备，包括存储器、处理器以及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述的多飞行器协同再入轨迹规划方法。
[0031]
另一方面，本发明还提供一种非暂态计算机可读存储介质，其上存储有计算机程序，所述计算机程序被执行时实现所述的多飞行器协同再入轨迹规划方法。
[0032]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0033]
本发明提供了一种多飞行器协同再入轨迹规划方法、系统、电子设备及介质，所述方法包括：建立飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型；根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型；对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型；根据期望的初始和终端条件生成可用的初始参考轨迹；根据所述处理后的轨迹规划问题模型和所述初始参考轨迹，采用序列凸优化方法确定多飞行器协同再入轨迹。本发明方法通过重新定义自变量对动力学模型进行改进，改进后的飞行器动力学模型能够直接控制末端时间，因此可以提高多飞行器协同再入轨迹规划效率；并且本发明方法确定的轨迹规划问题模型选取最小热流过载作为优化的目标函数，能够满足最小热流过载等更好的性能指标，从而提高了多飞行器协同再入轨迹规划精度。
附图说明
[0034]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获
得其他的附图。
[0035]
图1为本发明实施例提供的多飞行器协同再入轨迹规划方法的流程图；
[0036]
图2为本发明实施例提供的三个给定终端时间下的海拔、速度、弹道倾角和航向角随时间的变化曲线示意图；
[0037]
图3为本发明实施例提供的三个给定终端时间下的倾侧角指令曲线示意图；
[0038]
图4为本发明实施例提供的三个给定终端时间下的飞行轨迹的地面投影示意图；
[0039]
图5为本发明实施例提供的两飞行器海拔、速度、弹道倾角和航向角随时间的变化曲线示意图；
[0040]
图6为本发明实施例提供的两飞行器飞行轨迹的地面投影示意图；
[0041]
图7为本发明实施例提供的两飞行器目标函数的优化过程曲线示意图；
[0042]
图8为本发明实施例提供的两飞行器飞行轨迹仿真结果示意图；
[0043]
图9为本发明实施例提供的多飞行器协同再入轨迹规划系统的结构示意图。
具体实施方式
[0044]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0045]
本发明的目的是提供一种多飞行器协同再入轨迹规划方法、系统、电子设备及介质，能够解决多高超声速滑翔飞行器的协同再入轨迹规划问题，提高多飞行器协同再入轨迹规划效率和精度。
[0046]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0047]
图1为本发明实施例提供的多飞行器协同再入轨迹规划方法的流程图。参见图1，本发明一种多飞行器协同再入轨迹规划方法包括：
[0048]
步骤101：建立飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型。
[0049]
该步骤101首先建立了经典的高超声速滑翔飞行器(简称飞行器)再入动力学方程作为飞行器的动力学模型，然后为了实现末端时间自由且可控，改进了原方程，从而生成改进后的飞行器动力学模型。
[0050]
所述步骤101建立高超声速滑翔飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型，具体包括：
[0051]
步骤1.1：建立高超声速滑翔飞行器的动力学模型，所述动力学模型采用以时间t为自变量的运动学方程描述飞行器再入段无动力滑翔过程；假设地球为一个球体，考虑其自转，所述动力学模型的无量纲化方程如下：
[0052]
dr/dt＝vsinγ
[0053]
dθ/dt＝vcosγsinψ/(rcosφ)
[0054]
dφ/dt＝vcosγcosψ/r
[0055]
dv/dt＝-d-sinγ/r2+ω2rcosφ(sinγcosφ-cosγsinφcosψ)
[0056]
dγ/dt＝[lcosσ-cosγ/r2+v2cosγ/r+2ωvcosφsinψ+ω2rcosφ(cosγcosφ+sinγsinφcosψ)]/v
[0057]
dψ/dt＝[lsinσ/cosγ+v2cosγsinψtanφ/r-2ωv(tanγcosφcosψ-sinφ)+ω2rsinφcosφsinψ/cosγ]/v
[0058]
dσ/dt＝u
σ
[0059]
方程中从地心到飞行器的距离r被地球的平均半径r0＝6378140m无量纲化。状态量θ和φ分别是经度和纬度。v是相对于地球的飞行速度，被无量纲化，g0＝9.807m/s2是地球表面的重力加速度。γ、ψ分别是弹道倾角和速度航向角。常数ω是地球的自转角速度，其被无量纲化。自变量时间t被无量纲化。σ是倾侧角，为了保证飞行的稳定性，选取倾侧角的变化率u
σ
作为控制量，其被无量纲化。l和d分别是被g0无量纲化的升力和阻力加速度，它们分别具体表示为：
[0060][0061]
其中m为有量纲的飞行器质量，大气密度ρ＝ρ0exp(-h/hs)是飞行高度的函数，其中hs＝7200m，sa为飞行器参考面积，c
l
、cd分别为飞行器的升力系数和阻力系数，它们与飞行器的攻角和马赫数有关。本发明将攻角-速度剖面设计为常用的分段函数。
[0062]
步骤1.2：通过重新定义归一化时间τ∈[0,1]为自变量，对所述动力学模型进行改进，生成改进后的飞行器动力学模型其中x＝[r,θ,φ,v,γ,ψ,σ,t]
t
为飞行状态向量；r为地心到飞行器的距离，θ和φ分别表示经度和纬度，v是飞行器相对于地球的飞行速度，γ和ψ分别是弹道倾角和速度航向角，σ是倾侧角；u
σ
为倾侧角的变化率；时间控制量f(x,u
σ
)和f
ω
(x,u
σ
)分别表示所述改进后的飞行器动力学模型中不带有和带有地球自转角速度的项。
[0063]
本发明通过重新定义自变量对飞行器的动力学模型进行改进。假设飞行过程中的初末时间分别为t0和tf，重新定义归一化时间τ∈[0,1]为自变量，将原方程的时间区间映射到[0,1]上，并引入时间控制量u
t
＝t
f-t0，从而实际无量纲化时间t可以表示为：
[0064]
t＝t0+(t
f-t0)τ,τ∈[0,1]
[0065]
于是把原自变量无量纲化时间t作为新的状态量，得到其相对于归一化时间τ的微分方程为：
[0066][0067]
现在设飞行状态向量为x＝[r,θ,φ,v,γ,ψ,σ,t]
t
，把带有地球自转角速度的项分开，动力学微分方程组可以写成如下形式：
[0068][0069]
其中，f(x,u
σ
)代表不带有地球自转角速度的部分，f
ω
(x,u
σ
)代表带有地球自转角速度的部分，并且分别由f
1-f6和f
ω4-f
ω6
表示各等式方程中的项。则由上述式子得到以τ为自变量的运动方程作为改进后的飞行器动力学模型，可表示为：
[0070][0071]
其中u
t
为新增的时间控制量。
[0072]
步骤102：根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型。
[0073]
经典的飞行器再入路径约束可以表示为：
[0074][0075]
其中分别是飞行器能承受的热流密度过载n和动压q的无量纲化最大值。g
1-g3表示约束的函数表达式，kq是热流密度常数，ρ是大气密度，l和d分别是升力和阻力加速度。
[0076]
为了精确满足总再入飞行时间的协同约束，控制量u
t
需要设置为期望的飞行时间并全程保持。同时，飞行器的倾侧角及其变化率应当被限制在一定范围内以保证飞行器姿态的稳定性。
[0077][0078]
其中σ是倾侧角，σ
max
是倾侧角的限幅值；选取倾侧角的变化率u
σ
作为控制量，u
σmax
是倾侧角变化率的限幅值；u
t
是时间控制量，是期望的总飞行时间。
[0079]
优化过程对高超声速滑翔飞行器再入阶段的初末状态加以限制，主要针对其高度、经纬度、速度、弹道倾角和航向角，以实现多状态协同规划。
[0080]
[0081]
其中r0,θ0,φ0,v0,γ0,ψ0分别表示初始状态量，rf,θf,φf,vf,γf,ψf分别表示终端时刻状态量。v
fmin
，v
fmax
分别表示期望终端速度的最小值和最大值。
[0082]
本发明选取最小热流过载作为优化的目标函数，综合上述方程和不等式约束，带有飞行时间约束的协同再入轨迹规划可以描述为以下非凸最优控制问题：
[0083]
p0
[0084][0085]
subject to:
[0086]
式中，j表示优化函数，表示热流密度。状态向量为x＝[r,θ,φ,v,γ,ψ,σ,t]
t
，控制量是[u
t
,u
σ
]
t
，τ为归一化时间，f(x,u
σ
)和f
ω
(x,u
σ
)为系统的动力学方程，分别是飞行器能承受的热流密度过载n和动压q的无量纲化最大值。g
1-g3表示约束的函数表达式，kq是热流密度常数，ρ是大气密度，l和d分别是升力和阻力加速度。r是从地心到飞行器的距离，状态量θ和φ分别是经度和纬度。v是相对于地球的飞行速度，γ、ψ分别是弹道倾角和速度航向角。常数ω是地球的自转角速度，σ是倾侧角，σ
max
是倾侧角的限幅值；倾侧角的变化率u
σ
作为控制量，u
σmax
是倾侧角变化率的限幅值；u
t
是时间控制量，是期望的总飞行时间。r0,θ0,φ0,v0,γ0,ψ0分别表示初始状态量，rf,θf,φf,vf,γf,ψf分别表示终端时刻状态量。v
fmin
，v
fmax
分别表示期望终端速度的最小值和最大值。p0即为确定的轨迹规划问题模型。
[0087]
步骤103：对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型。
[0088]
步骤102确定的轨迹规划问题模型p0是一个高度非线性最优控制问题，本步骤103基于小扰动线性化方法，对非凸表达式进行近似处理；然后进行离散化，从而将其转换为可用序列凸优化方法求解的有限维二阶锥问题。
[0089]
所述步骤103对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型，具体包括：
[0090]
步骤3.1：基于小扰动线性化方法对所述轨迹规划问题模型进行线性化处理和凸
化处理，将所述轨迹规划问题模型p0转换为连续凸优化问题模型p1；
[0091]
动力学方程、过程约束和优化目标函数都是非凸形式。给定一个固定的历史状态可以根据小扰动线性化理论，通过在该点的一阶泰勒展开对它们进行线性化。
[0092]
在运动方程中，控制量与状态量是耦合的。因此需要借鉴多元函数的一阶泰勒展开进行线性化。用f(x,u)和f
ω
(x,u)来分别表示动力学方程中不带有和带有地球自转角速度的项，由于f
ω
(x,u)部分的值比较小，可以直接用固定参考状态近似：于是对动力学方程的线性化结果是：
[0093][0094]
其中
[0095]
非线性路径约束gj(r,v)都是与地心距r和速度v相关的函数，可以通过下式关于给定状态线性化：
[0096][0097]
其中
[0098]
另外，为了尽可能减小逼近误差，保证优化变量在给定参考点附近取值，引入信赖域约束，其中δ是信赖域的大小：
[0099][0100]
积分型目标函数也可以参考上式进行线性化。经过处理后，所述轨迹规划问题模型p0转换为连续凸优化问题模型p1。
[0101]
p1:
[0102]
subject to:
[0103]
式中，j表示优化函数，表示热流密度。状态向量为x＝[r,θ,φ,v,γ,ψ,σ,t]
t
，控
制量是u＝[u
t
,u
σ
]
t
，τ为归一化时间，用f(x,u)和f
ω
(x,u)来分别表示动力学方程中不带有和带有地球自转角速度的项，和带有地球自转角速度的项，g
1-g3表示约束的函数表达式，r是从地心到飞行器的距离，状态量θ和φ分别是经度和纬度。v是相对于地球的飞行速度，γ、ψ分别是弹道倾角和速度航向角。σ是倾侧角，σ
max
是倾侧角的限幅值；倾侧角的变化率u
σ
作为控制量，u
σmax
是倾侧角变化率的限幅值；u
t
是时间控制量，是期望的总飞行时间。r0,θ0,φ0,v0,γ0,ψ0分别表示初始状态量，rf,θf,φf,vf,γf,ψf分别表示终端时刻状态量。v
fmin
，v
fmax
分别表示期望终端速度的最小值和最大值。表示状态向量x的参考状态，δ是信赖域的大小。
[0104]
步骤3.2：对所述连续凸优化问题模型p1进行离散化处理，得到处理后的轨迹规划问题模型p2。
[0105]
为了求得连续最优控制问题p1的数值解，需要对其进行离散化处理。本发明采用梯形离散化的方法，首先把归一化时域分割为n个相等间隔，产生了n+1个自变量节点，每一段的步长是δτ＝1/n。所有离散的自变量节点表示为{τ0,τ1,τ2,...,τn}，其中τi＝τ
i-1
+δτ,i＝1,2,...,n。那么每个节点对应的状态量和控制量被离散为序列{x0,x1,x2,...,xn}和{u0,u1,u2,...,un}。基于上述序列，动力学方程可以按下式被数值积分：
[0106][0107]
其中，以k表示迭代次数，同时采用上一次的迭代结果作为本次求解的参考状态，即为即为即为δτ是离散步长，状态向量为x＝[r,θ,φ,v,γ,ψ,σ,t]
t
，控制量是u＝[u
t
,u
σ
]
t
，τ为归一化时间，用f(x,u)和f
ω
(x,u)来分别表示动力学方程中不带有和带有地球自转角速度的项，
[0108]
同时，其他的约束也分别被离散化，可表示为：
[0109][0110][0111]
[0112]
式中，带有k-1上标的变量都表示在上次迭代中的值，在这里作为参考值。g
1-g3表示约束的函数表达式，σ是倾侧角，σ
max
是倾侧角的限幅值；倾侧角的变化率u
σ
作为控制量，u
σmax
是倾侧角变化率的限幅值；u
t
是时间控制量，是期望的总飞行时间。r是从地心到飞行器的距离，v是相对于地球的飞行速度，表示状态向量x的参考状态，δ是信赖域的大小。
[0113]
飞行过程中随时间积累产生的热流负荷q也采用梯形法积分：
[0114][0115]
式中，δτ是离散步长，qi是每一个离散点的热流负荷，g1(r,v)表示热流密度的过程约束公式，r是从地心到飞行器的距离，v是相对于地球的飞行速度。
[0116]
于是p1被描述为离散的凸优化子问题，从而得到处理后的轨迹规划问题模型p2：
[0117]
p2
[0118][0119]
subject to:
[0120]
式中，j表示优化函数，qn是第n个离散点的热流负荷，i表示离散点的序号，状态向量为x＝[r,θ,φ,v,γ,ψ,σ,t]
t
，控制量是u＝[u
t
,u
σ
]
t
，δτ是离散步长，，δτ是离散步长，τ为归一化时间，用f(x,u)和f
ω
(x,u)来分别表示动力学方程中不带有
和带有地球自转角速度的项，g
1-g3表示约束的函数表达式，σ是倾侧角，σ
max
是倾侧角的限幅值；倾侧角的变化率u
σ
作为控制量，u
σmax
是倾侧角变化率的限幅值；u
t
是时间控制量，是期望的总飞行时间。r是从地心到飞行器的距离，v是相对于地球的飞行速度，状态量θ和φ分别是经度和纬度，γ、ψ分别是弹道倾角和速度航向角。r0,θ0,φ0,v0,γ0,ψ0分别表示初始状态量，rf,θf,φf,vf,γf,ψf分别表示终端时刻状态量。v
fmin
，v
fmax
分别表示期望终端速度的最小值和最大值。x
k-1
表示状态向量x的参考状态，δ是信赖域的大小。
[0121]
步骤104：根据期望的初始和终端条件生成可用的初始参考轨迹。
[0122]
所述步骤104根据期望的初始和终端条件生成可用的初始参考轨迹，具体包括：
[0123]
步骤4.1：根据期望的初始和终端条件确定控制量初始序列其中和分别表示i离散点的倾侧角变化率控制量和时间控制量。
[0124]
初始轨迹的规划方法以能量为自变量，定义了机动系数的概念来描述横纵向轨迹的耦合关系，然后在横侧向调整方向角偏差走廊并控制倾侧角反转以实现期望的侧向机动。但此时的状态集合中不包括飞行时间，因此根据飞行器动力学模型，按照下式得到该轨迹中各个节点所对应的时刻。注意这条数值积分的轨迹包括m+1个对应能量节点。
[0125][0126]
其中，i＝0,1,...,m表示各个离散点的序号，下标i和i+1表示在该时刻的状态量，{t0,t1,...,tm}表示时间序列，r是从地心到飞行器的距离，v是相对于地球的飞行速度，γ是弹道倾角。
[0127]
然后对得到的时间序列{t0,t1,...,tm}和对应倾侧角序列{σ0,σ1,...,σm}进行等时间间隔的线性插值。于是可以得出对应等时间间隔排列的倾侧角状态序列
[0128]
初始轨迹的飞行终端时间是tm，而求解过程的期望终端时间t
f*
可以在tm的附近任意取值，但不能有过大的调整，以免产生初始轨迹与其他约束之间的冲突，使得迭代求解失败。注意由于初始轨迹生成方法的限制，其终端状态可以与期望终端状态有一定的差值，但在允许范围内即可，理由同上。根据终端状态约束可以确定控制量初始序列接着，在考虑限幅约束的情况下，求得倾侧角变化率的初始序列：
[0129][0130]
其中，是对应等时间间隔排列的倾侧角状态序列元素，u
t,i
是i离散点的时间控
制量，表示i离散点的倾侧角变化率控制量，并作为初始迭代轨迹，δτ表示时间离散步长。
[0131]
步骤4.2：将所述控制量初始序列代入所述轨迹规划问题模型p0，通过四阶龙格库塔数值积分生成可用的初始参考轨迹。
[0132]
把上述步骤4.1得到的控制量初始序列代入动力学方程，通过四阶龙格库塔数值积分，重新得到初始轨迹的对应本文方程的状态量序列x
(0)
作为可用的初始参考轨迹。注意以上变量都需要提前无量纲化。
[0133]
步骤105：根据所述处理后的轨迹规划问题模型和所述初始参考轨迹，采用序列凸优化方法确定多飞行器协同再入轨迹。
[0134]
接下来介绍序列凸优化的求解过程。
[0135]
通过上述线性化处理，可以明确，给定参考轨迹的状态和控制量与真实轨迹之间的误差，对于近似过程的准确性有至关重要的影响。序列凸优化方法在轨迹规划中应用的具体过程就是，首先把某一条初始轨迹作为参考轨迹求解凸优化问题，再把这次得到的最优解作为下一次凸优化求解的参考轨迹，进行迭代。这一系列凸优化子问题将收敛至最终解，就可以作为原问题p1的解，过程可以如下描述：
[0136]
1)以k表示迭代次数。设k＝0，通过基于机动系数的三维轨迹快速生成方法，按初始状态约束，和末位置、高度、速度约束，生成初始轨迹。经过处理后得到第一条参考轨迹{x
(0)
,u
(0)
}。
[0137]
2)当k≥1时，把参考轨迹的状态变量和控制变量代入p1中，并求得该问题的解为{xk,uk}。
[0138]
3)检验下列序列收敛条件是否成立：
[0139]
sup|x
k-x
k-1
|≤ε,k≥1
[0140]
其中xk是第k次迭代的状态量结果，x
k-1
是第k-1次迭代的状态量结果。ε是序列收敛的容差阈值。若满足条件，则转至步骤4)；否则，令k＝k+1，重复步骤2)。
[0141]
4)得到最优轨迹{x*,u*}＝{xk,uk}。
[0142]
前述所有的x＝[r,θ,φ,v,γ,ψ,σ,t]
t
为状态向量，u＝[u
t
,u
σ
]
t
为控制量。
[0143]
下面采用仿真结果对本发明方法的有效性进行验证。
[0144]
利用matlab的建模工具箱yalmip建立了socp问题，并用mosek软件包进行问题的求解。仿真结果在采用intel-core i9-10900k 3.70ghz处理器的台式机上得到。
[0145]
仿真选用cav-h作为对象，飞行器参数设置：m＝907kg，sa＝0.484m2，飞行过程约束规定为n
max
＝3g，q
max
＝70kpa。飞行器攻角剖面设计为速度的分段线性函数：
[0146][0147]
其中选取参数α1＝20
°
，α2＝12
°
，v1＝5000m/s，v2＝2000m/s。飞行器的升力系数c
l
和阻力系数cd由气动数据插值得到。其他基本常数：ρ0＝1.752kg/m3，hs＝6700m，kq＝9.4369
×
10-5
。
[0148]
信赖域和迭代终止条件设置为：
[0149][0150][0151]
仿真初末状态约束见表1，h
0*
和h
f*
表示海拔高度，对倾侧角变化的限制为σ
max
＝80deg，u
σ,max
＝10deg/s。经过测试，初始轨迹的末端时间tm＝1502.55s，设定期望末端时间分别为1480s、1500s、1520s进行仿真。然后选取t
f*
＝1500s为例，展示其迭代过程，并与同条件下采用伪谱法求解的结果对比。
[0152]
表1初始状态和终端状态
[0153][0154]
图2为本发明实施例提供的三个给定终端时间下的海拔、速度、弹道倾角和航向角随时间的变化曲线示意图，其横坐标为时间(time)，纵坐标分别为三个给定终端时间(1480s、1500s、1520s)下的海拔(altitude)、速度(velocity)、弹道倾角(flight-path-angle)和航向角(heading-angle)。图3为本发明实施例提供的三个给定终端时间下的倾侧角指令曲线示意图，其横坐标为时间(time)，纵坐标为三种情况下的倾侧角(bankangle)指令曲线。图4为本发明实施例提供的三个给定终端时间下的飞行轨迹的地面投影示意图，其横坐标为经度(longitude)，纵坐标为维度(latitude)。图3和图4分别展示了三种情况下的倾侧角指令曲线和飞行轨迹的地面投影，可以看出飞行时间都被精确控制，同时各终端状态得到满足。并且当规定的飞行时间变长时，所规划的轨迹会增大纵向的机动，以延长飞行时间。每一次求解凸优化子问题的时间在0.5s至1.3s，这三次优化的迭代次数分别为6次、9次、5次，优化总时间分别为4.926s、7.865s、3.684s。
[0155]
仿真继续模拟了两架cav飞行器协同再入的情况。规定t
f*
＝1500s，cav1的初末条件参照表1，cav2相比来说只更改初始航向角ψ
0*
＝59deg。图5为本发明实施例提供的两飞行器海拔、速度、弹道倾角和航向角随时间的变化曲线示意图，图6为本发明实施例提供的两
飞行器飞行轨迹的地面投影示意图。图5和图6描述了协同再入飞行过程中，两飞行器的各状态量随时间的变化曲线以及轨迹地面投影，可以看出在两飞行器的飞行末端，总时间和各状态均保持一致，实现了协同再入规划，这为高超声速飞行器再入段编队的实现提供了方法基础。
[0156]
图7为本发明实施例提供的两飞行器目标函数的优化过程曲线示意图，其横坐标为步数(step number)，纵坐标为目标函数的值(value ofobjective function)。图8为本发明实施例提供的两飞行器飞行轨迹仿真结果示意图，其中validation表示验证曲线。图7展示了两飞行器对应的目标函数的优化过程。另外，为了验证本凸化方法的合理性，把优化的控制量结果，使用四阶龙格库塔方法代入原动力学方程数值积分，仿真步长为0.1s，结果展示在图8中。可以看出积分轨迹与优化轨迹非常接近，说明近似处理比较符合实际。
[0157]
综上，本发明提出的一种基于序列凸优化的多飞行器协同再入轨迹规划方法，可以在多约束条件下，规划多飞行器再入飞行时间相同的弹道轨迹，从而实现多状态协同的再入轨迹规划，提高了多飞行器协同再入轨迹规划效率和精度，并通过数值仿真验证了所提出的方法的有效性。
[0158]
此外，对应于上述提供的多飞行器协同再入轨迹规划方法，本发明还提供了一种多飞行器协同再入轨迹规划系统，参见图9，所述多飞行器协同再入轨迹规划系统包括：
[0159]
飞行器建模及模型改进模块901，用于建立飞行器的动力学模型，并通过重新定义自变量对所述动力学模型进行改进，生成改进后的飞行器动力学模型；
[0160]
轨迹规划问题模型建立模块902，用于根据所述改进后的飞行器动力学模型建立终端时间约束，附加飞行器的其他各项约束确定轨迹规划问题模型；
[0161]
轨迹规划问题模型处理模块903，用于对所述轨迹规划问题模型进行处理，得到处理后的轨迹规划问题模型；
[0162]
初始参考轨迹生成模块904，用于根据期望的初始和终端条件生成可用的初始参考轨迹；
[0163]
协同再入轨迹求解模块905，用于根据所述处理后的轨迹规划问题模型和所述初始参考轨迹，采用序列凸优化方法确定多飞行器协同再入轨迹。
[0164]
本发明实施例提供的多飞行器协同再入轨迹规划系统与上述实施例所述的多飞行器协同再入轨迹规划方法，其工作原理和有益效果类似，故此处不再详述，具体内容可参见上述方法实施例的介绍。
[0165]
本发明还提供一种电子设备，该电子设备可以包括：处理器、通信接口、存储器和通信总线。其中，处理器、通信接口、存储器通过总线完成相互间的通信。处理器可以调用存储器中的计算机程序，以执行所述的多飞行器协同再入轨迹规划方法。
[0166]
此外，上述的存储器中的计算机程序通过软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机、服务器或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器、随机存取存储器、磁碟或者光盘等各种可以存储程序代码的介质。
[0167]
进一步地，本发明还提供一种非暂态计算机可读存储介质，其上存储有计算机程序，所述计算机程序被执行时可以实现所述的多飞行器协同再入轨迹规划方法。
[0168]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
[0169]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。