基于随机响应面与混沌多项式的鲁棒轨迹优化

王培臣, 闫循良,*, 王 宽, 郑 雄

(1. 西北工业大学航天学院, 陕西 西安 710072; 2. 中国运载火箭技术研究院, 北京 100076)

0 引 言

作为飞行器设计的关键技术之一,轨迹优化逐渐成为高超声速等航空航天领域研究的热点问题[1-2]。然而,现有的轨迹优化研究通常仅注重标称情况下的轨迹性能优化提升,没有考虑真实情况下轨迹的鲁棒性和可靠性[3-6]。在诸如高超声速再入等实际飞行任务中,初始状态、动力学及环境参数等均存在诸多不确定性,势必使得飞行器沿标称轨迹飞行时产生较大的过程及终端偏差,进而增大约束违反的风险,降低飞行可靠性和终端任务精度。因此,为了提高优化轨迹的可靠性和抗干扰能力,部分学者逐步将研究重心聚焦于考虑不确定性的鲁棒轨迹优化[7-13],即在轨迹设计阶段事先考虑不确定性的影响,从而降低由轨迹规划迭代耗时和轨迹偏差修正带来的制导控制系统负担。

求解鲁棒轨迹优化问题的典型策略为:先将含有随机微分方程的随机最优控制问题转化为包含扩展确定性微分方程的确定最优控制问题,进而采用高效的轨迹优化算法对该高维问题进行数值求解[8],即包含不确定性量化传播和确定性轨迹优化两部分关键技术。典型的不确定性量化传播理论和方法包括蒙特卡罗(Monte Carlo,MC)方法[14]、线性协方差分析法[15-16]、无迹变换法[17]及混沌多项式展开(polynomial chaos expansion,PCE)方法[18-19]等。MC方法因简单易行而被广泛应用于不确定性量化及传播,但大量的数值模拟会导致统计结果收敛速度较慢[14];线性协方差分析法可对线性高斯过程模型中的不确定性进行量化传播,且具有较高的计算效率,但其不适用于高度非线性的轨迹优化问题[15];相比之下,PCE方法将随机变量表示为随机正交多项式的加权求和,可对具有任意分布类型的随机变量实现较为精确的逼近,同时能够以较低的计算代价达到与MC方法相当的精度,近年来在火星进入轨迹不确定性量化[18-21]、流体力学不确定性分析[22-25]中得到了广泛应用。PCE方法与系统模型结合的方式有两种,即嵌入式[17]与非嵌入式[19,23-25],嵌入式混沌多项式(intrusive polynomial chaos, IPC)将随机变量的混沌多项式展开代入动力学微分方程中,并利用伽辽金投影法将随机微分方程转换为一组以混沌多项式系数为状态量的高维确定性微分方程组,最终通过求解该微分方程组完成不确定性传播[21-22],具有较高的计算效率;非IPC(non-IPC, NIPC)则将原系统模型视为“黑箱”,进而根据模型的输入输出逼近混沌多项式系数,无需针对具体问题进行推导建模,因此能够处理一般非线性系统的不确定性量化传播问题,具有较好的通用性及工程应用价值。

而对于另一关键技术,即确定性轨迹优化而言,常见的方法主要包括间接法和直接法。间接法因对解的初值猜测较敏感、需要对最优性必要条件进行推导且无法处理复杂约束等原因,难以在复杂约束轨迹优化问题中得到广泛应用。近年来,以伪谱法[2]、凸优化[4-6]等为代表的直接法已被广泛应用于高超声速等航空航天领域的轨迹优化任务,并展现出了良好的数值计算优势和普适性等特点。其中,凸优化方法因具有多项式复杂度、收敛速度快,且其局部最优解理论上即为全局最优解等特点,已逐步被用于求解复杂约束下的轨迹优化问题以及考虑不确定性的高维状态扩展优化问题。

鉴于PCE方法与直接法的优势,诸多研究将二者进行有机结合,设计了多种基于PCE方法与直接法的鲁棒轨迹优化算法和求解策略。Fisher等[8]首次提出了数值求解随机最优控制问题的框架,该框架利用IPC将随机最优控制问题转化为高维确定性最优控制问题,并采用配点法对该问题进行求解。Li等[9]将NIPC与高斯伪谱法相结合,用于解决气动系数不确定的高超声速飞机的爬升鲁棒轨迹优化问题。Jiang等[10]将NIPC与伪谱法相结合,通过调用OpenCossan工具箱进行不确定量化传播,并采用并行计算技术提升高维确定性优化问题的求解效率,解决了初始状态及模型参数均不确定的火星再入鲁棒轨迹优化设计问题。Wang等[11]提出一种将IPC与凸优化结合的鲁棒轨迹优化方法,并与基于IPC与高斯伪谱法的方法进行对比,证明了所提方法在计算精度相当的情况下,计算效率显著提高。杨奔等[12]基于最小二乘法建立了NIPC与凸优化相结合的鲁棒轨迹优化模型,提升了气动参数不确定条件下轨迹的可靠性。

基于PCE和直接法的随机问题求解技术虽然已被逐步用于鲁棒轨迹优化,但仍存在计算效率低、通用性较差以及难于处理多维不确定性等问题,这些问题进而限制了其在复杂约束及强非线性鲁棒优化问题(如航空航天领域鲁棒轨迹优化)中的应用。首先,以高斯伪谱法为代表的传统直接法虽然解决了部分鲁棒轨迹优化问题,但其仍面临计算效率低、难于快速求解多维不确定问题等局限,而凸优化技术与PCE的结合仍然不够完善。其次,基于IPC的优化方法虽然具有较高的计算效率,但需要针对不同问题进行推导建模和代码改写,由于推导过程存在模型简化近似,只适用于处理多项式非线性问题[9];基于NIPC的优化方法虽然能够避免上述问题,但PCE系数求解大都基于求积、采样的方法,例如高斯求积(Gauss quadrature,GQ)方法[9-10]、稀疏网格求积法[13],而随着系统随机变量维数及非线性程度增大,这类方法的状态扩展维数会呈指数增长,使得状态离散节点数量显著增加,从而导致解决大规模NLP问题时非常困难,计算成本高昂且效率低下,尤其是当不确定因素维数n>3时,大量的内存消耗和计算成本会严重阻碍优化过程,甚至导致优化失败。

值得注意的是,除去上述基于求积的方法,NIPC系数也可采用随机响应面法(stochastic response surface method, SRSM)[23]进行计算。该方法基于代理模型思想,通过拟合系统输出与输入之间的关系获得混沌多项式模型,并结合典型抽样方法进行少量采样即可实现多项式系数计算,可以显著降低模型转化维数和规模,因此具有较高的精度和计算效率。目前,SRSM已被广泛应用于流体力学不确定性分析中[24-26],但尚未在鲁棒轨迹优化设计中广泛应用。

因此,本文针对航空航天领域鲁棒轨迹优化算法存在的问题,提出一种高效的鲁棒轨迹优化设计方法。通过构造基于NIPC与SRSM的不确定性量化传播模型,实现了不确定性随机问题的有效转化;随后,将该模型与凸优化方法相结合,设计了一种基于序列凸优化算法的求解策略,以实现对该高维确定性问题的快速求解。最终,以高超声速滑翔再入为例进行仿真,验证了本文方法的有效性。

1 鲁棒轨迹优化问题描述

一般的轨迹优化问题可以归结为如下确定性最优控制问题PD:

(1)

当系统参数向量w及状态初值x(t0)存在不确定性时,目标函数J、动力学函数f及约束函数均变为随机函数。因此,原确定性优化问题PD转变为随机优化问题。为了降低优化轨迹对不确定性的敏感度,即当存在不确定性时,预期的轨迹离散度在统计意义上尽可能小,通常需要在目标函数、约束函数中引入随机量的统计矩。因此,一般鲁棒轨迹优化问题PR可描述为

(2)

式中:下标μ和σ分别代表均值和标准差;kJ、kg为待设计权重系数;x0μ与x0σ分别为初始状态的均值和标准差;εf为状态终值标准差的上界。当采用上述模型进行轨迹优化时,需要引入不确定性量化传播技术计算相关统计量并处理随机动力学方程。

需要注意的是:①εf需要基于先验统计信息确定,若取值不当可能造成问题无解,为了提升优化算法的收敛性,可将xσ(tf)以罚函数的形式引入目标函数中;② 上述鲁棒轨迹优化模型同时考虑了轨迹的鲁棒性及可靠性,其中,在目标函数及端点约束中引入均值、标准差项可提升轨迹的鲁棒性,而通过引入考虑统计矩的过程约束模型,则可以保证以一定概率满足过程约束从而提升轨迹可靠性;③ 通过设计权重系数kJ、kg可实现轨迹的鲁棒性及可靠性之间的权衡,具体讨论见后文仿真部分。

下文将详细给出考虑初始状态x(t0)及系统参数向量w不确定性的鲁棒轨迹优化算法。首先,采用基于响应面法的PCE技术将原始鲁棒轨迹优化问题PR转化为状态扩展的高维确定性优化问题;随后,对该问题进行凸化、离散处理,并采用序列凸优化算法进行求解。

2 不确定量化传播建模与问题转化

本节将NIPC与SRSM相结合,同时引入拉丁超立方采样策略,构建了不确定性量化传播模型,从而将鲁棒轨迹优化问题转化为高维状态空间中的等价确定性优化问题。以下给出算法原理及实现步骤。

2.1 随机状态方程转化

假设随机状态向量x为n维,系统参数向量w为s维,则每个状态变量xi(t)及参数wk可展开为PCE多项式[20]:

(3)

(4)

式中:Δ=[Δ1,Δ2,…,Δd]为d维随机向量;Φj(Δ)为正交多项式函数,可通过与Δ分布类型相对应的一元正交基函数φj(Δl)(l=1,2,…,d)求张量积得到;xij(t)和wkj为PCE系数,展开项数P+1=[(p+d)!]/(p!d!),即由正交基函数φj(Δl)的阶数p和随机变量维数d决定;根据Askey法则[20],不同随机变量Δl的分布类型对应不同的正交基函数{φj(Δl)},具体对应形式如表1所示。

表1 正交基函数与随机变量分布类型对应关系Table 1 Correspondence relationship between distribution type of orthogonal basis function and random variables

PC(polynomial chaos,PC)系数wkj及xij(t0)可根据wk及xi(t0)的分布确定,而随机状态xi(t)的PC系数xij(t)未知,下面采用SRSM导出求解系数xij(t)的表达式。

将式(3)、式(4)代入式(2)的微分方程中,可得

(5)

采用拉丁超立方[24]进行ns次随机采样,并将随机样本点{Δ1,Δ2,…,Δns}代入方程(5),可得微分方程组:

(6)

式中:X与W分别表示随机状态x与参数w的PC系数向量,且有

X=[x10(t),x11(t),…,x1P(t),x20(t),x21(t),…,x2P(t),…,
xn0(t),xn1(t),…,xnP(t)]T

(7)

W=[w10,w11,…,w1P,w20,w21,…,w2P,…,ws0,ws1,…,wsP]T

(8)

(9)

设Ω(Δm)=[Φ0(Δm),Φ1(Δm),…,ΦP(Δm)],(m=1, 2,…,ns),定义

Θ(Δm)=In⊗Ω(Δm)

(10)

则Ψ=[Θ(Δ1)T,Θ(Δ2)T,…,Θ(Δns)T]T∈Rn·ns×n(P+1)。

若ns≥P+1,则可通过最小二乘回归法求解式的微分方程组,即有

(11)

式中:矩阵H=(ΨTΨ)-1ΨT∈Rn(P+1)×n·ns。为保证该模型的逼近精度,本文统一取采样点ns=2(P+1)[22],后文不再赘述。式(11)即为关于扩展系数X与控制量u的高维状态方程。

2.2 过程约束及目标函数的PC展开

由于非线性约束函数g(x(t),u(t),t,w)与目标函数J(x(t),u(t),t,w)均为随机状态变量x(t)的函数,故也可采用类似方法进行混沌多项式展开,具体如下:

(12)

(13)

式中:gj(x(t),u(t),t)与Jj(x(t),u(t),t)为对应函数的PC多项式系数。将随机样本点{Δ1,Δ2,…,Δns}代入方程(12)、方程(13)可得

(14)

(15)

(16)

(17)

(18)

(19)

类似地,有

(20)

(21)

(22)

2.3 统计量计算

根据文献[15],可将式(2)中随机状态x、过程约束函数g(x(t),u(t),t,w)及目标函数J(x(t),u(t),t,w)的均值和方差表示为PC多项式系数的函数,即

(23)

(24)

(25)

式中:〈·〉表示内积运算;gj,Jj(j=1,2,…,P)可通过式(21)、式(22)表示为X、u的非线性函数,进而可通过以上公式计算状态、过程约束与目标函数的统计量。

至此,已将鲁棒轨迹优化问题PR转化为关于扩展状态X及控制量u的高维确定性优化问题PED,进而可利用凸优化算法进行处理求解,限于篇幅,问题PED在此不再赘述。

3 基于凸优化的鲁棒优化求解策略

考虑到复杂轨迹优化问题大多是非凸的,如何将非凸问题进行凸化处理是应用凸优化技术求解最优轨迹的关键,并吸引了广大学者对此开展了深入研究,相关成果亦在航空航天领域最优控制及轨迹优化问题上得到了较为广泛的应用[4-6]。

本节引入逐次线性化[27]等技术对第2节中的高维确定性优化问题PED进行凸化处理,并将该问题离散处理为参数凸优化问题,进而建立相应的序列凸优化迭代算法,对该问题进行求解。序列凸优化方法通过外环迭代更新参考轨迹,内环求解凸问题逐步逼近原问题的最优解[27]。以下给出凸化、离散处理过程以及迭代算法的原理和实现步骤。

3.1 凸化处理

3.1.1 状态方程凸化

分析式(11)可知,H为常值矩阵,F(X,u,t,W)为关于扩展状态X的非线性函数,因此该式为非线性微分方程。下面采用逐次线性化方法对该式进行凸化处理。

将式(11)在参考轨迹{Xk,uk}附近进行一阶泰勒展开,可得

(26)

定义:

(27)

式中:∂F/∂X|(Xk,uk,t,W)与∂F/∂u|(Xk,uk,t,W)分别为参考轨迹上函数F(X,u,t,W)相对于状态X和控制量u的偏导数。其中,参考轨迹{Xk,uk}可通过迭代算法上一步优化结果得到,第一步优化采用的参考轨迹初值{X0,u0}可以X(t0)作为初始状态,以标称条件下的确定性优化结果作为控制量u(t)积分式(11)得到。

为保证线性化的精度,增加信赖域约束:

(28)

式中:δX和δu表示信赖域半径。

3.1.2 过程约束及目标函数凸化

类似地,考虑式(2)中的过程约束及目标函数均为关于扩展状态X及控制量u的非线性函数,仍采用逐次线性化方法对其进行凸化处理。

在参考轨迹{Xk,uk}附近进行一阶泰勒展开,可得:

gμ(X,u,t,W)+kggσ(X,u,t,W)=
GAX+GBu+GC

(29)

Jμ(X,u,t,W)+kJJσ(X,u,t,W)=
JAX+JBu+JC

(30)

定义:

(31)

(32)

3.1.3 状态不等式约束凸化

考虑式(2)中的状态不等式约束为关于扩展状态X的非线性函数,其凸化处理方法[11]如下。

由式(23)可知,随机状态变量xi(t)的方差可表示为

(33)

Q=LTL

(34)

状态变量xi的标准差可表示为

(35)

因此,式(2)中的状态不等式约束可无损凸化为如下二阶锥约束:

(36)

式中:ximax与ximin分别为第i个状态量对应的上下边界。

至此,问题PED转化为以下等价的高维确定性轨迹凸优化问题PED-Convex:

(37)

需要说明的是,问题PED-Convex的待优化状态为X,扩展为问题PR的(P+1)倍。

3.2 离散化处理

为了应用数值方法对问题PED-Convex进行求解,需要对状态量和控制量进行离散化处理。基于梯形法则,将自变量变化域[t0,tf]等间距离散为N个间隔,自变量离散为{t0,t1,…,tN},状态量离散为{X0,X1,…,XN},控制量离散为{u0,u1,…,uN}。

对式(37)的线性动力学方程进行离散化,可得:

(38)

类似地,对问题PED-Convex的其他约束及目标函数进行离散化处理,最终可得参数凸问题PED-disConvex:

(39)

该问题的待优化状态扩展为问题PR的(P+1)(N+1)倍。

3.3 序列凸优化求解策略

综上所述,非线性随机最优控制问题PR的解可通过迭代求解问题PED-disConvex进行逼近,并以当前获得轨迹不断更新参考轨迹,以提高解的收敛性和精度,即令

(40)

直到满足约束

(41)

式中:εX和εu为收敛误差限。

图1给出了基于序列凸优化的鲁棒轨迹优化算法原理框图。

图1 鲁棒轨迹优化算法原理Fig.1 Principle of robust trajectory optimization algorithm

4 仿真分析

本节首先通过定量分析对两种不确定性量化传播算法的状态扩展维数及计算效率进行比较,之后以典型的高超声速滑翔问题为例,开展多不确定性条件下的鲁棒轨迹快速优化仿真,以验证本文所设计算法的有效性。

4.1 不确定性量化传播算法定量对比分析

如前所述,当采用非嵌入式混沌多项式进行参数不确定性量化传播时,时常用到两种典型模型转化方法,即SRSM方法和GQ方法,以下通过定量分析对这两种方法的计算效率进行比较。

一般而言,对于具有p维随机变量的不确定问题,若PC多项式展开的阶数为p,根据文献[22]可知,GQ方法的状态扩展倍数为p,而SRSM的状态扩展倍数为p。表2给出了两种方法在不同p、d取值情况下的状态扩展情况。

表2 不同条件下GQ与SRSM的状态扩展性能对比Table 2 Comparison of state expansion between GQ and SRSM under different conditions

由表2可知,在p、d较小的情况下,两种方法的状态扩展倍数比较接近;随着p的增大,两种方法的状态扩展倍数均有所增大。相较于GQ方法,SRSM的增长速度较为缓慢;随着随机变量维数d的不断增大,GQ方法的扩展倍数呈指数增长,而SRSM增长较为缓慢。

可见,相较于GQ方法,本文采用的SRSM具有更高的计算效率,尤其在处理高维不确定性问题时具有显著优势,而GQ方法一般只能适用于d≤3的情况。此外,下文仿真部分亦给出了两种方法在精度与计算效率方面的进一步对比分析。

4.2 高超声速滑翔鲁棒轨迹优化仿真

高超声速滑翔飞行具有航程远、速度高、飞行环境复杂等典型特性,且面临多种复杂约束以及参数不确定性和干扰的影响,这使得其轨迹优化问题成为一个典型的快时变、非线性、强耦合、非凸的随机最优控制问题。文献[4-7]研究表明,标称条件下的再入滑翔轨迹优化具有较大的挑战性和应用价值,而考虑真实情况的多参数不确定性的鲁棒轨迹优化问题则具有更大的难度和挑战性。因此,本节以高超声速滑翔问题为例,开展多不确定性条件下的鲁棒轨迹快速优化仿真,以验证本文所设计算法的有效性。此外,由于文献[11]已将基于PC的高斯伪谱法与凸优化方法进行了对比,证明了基于PC的凸优化方法在计算精度相当的情况下,计算效率显著提高。因此,本文不再将以上两种方法进行比较,而是将所提基于随机响应面和混沌多项式的凸优化方法(简称为SRSM-PC-CO)与基于GQ及混沌多项式的凸优化方法(简称为GQ-PC-CO)进行比较,以体现本文算法在计算效率方面的优势。

4.2.1 问题描述

假设地球为旋转圆球,标称条件下的高超声速滑翔质心运动模型采用文献[5]的形式。针对目标点固定情况下的飞行时间最短滑翔问题,若结合任务需求预先选定攻角剖面,而仅将倾侧角作为待优化控制量,可将该轨迹优化对应的确定性最优控制问题描述为

(42)

再入滑翔飞行所面临的不确定性主要包括再入初始状态,以及大气密度、气动力系数等动力学参数的不确定性,具体可建模描述为

(43)

且有

(44)

式中:c为动力学参数向量;ρ为大气密度;CL,CD为升力、阻力系数;δ为服从特定分布且相互独立的随机变量;上划线“-”代表标称情况。

为了尽可能降低终端位置散布,将终端经纬度的标准差引入鲁棒轨迹优化问题的目标函数,则考虑初始状态及动力学参数不确定的滑翔鲁棒轨迹优化问题为

(45)

式中:rminμ(e),rminσ(e)为地心距下边界的均值与方差,在算法迭代求解过程中,可基于上一步参考轨迹并结合式(24)计算得到。上述问题符合一般鲁棒轨迹优化问题PR的形式,可通过本文所设计的鲁棒轨迹优化算法进行求解。

同时,对于上述考虑不确定性的滑翔轨迹鲁棒优化问题,初始状态不确定性为5,过程参数不确定性为3,若同时考虑这些不确定性,则随机变量维数d=8。考虑该问题具有较强非线性,设置混沌多项式阶数p=2,若采用GQ-PC-CO方法进行鲁棒轨迹优化,则需将待优化状态量扩展为标称情况的38(6 561倍),原始随机常微分方程随之扩展为6 561×5维耦合的确定性微分方程组。如考虑凸优化过程中的离散化处理,待优化参数量将会更大,以上高维扩展会导致计算效率下降甚至无法收敛等计算困难;而采用本文所提出SRSM-PC-CO方法进行求解时,待优化状态只需扩展(8+2)!/(8!2!)(45倍),相较而言,显著降低了待优化变量个数,因而具有良好的优化计算效率和收敛性。

4.2.2 仿真条件设置

表3 边界条件取值Table 3 Boundary condition values

表4 参数不确定性设置Table 4 Parameters uncertainty setting

所有仿真均在搭载Intel Core i7-8700 3.20 GHz Intel处理器的台式机完成,仿真环境为Matlab 2016b平台。基于自行开发的代码进行不确定性量化传播,同时基于CVX工具包[29]进行轨迹优化算法开发,并调用SDPT3求解器求解凸优化子问题PED-disConvex。

4.2.3 不确定性量化传播仿真分析

首先,基于标称条件进行确定性滑翔轨迹优化仿真,所得结果如图2所示。由图2(a)、图2(b)可以看出,所得速度-高度曲线连续光滑并且满足过程约束要求,热流密度和动压曲线在一段时间内靠近约束边界,具有超出约束边界的潜在风险;由图2(a)、图2(c)可以看出,优化结果与积分结果一致,且控制指令反转次数较少,优化结果具有较好的可行性。

图2 标称情况下滑翔轨迹优化结果Fig.2 Optimization results of glide trajectory under nominal conditions

仅考虑表4中的大气密度、升力和阻力系数3个动力学参数不确定性,基于上述优化结果开展GQ-PC、SRSM-PC方法的不确定性量化传播仿真,并与MC打靶结果进行对比,以分析不确定性对优化轨迹的影响,同时验证所给出不确定性量化算法的有效性。设置PC多项式阶数p=2,MC打靶次数为3 000次,部分仿真结果如图3～图5及表5所示。

图3 标称情况下优化控制量MC仿真结果Fig.3 MC simulation result of optimized control quantity under nominal condition

图4 不确定性量化传播的状态均值对比Fig.4 Comparison of state mean using uncertainty propagation

图5 不确定性量化传播的状态标准差对比Fig.5 Comparison of state standard deviation using uncertainty propagation

表5 不确定量化传播的终端状态统计量对比Table 5 Comparison of terminal state statistics using uncertainty propagation algorithm

由图3可知,标称情况下的轨迹优化结果易受到参数不确定性的影响,从而出现过程约束超出限制、终端精度下降等问题,即轨迹可靠性及鲁棒性均面临不同程度的降低。因此,需要进一步开展不确定条件下的鲁棒轨迹优化,以提升轨迹的抗干扰能力及可靠性,进而降低制导控制系统负担和任务失败风险。

图4和图5分别给出了不同方法对应的部分状态均值及标准差比较曲线,表5则进一步对比了终端状态的统计精度。

可以看出,3种方法所得统计结果较为一致,相较于MC方法,GQ-PC与SRSM-PC方法所计算的终端高度均值与标准差误差小于0.1 km,经纬度误差小于0.01°,说明GQ-PC与SRSM-PC方法均可逼近打靶结果,即可替代打靶方法进行不确定性量化传播。在计算精度相当的单次量化传播过程中,SRSM-PC方法仅需进行20次动力学方程积分,而GQ-PC方法需要进行27次积分,其计算量为SRSM-PC的1.35倍;若进一步增加不确定性维数并考虑优化算法的迭代计算过程,SRSM-PC方法在计算效率方面的优势将会更加明显。

4.2.4 两种算法的鲁棒轨迹优化仿真对比

采用本文所设计的SRSM-PC-CO算法进行鲁棒轨迹优化仿真,同时与现有的GQ-PC-CO优化方法及确定性优化(deterministic optimization,DO)方法进行对比,以验证本文算法的有效性及计算效率优势。设置权重系数kJ=kg=3,其余仿真条件同上,结果对比如图6所示。随后,利用图6(d)中的倾侧角指令分别进行MC打靶仿真,得到部分参数的统计结果,如表6所示。其中,Δθμf和Δφμf分别表示终端经度、纬度均值偏差,θσf和φσf分别表示终端经度、纬度标准差;tμf和tσf分别表示飞行时间的均值及标准差,Pe表示轨迹超出过程约束边界的概率,即事故发生率。

图6 GQ-PC-CO与SRSM-PC-CO的鲁棒轨迹优化结果Fig.6 Robust trajectory optimization results of GQ-PC-CO and SRSM-PC-CO

由图6(a)和图6(b)可以看出,两组鲁棒轨迹优化结果基本一致,且高度、热流密度和动压曲线相较于DO更远离过程约束边界,因而约束超出概率降低,优化轨迹的可靠性得到了明显提升。

由表6可知,两组鲁棒优化控制量的打靶统计结果较为一致,相较于DO,其终端经、纬度分布更加集中,均值偏差及标准差均有所降低,过程约束超出概率显著降低。以上结果表明,本文方法与GQ-PC-CO的优化结果较为一致,其轨迹的鲁棒性和可靠性均显著提升,验证了本文方法的有效性。

表6 两种方法优化结果的MC仿真结果对比Table 6 Comparison of MC simulation results of optimization results for two methods

此外,由图6(d)可知,相较于DO,鲁棒优化使得倾侧角指令在初始下降段基本维持在0°附近,随后继续以相对较小的幅值飞行,以增大纵向平面内的法向爬升过载,进而提升前段轨迹高度,最终降低热流密度以遵守过程约束,而为了完成航向调整,在飞行中后期则需要采用相对更大的倾侧角幅值实现快速转弯,如图6(c)所示。

表7给出了两种算法的数值处理性能对比。可见,在相同条件下,本文算法对应的待优化变量数目显著降低,在迭代步数相当情况下的计算耗时降低了约62%,进一步说明了SRSM-PC-CO鲁棒轨迹优化算法在计算效率方面具有显著优势。

表7 两种优化算法的数值计算性能对比Table 7 Numerical computation performance comparison of two optimization algorithms

4.2.5 基于权重参数调整的轨迹优化可靠性与鲁棒性能分析

为了进一步探讨所提出轨迹优化算法在可靠性和鲁棒性方面的权衡能力,以更好地满足不同任务的需求侧重,下文开展基于权重参数调整的鲁棒优化仿真与分析。由前文鲁棒优化建模可知,调整权重参数kg、kJ可以满足不同的可靠性和鲁棒性需求。若令kg=3,kJ=0,则优化模型演化为基于可靠性的优化[10](reliability-based optimization, RBO);若令kg=kJ=3,则问题演化为鲁棒优化(robust optimization, RO)。

以表4中的3个动力学参数及5个初始状态作为不确定参数,采用SRSM-PC-CO算法进行仿真,同时将其与DO对比,以探究参数调整所带来的轨迹特性改变。表8给出了不同条件下的权重系数取值,其他仿真条件与前文一致。仿真结果如图7所示。

表8 权重系数设置Table 8 Weight factors setting

图7 DO、RBO与RO的结果参数对比Fig.7 Parameter comparison of results of DO, RBO, and RO

由图7(a)、图7(b)可知,相较于图3(b)的确定性优化结果,RO与RBO对应的高度、热流密度和动压曲线更远离约束边界,即优化轨迹的可靠性得到了明显提升;与RBO相比,RO对应的高度在第一个波峰处和末段更高,说明这种轨迹形式有利于降低随机不确定性的影响,提升终端位置和飞行时间指标的鲁棒性。由图7(c)可知,与前文分析类似,为了降低高度曲线在第一个波谷处的热流峰值,即提升轨迹可靠性,RBO与RO对应初始下降段的倾侧角幅值均较小,且大小基本一致,这表明RBO与RO对热流峰值的可靠性相当;在波谷之后,RO继续保持更小的倾侧角幅值,以保持更大的飞行高度;而在飞行后段,RO则采用了更大的倾侧角幅值,以实现航向快速调整,这一现象在图7(d)的地面航迹曲线结果中亦得到了进一步印证。事实上,轨迹可靠性和鲁棒性的提升对飞行器在不同飞行阶段的控制能力提出了更高要求,且没有统一的规律可循,需要针对不同的问题进行专门分析。

利用图7(c)中的倾侧角优化结果进行MC打靶仿真,得到部分参数的统计结果,如图8及表9所示,其中,图8给出了打靶结果中部分参数的近似概率密度函数。

图8 DO、RBO与RO的MC仿真结果概率密度函数Fig.8 Probability density functions of DO, RBO, and RO’s MC simulation result

表9 DO、RBO与RO的MC仿真数据对比Table 9 Data comparison of MC simulation of DO, RBO and RO

由图8(a)、图8(b)及表9可知,RBO与RO的终端经、纬度统计精度基本一致,且相较于DO有明显改善,而由于考虑了终端性能和优化指标的鲁棒性,RO对应的终端经、纬度和飞行时间散布在三者中最小。由图8(c)～图8(e)及表9可知,相较于DO结果,RBO与RO的热流、动压、过载峰值分布均远离约束边界,约束超出的概率由DO的67.4%分别降低到0.90%和0.88%,可见两者的可靠性均有所提升。以上结果表明,与RBO结果相比,RO在保证轨迹可靠性的同时,对不确定性的敏感度最低,具有更好的鲁棒性。由图8(f)及表9可知,相较于DO结果,RBO与RO对应的飞行时间均值较大,这是由于不确定性的引入和可靠性约束的施加会造成优化求解的可行域减小,迫使轨迹高度提升,进而导致飞行时间增大。而相较于RBO,RO的飞行时间均值更大,这是由于RO的目标函数综合考虑了飞行时间的均值及标准差,而RBO的目标函数仅考虑了飞行时间的均值。因此,轨迹鲁棒性与可靠性的提升是以牺牲一定最优性为代价的。在实际应用时,可根据设计需求调整权重系数,以满足轨迹的可靠性、鲁棒性和最优性综合需求。

表10对比了3种仿真条件下的计算耗时及待优化变量数。由表10可以看出,即使同时考虑初始状态及参数不确定性,RBO及RO的优化耗时均在可接受范围内,若采用并行计算,优化耗时可进一步缩短。

表10 DO、RBO与RO的数值计算性能对比Table 10 Numerical computation performance comparison of DO, RBO and RO

5 结 论

针对存在多参数不确定性的轨迹优化问题,本文研究了一种通用、高效的鲁棒轨迹优化方法,并以高超声速再入滑翔轨迹优化为例进行了仿真分析与对比验证,相关研究结论如下:

(1) 与确定性轨迹优化和传统鲁棒轨迹优化算法相比,所设计鲁棒轨迹优化方法可以获得兼具可靠性和鲁棒性的优化轨迹,且具有相当的精度和显著的计算效率优势。

(2) 轨迹鲁棒性和可靠性的提升是以牺牲一定最优性为代价的,在实际应用时,可根据设计需求调整权重系数,以综合满足轨迹的可靠性、鲁棒性和最优性需求。

(3) 所构建的基于随机响应面法及NIPC的不确定性量化传播模型,可实现初始状态和过程动力学等多参数不确定随机问题的有效转化,且具有状态扩展维数少、计算效率高的优势。

(4) 与基于IPC的优化方法相比,本文方法将原系统模型视为“黑箱”,无需针对具体问题进行推导建模,因而具有较高的工程应用价值。