采用滚动伪谱优化的组合动力飞行器上升段制导方法

唐湘佶,李兆亭,张洪波

(国防科技大学 空天科学学院,湖南 长沙 410073)

组合动力飞行器是未来空天飞行器发展的重要方向,包括吸气式临近空间高超声速飞行器、水平起降空天往返飞行器[1]等。本文针对可水平起降的涡轮基组合动力(TBCC)飞行器,飞行任务场景设定为水平起飞,加速爬升后在临近空间进行高超声速巡航。该场景适用于高超声速飞机、两级入轨空天飞行器一子级[2]等飞行器。

上述飞行场景的任务需求和组合动力飞行器上升段的飞行特性要求制导系统具有较好的自适应能力。上升过程飞行跨度大,具有宽速域多模态气动/推进强耦合特性[2];动力系统对控制及其动态过程的约束严格,性能受状态与环境影响大;不同速度区间内动力形式不同,系统特性多变,非线性强;飞行器模型与环境存在不确定性,飞行过程中误差累积效应显著。复杂多变的系统和飞行环境加上两者间的强耦合关系导致制导难度大,需要制导算法具有较好的精度以及自适应性。

上升制导中较为常见的标称轨迹跟踪制导对于较大偏差和扰动的纠偏能力有限,对于复杂多变的系统模型,制导律设计难度大,且难以满足对自适应能力的需求。提高制导自适应能力的技术途径之一是基于在线解算最优控制的制导技术,这一方法在工程中已被应用,例如运载火箭上升段的迭代制导[3]等。随着计算机技术的发展以及数值计算方法的深入应用,一系列基于轨迹优化技术的制导新方法被提出[4],其中最优控制问题在线求解方法包括间接法[5-6]、直接法[4,7]和混合法[8]。

直接法中的伪谱法具有快速求解最优控制问题的潜力,可用于最优控制实时计算,实现最优反馈控制[9]。国内外对伪谱法在最优反馈控制问题及制导方面的应用进行了深入研究。ROSS等[10]将基于伪谱法的实时最优反馈控制应用于卫星姿态控制。BOLLINO等[11]研究了基于伪谱法的X-33飞行器再入段最优非线性反馈制导。张友安等[12]研究了伪谱法在高超声速飞行器再入制导中的应用,理论上分析了闭环控制系统的有界稳定。闫循良等[13]将改进的伪谱反馈控制应用于运载火箭的远程变轨制导,通过状态量缩减降低了制导问题规模,保证了制导算法实时性。许东欢等[14]针对空天飞行器研究了基于伪谱法的闭环制导算法可行性。张志国等[4]研究了伪谱法在大气层外运载火箭上升段制导律中的应用,与迭代制导的对比结果验证了该方法的制导精度。上述研究表明了经典伪谱法可用于制导,但针对的问题约束相对较少,控制指令计算难度相对较低,部分研究中未探讨方法对误差的适应性及计算实时性。

本文针对涡轮基组合动力飞行器的上升过程,研究了采用滚动伪谱优化的制导律。应用团队自研伪谱法求解器SPTOS,利用分段伪谱法在线求解最优控制。分段伪谱法是对经典伪谱法的改进,采用分段拟合提高收敛速率,该方法尤其适用于变量变化不够平滑、有间断点等情形[15]。制导过程以双层滚动优化框架为基础,引入离散网格的同步滚动更新改进求解效率,制导周期根据实时状态偏差自适应调整。介绍了仿真及制导过程中所用动力学模型和相关约束,及采用滚动伪谱优化的制导方法,展示了仿真结果并进行了讨论。

1 飞行器动力学模型及约束

针对制导过程,主要研究飞行器的质心运动,而加速爬升的运动过程主要在纵平面内。基于上述分析,建立飞行器质心运动方程及相关约束条件。

1.1 质心运动方程

仅考虑纵平面内的运动,将状态量取为x=(HRvθm)T,分别表示高度,航程,速度,速度倾角以及飞行器质量。控制量取为u=(φtα)T,分别表示节流系数与攻角。采用圆球假设并忽略地球自转,建立如下飞行器质心运动方程:

(1)

式中:re为地球半径;ρ为大气密度;FP,FL,FD分别为发动机推力、升力与阻力;KS为空气/燃料流量转换系数;Sc为进气口截面积;CA(Ma,α)为等效进气面积系数,是马赫数和攻角的函数,通过插值计算。

气动力的计算公式如下:

(2)

式中:Sref为飞行器参考面积;CL(Ma,α),CD(Ma,α)分别为升力系数和阻力系数,与马赫数和攻角有关,通过插值计算。

本文考虑的组合动力系统包括3种工作模态:在Ma<2时是涡轮动力,Ma=2～6为冲压模态,Ma>6时是超燃冲压动力。推力可统一描述成如下形式:

(3)

式中:g0为水平面重力加速度;Isp(Ma,φt)为发动机比冲,与马赫数和节流系数有关,通过插值计算。上述飞行器模型参数及插值数据参考文献[16]。

组合动力系统的特性如图1所示。主要性能参数受状态和控制的影响大,由图1(a)可见,随着速度的增加,比冲整体减小;不同动力形式下比冲受速度的影响不同;节流系数提高,比冲增加。如图1(b)所示,在低速段,进气效率受攻角影响小;跨声速段的进气效率有所下降;在冲压动力段,气体压缩过程使进气效率不断提升,由于需要借助机体实现来流的减速增压,高速推进时的进气效率受攻角影响明显。

图1 组合动力系统特性

1.2 约束条件

鉴于飞行任务要求及飞行器自身特性,实际飞行过程中存在诸多约束条件,本文考虑3类典型约束:状态量约束、控制量约束及过程约束。

①状态量约束。状态量在初始时刻t0及末时刻tf的约束形式为

H(t0)=H0,v(t0)=v0,θ(t0)=θ0,m(t0)=m0
H(tf)=Hf,v(tf)=vf,θ(tf)∈[θmin,θmax]

②控制量约束。节流系数、攻角及其变化率均存在上下界,约束形式为

式中:KH为热流密度常数。

上述符号中,下标min表示下界,max表示上界;下标0表示初始时刻的值,f表示终端时刻的值。

2 采用滚动伪谱优化的闭路制导方法

本文制导方法的核心为数值最优制导律。在双层滚动优化框架下,通过外层滚动优化使剩余时间内制导指令随状态反馈更新,实现闭路制导过程,并通过制导周期的在线调整进一步提高方法的自适应能力;内层利用分段伪谱法实时求解最优控制,在内层迭代过程中引入网格更新策略以提高计算实时性。

2.1 基于双层滚动优化的制导策略

制导策略如图2所示。

图2 基于双层滚动优化的制导策略

为提高初次指令解算效率,将分段伪谱法离线规划的最优轨迹作为首次指令更新时的参考。制导过程中单次指令解算的具体步骤如下:

①根据实时状态反馈x(ti)建立多约束最优控制问题Pi(x(ti)),该问题的初始状态与时刻即为状态反馈x(ti)及相应的采样时刻ti;

②通过分段伪谱法将连续问题离散化为非线性规划(NLP)问题,离散时的初始离散网格根据前次解算的离散网格Mi-1进行更新;

③对前次优化结果Ui-1,Xi-1进行插值,获得新的离散点处的迭代初值,然后求解NLP问题;

⑤当出现求解结果不收敛、求解时间过长等情况时,在相应时刻不进行制导指令更新,沿用前次计算结果,直至成功在线求解后再次更新。

制导周期Tg采用自适应调整方法,由实时状态偏差Δx在线控制。自适应周期制导过程如图3所示,图中,tc为在线解算时间。制导过程中按一定采样周期实时获取飞行状态,并计算与参考状态的偏差,设定偏差容限ΔxLT,在Δx(ti)≥ΔxLT的一系列采样时刻ti进行制导指令的在线解算与更新。实际制导周期Tg受偏差累积特性及在线解算时间tc影响,为可变周期,当解算时间非常短时该部分影响可忽略。

图3 自适应周期制导过程

本文以高度偏差ΔH作为周期自适应调整依据,为适应飞行过程中不同的飞行状态与不同阶段的制导需求,建立偏差容限函数,形式为

式中:Hd为误差容限切换高度。

制导周期的自适应调整能使指令更新频率与实际状态偏差累积特性相匹配。偏差容限的作用一方面在于避免误差累积过大导致制导修正的难度过高;另一方面,当状态偏差较小时,更新前后的控制量差别相对较小,指令解算速度更快。因此设定合适的偏差容限能够较好地平衡计算效率与制导精度。

2.2 分段伪谱法求解最优控制问题

多约束非线性最优控制问题具有通用的标准形式,包括性能指标和各类约束。本文以燃料最省作为优化目标,等价于使终端质量最大,问题形式如下:

如图4所示,分段伪谱法的基本原理是将连续问题直接离散化,利用数值优化算法求解;之后评估离散解对原问题解的逼近程度,若达不到精度要求,则进行离散网格的细化,通过构造序列NLP问题迭代求解;否则直接输出结果。

图4 分段伪谱法求解流程

本文选择Legendre-Gauss-Radau(LGR)点作为离散点,其分布由Legendre正交多项式的根决定。LGR点的定义域为[-1,1),便于描述各分段点处连续条件。LGR点与各区间内时间具有如下关系:

离散化的关键是动力学微分约束的转换,其余约束可直接应用于离散变量集{X,U}。动力学约束的转换关键在于对微分算子的有限逼近。各区间归一化后,分别用拉格朗日插值多项式近似原问题的状态量,LGR点数目Nk确定时,插值基函数形式固定。计算公式如下:

将决策变量取为{X,U}及终端时刻tf,则原连续问题转换成如下NLP问题:

式中:l=1,2,…,Nk;k=1,2,…,K。

本文运用开源求解器IPOPT求解上述大规模NLP问题,其基本原理为内点法,给定解的初始参考值,经过迭代可收敛至NLP问题最优解。序列NLP问题构造过程中的数值解评估方法及离散网格细化方法参考文献[15]。

2.3 离散网格滚动更新

离散网格关系到NLP问题的规模及问题特性,是求解效率的主要影响因素之一。与分段伪谱法的网格细化不同,滚动更新针对初始离散网格。建立合适的初始离散网格能降低NLP问题难度,还能避免构造序列NLP问题,提高求解效率。

离散网格M的滚动更新过程为对3个决定参数(分段数K、各区间的占比hk与离散点数Nk)的迭代计算。根据文献[17]中的收敛理论,对于分段伪谱法,当解的逼近程度一定时,区间长度与段内离散点数近似具有等比关系,区间越短,达到一定逼近程度所需要的插值多项式阶数越低,离散点个数越少。在理想情况下,滚动优化过程中每次在线求解的结果是对上一次计算结果的截取。为相应调整NLP问题的规模,离散网格也应进行截取。

图5 离散网格更新示意图

在继承离散网格Mi-1的基础上,根据当前时刻t′0=t0+Δt,对初始离散网格的3个参数进行调整,方法如下:

①以前次成功求解时的离散网格Mi-1为参考,判断当前时刻t′0所处区间顺序z;

3 仿真分析

仿真中初始状态与期望终端状态如表1所示,各变量含义如前所述。

表1 初始状态与期望终端状态

表2 变量范围约束

采样周期取0.1 s,偏差容限函数参数ΔHmax1=30 m,ΔHmax2=10 m,Hd=43 km。

仿真硬件条件为Intel Core i7-4712MQ、2.3 GHz的PC。采用4阶定步长Runge-Kutta方法进行数值积分,步长取0.05 s。

取大气密度偏差为10%,升力系数偏差为10%,阻力系数偏差为10%,推力偏差为2%。组合极值拉偏仿真的结果如图6和图7所示,控制量及其变化率均在约束范围内,终端速度倾角为0.07°,也满足约束。

图6 速度和高度随时间变化曲线

图7 攻角和节流系数随时间变化曲线

由图可见,加速到期望速度的时间相比于高度爬升过程更短,这是由于随高度增加,大气密度减小,耗油率增加,因此为节省燃料,优先完成加速过程,这体现出在误差条件下,迭代更新后的制导指令仍保持最优控制的原则。虽然指令的整体动态变化显著,但一个周期内指令是根据离散解进行插值得到,在小范围内能够保证控制量相对平缓,易于控制系统实现。

动压与热流密度均在上限约束范围内,计算结果如图8所示。仿真结果表明:制导算法具有可行性;采用滚动伪谱优化的制导方法能够在满足基本约束的条件下完成上升制导任务。

图8 热流密度与动压随时间变化曲线

指令更新结果与相应参考状态的对比结果如图9和图10所示。以一个制导周期为例,为突出偏差影响,在高度-速度域内进行参考状态量的对比。

图9 参考状态对比

由图9可见,制导过程中,剩余时间内制导指令随状态反馈自适应修正,修正后的控制量对应的参考状态相应改变。

指令迭代更新对比如图10所示。外层滚动优化输入状态反馈输出修正后的控制,状态反馈隐含了误差信息,为保持最省燃料爬升,剩余时间内攻角指令不断调整。

以2个制导周期内的攻角指令为例,在引入状态反馈后,剩余时间内控制发生改变;仅截取一个制导周期内的指令作用于飞行器,该周期内的指令变化较为规律。

自适应制导周期如图11所示,在线求解次数为34次,每次求解均收敛至最优解。平均周期约为4 s,最小周期为0.6 s,此时高度超过Hd,指令更新频率高。

图11 自适应制导周期

在线计算耗时如图12所示。由图12可看出,平均在线计算时间约为270 ms,大部分时间在200 ms左右,满足一定的实时性要求。在滚动优化的内层迭代过程中,离散网格同步滚动更新,改善了NLP问题的特性与规模;基于前次优化结果的迭代使优化算法收敛更快,提高了在线求解实时性。

图12 在线计算时间

本文仅从计算方法的角度寻求在线快速求解。在实际应用中,结合硬件设备优化计算程序,可进一步提高求解效率。

为进一步验证本文方法对误差和飞行器特性的适应性,进行了500次的蒙特卡洛打靶仿真。各误差项的3σ值如表3所示。

表3 蒙特卡洛打靶仿真误差项

终端高度误差与速度误差的统计结果如表4所示。由表可见,状态偏差均值为负,在最省燃料爬升过程中推进效率接近最优,在误差条件下容易出现飞行器爬升能力趋近饱和,控制修正能力有限,使加速爬升的末端高度状态低于期望值。

表4 打靶仿真统计结果

随着周期自适应调整,任务末段更新频率的提高保证了一定的终端制导精度。终端高度与速度精度整体均在可接受范围内,在不同误差条件下,制导方法均能够适应复杂多变的组合动力系统特性。

4 结论

本文研究了采用滚动伪谱优化的闭路制导方法。在双层滚动优化框架下,根据实时状态反馈建立多约束最优控制问题,加入离散网格同步更新,运用分段伪谱法进行快速计算,并在线自适应调整制导周期,从而实现制导指令的迭代更新。对制导过程进行数值仿真,主要结论如下:

①制导方法的精度满足任务需求,在线求解过程始终满足优化指标,具有最优反馈性质。根据状态反馈进行指令的迭代更新,能够抑制模型参数不确定和外界干扰引起的误差影响,使飞行器达到期望的飞行状态,方法具有一定的鲁棒性。

②引入离散解迭代以及离散网格滚动更新,提高了分段伪谱法在线求解多约束最优控制问题的效率,使该方法满足一定实时性要求。

③制导周期的自适应调整能够较好地平衡制导效率与制导精度,提高了制导方法自适应性。