基于伪谱凸优化和L1罚函数的弹道规划方法研究

王庆海,陈 琦,王中原,尹秋霖

(南京理工大学 能源与动力工程学院,江苏 南京 210094)

轨迹优化方法按照是否得到解析解被分为2大类:解析法和数值解法。解析法基于最优控制原理,针对实际问题模型,根据极值条件和边界条件推导出解析解。史金光等利用庞特里亚金极小值原理,设计了滑翔制导炮弹以最大射程为指标的最优控制参数。解析法对于简单系统较为有效,但难以求解复杂的多约束非线性优化问题。数值解法通过某种离散方法,对时间、状态和控制量进行离散,将非线性最优控制问题转化为有限参数规划问题,然后采用合适的参数规划算法求出最优解。直接打靶法是数值解法中最常用的一种。涂良辉等采用直接打靶法规划再入飞行器无动力跳跃轨迹。李瑜采用直接打靶法,对助推-滑翔导弹的最大射程和可达区域进行了分析。直接打靶法只对控制量进行离散,其优化参数少,便于理解,但是直接打靶法需要使用高阶积分算法求解状态变量,计算量大,效率不高。根据FRABIEN的研究,直接打靶法对初始猜测敏感度很高。伪谱法是另一种最常用的数值解法。陈琦等针对滑翔制导炮弹,采用Gauss伪谱法离散弹道,以最短飞行时间为性能指标,利用序列二次规划算法(SQP)规划了滑翔段弹道。黄诘等通过Radau伪谱法离散导弹轨迹,分别取最大终端速度和最大落角为性能指标,采用SQP方法规划出最优攻击轨迹。张科南等利用SQP方法,为高超声速飞行器规划了多约束突防弹道。伪谱离散与序列二次规划法结合的伪谱法在求解非线性最优化问题时十分有效,但其求解速度较慢。造成这一问题的主要原因是序列二次规划方法在每一次迭代中不仅要更新Jacobian矩阵,还需要更新Hessian矩阵。

凸问题在数学上已经被证明,使用内点法可以保证在多项式时间内获得全局最优解。由于这一固有的优点,近年来,利用凸优化求解轨迹优化问题成为研究热点。陈嘉澍等采用Chebyshev伪谱法离散高超声速飞行器动态方程,结合无损凸化方法,将原问题转化为标准二阶锥规划问题求解。ACIKMESE等以最小燃料消耗为指标,利用序列凸优化方法规划最优火星着陆轨迹。同样以火星着陆为研究背景,文献[13-14]以最小着陆误差为性能指标,利用凸优化算法规划着陆轨迹。王金波等利用翻转Radau伪谱法离散一级火箭回收轨迹,将空心非凸推力约束转换为二阶锥约束,以燃料消耗量为性能指标,通过序列凸优化方法得到最优回收轨迹。文献[16]采用梯形离散方法离散连续变量,利用序列凸优化方法规划了空地导弹的最大存速轨迹,但精度较差。ROH等采用梯形离散方法,离散导弹简化三自由度模型,将动态方程线性凸化,采用带L1惩罚的序列凸优化算法(LPSCP)实现实时弹道规划。

文献[9]提出的(LPSCP)算法信赖域宽度较大,且恒定;算法中L1惩罚系数初值较大,且呈指数增长。该算法在解决简单轨迹规划问题时表现良好,但面对较复杂的弹道规划问题时表现较差。较大的信赖域宽度造成线性近似精度低;指数增长的惩罚系数会在迭代过程中引起浮点数计算误差,尤其制导炮弹轨迹优化问题中包含大量小数量级决策变量,严重时会导致无法求解。

基于上述考虑,本文提出一种改进的L1惩罚序列凸规划算法(ILPSCP),算法中采用指数衰减策略对相对信赖域宽度进行更新,并采用指数增长带上界的L1惩罚系数。仿真结果显示:与传统L1惩罚序列凸规划算法(LPSCP)相比,本文提出的ILPSCP算法收敛速度更快,稳定性更强。

1 建模

本文以控制能量最优为性能指标,规划制导炮弹打击固定目标的滑翔弹道。为实现一定的打击效能,要求落角不大于给定上界值,末速度应不小于给定下界值。

1.1 制导炮弹纵向平面内滑翔段弹道模型

本文对制导炮弹纵向二维平面内方案弹道进行规划,对应的动态方程为

(1)

式中:为速度,为弹道倾角,为阻力,为升力。阻力和升力的计算公式为

式中:为特征面积,0为零升阻力系数,为诱导阻力系数,为升力系数导数,为攻角。

1.2 模型约束

假设制导炮弹在滑翔起控点速度为,起控点弹道倾角为,起控点坐标为(,),目标点坐标为(,)。为实现较好的打击效能,末速度约束下界为,落角约束上界为。

于是建立初始状态约束为

(2)

末状态约束:

(3)

考虑弹体稳定性,攻角不能过大,建立控制量约束:

≤≤

(4)

1.3 最优化模型

参考文献[5],以滑翔能量消耗最低为性能指标,对应二维方案弹道的性能指标函数为

需要注意的是,初始时间和终止时间是可变的,也是问题的决策变量。

综上,制导炮弹滑翔轨迹最优化问题模型为

2 伪谱离散和凸化

2.1 模型一般化

为了方便后续公式推导,本文先建立制导炮弹滑翔段控制能量最优轨迹优化模型对应的一般化模型。一般化模型中状态变量、控制变量以及状态函数都以向量形式表示。

状态向量∈×1,控制向量=×1,本文中=4,=1。状态向量:

=()

控制向量为

=

于是动态方程可以写为

(5)

式中:∈×1为状态函数向量。

目标函数可以写为

(6)

式中:()为控制向量的Euclidean范数。状态向量、控制向量和时间变量的上下界约束可以写为

≤≤

(7)

≤≤

(8)

≤≤

(9)

式中:L表示下边界约束,U表示上边界约束。

综上,制导炮弹滑翔段控制能量最优轨迹规划模型对应的一般化最优化模型为

2.2 Radau伪谱离散

伪谱法采用全局多项式插值,相对于传统的等距梯形离散方法,伪谱离散可以获得更高的精度。Radau伪谱法使用阶Radau勒让德多项式的根为节点,其节点是非对称、非等间距的,采用它的根作为插值节点可以避免等距高阶多项式插值经常出现的龙格现象。相比Gauss伪谱离散,Radau伪谱离散没有附加的积分约束,更为简洁。Radau多项式是阶勒让德多项式和-1阶勒让德多项式-1的差:

()=()--1(),∈[-1,1]

阶Radau多项式有个根,这个根包含=-1,<1。可以看出,配置点包含初始点,但不包含终止点。

定义∈[-1,1] 为伪谱时间,∈[,]为物理时间。伪谱时间和物理时间之间存在如下映射关系:

(10)

Radau伪谱法的节点包含阶Radau多项式的个根和一个附加终止点,共计+1个节点。

定义为第个节点对应的拉格朗日多项式函数的基函数,写作:

于是,状态向量的第个分量()可以用经过这+1个节点的拉格朗日多项式函数来近似。

(11)

对式(11)两边求导,可得:

于是第个状态变量在第个配置点处微分的值为

(12)

式中:∈(+1)×为状态矩阵,为第个状态变量在第个离散节点的状态。在这里定义控制矩阵为∈×,为第个控制变量在第个离散节点的控制量。∈×(+1)被称为微分矩阵,其每一行对应着相应配置点,每一列对应着相应拉格朗日基函数。

对式(10)两边求导,并代入式(5)可得:

定义伪谱域函数:

即:

离散形式为

(13)

式中:∈×。

=(,,,)

式中矩阵下标对应第个离散节点。结合式(12)和式(13)可得:

=(,,,)

(14)

由式(14)可知,本文研究的最优化问题的决策变量包含:状态变量、控制量、起始时间和结束时间。为了方便数值求解,需要根据标准数值规划模型,将所有决策变量转换成列向量。

定义决策变量列向量为

=((vec())(vec()))

式中:vec表示矩阵列变换操作,∈×1,为决策变量总数。

=(+1)++2

式中:(+1)对应状态决策变量总数,对应控制决策变量总数,常数2对应初末时间,为控制变量数目。

定义状态函数列向量为∈×1:

()=vec(())

于是式(14)可以改写为

=()

(15)

=(diag(,)()×(+2))

Radau伪谱法的离散求积公式为

式中:为被积函数在第个配置点的值;为对应配置点的积分权系数。于是式(6)的离散形式为

式中:表示第个控制变量。引入辅助变量建立二阶锥约束:

(16)

(17)

2.3 线性凸化

(18)

=

(19)

式中:∈(·)×,∈(·)×1。

弹道规划问题中,决策变量分为状态变量、控制变量和时间变量3个大类。每一个大类之间数量级相差较大,而大类之中决策变量的数量级也相差悬殊。本文采用比例缩放的思想,指定不同决策变量的缩放比例系数。比例缩放系数定义为决策变量上边界值与下边界值的差,即:

设定统一的相对信赖域宽度,从而实现简单地调整一个相对信赖域宽度,完成对所有决策变量信赖域宽度的调整。

相对信赖域约束定义为

(20)

式中:相对信赖域宽度为一个足够小的数。

同理,定义相对误差为

式(18)、式(20)的物理意义为:新轨迹点在原轨迹点附近,由于轨迹由轨迹点连接而成,即新轨迹在原轨迹附近,从变分学角度看,新轨迹在原轨迹的邻域内。

式(17)为线性目标函数,式(19)为线性等式约束方程,式(16)和式(20)为锥约束,于是式(16)～式(17)、式(19)～式(20)组成了标准凸优化模型,称该凸优化模型为原子问题:

3 改进L1惩罚序列凸规划算法

由文献[9]可知,采用普通序列凸规划算法(SCP)求解弹道规划问题时,其对初始预测十分敏感。如果不能提供合适的初始预测,子问题的解空间将为空,无法求解。合适的初始预测依靠丰富的经验,当问题较为复杂时,经验预测也难以适用。

L1惩罚序列凸规划算法(LPSCP)成功解决了序列凸规划算法(SCP)对初始预测敏感这个问题。LPSCP算法通过引入松弛辅助变量,,,,在原问题的基础上建立L1惩罚子问题。子问题定义为

式中:∈(·)×1和∈(·)×1为一组用于松弛线性等式约束方程的辅助变量;∈×1和∈×1为一组用于松弛信赖域不等式约束的辅助变量;为违反等式约束的惩罚系数,为违反信赖域不等式约束的惩罚系数;=1,2,…,。

2个惩罚系数的获得方式为

式中:为违反原子问题约束时,随着迭代次数呈指数增长的总惩罚系数,增长速率由系数确定。∈[0,1]为一个权系数,代表对2种约束惩罚权重。越大,对违反动态方程等式约束的惩罚越大,对违反信赖域不等式约束惩罚越少,反之同理。当=1时,代表只对动态约束采取L1惩罚;当=0时,代表只对信赖域约束采取L1惩罚。的具体取值根据实际情况而定,本文研究的制导炮弹轨迹优化问题动态约束直接影响初始预测的敏感性,相比之下更为重要,所以惩罚权系数取值最好不小于0.5。传统L1惩罚序列凸规划(LPSCP)算法流程如图1所示。

图1 传统LPSCP算法流程图

从式(18)、式(20)可以看出,信赖域宽度不仅影响线性近似的精度,还将决定每一次迭代优化的最大收敛步长。相对信赖域宽度越小,则最大收敛步长越小,但线性近似精度越高;相对信赖域宽度越大,则迭代最大收敛步长越大,但线性近似精度越低。传统序列凸规划算法使用固定信赖域宽度,如果相对信赖域宽度太大,则无法满足线性近似条件,会导致算法失效。为了兼具收敛速度和近似精度,改进的L1惩罚序列凸规划(ILPSCP)算法采用指数衰减策略更新相对信赖域宽度:

=max(×,)

式中:为初始相对信赖域宽度,由于采用了指数衰减策略,可以取得适当大一些,有利于前期迭代快速收敛;为衰减系数,越大相对信赖域宽度衰减越慢,的取值一般在08～098之间;为迭代次数;为相对信赖域宽度下界。算法的终止条件为前后两次迭代决策变量之间的相对差值不大于相对误差容错系数。信赖域宽度设置下界的原因在于相对信赖域宽度必须大于等于相对误差容错系数,即:

≥

否则当<,会导致解的精度并未达到预期,迭代过程便提前结束。

理论上,当总惩罚系数足够大时,子问题和子问题求得的解相同。然而在研究一个数值算法时,有必要考虑计算机浮点数计算精度。计算机运算中,数量级差距太大的数进行运算时可能会出现误差,这对于决策变量本身数量级差距较大的弹道规划问题是不利的,惩罚系数太大时,将无法求解。如图2所示,当惩罚系数达到10时,求得的所有的决策变量都在0附近,显然不是期望的弹道,是无效解。

图2 惩罚系数cex=107 时的仿真结果

综合考虑,为惩罚系数设定一个上界,使惩罚系数满足:

=min(×,)

改进的L1惩罚序列凸规划(ILPSCP)算法流程如图3所示。

图3 ILPSCP算法流程图

4 仿真和分析

本节以制导炮弹纵向平面内滑翔段弹道规划问题为例,分别采用传统L1惩罚序列凸规划算法(LPSCP)和本文提出的改进L1惩罚序列凸规划算法(ILPSCP)进行仿真对比。为验证本文提出的ILPSCP算法的准确性,使用非线性问题最优化通用工具箱GPOPS2进行仿真,作为参照组。

4.1 仿真参数设置

制导炮弹参数如表1所示。表2为各决策变量的上下边界值。

表1 制导炮弹参数

表2 决策变量的上下边界值

直接打靶法对初始预测敏感,需要遴选较好的初始预测才能求解,然而本文提出的ILPSCP算法具备强大的冷启动能力。正常情况下,弹道的始末端点坐标预测值会选择制导炮弹滑翔起点坐标和目标的坐标,这种初始预测经过拉格朗日插值得到的初始预测弹道轨迹比较接近期望的最优弹道轨迹,有利于算法的求解。表3为本文仿真使用的始末变量预测值,下标“0”对应变量起始值,下标“f”对应变量末值,下标“g”表示预测值。为了验证ILPSCP算法的冷启动能力,表3提供的始末端点坐标预测值极其恶劣,是2个重合的点,并且这个重合点即不是制导炮弹滑翔起点坐标,也不是目标坐标。

表3 始末变量初始预测值

表4中为弹道规划实例参数,主要包含各状态变量的初始点约束、终点约束以及规划算法中需要的各参数。为保证制导炮弹在飞行末段具备一定的机动能力,设置飞行速度下界为200 m/s;考虑到制导炮弹飞行稳定性;攻角绝对值最大值为π/12;其他参数见表4。

表4 实例仿真参数表

4.2 仿真结果与分析

仿真结果如图4～图9所示。图4为飞行轨迹对比图;图5为飞行速度对比图;图6为飞行过程中弹道倾角对比图;图7为攻角曲线对比图,也是该模型的控制曲线;图8为本文提出的ILPSCP算法和传统LPSCP算法的目标函数收敛曲线对比图;图9为GPOPS2和ILPSCP算法解的节点分布图。

图4 飞行轨迹对比图

图5 速度曲线对比图

图6 弹道倾角曲线对比图

图7 攻角曲线对比图

图8 目标函数收敛曲线图

图9 离散节点分布对比图

图4～图7中,给出了传统LPSCP算法第199次和第200次迭代的结果,可以看出,这两组曲线都没能和参照组曲线重合,即没有收敛到最优解。结合图8可以看出,LPSCP算法对应的目标函数,仅需20次迭代便趋于一个值附近,并在该值附近震荡。但在经过50次、100次甚至150次迭代之后,该目标函数曲线仍然持续震荡,难以收敛,最后在0.607和0.690之间持续跳变。可以预计,即使进行更多迭代仍然难以收敛。形成鲜明对比的是,图8中ILPSCP算法对应的目标函数曲线在经过30次迭代后便基本稳定,之后变化量极小。

由图4～图7和图9可以看出,ILPSCP算法的节点和GPOPS2的节点虽然不同,但ILPSCP算法规划出的各变量曲线与参照组GPOPS2高度重合,这证明本文提出的ILPSCP算法是有效的。

本文的弹道规划模型的性能指标为攻角平方积分,可以发现,图7中LPSCP算法的2条攻角曲线的平方积分比其他2种算法的更小,从图8中也可以看出,LPSCP算法的性能指标小于ILPSCP算法,这与本文描述的ILPSCP算法优于LPSCP算法这一观点不符。出现这种现象的原因是:LPSCP算法使用固定相对信赖域=03,较大的信赖域宽度对于简单问题影响不大,但对于本文较复杂的弹道规划问题,会造成线性近似假设条件不成立,近似误差大。过大的线性近似误差导致近似模型无法反映制导炮弹物理模型,所以,对于制导炮弹模型,图7中LPSCP算法的攻角曲线和图8中LPSCP算法的目标函数曲线都是无效的。过大的线性近似误差,也是导致LPSCP算法无法收敛的原因。采用=03的LPSCP算法在2个无效解之间跳变,而=03的ILPSCP算法收敛到最优解,显然ILPSCP算法比LPSCP算法更适合解复杂的弹道规划问题。

LPSCP算法在=03时便因为信赖域宽度过大无法求解,而ILPSCP算法在更大的初始相对信赖域宽度情况下,仍然能稳定求解。增大后,ILPSCP算法收敛到指定精度(最优解与=03时相同)所需迭代次数,如表5所示。随着初始相对信赖域的增大,所需迭代次数增加,但相比信赖域宽度的增长倍数,其增长幅度较小。相比于传统的LPSCP算法,即使设置了更大的初始相对信赖域宽度,ILPSCP算法仍然能快速稳定收敛,这进一步证明了ILPSCP算法的优越性。

表5 不同δ0下ILPSCP算法迭代次数

5 结束语

本文推导了时间自由且以控制能量最优为指标的弹道规划一般化模型,针对该模型采用伪谱离散和线性凸化完成凸优化建模。针对传统L1序列凸规划算法(LPSCP)在解决复杂的初末时间自由、非线性强的弹道规划问题时存在稳定性较差等问题,本文提出了一种改进的L1伪谱凸规划算法(ILPSCP)。ILPSCP算法中引入指数衰减的相对信赖域宽度和带上界的L1惩罚系数。仿真结果证明,ILPSCP算法收敛速度快,精度高,冷启动能力强,很好地解决了传统L1序列凸规划算法(LPSCP)存在的问题。