一种基于SRSM和NIPC的飞行器鲁棒轨迹优化方法

一种基于srsm和nipc的飞行器鲁棒轨迹优化方法
技术领域
1.本发明属于飞行器技术领域，具体涉及一种基于srsm和nipc的飞行器鲁棒轨迹优化方法。


背景技术：

2.作为飞行器设计的关键技术之一，轨迹优化逐渐成为高超声速等航空航天领域研究的热点问题。然而，现有的轨迹优化研究通常仅注重标称情况下的轨迹性能优化提升，而没有考虑真实情况下轨迹的鲁棒性和可靠性。在实际飞行任务中，初始状态、动力学及环境参数等均存在诸多不确定性，势必使得沿标称轨迹飞行时产生较大的过程及终端偏差，进而增大约束违反的风险，降低飞行可靠性和终端任务精度。因此，为了提高优化轨迹的可靠性和抗干扰能力，部分学者逐步将研究重心聚焦于考虑不确定性的鲁棒轨迹优化，即在轨迹设计阶段事先考虑不确定性的影响，从而降低轨迹规划迭代耗时和轨迹偏差修正所带来的制导控制系统负担。
3.求解鲁棒轨迹优化问题的典型策略为：先将含有随机微分方程的随机最优控制问题转化为包含扩展确定性微分方程的确定最优控制问题，进而采用高效的轨迹优化算法对该高维问题进行数。值求解，即包含不确定性量化传播和确定性轨迹优化两部分关键技术。典型的不确定性量化传播理论和方法包括蒙特卡洛方法(monte carlo，mc)、线性协方差分析法、无迹变换法及混沌多项式展开法(polynomial chaos，pc)等。蒙特卡洛方法因其简单易行而被广泛应用于不确定性量化及传播，但大量的数值模拟会导致统计结果收敛速度较慢；线性协方差分析法可对线性高斯过程模型中的不确定性进行量化传播，且具有较高的计算效率，但其不适用于高度非线性的轨迹优化问题；相比之下，pce将随机变量表示为随机正交多项式的加权求和，可对具有任意分布类型的随机变量实现较为精确的逼近，同时能够以较低的计算代价达到与mc方法相当的精度，近年来在火星进入轨迹不确定性量化、流体力学不确定性分析中得到了广泛应用。pce方法与系统模型结合的方式有两种，即嵌入式与非嵌入式，嵌入式混沌多项式(intrusive polynomial chaos，ipc)将随机变量的混沌多项式展开代入动力学微分方程中，并利用伽辽金投影法将随机微分方程转换为一组以混沌多项式系数为状态量的高维确定性微分方程组，最终通过求解该微分方程组完成不确定性传播，具有较高的计算效率；非嵌入式混沌多项式(non-intrusive polynomial chaos，nipc)则将原系统模型视为“黑箱”，进而根据模型的输入输出逼近混沌多项式系数，因此能够处理一般非线性系统的不确定性量化传播问题，且具有较好的通用性及工程应用价值。
4.而对于另一关键技术即确定性轨迹优化而言，常见的方法主要包括间接法和直接法。间接法因对解的初值猜测较敏感，需要对最优性必要条件进行推导且无法处理复杂约束等原因，难以在复杂约束轨迹优化问题中得到广泛应用。近年来，以伪谱法、凸优化等为代表的直接法已被广泛应用于高超声速等航空航天领域轨迹优化任务，并展现了良好的数值计算优势和普适性等特点。其中，凸优化方法因具有多项式复杂度、收敛速度快，且其局部最优解理论上即为全局最优解等特点，已逐步被用于求解复杂约束下的轨迹优化问题以
及考虑不确定性的高维状态扩展优化问题。
5.鉴于pce与直接法的优势，诸多研究将二者进行有机结合，设计了多种基于pce与直接法的鲁棒轨迹优化算法和求解策略。fisher等首次提出了数值求解随机最优控制问题的框架，该框架利用ipc将随机最优控制问题转化为高维确定性最优控制问题，并采用配点法对该问题进行求解。li等将nipc与高斯伪谱法相结合，用于解决气动系数不确定的高超声速飞机的爬升鲁棒轨迹优化问题。jiang等将nipc与hp伪谱法相结合，通过调用opencossan工具箱进行不确定量化传播，并采用并行计算技术提升高维确定性优化问题的求解效率，解决了初始状态及模型参数均不确定的火星再入鲁棒轨迹优化设计问题。wang等提出了一种ipc与凸优化结合的鲁棒轨迹优化方法，并与基于ipc与高斯伪谱法的方法进行对比，证明了所提方法在计算精度相当的情况下，计算效率显著提高。杨奔等基于最小二乘法建立了nipc与凸优化相结合的再入轨迹鲁棒优化模型，提升了气动参数不确定条件下轨迹的可靠性。
6.基于pce和直接法的随机问题求解技术虽然已被逐步用于鲁棒轨迹优化，但仍存在计算效率低、通用性较差以及难于处理多维不确定性等问题，因而限制了其在复杂约束及强非线性鲁棒优化问题(如航空航天领域鲁棒轨迹优化)中的应用。首先，以高斯伪谱法为代表的传统直接法虽然解决了部分鲁棒轨迹优化问题，但其仍面临计算效率低、难于快速求解多维不确定问题等局限性，而凸优化技术与pce的结合仍然不够完善。其次，基于ipc的优化方法虽然具有较高的计算效率，但需要针对不同问题进行推导建模和代码改写，由于推导过程存在模型简化近似，故其通用性较差，且只适用于处理多项式非线性问题；基于nipc的优化方法虽然能够避免上述问题，但pce系数求解大都基于求积、采样的方法，例如高斯求积方法(gauss quadrature，gq)、稀疏网格求积法，而随着系统随机变量维数及非线性程度增大，这类方法的状态扩展维数会呈指数增长，使得状态离散节点数量显著增加，从而导致解决大规模nlp问题时非常困难，计算成本高昂且效率低下，尤其是当不确定因素维数n＞3时，大量的内存消耗和计算成本会严重阻碍优化过程，甚至导致优化失败。
7.值得注意的是，除去上述基于求积的方法，nipc系数也可采用随机响应面法(stochastic response surface method，srsm)进行计算。该方法基于代理模型思想，通过拟合系统输出与输入之间的关系获得混沌多项式模型，并结合典型抽样方法进行少量采样即可实现多项式系数计算，可以显著降低模型转化维数和规模，因此具有较高的精度和计算效率。目前，srsm已广泛应用于流体力学不确定性分析中，但尚未在鲁棒轨迹优化设计中广泛应用。


技术实现要素：

8.为了解决上述技术问题，本发明提出了一种基于srsm和nipc的飞行器鲁棒轨迹优化方法。
9.为了实现上述目的，本发明提供如下技术方案：
10.一种基于srsm和nipc的飞行器鲁棒轨迹优化方法，包括以下步骤：
11.构建参数不确定条件下的飞行器鲁棒轨迹优化问题pr；
12.基于非嵌入式混沌多项式nipc与随机响应面法srsm构建不确定性量化传播模型，将飞行器鲁棒轨迹优化问题pr转化为高维确定性最优控制问题p
ed
；
13.利用逐次线性化方法将高维确定性最优控制问题p
ed
凸化为凸优化问题p
ed-convex
；
14.将凸优化问题p
ed-convex
进行离散化处理，得到参数凸优化问题p
ed-disconvex
，并构建基于序列凸优化的优化求解算法，通过迭代求解凸优化问题p
ed-disconvex
逼近飞行器鲁棒轨迹优化问题pr的解；
15.将飞行器的总体参数、约束条件参数、不确定参数及权重系数作为参数输入，通过调用基于序列凸优化的优化求解算法进行优化求解，得到最优控制量，完成飞行器轨迹优化。
16.优选地，所述飞行器鲁棒轨迹优化问题pr的模型为：
[0017][0018]
式(1)中，x(t)为状态向量，u(t)为控制向量，w为系统随机参数向量；u
min
与u
max
为控制量的上下边界；j为目标函数，f为动力学函数，g为过程约束函数，下标μ和σ代表对应变量的均值和标准差，下标0和f分别表示初始与终端状态；kj、kg为待设计权重系数，x
0μ
与x
0σ
分别为初始状态的均值和标准差，εf为状态终值标准差的上界，为状态向量的随机微分向量。
[0019]
优选地，所述基于非嵌入式混沌多项式nipc与随机响应面法srsm构建不确定性量化传播模型，将飞行器鲁棒轨迹优化问题pr转化为高维确定性最优控制问题p
ed
，具体为：
[0020]
通过引入拉丁超立方采样策略，并基于随机响应面法srsm将含不确定参数的状态向量的随机微分向量扩展为关于pc多项式系数的高维常微分方程，并通过pc多项式系数计算问题pr中状态、约束、目标函数的均值及标准差，进而将飞行器鲁棒轨迹优化问题pr转化为高维确定性最优控制问题p
ed
。
[0021]
优选地，所述将飞行器鲁棒轨迹优化问题pr转化为高维确定性最优控制问题p
ed
，具体包括以下步骤：
[0022]
(1)状态向量的随机微分向量的pc展开
[0023]
假设随机状态向量x为n维，系统参数向量w为s维，则每个状态变量xi(t)及参数wk展开为pc多项式：
[0024]
[0025][0026]
其中，δ＝{δ1,
…
,δd}为d维随机向量，φj(δ)为正交多项式函数，通过与δ分布类型相对应的一元正交基函数φj(δ
l
),l＝1,...,d求张量积得到；x
ij
(t)和w
kj
为pc系数，展开项数p+1＝[(p+d)！]/(p！d！)，即由正交基函数φj(δ)的阶数p和随机变量维数d决定；
[0027]
pc系数w
kj
及x
ij
(t0)根据wk及xi(t0)的分布确定，采用随机响应面法导出求解系数x
ij
(t)的表达式；
[0028]
将式(2)、(3)代入式(1)的微分方程中得：
[0029][0030]
其中，i＝1,2,...,n，k＝1,2,...,s；
[0031]
采用拉丁超立方方法进行ns次随机采样，并将随机样本点代入方程(4)中，得微分方程组：
[0032][0033]
式(5)中，x与w分别表示随机状态向量x与参数w的pc系数向量，f(x,u,t,w)为关于扩展状态x的非线性函数，且有：
[0034][0035]
w＝[w
10
,
…
,w
1p
,w
20
,
…
,w
2p
,
…
,w
s0
,
…
,w
sp
]
t (7)
[0036][0037]
设ω(δm)＝[φ0(δm),φ1(δm),
…
,φ
p
(δm)]，(m＝1,2,...,ns)，定义：
[0038][0039]
则
[0040]
若ns≥p+1，则通过最小二乘回归法求解式(5)的微分方程组，即有：
[0041]
[0042]
式(10)中，常值矩阵
[0043]
采样点ns＝2(p+1)，通过积分式(10)即获得每一时刻的pc系数向量x(t)；
[0044]
(2)过程约束及目标函数的pc展开
[0045]
由于过程约束函数g(x(t),u(t),t,w)与目标函数j(x(t),u(t),t,w)均为状态向量x(t)的函数，故采用上述进行混沌多项式展开，具体如下：
[0046][0047][0048]
式(11)和(12)中，gj(x(t),u(t),t)与jj(x(t),u(t),t)为对应函数的pc多项式系数；
[0049]
将随机样本点代入方程(11)、(12)得：
[0050][0051][0052]
式(13)和(14)中，与分别表示函数g(x(t),u(t),t,w)与j(x(t),u(t),t,w)的pc系数向量，且有
[0053][0054][0055]
[0056][0057]
类似地，有：
[0058]
ψ＝[ω(δ1)
t
,ω(δ2)
t
,
…
,ω(δ
ns
)
t
]
t (19)
[0059]
同理，通过最小二乘回归法求解系数向量与
[0060][0061][0062]
式(20)和(21)中，矩阵函数g(x,u,t,w)与j(x,u,t,w)参考式(17)、(18)确定；
[0063]
(3)统计量计算
[0064]
将式(1)中随机状态向量x、过程约束函数g(x(t),u(t),t,w)及目标函数j(x(t),u(t),t,w)的均值和方差表示为pc多项式系数的函数；
[0065][0066][0067][0068]
其中，＜,＞表示内积运算；gj、jj(j＝1,2,...,p)通过式(20)、(21)表示为x、u的非线性函数，进而通过式(22)、(23)、(24)计算状态、过程约束与目标函数的统计量。
[0069]
优选地，所述利用逐次线性化方法将高维确定性最优控制问题p
ed
凸化为凸优化问题p
ed-convex
，具体为将不确定性量化传播模型与凸优化方法相结合，通过逐次线性化方法对高维确定性最优控制问题p
ed
中的高维确定性微分方程、非凸约束、目标函数进行凸化处理，将高维确定性最优控制问题p
ed
转化为凸优化问题p
ed-convex
。
[0070]
优选地，所述将高维确定性最优控制问题p
ed
凸化为凸优化问题p
ed-convex
，具体包括以下步骤：
[0071]
(1)凸化处理
[0072]
采用逐次线性化方法对式(10)进行凸化处理；
[0073]
将式(10)在参考轨迹{xk,uk}附近进行一阶泰勒展开，得：
[0074][0075]
定义：
[0076][0077]
其中，与分别为参考轨迹上函数f(x,u,t,w)相对于状态x和控制量u的偏导数；
[0078]
为保证线性化的精度，增加信赖域约束：
[0079][0080]
式(27)中，δ
x
、δu表示信赖域半径；
[0081]
考虑式(1)中的过程约束及目标函数均为关于扩展状态x及控制量u的非线性函数，仍采用逐次线性化方法对其进行凸化处理；
[0082]
在参考轨迹{xk,uk}附近进行一阶泰勒展开，得：
[0083][0084][0085]
定义：
[0086]
[0087][0088]
至此，问题p
ed
转化为以下等价的高维确定性轨迹凸优化问题p
ed-convex
：
[0089][0090]
优选地，所述将凸优化问题p
ed-convex
进行离散化处理，得到参数凸优化问题p
ed-disconvex
，具体为：
[0091]
基于梯形法则对扩展状态量、控制量进行离散化处理，将凸优化问题p
ed-convex
转化为参数凸优化问题p
ed-disconvex
，进而构建基于序列凸优化的优化求解策略，通过迭代求解问题p
ed-disconvex
逼近鲁棒优化问题pr的解。
[0092]
优选地，所述基于梯形法则对扩展状态量、控制量进行离散化处理，将凸优化问题p
ed-convex
转化为参数凸优化问题p
ed-disconvex
，进而构建基于序列凸优化的优化求解策略，通过迭代求解问题p
ed-disconvex
逼近鲁棒优化问题pr的解，具体包括以下步骤：
[0093]
(1)离散化处理
[0094]
基于梯形法则，将自变量变化域[t0,tf]等间距离散为n个间隔，自变量离散为{t0,t1,
…
,tn}，状态量离散为{x0,x1,
…
,xn}，控制量离散为{u0,u1,
…
,un}；
[0095]
对式(25)的线性微分方程进行离散化，得：
[0096][0097]
其中，t＝(t
f-t0)/n；
[0098]
对问题p
ed-convex
的其他约束及目标函数进行离散化处理，最终得问题p
ed-disconvex
：
[0099][0100]
(2)序列凸优化求解策略
[0101]
飞行器鲁棒轨迹优化问题pr的解通过迭代求解问题p
ed-disconvex
进行逼近，并以当前获得轨迹来不断更新参考轨迹，即令：
[0102][0103]
直到满足约束：
[0104][0105]
式(36)中，ε
x
、εu为收敛误差限。
[0106]
优选地，所述将飞行器的总体参数、约束条件参数、不确定参数及权重系数作为参数输入，通过调用基于序列凸优化的优化求解算法进行优化求解，得到最优控制量，完成飞行器轨迹优化，具体包括以下步骤：
[0107]
设置飞行器总体参数，构建飞行器质心运动模型；
[0108]
设置不确定参数w服从的分布类型及标准差，进而确定pc多项式的基函数表达式、随机变量维数d、多项式阶数p及拉丁超立方采样点；
[0109]
设置飞行轨迹初始、终端状态约束、控制量约束、过程约束、目标函数及权重系数，进而构建具体的轨迹优化问题；
[0110]
通过调用优化求解算法对轨迹优化问题进行优化求解，得到最优控制量，完成飞行器轨迹优化。
[0111]
优选地，所述飞行器总体参数包括质量、参考面积、气动系数。
[0112]
本发明提供的基于srsm和nipc的飞行器鲁棒轨迹优化方法具有以下有益效果：
[0113]
通过构造基于nipc与srsm的不确定性量化传播模型，实现了不确定性随机问题的有效转化；随后，将该模型与凸优化方法相结合，设计了一种基于序列凸优化算法的求解策略，以实现对该高维确定性问题的快速求解，可以获得兼具可靠性和鲁棒性的优化轨迹，且具有相当的精度和显著的计算效率优势。基于所构建的基于随机响应面法及非嵌入式混沌多项式的不确定性量化传播模型，可实现初始状态和过程动力学等多参数不确定随机问题
的有效转化，且具有状态扩展维数少、计算效率高的优势。
附图说明
[0114]
为了更清楚地说明本发明实施例及其设计方案，下面将对本实施例所需的附图作简单地介绍。下面描述中的附图仅仅是本发明的部分实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0115]
图1为本发明实施例1的基于srsm和nipc的飞行器鲁棒轨迹优化方法的原理框图；
[0116]
图2为标称情况下滑翔轨迹优化结果；
[0117]
图3为标称情况下优化控制量打靶结果；
[0118]
图4为不确定性量化传播的状态均值对比结果；
[0119]
图5为不确定性量化传播的状态标准差对比结果；
[0120]
图6为gq-pc-co与srsm-pc-co的鲁棒轨迹优化结果对比结果；
[0121]
图7为do、rbo与ro的结果参数对比结果；
[0122]
图8为do、rbo与ro打靶结果概率密度函数曲线。
具体实施方式
[0123]
为了使本领域技术人员更好的理解本发明的技术方案并能予以实施，下面结合附图和具体实施例对本发明进行详细说明。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
[0124]
实施例1
[0125]
本发明针对航空航天领域鲁棒轨迹优化算法存在的问题，提出了一种通用、高效的鲁棒轨迹优化设计方法。通过构造基于nipc与srsm的不确定性量化传播模型，实现了不确定性随机问题的有效转化；随后，将该模型与凸优化方法相结合，设计了一种基于序列凸优化算法的求解策略，以实现对该高维确定性问题的快速求解。最终，以高超声速滑翔再入为例进行了仿真，验证了本发明方法的有效性。
[0126]
基于上述原理，实施例提供了一种基于srsm和nipc的飞行器鲁棒轨迹优化方法，具体包括：
[0127]
步骤1、构建参数不确定条件下的飞行器鲁棒轨迹优化问题pr，该问题为随机优化问题。
[0128]
一般的轨迹优化问题可以归结为如下确定性最优控制问题pd：
[0129][0130]
式中，x、u分别为状态、控制向量，w为系统参数向量；为系统
动力学方程；g(x(t),u(t),t,w)≤0为过程约束，且包含了状态边界约束；u
min
与u
max
为控制量的上下边界；x(t0)＝x0，x(tf)＝xf为端点约束。
[0131]
当系统参数向量w及状态初值x(t0)存在不确定性时，目标函数j、动力学函数f及约束函数均变为随机函数。因此，原确定性优化问题pd转变为随机优化问题。为了降低优化轨迹对不确定性的敏感度，即当存在不确定性时，预期的轨迹离散度在统计意义上尽可能小，通常需要在目标函数、约束函数中引入随机量的统计矩。因此，一般飞行器鲁棒轨迹优化问题pr可描述为：
[0132][0133]
式(1)中，x(t)为状态向量，u(t)为控制向量，w为系统随机参数向量；u
min
与u
max
为控制量的上下边界；j为目标函数，f为动力学函数，g为过程约束函数，下标μ和σ代表对应变量的均值和标准差，下标0和f分别表示初始与终端状态；kj、kg为待设计权重系数，x
0μ
与x
0σ
分别为初始状态的均值和标准差，εf为状态终值标准差的上界，为状态向量的随机微分向量。当采用上述模型进行轨迹优化时，需要引入不确定性量化传播技术计算相关统计量并处理随机动力学方程。
[0134]
需要注意：(1)εf需要基于先验统计信息确定，若取值不当可能造成问题无解，为了提升优化算法的收敛性，可将x
σ
(tf)以罚函数的形式引入目标函数中；(2)上述鲁棒轨迹优化模型同时考虑了轨迹的鲁棒性及可靠性，其中，在目标函数及端点约束中引入均值、标准差项可提升轨迹的鲁棒性，而通过引入考虑统计矩的过程约束模型，则可以保证以一定概率满足过程约束从而提升轨迹可靠性；(3)通过设计权重系数kj、kg可实现轨迹的鲁棒性及可靠性之间的权衡。
[0135]
如图1所示，下面将详细给出考虑初始状态x(t0)及系统参数向量w不确定性的鲁棒轨迹优化算法。首先，采用基于响应面法的混沌多项式展开技术将原始飞行器鲁棒轨迹优化问题pr转化为状态扩展的高维确定性优化问题；随后，对该问题进行凸化、离散处理，并采用序列凸优化算法进行求解。
[0136]
步骤2、基于非嵌入式混沌多项式nipc与随机响应面法srsm将飞行器鲁棒轨迹优化问题pr转化为高维确定性最优控制问题p
ed
，构建不确定性量化传播模型。
[0137]
具体为通过引入拉丁超立方采样策略，并基于随机响应面法srsm将含不确定参数的状态向量x(t)扩展为关于pc多项式系数的高维常微分方程，并通过pc多项式系数计算问题pr中状态、约束、目标函数的均值及标准差，进而将飞行器鲁棒轨迹优化问题pr转化为等
效的高维确定性最优控制问题p
ed
。
[0138]
下面将nipc与srsm相结合，同时引入拉丁超立方采样策略，建立了不确定性量化传播模型，从而将飞行器鲁棒轨迹优化问题转化为高维状态空间中的等价确定性优化问题。
[0139]
(1)状态x的pc展开
[0140]
假设随机状态向量x为n维，系统参数向量w为s维，则每个状态变量xi(t)及参数wk可展开为pc多项式：
[0141][0142][0143]
其中，δ＝{δ1,
…
,δd}为d维随机向量，φj(δ)为正交多项式函数，可通过与δ分布类型相对应的一元正交基函数φj(δ
l
),l＝1,...,d求张量积得到；x
ij
(t)和w
kj
为pc系数，展开项数p+1＝[(p+d)！]/(p！d！)，即由正交基函数φj(δ)的阶数p和随机变量维数d决定；正交基函数{φj(δ
l
)}与随机变量δ
l
的关系参照。
[0144]
pc系数w
kj
及x
ij
(t0)可根据wk及xi(t0)的分布确定，而随机状态向量xi(t)的pc系数x
ij
(t)未知，下面采用随机响应面法导出求解系数x
ij
(t)的表达式。
[0145]
将式(2)、(3)代入式(1)的微分方程中得：
[0146][0147]
其中，i＝1,2,...,n，k＝1,2,...,s。
[0148]
采用拉丁超立方方法进行ns次随机采样，并将随机样本点代入上述方程(4)，可得微分方程组：
[0149][0150]
式(5)中，x与w分别表示随机状态向量x与参数w的pc系数向量，且有：
[0151][0152]
w＝[w
10
,
…
,w
1p
,w
20
,
…
,w
2p
,
…
,w
s0
,
…
,w
sp
]
t (7)
[0153][0154]
设ω(δm)＝[φ0(δm),φ1(δm),
…
,φ
p
(δm)]，(m＝1,2,...,ns)，定义：
[0155][0156]
则
[0157]
若ns≥p+1，则可通过最小二乘回归法求解式(5)的微分方程组，即有：
[0158][0159]
式(10)中，矩阵为保证该模型的逼近精度，本发明统一取采样点ns＝2(p+1)，后文不再赘述。通过积分上式即可获得每一时刻的pc系数向量x(t)。
[0160]
(2)过程约束及目标函数的pc展开
[0161]
由于非线性约束函数g(x(t),u(t),t,w)与目标函数j(x(t),u(t),t,w)均为状态向量x(t)的函数，故也可采用类似方法进行混沌多项式展开，具体如下：
[0162][0163][0164]
式(11)和(12)中，gj(x(t),u(t),t)与jj(x(t),u(t),t)为对应函数的pc多项式系数。
[0165]
将随机样本点代入方程(11)、(12)可得：
[0166][0167][0168]
式(13)和(14)中，与分别表示函数g(x(t),u(t),t,w)与j(x(t),u(t),t,w)的pc系数向量，且有
[0169][0170]
[0171][0172][0173]
类似地，有
[0174][0175]
同理，可通过最小二乘回归法求解系数向量与
[0176][0177][0178]
式(20)和(21)中，矩阵函数g(x,u,t,w)与j(x,u,t,w)参考式(17)、(18)确定。
[0179]
(3)统计量计算
[0180]
将式(1)中随机状态向量x、过程约束函数g(x(t),u(t),t,w)及目标函数j(x(t),u(t),t,w)的均值和方差表示为pc多项式系数的函数；
[0181][0182][0183][0184]
其中，＜,＞表示内积运算；gj、jj(j＝1,2,...,p)可通过式(20)、(21)表示为x、u的
非线性函数，进而可通过以上公式(22)、(23)、(24)计算状态、过程约束与目标函数的统计量。
[0185]
至此，已将飞行器鲁棒轨迹优化问题pr转化为关于扩展状态x及控制量u的高维确定性优化问题p
ed
，进而可利用凸优化算法进行处理求解，限于篇幅，问题p
ed
在此不再赘述。
[0186]
如前所述，采用非嵌入式混沌多项式nipc进行参数不确定性量化传播时，时常用到两种典型模型转化方法，即随机响应面法(srsm)和高斯求积法(gq)，以下通过定量分析对这两种方法的计算效率进行比较。
[0187]
一般而言，对于具有d维随机变量的不确定问题，若pc多项式展开的阶数为p，gq法的状态扩展倍数为(p+1)d，而srsm方法的状态扩展倍数为2(p+d)！/(p！d！)。表1给出了两种方法在不同p、d取值情况下的状态扩展情况。
[0188]
表1不同条件下gq与srsm的状态扩展性能对比
[0189][0190]
由表1可知，在p、d较小的情况下，两种方法的状态扩展倍数比较接近；随着p的增大，两种方法的状态扩展倍数均有所增大，相较于gq方法，srsm方法的增长速度较为缓慢；随着随机变量维数d的不断增大，gq方法的扩展倍数呈指数增长，而srsm方法增长较为缓慢。可见，相较于gq方法，本发明采用的srsm方法具有更高的计算效率，尤其在处理高维不确定性问题时具有显著优势，而gq方法一般只能适用于d≤3的情况。此外，后文仿真部分亦给出了两种方法在精度与计算效率方面的进一步对比分析。
[0191]
步骤3、利用凸优化方法将不确定性量化传播模型中的高维确定性最优控制问题p
ed
凸化为凸优化问题p
ed-convex
。
[0192]
具体为将不确定性量化传播模型与凸优化方法相结合，通过逐次线性化方法对高维确定性最优控制问题p
ed
中的高维确定性微分方程、非凸约束、目标函数进行凸化处理，将高维确定性最优控制问题p
ed
转化为凸优化问题p
ed-convex
。
[0193]
步骤4、将凸优化问题p
ed-convex
进行离散化处理，得到参数凸优化问题p
ed-disconvex
，并基于参数凸优化问题p
ed-disconvex
构建序列凸优化算法，通过迭代求解问题p
ed-disconvex
逼近鲁棒优化问题pr的解。
[0194]
具体为基于梯形法则对扩展状态量、控制量进行离散化处理，将凸优化问题p
ed-convex
转化为参数凸优化问题p
ed-disconvex
，进而构建基于序列凸优化的优化求解策略，通过迭代求解问题p
ed-disconvex
逼近鲁棒优化问题pr的解。
[0195]
考虑到复杂轨迹优化问题大多是非凸的，如何将非凸优化问题进行凸化处理是应用凸优化技术求解最优轨迹的关键，并吸引广大学者对此开展了深入研究，相关成果亦在航空航天领域最优控制及轨迹优化问题得到了较为广泛的应用。
[0196]
这里引入逐次线性化等技术对上节中的高维确定性优化问题p
ed
进行凸化处理，并将该问题离散处理为参数凸优化问题，进而建立相应的序列凸优化迭代算法对该问题进行求解。序列凸优化方法通过外环迭代更新参考轨迹，内环求解凸优化问题逐步逼近原问题
的最优解。以下给出凸化、离散处理过程以及迭代算法原理和实现步骤。
[0197]
(1)凸化处理
[0198]
分析式可知，h为常值矩阵，f(x,u,t,w)为关于扩展状态x的非线性函数，因此该式为非线性微分方程。下面采用逐次线性化方法对式(10)进行凸化处理。
[0199]
将式(10)在参考轨迹{xk,uk}附近进行一阶泰勒展开，可得：
[0200][0201]
定义
[0202][0203]
其中，与分别为参考轨迹上函数f(x,u,t,w)相对于状态x和控制量u的偏导数。
[0204]
为保证线性化的精度，增加信赖域约束：
[0205][0206]
式中，δ
x
、δu表示信赖域半径。
[0207]
类似地，考虑式(1)中的过程约束及目标函数均为关于扩展状态x及控制量u的非线性函数，仍采用逐次线性化方法对其进行凸化处理。
[0208]
在参考轨迹{xk,uk}附近进行一阶泰勒展开，可得：
[0209][0210][0211]
定义：
[0212]
[0213][0214]
其中，相关偏导数的定义与前文类似，此处不再赘述。
[0215]
至此，问题p
ed
转化为以下等价的高维确定性轨迹凸优化问题p
ed-convex
：
[0216][0217]
需要说明的是，问题p
ed-convex
的待优化状态为x，扩展为问题pr的p+1倍。
[0218]
(2)离散化处理
[0219]
为了应用数值方法对问题p
ed-convex
进行求解，需要对状态量和控制量进行离散化处理。基于梯形法则，将自变量变化域[t0,tf]等间距离散为n个间隔，自变量离散为{t0,t1,
…
,tn}，状态量离散为{x0,x1,
…
,xn}，控制量离散为{u0,u1,
…
,un}。
[0220]
对式(25)的线性动力学方程进行离散化，可得：
[0221][0222]
其中，t＝(t
f-t0)/n。
[0223]
类似地，对问题p
ed-convex
的其他约束及目标函数进行离散化处理，最终可得问题p
ed-disconvex
：
[0224][0225]
至此，p
ed-convex
就转化为参数凸优化问题，可通过经典内点法求解。
[0226]
(3)序列凸优化求解策略
[0227]
综上所述，飞行器鲁棒轨迹优化问题pr、即非线性随机最优控制问题的解可通过迭代求解问题p
ed-disconvex
进行逼近，并以当前获得轨迹来不断更新参考轨迹，以提高解的收敛性和精度，即令：
[0228][0229]
直到满足约束：
[0230][0231]
式(36)中，ε
x
、εu为收敛误差限。
[0232]
步骤5、将飞行器总体参数、约束条件参数、不确定参数及权重系数作为输入，通过序列凸优化算法进行优化求解，输出轨迹优化的控制量，具体包括以下步骤：
[0233]
(1)设置飞行器总体参数，构建飞行器质心运动模型；
[0234]
(2)设置初始状态及环境参数不确定性服从的分布类型及标准差，进而确定pc多项式的基函数表达式、随机变量维数d、多项式阶数p及拉丁超立方采样点。
[0235]
(3)设置飞行轨迹初始、终端状态约束、控制量约束、过程约束、目标函数及权重系数。
[0236]
(5)将参考轨迹代入不确定量化传播模型中，得到扩展状态的初值及控制量初值以构建参数凸优化问题p
ed-disconvex
。
[0237]
(5)基于序列凸优化求解策略，通过迭代求解问题p
ed-disconvex
逼近原鲁棒优化问题pr的解，输出鲁棒轨迹优化控制量u
l
。
[0238]
下面以高超声速滑翔鲁棒轨迹优化为例，来对本实施例提供的优化算法进行验证。
[0239]
高超声速滑翔飞行具有航程远、速度高、飞行环境复杂等典型特性，且面临多种复
杂约束以及参数不确定性和干扰的影响，这使得其轨迹优化问题成为一个典型的快时变、非线性、强耦合、非凸的随机最优控制问题。研究表明，标称条件下的再入滑翔轨迹优化具有较大的挑战性和应用价值，而考虑真实情况下多参数不确定性的鲁棒轨迹优化问题则具有更大的难度和挑战性。因此，下面以典型的高超声速滑翔飞行为例，开展多不确定性条件下的鲁棒轨迹快速优化仿真，以验证本发明所设计算法的有效性；同时，给出了基于高斯求积策略及混沌多项式的凸优化方法(记为gq-pc-co)的仿真结果，通过对比以体现本发明算法在计算效率方面的优势。
[0240]
第一、高超声速滑翔鲁棒轨迹优化问题描述
[0241]
假设地球为旋转圆球，标称条件下的高超声速滑翔质心运动模型采用现有的形式。针对目标点固定情况下的飞行时间最短滑翔问题，若结任务需求预先选定攻角剖面，而仅将倾侧角作为待优化控制量，可将该轨迹优化对应的确定性最优控制问题描述为：
[0242][0243]
式中：自变量为无量纲能量e＝1/r-v2/2，无量纲状态量x＝[r,θ,φ,γ,ψ]
t
；r为无量纲地心距，v为无量纲速度，θ为地心经度，φ为地心纬度，γ为当地速度倾角，ψ为航迹偏航角，υ为倾侧角，r
min
为考虑过程约束的地心距下边界，具体可根据最大热流密度最大动压q
max
、最大过载n
max
计算。
[0244]
再入滑翔飞行所面临的不确定性主要包括再入初始状态，以及大气密度、气动力系数等动力学参数的不确定性，具体可建模描述为：
[0245][0246]
且有：
[0247][0248]
c＝[ρ,c
l
,cd]
t
ꢀꢀꢀꢀꢀ
(40)
[0249][0250]
其中，c为动力学参数向量，ρ为大气密度，c
l
、cd为升力、阻力系数，σ为标准差，δ为标准差与标称值之比，上划线
“‑”
代表标称情况。随机变量δ之间相互独立且服从标准正态分布。
[0251]
为了尽可能降低终端落点散布，将终端经纬度的标准差引入鲁棒轨迹优化问题的目标函数中，即有：
[0252]
j＝μ(tf)+kj[σ(tf)+σ(θf)+σ(φf)] (42)
[0253]
对于上述考虑不确定性的滑翔轨迹鲁棒优化问题而言，初始状态不确定性为5，过程参数不确定性为3，若同时考虑这些不确定性，则随机变量维数d＝8。设置混沌多项式阶数p＝2，若采用gq-pc-co方法进行鲁棒轨迹优化，则需将待优化状态量扩展为标称情况的38＝6561倍，原始随机常微分方程随之扩展为6561
×
5维耦合的确定性微分方程组，如考虑凸优化过程中的离散化处理，待优化参数量将会更大，以上高维扩展会导致计算效率下降甚至无法收敛等计算困难；而采用本发明所提出的基于srsm和nipc的凸优化(srsm-pc-co)方法进行求解时，待优化状态只需扩展(8+2)！/(8！2！)＝45倍，相较而言，显著降低了待优化变量个数，因而具有良好的优化计算效率和收敛性。
[0254]
仿真分析
[0255]
(1)仿真条件设置
[0256]
采用美国的cav-h飞行器进行仿真，标称边界条件设置如表2所示，参数不确定性设置如表3所示。主要约束设置为：q
max
＝100kpa，n
max
＝3，υ
max
＝80
°
。
[0257]
综合考虑航程能力和防热需求，将攻角剖面预设如下
[0258][0259]
表2边界条件取值
[0260][0261]
表3参数不确定性设置
[0262][0263]
仿真中采用固定数目离散点，取值n＝300，序列凸优化算法的信赖域半径和收敛误差限设置如下：
[0264][0265]
所有仿真均在搭载intel corei7-87003.20ghz intel处理器的台式机完成，仿真环境为matlab2016b平台。基于自行开发的代码进行不确定性量化传播，同时基于cvx工具包进行轨迹优化算法开发，并调用sdpt3求解器求解凸优化子问题p
ed-disconvex
。
[0266]
第二、不确定性量化传播仿真分析
[0267]
首先，基于标称条件进行确定性滑翔轨迹优化仿真，所得结果如图2所示，(a)h-v剖面，(b)过程约束，(c)倾侧角指令。由图2(a)、(b)可以看出，所得速度-高度曲线连续光滑并且满足过程约束要求，热流密度和动压曲线在一段时间里靠近约束边界，具有超出约束边界的潜在风险；由图2(a)、(c)可以看出，优化结果与积分结果一致，且控制指令反转次数较少，优化结果具有较好的可行性。
[0268]
仅考虑大气密度、升力和阻力系数3个动力学参数不确定性，基于上述优化结果开展gq-pc、srsm-pc方法的不确定性量化传播仿真，并与蒙特卡洛打靶结果进行对比，以分析不确定性对优化轨迹的影响，同时验证所给出不确定性量化算法的有效性。设置pc多项式阶数p＝2，蒙特卡洛打靶次数为3000次，部分仿真结果如图3～5所示。
[0269]
图3(a)终端经纬度散布，(b)热流约束，由图3可知，标称情况下的轨迹优化结果易受到参数不确定性的影响，从而出现过程约束超出限制、终端精度下降等问题，即轨迹可靠性及鲁棒性均面临不同程度的降低。因此，需要进一步开展不确定条件下的鲁棒轨迹优化，以提升轨迹的抗干扰能力及可靠性，进而降低制导控制系统负担和任务失败风险。
[0270]
图4～5分别给出了不同方法对应的部分状态均值及标准差比较曲线。图4(a)高度和速度偏角均值，(b)经度和纬度均值；图5(a)高度和速度偏角标准差，(b)经度和纬度标准差。可以看出，三种方法所得统计结果较为一致，说明gq-pc与srsm-pc方法均可逼近打靶结果，即可替代打靶方法进行不确定性量化传播。在计算精度相当的单次量化传播过程中，srsm-pc方法仅需进行20次动力学方程积分，而gq-pc方法需要进行27次积分，其计算量为srsm-pc的1.35倍；若进一步增加不确定性维数并考虑优化算法的迭代计算过程，srsm-pc方法在计算效率方面的优势将会更加明显。
[0271]
第三、两种算法的鲁棒轨迹优化仿真对比
[0272]
为验证所设计算法在鲁棒轨迹优化方面的能力，下面采用上面提出的srsm-pc-co算法进行鲁棒轨迹优化仿真，同时与现有的gq-pc-co优化方法及确定性优化方法(do)进行对比，以验证本发明算法的有效性。设置权重系数kj＝kg＝3，其余仿真条件同上，结果对比如图6所示，图6(a)高度剖面，(b)过程约束，(c)地面航迹，(d)倾侧角指令。随后，利用图6(d)中的倾侧角指令分别进行蒙特卡洛打靶仿真，得到部分参数的统计结果表4所示。其中，δθ
μf
、δφ
μf
分别表示终端经度、纬度均值偏差，θ
σf
、φ
σf
分别表示终端经度、纬度标准差；t
μf
、t
σf
分别表示飞行时间的均值及标准差，pe表示轨迹超出过程约束边界的概率，即事故发生率。
[0273]
表4两种方法优化结果的打靶数据对比
[0274][0275]
由图6(a)～(b)可以看出，两组鲁棒轨迹优化结果基本一致，且高度、热流密度和动压曲线相较于do更远离过程约束边界，优化轨迹的可靠性得到了明显提升。由表5可知，两组鲁棒优化控制量的打靶统计结果较为一致，相较于do，其终端经纬度分布更加集中，均值偏差及标准差均有所降低，过程约束超出概率显著降低。以上结果表明，本发明方法与gq-pc-co的优化结果较为一致，其轨迹的鲁棒性和可靠性均显著提升，验证了本发明方法的有效性。
[0276]
此外，由图6(d)可知，相较于do，鲁棒优化使得倾侧角指令在初始下降段基本维持在0
°
附近，随后继续以相对较小的幅值飞行，以增大纵向平面内的法向爬升过载，进而提升前段轨迹高度，最终降低热流密度以遵守过程约束，而为了完成航向调整，飞行中后期则需要采用相对更大的倾侧角幅值实现快速转弯，如图6(c)所示。
[0277]
表5给出了两种算法的数值处理性能对比。可见，在相同条件下，本发明算法对应的待优化变量数目显著降低，迭代步数相当情况下的计算耗时降低了约62％，进一步说明了srsm-pc-co鲁棒轨迹优化算法在计算效率方面具有显著优势。
[0278]
表5两种优化算法的数值计算性能对比
[0279][0280]
第四、基于权重参数调整的轨迹优化可靠性与鲁棒性能分析
[0281]
为了进一步探讨所提出轨迹优化算法在可靠性和鲁棒性方面的权衡能力，以更好的满足不同任务的需求侧重，下文开展基于权重参数调整的鲁棒优化仿真与分析。
[0282]
由前文鲁棒优化建模可知，调整权重参数kg、kj可以满足不同的可靠性和鲁棒性需求。若令kg＝3，kj＝0，则优化模型演化为基于可靠性的轨迹优化(reliability-based optimization，rbo)；若令kg＝kj＝3，则问题演化为鲁棒轨迹优化(robust optimization，ro)。综合考虑动力学参数及初始状态不确定性，采用srsm-pc-co算法进行仿真，同时将其与确定性优化结果(do)对比，以探究参数调整所带来的轨迹特性改变。表6给出了不同条件下的权重系数取值，其他仿真条件与前文一致。仿真结果如图7所示，图7(a)高度剖面，(b)过程约束，(c)倾侧角指令，(d)地面航迹。
[0283]
表6权重系数设置
[0284][0285]
由图7(a)、(b)可知，ro与rbo对应的高度、热流密度和动压曲线相较于do更远离过程约束边界，即优化轨迹的可靠性得到了明显提升；与rbo相比，ro对应的高度在第一个波峰处和末段更高，说明这种轨迹形式有利于降低随机不确定性的影响，提升终端位置和飞行时间指标的鲁棒性。由图7(c)可知，与前文分析类似，为了降低高度曲线第一个波谷处的热流峰值即提升轨迹可靠性，rbo与ro对应初始下降段的倾侧角幅值均较小，且大小基本一致，表明rbo与ro对热流峰值的可靠性相当；在波谷之后，ro继续保持更小的倾侧角幅值以保持更大的飞行高度；而在飞行后段，ro则采用了更大的倾侧角幅值以实现航向快速调整，这一现象在图7(d)的地面航迹曲线结果中亦得到了进一步印证。事实上，轨迹可靠性和鲁棒性的提升对飞行器不同飞行阶段的控制能力提出了更高要求，且没有统一的规律可循，需要针对不同的问题进行专门分析。
[0286]
利用图7(c)中的倾侧角优化结果进行蒙特卡洛打靶仿真，得到部分参数的统计结果，如图8及表7所示，(a)终端经度，(b)终端纬度，(c)热流密度峰值，(d)动压峰值，(e)过载峰值，(f)飞行时间；其中，图8给出了打靶结果中部分参数的近似概率密度函数。由图8(a)、(b)及表7可知，rbo与ro的终端经纬度统计精度基本一致，且相较于do有明显改善，而由于考虑了终端性能和优化指标的鲁棒性，ro对应的终端经纬度和飞行时间散布在三者中最小。由图8(c)～(e)及表7可知，rbo与ro的热流、动压、过载峰值分布均远离约束边界，约束超出的概率由do的67.4％分别降低到0.90％和0.88％，可见两者的可靠性均有所提升。以上结果表明，与rbo结果相比，ro在保证轨迹可靠性的同时，对不确定性的敏感度最低，具有更好的鲁棒性。由图8(f)及表7可知，相较于do结果，rbo与ro对应的飞行时间均值较大，这是由于不确定性的引入和可靠性约束的施加会造成优化求解的可行域减小，迫使轨迹高度提升进而导致飞行时间增大。而相较于rbo，ro的飞行时间均值更大，这是由于ro的目标函数中综合考虑了飞行时间的均值及标准差，而rbo的目标函数中仅考虑了飞行时间的均值。因此，轨迹鲁棒性与可靠性的提升是以牺牲一定最优性为代价的，实际应用时，可根据设计需求调整权重系数，以满足轨迹的可靠性、鲁棒性和最优性综合需求。
[0287]
表7do、rbo与ro打靶数据对比
[0288][0289]
表8对比了三种仿真条件下的计算耗时及待优化变量数，可以看出，即使同时考虑初始状态及参数不确定性，rbo及ro的优化耗时均在可接受范围内，若采用并行计算，优化耗时可进一步缩短。
[0290]
表8do、rbo与ro的数值计算性能对比
[0291][0292]
综上所示，针对存在多参数不确定性的轨迹优化问题，本发明研究了一种通用、高效的鲁棒轨迹优化方法，并以高超声速再入滑翔轨迹优化为例进行了仿真分析与对比验证，得到：
[0293]
(1)与确定性轨迹优化和传统鲁棒轨迹优化算法相比，所设计鲁棒轨迹优化方法可以获得兼具可靠性和鲁棒性的优化轨迹，且具有相当的精度和显著的计算效率优势。
[0294]
(2)轨迹鲁棒性和可靠性的提升是以牺牲一定最优性为代价的，实际应用时，可根据设计需求调整权重系数，以综合满足轨迹的可靠性、鲁棒性和最优性需求。
[0295]
(3)基于所构建的基于随机响应面法及非嵌入式混沌多项式的不确定性量化传播模型，可实现初始状态和过程动力学等多参数不确定随机问题的有效转化，且具有状态扩展维数少、计算效率高的优势。
[0296]
(4)与基于嵌入式混沌多项式的优化方法相比，本发明方法将原系统模型视为“黑箱”，无需针对不同问题进行推导建模，因而具有较高的通用性和工程应用价值；此外，本发明方法不仅适用于一般的鲁棒轨迹优化问题，也可推广到一般的非线性鲁棒动态优化问题。
[0297]
以上所述实施例仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换，均属于本发明的保护范围。