基于可变信赖域的飞行器轨迹快速优化

2023-07-29王永海赵雅心赵亮博

计算机仿真 2023年6期

李瑾,王永海,赵雅心,赵亮博

(1. 北京航天长征飞行器研究所,北京 100076;2. 中国运载火箭技术研究院,北京 100076)

1 引言

升力式飞行器返回大气层过程中,可充分利用其气动性能和机动能力,完成大范围、长距离、高马赫数下的飞行,可对敌防御系统进行有效突破。在飞行过程中,其会面临复杂的外部气动环境变化:非线性的驻点热流条件、飞行器过载以及动压限制等,同时也需考虑飞行路径约束(如禁飞区)。因此针对飞行器再入段的轨迹优化近年来成为研究的热点,并已经形成诸多有效的方法,如直接法、间接法、精确解法和现代启发式算法等。由于伪谱法[1]和凸优化方法作为直接法中两个比较高效的方法,其得到广泛使用。目前主要针对提高优化效率、缩短规划时间对上述方法展开大量研究和改进,以实现可应用于未来弹上实时轨迹优化的解决方案[2]。

针对伪谱法,其主要包括Gauss、Legendre和Redau三类伪谱法,并且目前大部分解决方案是对上述三种方法的适应性改进[3]。但伪谱法的本质是一种离散和转化方法,其最终计算性能指标需与求解方案相结合,如大多数情况下伪谱法与序列二次规划(SQP)一同使用以完成优化计算,考虑到SQP问题中存在的对初值灵敏度、收敛速度等方面的不足,其仍需要进一步对算法进行改进[3]。而现代启发式算法如遗传算法、粒子群算法等虽然在迭代速度方面有着独特的优势,但对再入过程的复杂约束条件会出现精度不高,计算后期容易陷入局部极小值等问题。

而凸优化作为近年来逐渐兴起的优化算法,其对于初值不敏感,并且当问题满足凸优化的模型时,具有局部最优解为全局最优解的优势,可有效避免陷入局部极小值的问题。并且随着计算方法的提升,大规模凸优化问题能够在短时间内求得最优解[4],该方法也逐渐展开应用。其率先在宇航领域展开应用,美国JPL实验室的Acikmese[5]等人率先将凸优化方法应用于火星着陆器动力下降段制导任务中,并取得很好的效果;刘新福和Lu等采用凸优化方法对飞行器不同场景下轨迹优化方法进行了详细的分析。在可重复使用火箭回收[6]、空间轨道转移[7]、地外天体着陆以及可重复使用飞行器轨迹优化场景中多次应用该方法[8-10],通过对非凸问题进行转化,得到凸优化计算的模型,高效快速完成求解计算。周祥[11]等提出指令反解方法嵌套在凸优化模型中,使得计算精度得到提高,完成再入三维轨迹剖面的快速规划。王振博[12]采用序列凸优化方法,通过增加倾侧角变化率作为新的控制输入,对高非线性条件下的行星再入最优控制问题进行了快速求解。杨奔[13]针对凸优化方法解决再入飞行器轨迹优化过程中出现的“伪不可行”问题采用b样条曲线离散控制量并基于回溯直线搜索方法,提高数值仿真的稳定性和光滑性。王劲博[14]、马林[15]等结合序列凸优化和无损凸化方法展开了系列研究并给出相应理论证明,应用于可重复使用火箭轨迹优化和动力软着陆问题中。

本文以升力式飞行器三自由度归一化再入运动模型为基础,以能量e为自变量,减少运动方程数目,同时避免终端时间自由问题。采用基于可变信赖域的序列凸化和松弛方法对约束条件进行凸化处理。通过自动调节信赖域范围和权值函数大小,以减少迭代中的振荡问题。最后将凸问题离散化处理后,进行数值仿真计算。通过分析优化结果,验证算法的可行性和计算效率。

2 轨迹优化算法描述

2.1 再入运动模型

升力式飞行器再入段没有发动机推力作用,仅在重力和空气动力作用力下完成飞行。为简化优化数目,选取能量为自变量,考虑地球自转角速度的影响,建立再入段的三自由度无量纲运动模型如下

(1)

(2)

上式中需要注意的,速度V无量纲,密度ρ为有量纲的气动密度,是关于r的函数ρ=ρ0exp(-(r-1)R0/hs),hs表示高度系数,参考面积Sref和质量m均为有量纲的物理量。气动升力系数CL和阻力系数CD本文假设只与飞行马赫数以及攻角有关。

由地球自转引起的哥式加速度项C1、C2以及牵连加速度项C3、C4如下所示。

(3)

2.2 约束条件

终端等式约束为高度和经纬度,飞行航迹角给定在一个范围内。同理航向角由于滑翔段后还将有其它飞行阶段(如能量管理阶段或者是末制导阶段),其也需限制在一定范围内。不失一般性,这里假设终端约束条件如下所示

(4)

其飞行过程中,需满足热流、动压和过载要求,其无量纲形式如下所示

(5)

除上述约束外,控制量倾侧角约束如下

(6)

2.3 优化问题描述

综合考虑上述升力式飞行器运动模型和相关约束条件,考虑时间最短优化指标,其轨迹优化问题(P1)可以表示为

subto(1) (4) (5) (6)

(7)

3 凸化与离散

由于升力式飞行器轨迹优化问题是一高度非线性的最优控制问题,其中的性能指标、状态方程以及约束条件大部分难以直接利用凸优化方法进行求解。因此需对其进行凸化处理,本文主要采用序列凸优化方法(SCP)和松弛方法进行凸化,通过求解一系列凸优化的子问题,从而达到逐步逼近原有问题的解。

3.1 运动方程和控制量凸化

根据上文选取确定的状态量和控制量,将运动方程组进行线性化处理,改写为以下形式:

(8)

其中

b(xk)=f0(xk)-A(xk)xk+fΩ(xk)

A中各项系数可通过推导计算求得。式中的xk为参考轨迹的参考点的状态量。并注意矩阵A中对矢径r求导的相关项中,阻力加速度D对r的导数,Dr=-DR0/hs。由于线性化展开发生在xk处,需引入信赖域约束以确保其线性化的有效性,即

|x-xk|≤δ

(9)

针对控制量u的约束条件(6)中,第一个不等式约束满足凸化要求,但第二条约束条件是非凸的二次等式约束,拟对其采用松弛方法将其转化为不等式约束

(10)

3.2 约束和指标凸化

针对过程约束,将其转化为关于矢径r的对数形式的约束,并进行整合,如下所示

(11)

将上述约束进行整理改写为

(12)

针对飞行过程中,还将有飞行路径约束(如禁飞区)和其它约束,这部分将在后文中进行添加以验证算法能力。

针对性能指标凸化问题,通过辅助变量将待积函数转化为线性形式,考虑优化初期经纬度等式约束较为严格,难以得到全局最优解,以罚函数形式将其加入性能指标当中。同时为确保控制量约束松弛带来的最优解与原问题保持一致,增加正则项,因此性能指标如下所示

(13)

其中,Δθ,Δφ为待优化变量,cθ,cφ为惩罚项因子,迭代过程进行适应性调整。

3.3 最优问题离散化

为了进行后续数值求解,需要将上述最优问题进行离散化处理,将能量区间[e0,ef]进行等距划分,得到离散的自变量E={e0,e1,e2,…,eN},其中N为离散参数。对应的离散状态量和离散控制量分别如下所示

X={x0,x1,x2,…,xN},U={u0,u1,u2,…,uN}

(14)

其中,xi=[ri,θi,φi,γi,ψi]T,ui=[u1,u2]T。根据上述线性化的方程以及梯形积分公式有

i=1,2,…,N

(15)

合并之后得

(16)

Hi-1,Ki,Gi相关参数由式(15)易知。将上式整理成矩阵形式Mz=Y,其选取优化变量z为z=[x0,x1,x2,…,xN,u0,u1,u2,…,uN,Δθ,Δφ],M矩阵和Y向量如下所示

同样,其性能指标也需要离散化处理,主要为积分项进行离散化处理。并将其转化为J=Lz的形式。如下所示

(17)

通过上述离散化处理,凸问题P1的离散形式如P2所示,P2问题已经转化为标准的具有线性性能指标函数,可以进行凸优化求解。

(18)

3.4 改进信赖域算法

针对序列凸优化算法,其旨在通过序列线性化方式将非凸问题转化,但该过程中需要满足线性化要求的信赖域指标。固定的信赖域会在前期优化中会出现较大的振荡现象。为优化该过程,设计了基于可变信赖域和权值函数的方法如下所示,使迭代过程中更加接近优化值。

δ=-Kδ0

K=-ae-bk

(19)

其中,k为迭代次数,a,b分别为设置参数,可根据需求选取合适数值,本文中设置a=1.19,b=-0.36。同理性能指标中的惩罚函数因子cθ和cφ也按照上述方法进行设置。

程序流程图如下所示。

图1 优化流程图

4 数值仿真

4.1 问题描述

本文中,拟采用美国某一可重复使用飞行器(RLV)的飞行模型[16]。其质量为m=104305kg,参考面积为S=391.22m2,其气动模型拟合函数如下所示

(20)

图2 气动力系数与速度之间的关系

攻角方案采用给定的分段函数式剖面,其具体形式如下所示

(21)

4.2 计算结果分析

本文的计算仿真采用MATLAB平台下的Mosek求解器进行求解凸优化问题,其可大大提高计算的效率。运行计算机搭载IntelCorei5-6500 3.20GHz处理器。每次优化CPU的平均时间为0.174s。根据飞行结果,其飞行时间为1374.05秒,飞行距离为5824.06千米,下面对部分飞行参数结果进行展示分析。

图3中为飞行高度的变化情况,结合下图8的飞行轨迹图可知其优化曲线满足约束要求,并且在逐步迭代过程中曲线逐渐趋于最优解。分析下述飞行终点与设置值大小如下表2所示,可知终端高度和航向角满足等式约束,而经纬度为满足优化初期结果正确性,于性能指标中增加了惩罚因子,使得终端纬度与目标值有一定误差。速度倾角为给定的区间范围,其数值满足终端条件。

表2 规划终端结果表

图3 迭代飞行轨迹变化情况

由飞行控制量倾侧角的变化情况可知,其数值变化主要集中于10度上下的区间范围,由于飞行过载的限制,整个飞行轨迹为一个大半径的圆弧,其形式与倾侧角变化过程正好对应。除在接近终端约束过程中需要调节速度倾角和终端位置坐标,因此倾侧角会出现一个幅度较大的变化,其满足飞行要求。如图4。

图4 倾侧角变化情况

由图5,6所示,速度倾角在滑翔飞行过程中始终保持在0度附近,满足滑翔飞行要求;而航向角由于初始值和终端值的设定,其需要完成一定范围的角度变化。并且同飞行高度变化曲线和控制量变化结合来看,其速度倾角在滑翔飞行段由于高度变化存在小幅度的跳动现象,并且在飞行终端,为满足终端硬性约束条件,倾侧角和速度倾角同步发生一定的振荡。飞行高度与速度的变化曲线,如图8所示,其过程满足速度变化要求。