热率约束下高超声速飞行器Gauss时间网格参数化轨迹规划

刘平 ,刘航 ,仇国庆 ,刘兴高

(1.重庆邮电大学自动化学院,重庆 400065;2.浙江大学控制科学与工程学院,浙江杭州 310027)

1 引言

高超声速飞行器(hypersonic vehicle,HV)具有极高的飞行速度,在军事和民用领域具有广阔应用前景[1].近年来,HV再入轨迹优化研究受到国内外研究者广泛关注[2–3].考虑到HV在再入段一般无动力源,完全依靠气动力实现远距离快速精准打击,所以再入段轨迹规划结果对飞行器能否准确安全完成飞行任务至关重要[4–6].然而,HV在再入段做滑翔飞行时需要进入大气层,飞行环境复杂且速度高(马赫数>5),此时在气动加热产生的热量和内部应力等共同作用下机体表面温度会长时间在1000◦C以上[7],这给HV热防护系统(thermal protection system,TPS)的结构设计方面带来了巨大的挑战,所以热防护问题也成为HV结构设计的主要研究内容之一[8].但是,HV热防护系统也不能无限增加防护指标,如果HV 热防护系统的结构设计得过于复杂,虽然能够有效增强飞行器防热,但会增加飞行器重量和体积,对隔热材料要求也会更为苛刻,进而增加飞行器设计和维护成本.

文献[9]指出,考虑到HV做无动力滑翔飞行,通过改变攻角的大小来改变气动力产生的升降力,因而攻角的改变也能影响到热防护系统的作用.因此,从飞行器控制角度出发,通过对攻角参数的优化,选取最优攻角参数给TPS设计提供了新的思路.近些年,国内外的诸多研究者[10–11]也基于此思路将HV飞行过程中热载(即在整个飞行时间内对热率进行积分)作为优化性能指标,基于优化理论减少了飞行时间内由气动加热产生的热量和热防护系统的重量.但是由于HV需要满足的热率路径约束一般按最大热率设计,传统优化方法还是难以降低热保护系统的设计难度和材料成本;另一方面,HV的热防护系统主要作用于再入段,而HV再入段飞行要求在大气层外缘的高超声速飞行器安全可靠的到达稠密大气层中的某一个位置,并且满足路径约束、终端约束和控制变量约束等条件[12],再入段轨迹规划作为实现HV再入段飞行任务的重要手段,对热防护系统选择和HV飞行任务的成功实施均具有重要影响.

本质上来讲,HV再入轨迹规划是含有复杂约束条件的非线性最优控制问题,一般通过优化控制向量,使飞行器从初始状态开始,在满足一系列复杂约束条件后达到预定终端状态并且实现某个性能指标最优[13].因而,最优控制算法是HV再入轨迹规划的有效方法[14].当前,最优控制问题的求解方法包括数值法和启发式算法等[15–16].数值法主要分为间接法和直接法[17],其中间接法是基于庞特里亚金极小值原理将最优控制问题转化为两点边值问题进行理论求解,进而得到比较精确的解[4];直接法则是将最优控制问题(optimal control problem,OCP)离散化为非线性规划(nonlinear programming,NLP)问题再进行数值优化方法求解,其主要分为伪谱法(Gauss pseudospectral method,GPM)和控制向量参数化方法(control variable parameterization,CVP),直接法的优点在于能够有效降低优化问题的复杂程度,具有较大的收敛域,更易数值求解实现.其中:伪谱法主要思想是通过采用多种正交多项式对状态变量和控制向量离散逼近,然后将OCP近似转化为NLP问题进行求解.近年来,该类法也被广泛应用于HV的轨迹优化问题求解,比如文献[10]基于传统Gauss伪谱方法开展了HV再入轨迹快速优化研究;文献[18]提出了分段区间配点网格重构的自适应伪谱法,提升了HV优化问题求解效率;文献[19]则在传统Gauss伪谱法基础上引入了分段优化思路适应不同飞行任务.伪谱法具有收敛域范围广,收敛速度快的优点,但文献[14]也指出其存在生成的NLP问题规模较大,随着求解问题的复杂化增加求解难度而不易找到满意的优化解的缺陷.

相比于伪谱法,CVP方法则只离散控制向量,而微分方程保留在转化后的NLP问题中,所以产生的NLP问题规模较小,同时得到的状态轨迹也更加精确,但求解NLP问题需要多次求解微分方程组,控制变量离散时间节点选取对优化结果有影响.当前,CVP法也在HV轨迹优化问题求解上得到了广泛应用,其离散网格选择方式与路径约束处理成为了研究热点方向[20,22].由于HV具有复杂的约束条件,需要在整个再入过程满足热率、动压和过载等路径约束[14],相比其他优化问题,处理起来具有一定难度.为此,Teo等[23]提出了增强时间转换策略(enhanced control variable parameterization,ECVP),直接将时间网格的长度作为待优化参数,使原来时间网格的节点映射到新的时间尺度上;Zhang等[24]使用了经典的CVP方法和ECVP方法,用于HV轨迹优化求解;Liu等[25]基于CVP法将光滑化惩罚函数和终端约束处理用于轨迹优化的求解.然而,为了获得更加精确的优化结果,还需要对时间网格进行优化,使其在较少时间网格数量下划得更加接近实际控制曲线[26].因此,研究者们也提出了各自的时间网格优化策略[25,27–28].

HV再入轨迹优化的航程距离能够很好反映飞行器在再入过程的机动性能指标,为了进一步提升HV轨迹优化问题的求解效果,同时得到精确的状态轨迹约束满足曲线,本文以CVP方法为基础,借鉴伪谱法离散配点方式,针对HV轨迹优化问题的约束处理和控制变量离散时间网格选择提出一种非均匀Gauss离散时间网格参数化CVP方法.首先,将求解的轨迹优化问题数学模型化,采用光滑化惩罚函数法将路径约束转化到目标函数中;其次,求解Legendre多项式的零点,通过时间转换公式得到控制时域内的高斯了离散时间节点;然后采用终端约束惩罚方法将优化问题转化为无约束问题;最后,对热率约束下HV轨迹优化的最大航程问题进行数值仿真,并与先前求解方法和文献方法对比以验证提出方法的正确和有效性.在此基础上,本文进一步结合HV再入飞行任务要求,开展不同热率约束下HV最优轨迹优化算法研究,通过飞行器参数最优选择,在满足HV飞行约束的情况下分析降低热载约束限值对飞行目标性能的影响,以此为HV的TPS设计选择提供理论依据.

本文的结构安排包括5个部分:第1节为引言部分;第2节是HV轨迹优化问题数学建模;第3节描述了轨迹优化问题路径约束处理方法和提出的基于Gauss时间网格离散CVP方法,并描述了提出方法的求解结构和算法步骤;第4节进行了CAV–HV实例数值仿真求解和结果分析;第5节对本文工作进行了总结.

2 HV轨迹优化问题数学模型

2.1 HV运动方程模型

结合文献[4,25]研究,将HV对象在再入飞行过程中假设为一个质点,同时考虑地球非旋转.因而,得到描述HV质心运动的微分方程模型为

式中:h(t)为飞行器的海拔高度;V(t)为飞行速度;γ(t)为航迹倾角;r(t)为航程;Re为地球半径;m为飞行器的质量;g为重力常数;L和D为飞行器的升力和阻力,分别表示为

其中:Sa为飞行器的受力参考面积;ρ(h)表示大气密度随高度的变化关系,表达式为

CL和CD是升力L和阻力D的系数,分别是关于攻角α的函数

攻角α为控制变量,其控制边界为

同时,在HV再入过程中,飞行环境恶劣且受到自身结构限制,因此HV在整个飞行过程中需要严格满足热率、动压和过载等路径约束,分别表示为

此外,对于HV的再入轨迹规划,通常需要从预先设定的初始时刻状态到达指定的终端时刻状态,因而HV再入轨迹规划过程还需要满足初始和终端约束条件表示为

2.2 轨迹优化问题模型

结合模型分析,定义动态系统状态向量为x(t):=[h(t)V(t)γ(t)r(t)]T,控制向量为u(t):=α,则HV质心运动微分方程可以简化表示成如下形式:

式中:f:R×Rn×Rm是关于时间t∈R、状态向量x∈Rn和控制向量u∈Rm的非线性函数,n和m分别表示状态向量和控制向量的维度.结合式(11)–(12),系统的初始和终端状态值表示为

同时,结合通用路径约束形式将路径约束(8)–(10)写成一般形式为

控制边界条件表示为

在HV轨迹优化过程中,通常存在各种各样的优化目标,如飞行时间最短、能量消耗最少等,这些性能指标通常用目标函数表示.本文采用Bolza形式的目标函数进行表示,其包含终值项(Mayer形式)和积分项(Lagrange形式),具体如下式所示:

其中:Φ0(·)为终值项,L0(·)为积分项,Φ0和L0为连续可导函数.最终,HV再入轨迹优化问题可以描述为如问题(P1)所示的最优控制问题.

问题1(P1) 给定连续动态方程(13)和状态初始条件(14),通过优化控制向量u(t)使得其在整个控制时域t内满足路径约束条件(16)和边界约束(17),在到达指定终端约束(15)的同时,目标函数(18)为最优.

3 非均匀Gauss网格参数化轨迹优化方法

3.1 路径约束处理

由于飞行器在飞行任务中面临异常恶劣的飞行环境,故在整个飞行时间内要严格满足热率、动压和过载路径约束,这给轨迹优化的求解带来了难度[25].本文在先前研究发现[29],采用惩罚函数法可以很好地处理最优控制问题路径约束.因此,本文采用l1光滑化近似函数将HV轨迹优化问题的约束近似化处理,然后采用惩罚函数将路径约束增广进目标函数中进行求解.

首先,运用max函数对式(16)中每个不等式路径约束Gr(x(t),u(t),t),r∈∆={1,2,3}进行表示

由于式(19)为非光滑化函数,因而运用光滑化近似函数进行处理为

其中:Pr是关于Gr[x(t),u(t),t]和ε的函数,ε(0<ε<1)是光滑化因子.分析后可知,当ε取值较小时,式(19)可以光滑近似化为

式(21)可由以下定理1保证:

定理1设光滑化因子ε>0,对光滑近似函数pr{Gr[x(t),u(t),t],ε},有

限于篇幅,定理1的证明请参阅文献[29].

进一步,为了处理光滑化后的不等式约束,引入附加状态变量xn+1,令

由xn+1(t0)=0,可得

将式(24)通过引入惩罚因子ρ增广进式(18),问题(P1)的目标函数由式(18)转化成为如下式子:

经过光滑化近似和附加状态变量转换后,不等式约束项将通过新的状态变量转换进微分方程中,一方面减少了不等式约束个数;另一方面路径约束满足情况也可以通过新的状态轨迹精确表示.

注1本文光滑化函数中光滑化因子ε的选取公式参考文献[29]推荐值为ε=r为路径约束的个数,ς为求解NLP问题时设置的求解精度;为了保证路径约束的近似化后的求解精度,惩罚因子推荐选取为ρ≥ρ∗(ρ∗=10/ς).

3.2 非均匀Gauss控制向量参数化

通过对文献[25]得到的HV优化轨迹分析后发现,其控制曲线轨迹呈两端变化剧烈而中间变化平缓的趋势,具有类似于Gauss分布的特性.因此,本文提出非均匀高斯离散方法进行控制向量参数化.引入新的时间变量τ将控制时域[t0,tf]采用式(26)进行时间尺度变换

从而时间区间从[t0,tf]转化到[−1,1].进一步采用Legendre多项式在[−1,1]区间进行离散配点选择,其多项式公式为

由Legendre多项式正交多项式的性质可得,Pn(x)所有的零点都是位于开区间(−1,1)内的单重实根.结合文献[30]的分析和证明可知,当αn和βn的取值分别为αn=0,βn=时,如下矩阵:

的特征值即为Legendre多项式的根.

因而,通过对N次多项式(27)求解,便可以得到开区间(−1,1)内的N个离散点,记作τi(i=1,2,···,N),该离散节点符合Gauss 分布,相邻时间节点的长度不相等,即τk−τk−1τk+1−τk,k=1,2,···,N−1.

结合得到的N个非均匀离散时间节点,则可采用控制向量参数化将转换后的控制时域[−1,1]分成N个时间网格,每个时间网格的节点记为

t0<···

在每个时间段[tk−1,tk]内,将连续的控制变量近似为一组参数函数,控制变量在每个时间段上近似为

此外,为了减少约束个数,采用终端约束转换法将终端约束增广进目标函数中,则问题(31)的目标函数被进一步转化成为如下形式:

其中:ηi,(i=1,2,3)为惩罚系数,其参数取值选择可参见文献[25].最终,问题(P1)被转化为如下NLP问题.

问题2(P2) 给定连续动态方程(13)和状态初始条件(14),通过在整个控制时域[−1,1]内优化控制向量参数σ,使目标函数(32)最大.

经过控制变量参数化后,问题(P2)对于问题(P1)的收敛性可以由定理2保证,篇幅所限,其证明过程读者可以参阅文献[29,31].

3.3 非均匀Gauss控制向量参数化优化方法步骤

基于上述非均匀Gauss时间网格控制向量参数化,给出本文提出方法(简称为CVP–G方法)的实现过程,其算法流程如图1所示,算法步骤如下.

图1 基于非均匀Gauss控制向量参数化优化算法流程图Fig.1 Flow chart of non-uniform Gauss control variable parameterization optimization approach

步骤1输入离散个数N,初始优化参数u0,tf;设置优化问题的求解精度,微分方程求解精度,路径约束惩罚因子ρ,光滑化因子ε,终端状态约束惩罚系数ηi,(i=1,2,3),状态初始值x0,控制边界umin和umax;

步骤2采用光滑化惩罚函数法处理路径约束;

步骤3计算N次Legendre多项式,得到N个非均匀时间网格;

步骤4使用非均匀Gauss时间网格进行控制向量参数化,将问题(P1)离散化为NLP问题(31);

步骤5终端状态惩罚函数处理,得到目标函数(32);

步骤6计算问题(P2)的目标函数关于优化参数的梯度信息;

步骤7结合梯度信息使用SQP算法求解NLP问题;

步骤8输出优化后控制轨迹.

4 仿真校验

以美国波音公司提出的通用航空器(common aero vehicle,CAV)HV模型为测试对象[32],选择再入段HV纵向再入航程作为优化目标函数进行数值优化仿真,因此问题(P1)中相应的目标函数为

CAV高超声速飞行器的各项气动参数取值如表1所示,轨迹优化的初始值和终端约束值见表2.本文的数值仿真实验在配置为1.6 GHz Intel Core i5-8286U处理器和16 GB 2666 MHz Samsung DDR4内存的个人电脑MATLAB 2018a平台上进行.

表1 CAV高超声速飞行器的各项气动参数取值Table 1 Aerodynamic parameters of CAV hypersonic vehicle

表2 CAV–HV轨迹优化边界值约束Table 2 Boundaries of CAV–HV trajectory optimization

在CVP–G方法数值测试过程中,使用fmincon求解器在SQP算法下求解NLP问题,优化误差精度设置为10−4,微分方程求解精度为10−6,所采用的惩罚因子为ρ=10,3个终端约束的惩罚系数分别设定为η1=100,η1=100,η1=10−6.同时,采用文献[25]控制变量离散化个数,设置离散时间节点数为N=30.此外,为了验证本文提出方法求解的有效性,选取Zhang等[24]使用的经典CVP方法和增强型CVP(enhanced CVP,ECVP)方法,Liu等[25]提出的CVP–B–P和CVP–S–P方法在同样热率约束条件下与本文方法进行对比测试.

表3列出了本文方法求解结果与文献方法的优化结果对比.在热率约束为800,000 W/m2时,本文方法得到8160.0 km的纵向再入航程和2088.7 s的飞行时间,与Zhang等[24]使用经典CVP方法和ECVP方法求解得到的6831.0 km和7265.8 km最大航程,Liu等[25]提出的CVP–B–P和CVP–S–P方法得到的7724.8 km和7839.9 km优化结果相比,本文的CVP–G方法求解得到的最大航程和飞行时间均得到提升,比文献[24](EC-VP)和文献[25](CVP–S–P)得到的航程分别提升894.2 km和320.1 km,表明了所提出的非均匀时间节点离散对HV轨迹优化的有效性.

表3 不同方法CAV–HV轨迹优化结果Table 3 Test results of different CAV–HV trajectory optimization methods

进一步,图2给出了CVP–G和CVP–S–P方法在相同离散时间节点数下仿真计算得到的控制曲线,其中CVP–G方法为非均匀高斯离散时间网格,CVP–S–P方法为均匀时间网格.由图可知,本文提出方法得到的控制曲线在中间段振荡较小,在两端的控制也较CVP–S–P方法更为精细.在此控制下,本文的CVP–G方法得到的航程也比CVP–S–P方法远,如图3所示.由于CVP–G和CVP–S–P方法均采用了光滑化处理策略,其热率约束满足情况如图4所示,由图可知,两种方法均能很好满足热率约束限值,显示出本文方法的作用.行高斯配点离散处理,在相同热率约束下对此HV轨迹优化问题进行求解并获得了8237.8 km的最大航程距离,相比于本文优化结果,文献[33]得到的最大航程虽然提高了77.8 km,但其热率路径约束在89∼109 s时违反上限值;文献[4]针对GPM方法在非离散点出现违反约束限值的情况,对GPM方法提出了改进,改进方法在第3次迭代后热率约束得到满足.与文献结果对比后发现,本文只需采用一次光滑化和附加状态转换便可对不等式约束进行很好处理并得到精确的约束满足.

图2 CVP–G和CVP–S–P的控制曲线(N=30)Fig.2 Control curves of CVP–G and CVP–S–P methods(N=30)此外,文献[33]采用伪谱法对不等式路径约束进

图3 CVP–G和CVP–S–P的航程曲线(N=30)Fig.3 Downrange curves of CVP–G and CVP–S–P methods(N=30)

图4 CVP–G和CVP–S–P的热率差值变化曲线(N=30)Fig.4 Heating rate deviation curves of CVP–G and CVP–S–P methods(N=30)

为了进一步分析不同热率约束限值对CAV–HV航程规划的影响,以此在HV 热防护系统设计与轨迹规划中找到最佳组合.采用本文CVP–G 方法对最大热率约束从800,000 W/m2以20,000 W/m2递减至680,000 W/m2下的6种约束条件分别进行轨迹优化求解,以分析热率约束与轨迹规划之间的关系.表4给出了本文CVP–G方法在6 组不同最大热率约束下的CAV–HV轨迹优化结果,包括最大飞行航程和飞行时间两个指标.由表4可知,在6 组测试中,本文CVP–G方法均成功求解,进一步表明所提出方法的有效性.分析表4可知,随着热率约束限值的下降,CAV–HV规划得到的最大航程和飞行时间均呈减少趋势,这表明热率约束限值对HV航程具有显著影响.

表4 不同热率约束下CVP–G方法优化结果Table 4 Optimization results of CVP–G method with different heating rate constraints

为深度分析热率约束对航程的影响程度,图5给出了航程随热率约束变化的趋势图,由图可知,热率约束降低对航程有减少作用,但其降低作用较为轻微.更进一步,图6 给出了当热率从800,000 W/m2降为680,000 W/m2时,热率和航程变化的百分比趋势,分析可知,当热率约束最大降低15%时,对应的最大航程从8160 km降到7902.5 km,航程损失仅3.16%,由此表明单纯增大HV热防护系统的热率防护上限极值对增加HV的航程作用有限.同时,通过与表3种数据对比后可知,本文CVP–G方法在680,000 W/m2热率约束下得到的最大航程为7902.5 km,优于CVP–S–P方法在800,000 W/m2热率约束下的最大航程值,表明本文方法在得到更优航程规划的同时,对于HV热防护系统的热率防护限值设计具有重要的指导作用,将为高超声速飞行器的热率监测设备工作范围设计提供支撑.

图5 CVP–G规划航程随热率约束变化趋势图Fig.5 Trend of CVP–G optimization range variation with different heating rate constraints

图6 CVP–G规划航程随热率约束变化百分比趋势图Fig.6 Percent trend of CVP–G optimization range variation with different heating rate constraints

为了展示不同热率约束下CVP–G方法优化得到的控制曲线和状态曲线变化,为HV热防护系统设计提供参考,图7给出了6种热率约束下的CAV–HV最优攻角控制曲线,图8–11展示了相对应的高度、速度、航迹倾角、航程状态变化曲线,图12给出了不同热率约束限值下的热率差值变化曲线.由图可知,在不同最大热率约束下,控制曲线具有同样的趋势,但在60∼100 s时攻角控制变化最为剧烈,这与图12中的热率差值在60∼100 s剧烈变化相对应,表明攻角的迅速变化能够减少由气动力产生的热量,从而满足热率约束;同时,分析图8–11状态曲线可知,随着热率约束限值的降低,HV存在的振荡现象在一定程度上呈减弱趋势,表明降低热率约束限值能够加强飞行器飞行弹道的稳定性.综合以上分析可以发现,采用本文提出的CVP–G优化算法可以通过对HV飞行器攻角的优化控制使飞行器在较低的热率约束限值下得到航程规划参考,将为不同航程需求下飞行器热保护系统设计提供理论算法支撑.同时,通过分析不同热率约束下的航迹规划还可以为相关监测的仪器仪表设备所需要承受的工作环境提供理论计算依据.

图7 不同热率约束下攻角控制曲线(CVP–G方法)Fig.7 Angle of attack control curves under different heating rate constraints(CVP–G method)

图8 不同热率约束下CAV–HV高度曲线(CVP–G方法)Fig.8 Altitude curves under different heating rate constraints(CVP–G method)

图9 不同热率约束下CAV–HV速度曲线(CVP–G方法)Fig.9 Velocity curves under different heating rate constraints(CVP–G method)

图10 不同热率约束下CAV–HV航迹倾角曲线(CVP–G方法)Fig.10 Path angle curves under different heating rate constraints(CVP–G method)5 结论

图11 不同热率约束下CAV–HV航程曲线(CVP–G方法)Fig.11 Downrange curves under different heating rate constraints(CVP–G method)

图12 不同热率约束限值下热率差值变化曲线(CVP–G方法)Fig.12 Heating rate deviation curves under different heating rate constraints(CVP–G method)

本文提出一种用于热率约束下HV再入航程轨迹优化问题计算的非均匀Gauss离散时间网格控制变量参数化算法.该方法采用了光滑化近似和附加状态变量转换相结合的不等式路径约束处理方式;在CVP算法框架下,提出了基于Gauss时间网格的控制变量参数化策略,可以有效改善HV攻角的控制精度,提升HV再入航程.通过对CAV–HV通用航空器的轨迹优化仿真测试发现,与同等离散参数下的均匀离散网格CVP–S–P方法相比,改进方法得到的再入航程距离更远,最大航程增加320.1 km,提升4.1%,显示了所提出方法对HV轨迹优化的有效性和实用价值.此外,基于本文算法,还分析了不同热率约束限值对航程规划的影响,仿真结果显示当热率限值降低15%时,最大航程损失仅3.16%,表明单纯增大HV热防护系统的热率防护上限极值对于增加HV的航程作用有限,这将为高超声速飞行器TPS的设计提供重要的理论指导.下一步,将联合具有实验平台的研究院所,共同开展不同热率约束下HV航迹规划实验验证工作,并在热率约束限值与最大航程规划中寻求最佳组合,为HV的设计提供理论实验支撑.