不确定条件下采用协方差描述函数法的再入轨迹鲁棒优化

刘 旭,李 响,张后军,郭宇恒,王晓鹏

(1.北京理工大学宇航学院,北京 100081；2.北京机电工程研究所,北京 100083；3.上海机电工程研究所,上海 200233)

0 引 言

传统的再入轨迹优化不考虑不确定性因素,即默认优化模型中所有参数都是准确的[1-3],而实际飞行过程中不确定性因素客观存在,如飞行器气动参数误差、大气参数误差等,这些误差的影响一般通过后续的控制系统设计消除。近年来,有学者在飞行轨迹优化阶段就将这些不确定性因素考虑进去,提出了飞行轨迹的鲁棒优化(Robust design optimization,RDO)[4]与可靠性优化(Reliability-based design optimization,RBDO)[5]。RDO在最小化目标函数期望的同时,尽量使目标函数对不确定性因素不敏感；RBDO着重于在不确定性因素影响下将失效概率限制在某个值以内。通过RDO或RBDO优化后得到的飞行轨迹对不确定性因素的敏感度降低,自身就有一定抗干扰的能力。

文献[4]指出在考虑不确定性的情况下,相较于飞行器静态鲁棒优化(如气动优化、结构优化等),飞行器动态鲁棒优化(如飞行轨迹优化)的研究比较有限,因此提出一种不确定条件下的飞行轨迹鲁棒优化方法；文献[6]指出飞机轨迹规划和预测中的不确定性给未来空中交通管理系统带来了重大挑战,因此提出一种考虑气象不确定性的基于最优控制的飞机轨迹规划方法；文献[7]注意到不确定性(尤其是风场不确定性)对于无人机动态滑翔轨迹的影响,在量化不确定性后提出一种动态滑翔轨迹鲁棒优化方法；火星大气密度具有的强不确定性是火星进入段飞行任务的主要干扰源,对此文献[8]提出一种火星大气进入段的轨迹鲁棒优化,得到的最优轨迹在不确定条件下的鲁棒性明显优于确定性优化,而文献[5]则是提出用序列可靠性优化(Reliability-based sequential optimization,RBSO)来降低参数不确定性对火星进入段轨迹的影响,同样得到优于确定性优化的结果；文献[9]应用RBDO优化飞行器在不规则小行星软着陆的轨迹,将不确定动态模型的轨迹优化转化为带有可靠性约束的参数优化,并应用序列优化和可靠性评估(Sequential optimization and reliability assessment,SORA)流程求解；文献[10]利用脱敏最优控制(Desensitized optimal control,DOC),将标称轨迹设计和标称轨迹跟踪结合在一起,得到的火星下降段脱敏轨迹对不确定性和扰动的敏感性大幅降低。

在飞行轨迹鲁棒优化或可靠性优化的实际计算中,不确定性传播(Uncertainty propagation,UP),即如何根据输入不确定量计算得到输出不确定量的统计特征(一般指均值和方差)是关键步骤。传统蒙特卡洛(Monte Carlo,MC)仿真法易于操作,对各种复杂问题均有很强的适应性,得到了较广泛的应用,也常用于检验其他方法的准确性。但MC仿真法的精度依赖于大量样本点上的仿真次数,为得到可靠结果需要的高昂计算代价限制了它在轨迹鲁棒优化中的应用。混沌多项式展开(Polynomial chaos expansion,PCE)是一种基于Askey正交多项式随机展开的不确定性传播分析法,能高效地计算输出不确定性变量的统计特征[11],近年来在飞行轨迹鲁棒优化中得到了较广泛的应用[4,6-8]。基于可靠性的飞行轨迹优化[5,9]多采用最可能失效点(Most probable point,MPP)法进行不确定性分析,包括一次可靠度法(First-order reliability method,FORM)[12]和二次可靠度法(Second-order reliability method,SORM)[13],其基本思路是在MPP点进行局部展开来近似极限状态函数,并以此计算失效概率。文献[14]引入人工神经网络模型,用训练后的神经网络作为复杂系统的替代模型,改善单独使用MC法的不足。文献[15] 通过Kriging代理模型来构建目标变量和不确定性输入参数的响应关系,以此来分析目标变量的不确定性信息。

除了前述PCE方法与基于MPP的方法,协方差分析描述函数技术(Covariance analysis describing function technique,CADET)是另外一种高效的随机系统统计特征计算方法。该方法最早由Gelb提出[16],先采用统计线性化方法将非线性系统转化为近似线性系统,再用协方差分析法导出系统随机状态量均值和协方差的传播方程,只需一次求解就能确定特定随机系统的统计特征,在飞行器制导系统性能统计分析中得到了广泛应用。文献[17]应用CADET进行单发武器的射击效能分析；文献[18] 基于CADET建立了制导炸弹半实物仿真系统的误差传播理论,分析了仿真姿态角和姿态角速度的误差对仿真结果的影响；文献[19]基于闭环控制系统下线性化的相对运动动力学模型,采用CADET对定义航天器误差椭球的协方差矩阵进行分析；文献[20] 针对火星探测捕获的误差传播,利用CADET分析单一误差和综合误差对不同捕获策略入轨精度的影响。现有的文献中,CADET主要用于制导精度分析,本研究将CADET引入再入轨迹鲁棒优化中,实现了再入轨迹鲁棒优化的高效计算。

1 飞行轨迹鲁棒优化模型

不失一般性,先以传统确定性条件下Bolza型飞行轨迹优化问题构建飞行轨迹鲁棒优化模型,可以描述为：将u(t)作为控制量,使得

(1a)

(1b)

C[x(t),u(t),t]≤0

(1c)

H[x(t0),u(t0),x(tf),u(tf),tf]=0

(1d)

式中：x(t)是状态向量；u(t)是控制向量；f是动态约束向量函数；C是过程约束向量函数；H是端点约束向量函数。

现在考虑不确定性因素的影响,将随机量s引入式(1),得到：将u(t)作为控制量,使得

u(t),t]dt

(2a)

(2b)

C[x(t,s),u(t),t]≤0

(2c)

H[x(t0,s),u(t0),x(tf,s),u(tf),tf]=0(2d)

由于引入了随机量s,式(1)中的状态量x(t)由确定量转变为随机量x(t,s)。相应地,式(1)中确定性的目标函数式(1a)、动力学方程式(1b)、过程约束式(1c)和端点约束式(1d)分别转变为随机条件下的式(2a)、(2b)、(2c)和(2d)。

由于含有随机量s,式(2)中的目标函数式(2a)和约束条件式(2c-2d)与传统轨迹优化中的目标函数式(1a)和约束条件式(1c-1d)有本质的不同。现有的轨迹优化方法都是针对如式(1)所示的确定性问题发展起来的,不能直接处理式(2)所示的含有随机量的优化问题。为此,引入飞行器气动和结构优化中鲁棒优化的理念[21-22],将传统确定条件下飞行轨迹优化模型转化为如式(3)所示不确定条件下的飞行轨迹鲁棒优化模型：将u(t)作为控制量,使得

min μ(J)+kJ·σ(J)

(3a)

(3b)

μ(C[x(t,s),u(t),t])+kC·σ(C[x(t,s),

u(t),t])≤0

(3c)

μ(H[x(t0,s),u(t0),x(tf,s),u(tf),tf])=0

(3d)

σ(H[x(t0,s),u(t0),x(tf,s),u(tf),tf])≤εH

(3e)

式中：μ(·)和σ(·)分别表示随机量的均值和标准差；kJ和kC是权重系数；εH是随机端点约束标准差的容许值。

式(3)所表达的具体含义如下：

式(3a)与目标函数式(1a)对应,是随机目标函数的均值μ(J)和标准差σ(J)的加权和,通过最小化随机目标函数的均值来表达确定性目标函数式(1a)的设计意图,同时通过最小化随机目标函数的标准差来限制它的离散程度。目标函数的最优性和鲁棒性可以通过选取权重系数达到设计要求,理论上kJ值越大,目标函数的离散程度越小。

式(3c)与过程约束式(1c)对应,式(1c)是不等式约束,在鲁棒优化中,不等式约束通过随机约束的均值μ(C)和标准差σ(C)的加权和来表达,kC越大,所得最优轨迹的鲁棒性越高,满足过程约束的概率越大,但求解难度也会相应增加。

式(3d)、(3e)与端点约束式(1d)对应,式(1d)是等式约束,在鲁棒优化中,等式约束通过随机约束的均值式(3d)和标准差式(3e)两组约束来表达,这与不等式约束是不同的。均值约束要求等于原问题中的值,而标准差约束要求小于εH,εH越小,鲁棒性越高。

总体来看,鲁棒优化模型式(3)通过要求随机优化模型式(2)中随机量的均值满足确定性优化模型式(1)中相应的约束,来表达原有的设计意图；同时通过约束随机优化式(2)中随机量的标准差,限制其离散程度,提高最优轨迹的鲁棒性。这样从统计角度来看,式(3)得到的最优轨迹在随机扰动的影响下偏离设计值的程度要小于式(1),从而达到鲁棒优化设计的目标。

值得注意的是,鲁棒优化模型中式(3b)目前仍是含有随机量s的动力学方程,现有轨迹优化方法(如伪谱法等)难以直接处理这类问题,在下节中将讨论如何利用协方差分析描述函数技术(CADET)将含有随机量的鲁棒优化问题式(3)转化为等效的确定性优化问题,从而可以采用伪谱法等现有轨迹优化方法求解。

2 采用CADET重构轨迹鲁棒优化模型

对于式(3)所示的鲁棒优化模型,如何高效计算式(3a,3c-3e)中的统计特征(均值μ(·)和标准差σ(·))是关键步骤。蒙特卡洛法程序结构简单,易于实现,得到了广泛应用。根据大数定律,样本数N→∞时,蒙特卡洛法的结果以概率1收敛于真实值。一般认为,样本点数量达到104时,蒙特卡洛法的结果才是稳定收敛和可靠的[23],工程上通常要成千上万样本点才能得到满足精度要求的结果。对单次随机过程的统计特征分析而言,蒙特卡洛法是可以接受的；而轨迹优化计算是一个迭代过程,在每一迭代步都需要计算均值μ(·)和标准差σ(·)的值,在这种情况下,蒙特卡洛法将导致巨大的计算量,在实际计算中是不可行的。根据上述情况,本文采用CADET计算随机系统的统计特征,并将含有随机量的式(3)转化为等效的确定性优化模型。

2.1 CADET方法

一般连续时域非线性随机系统可以表示为：

(4)

式中：x(t,s)是n维随机状态向量；w(t)是m维强迫函数向量,包含作用于系统的随机扰动s和控制输入u；f(x,t)是x(t,s)的非线性矢量函数；G(t)是n×m的连续函数矩阵。

协方差分析描述函数法,就是先利用描述函数理论对非线性系统进行统计线性化,然后对线性化后的系统进行协方差分析。所谓统计线性化,是根据状态量x(t,s)的概率密度函数形式,在较大的x(t,s)变化范围内,用一个线性函数来逼近非线性函数f(x,t)。其本质是对系统的每一个非线性用一个或几个近似的、对输入幅度灵敏的线性增益来代替[24]。这种方法的优点是不要求f(x,t)连续可微,使许多不能用真线性化进行线性化处理的实际函数有了线性化的可能[25]。

(5)

(6)

(7)

式中：b(t)是w(t)的确定性分量,Q(t)是w(t)随机分量的谱密度矩阵。只需确定初始条件m(0)和p(0),利用式(7)进行一次计算就能获得任意时刻系统状态量的统计特征,而且精度与几千次的蒙特卡洛仿真相当,极大程度上节省了计算时间与成本。

由式(6)可知,描述函数的计算依赖于状态量的概率密度函数,一般来说十分复杂,甚至不能求得解析式。工程上常假定x(t,s)服从联合正态分布,从而描述函数的计算可简化为

(8)

这里假设随机状态量x(t,s)服从正态分布,只有对高斯输入的线性系统才是严格成立的,而对非高斯输入的线性和非线性系统往往是近似的。对于非线性系统,尽管输入是高斯的,输出一般来说是非高斯的。由中心极限定理可知：当经过低通滤波器后,随机过程趋于高斯过程。因此,当信号通过系统传递时,可以借助系统的线性部分来保证非高斯的非线性输出结果近似是正态的[24]。前述高斯假设的有效性已经得到了广泛的研究和检验。

在附录A中,给出了采用CADET方法得到的再入轨迹非线性动力学模型的统计线性化模型。

2.2 轨迹鲁棒优化的等效确定性优化模型

式(3)所示的轨迹鲁棒优化模型尽管表达了均值、标准差等引入随机量后的性能信息,但式(3b)仍是关于x(t,s)的随机微分方程,含有随机量s,且其余各式中关于均值和标准差的求解方法也没有给出,式(3)是无法用于实际计算的。现在应用CADET,将式(3)中含有的随机变量s消去,得到等效的确定性优化模型,具体步骤为将式(3b)转化为关于x(t,s)均值和协方差的常微分方程,同时将式(3)中目标函数、过程约束和端点约束表示为关于x(t,s)均值和协方差的函数,从而得到等效的确定性轨迹鲁棒优化模型：将u(t)作为控制量,使得

minJμ(m(t),p(t),u(t),t)+kJ·Jσ(m(t),

p(t),u(t),t)

(9a)

(9b)

(9c)

Cμ(m(t),p(t),u(t),t)+kC·Cσ(m(t),

p(t),u(t),t)≤0

(9d)

Hμ(m(t),p(t),u(t),t)=0

(9e)

Hσ(m(t),p(t),u(t),t)≤εH

(9f)

式中：m(t)和p(t)分别是随机状态量x(t,s)的均值和协方差,根据CADET中的式(7)计算；(·)μ和(·)σ分别表示用于计算(·)的均值和标准差的函数。

对比式(9)和式(3),可以发现,采用CADET后,含有随机量s的鲁棒优化模型式(3)转化为不含随机变量的等效确定性优化模型式(9),应用现有的轨迹优化算法完全可以求解。换言之,通过求解确定性优化问题式(9),就能得到不确定条件下,鲁棒优化意义下的最优控制量和最优轨迹。

3 算例与分析

3.1 确定性再入轨迹优化模型

文献[26]中的再入轨迹优化问题是一例经典的轨迹优化问题,希望在满足终端和过程约束的条件下,最大化终点的横向距离(即纬度),在确定性条件下可以描述为：将α(t)和σ(t)作为控制量,使得

maxJ=φ(tf)

(10a)

(10b)

(10c)

(10d)

式中：r是到地球球心的距离；θ是经度；φ是纬度；v是速度；γ是航迹角；ψ是航向角；α是攻角；σ是倾侧角；m是质量；μ是地心引力常数；Re是地球半径；D是阻力加速度；L是升力加速度；q是热流密度。式(10c)表示端点约束,式(10d)表示过程约束。关于该模型中各参数的详细说明可参见文献[26],这里不再详述。

进一步,阻力加速度和升力加速度表示为：

(11)

式中：CD和CL分别是阻力系数和升力系数；ρ是大气密度S是特征面积；h是高度；H是密度尺度高；ρ0是海平面大气密度。

在传统飞行轨迹优化中,气动参数CD·和CL·视为确定值,根据实验得到的气动数据求得。

3.2 再入轨迹鲁棒优化模型

由于制造工艺误差、测量误差等各方面因素,将气动参数CD·和CL·视为不确定性量更为合理。本文中假设：

(12)

式中：(·)N表示(·)的参考值；s1,s2表示随机变量。本文仿真实验中,s1～N(0.01,0.012),s2～N(-0.01,0.012)且s1,s2相互独立。需要说明的是,这里s1,s2的均值不为0,目的是使干扰条件更苛刻,以显示方法的有效性。关于s1,s2的概率分布规律等因素本身也是值得研究的问题,本研究着眼于再入轨迹鲁棒优化方法,因此在这里对此不作过多讨论。

引入随机变量s1,s2后,确定性轨迹优化问题式(10)转化为如式(2)所示的随机轨迹优化问题,在此基础上进一步转化为式(3)所示的轨迹鲁棒优化问题,具体如下：

式(10)中的目标函数式(10a)变为：

maxJ=μ(φ(tf,s))-kJ·σ(φ(tf,s))

(13)

动力学方程式(10b)转化为含随机变量的动力学方程：

(14)

式中：x(t,s)表示随机状态向量；f表示式(10b)的右端项函数。

相应地,式(10)中的确定性端点约束式(10c)和过程约束式(10d)转化为用均值和标准差描述的统计意义上的约束。其中,等式约束式(10c)转化为由均值和标准差表示的两组约束式(15a)和式(15b),不等式约束(10d)转化为式(15c)：

(15a)

(15b)

(15c)

需要说明：式(15b)仅对原确定性问题中有终端约束的状态变量施加了标准差约束；所考虑的气动不确定性不会影响到初始状态,因此式(15b)中没有设置对初值标准差的约束；此外,用来限制离散程度的标准差容许值和权重系数由设计者根据实际情况设定,本研究中设置如下：εhf=4500 m,εvf=250 m/s,εγf=3°,kJ=3,kCr=3,kCθ=3,kCv=3,kCγ=3,kCq=3。

综上,再入轨迹鲁棒优化模型(RO)表示为：将α(t)和σ(t)作为控制量,使得

(16)

至此,已构建起2个再入轨迹优化问题,即传统确定性优化问题DO(式(10))以及同时考虑均值和标准差约束的鲁棒优化问题RO(式(16))。在下节中,将分别求解这两个优化问题,并对结果进行对比分析。

3.3 计算结果及讨论

鲁棒优化模型式(16)含有式(14),而式(14)是含有随机变量s的随机飞行动力学方程,在实际计算求解式(16)时,需要采用CADET方法将式(14)转化为形如式(7)的关于均值和协方差的确定性动力学方程,具体的转化结果可参考附录A。

3.3.1最优控制量的对比分析

采用伪谱法求解这两个轨迹优化问题(关于伪谱法的相关论述这里不再展开,可参看相关文献[27-29]),得到传统优化DO的最优控制量：αDO(t)和σDO(t),鲁棒优化RO的最优控制量：αRO(t)和σRO(t),如图1所示,为了更清楚地显示,两个小方框给出了相应时间段的放大图。

图1 最优控制量(方框内为放大图)Fig.1 Optimal control

观察图1中控制量的时间历程曲线,可得如下信息：

鲁棒优化最优控制量αRO(t)、σRO(t),分别与对应的确定性优化的最优控制量αDO(t)、σDO(t)的时间历程趋势基本一致,但αDO(t)和σDO(t)相对平缓,而αRO(t)和σRO(t)的波动较大,且越接近终端时刻波动越明显。

αRO(t)和σRO(t)所需的时间比αDO(t)和σDO(t)的更长(RO为2289.4 s,DO为2200.1 s)。

上述现象表明,鲁棒优化为了降低不确定性,比传统确定性优化付出了更多的控制成本,这是符合客观实际的。此外,需要说明的是,图1显示攻角在末端变化较快,这是图形比例原因导致的(总的飞行时间较长),攻角实际最大变化率绝对值为3.43°/s,工程上是完全可以实现的。

3.3.2目标和约束的鲁棒性对比分析

为验证两种优化模型得到的最优控制量在不确定因素影响下的效果,将αD0(t)和σD0(t)、αRO(t)和σRO(t)分别代入式(14)所示的含随机变量的动力学方程中,进行10000次的MC仿真,并对仿真结果进行对比分析,主要从三个方面进行：目标函数的鲁棒性对比；终端约束的鲁棒性对比；过程约束的鲁棒性对比。

首先给出相应的统计结果：

目标函数终端时刻纬度的统计特性及散布情况如表1和图2所示；终端约束终端时刻高度、速度和航迹角的统计特性及散布情况如表2、表3和图3所示；过程约束热流密度约束的情况如图4和表4所示。图2和图3中横坐标表示MC仿真次数,对图表的解释及分析见后。

表1 终端时刻纬度的统计特性Table 1 Statistical property of latitude at terminal time

图2 终端时刻纬度的散布情况Fig.2 Dispersion of latitude at terminal time

表2 终端时刻高度、速度和航迹角的均值和相对误差Table 2 Mean and relative error of height,speed and flight path angle at terminal time

表3 终端时刻高度、速度和航迹角的标准差Table 3 Standard deviation of height,speed and flight path angle at terminal time

图3 终端时刻高度、速度和航迹角的散布情况Fig.3 Dispersion of height,speed and flight path angle at terminal time

图4 热流密度的散布情况Fig.4 Dispersion of heat flux density

通过分析上述图2到图4中的结果及表1到表4中的数据,可以得到以下结论：

表4 热流密度约束的满足情况Table 4 Satisfaction of heat flux density constraint

1)目标函数结果对比

表1中第一、二行数据显示,对于目标函数φ(tf)的均值,DO大于RO,而对于φ(tf)的标准差,RO小于DO。

φ(tf)的均值越大而标准差越小越好,因为极大化φ(tf)是设计目标,而较小的标准差则意味着抗干扰能力较强。因此上述现象表明,RO为提高系统对不确定性因素的鲁棒性,牺牲掉一部分φ(tf)的均值,这是符合设计预期的。

为了显示出均值和标准差的综合效果,表1中第三行数据显示了均值和标准差的加权和,kJ取3时,RO表现优于DO。

图2显示了气动不确定条件下,终端时刻纬度的散布情况,可以直观地看出RO中的样本点分布更为紧密,表明RO的目标函数对不确定性因素的敏感程度更低。

2)终端约束结果对比

表2中数据显示,在气动不确定因素影响下,DO的终端高度、速度和航迹角的均值与约束要求的值有较大偏差,最大达到51.40%。而RO的终端状态均值与终端约束要求的值基本吻合。这说明,在气动不确定因素影响下RO能使终端状态以约束要求值为中心散布,而DO的终端状态散布偏离了约束要求值。

表3中数据显示,RO中的高度和航迹角的标准差小于DO中对应量的标准差,这说明,在气动不确定影响下,RO中的高度和航迹角的分布更紧密。此外,尽管DO中速度的标准差小于RO中速度的标准差,但由于DO中速度的均值偏离了要求的均值,较小的标准差并未带来实际的收益,这在图3中得到了体现,具体如下。

图3表明了气动不确定条件下,高度、速度和航迹角终端状态的散布情况,红实线是设计要求,表示均值,上、下红虚线表示均值加减标准差,则落入上、下虚线之间部分的样本点越多,说明抗干扰能力越强,鲁棒性越好。DO的10000个样本点中分别有5937(高度)、6599(速度)和6310(航迹角)个在这个范围内,占比为59.37%、65.99%和63.10%；而RO的10000个样本点中,分别有6883(高度)、6766(速度)和8267(航迹角)个在这个范围内,占比为68.83%、67.66%和82.67%。显然,RO所得到的结果比DO的鲁棒性更好。

3)过程约束结果对比

本算例中共有7个过程约束,其中热流密度的过程约束是主动约束,因此这里以热流密度过程约束来讨论不确定因素的影响。

图4a和图4b分别显示了DO和RO中q(t)时间历程样本的分布情况,红线代表飞行过程中热流密度的约束上限。对比图4a和图4b,直观地显示出DO中大量热流密度时间历程样本超出了约束限制,而在RO中,绝大多数样本都是满足约束的。表4中给出了具体值,在DO的10000条q(t)时间历程样本中,有17条完全满足约束,占比为0.17%；而在RO的10000条q(t)时间历程样本中,有9933条完全满足约束,占比为99.33%。上述结果表明,在气动不确定条件下,RO中的热流密度过程约束有极强的抗干扰能力,而DO中的热流密度过程约束的抗干扰能力很弱。

总的来说,在不确定性因素的影响下,再入轨迹鲁棒优化的最优轨迹相较于确定性优化的最优轨迹表现出更强的抗干扰能力,而对于不同的物理性能指标,抗干扰的能力有所不同,在本算例中,对热流密度的抗干扰能力最强。

4 总结与展望

传统的再入轨迹优化是在确定性条件下展开的,本研究将不确定性参数引入再入轨迹优化中,构建了再入轨迹鲁棒优化模型,使得到的飞行轨迹具有一定的抗干扰能力。将CADET作为不确定性分析手段,引入再入轨迹鲁棒优化的实际计算中,实现了再入轨迹鲁棒优化的高效计算。仿真结果表明,相比于确定性轨迹优化,本文方法得到的最优轨迹在不确定性因素的干扰下具有更强的鲁棒性。

本研究的着眼点是再入轨迹鲁棒优化方法,而轨迹优化中的不确定性来源众多,将来还应对输入不确定性的合理建模展开深入研究。此外,协方差分析法的精度会随着系统非线性程度的增加而降低,对于非高斯输入的有效性也有待考察,因此准确高效的不确定性传播方法依然是未来关注的重点。