基于偽譜凸優(yōu)化和L1罰函數(shù)的彈道規(guī)劃方法研究

王慶海,陳 琦,王中原,尹秋霖

(南京理工大學 能源與動力工程學院,江蘇 南京 210094)

軌跡優(yōu)化方法按照是否得到解析解被分為2大類:解析法和數(shù)值解法。解析法基于最優(yōu)控制原理,針對實際問題模型,根據(jù)極值條件和邊界條件推導出解析解。史金光等利用龐特里亞金極小值原理,設(shè)計了滑翔制導炮彈以最大射程為指標的最優(yōu)控制參數(shù)。解析法對于簡單系統(tǒng)較為有效,但難以求解復雜的多約束非線性優(yōu)化問題。數(shù)值解法通過某種離散方法,對時間、狀態(tài)和控制量進行離散,將非線性最優(yōu)控制問題轉(zhuǎn)化為有限參數(shù)規(guī)劃問題,然后采用合適的參數(shù)規(guī)劃算法求出最優(yōu)解。直接打靶法是數(shù)值解法中最常用的一種。涂良輝等采用直接打靶法規(guī)劃再入飛行器無動力跳躍軌跡。李瑜采用直接打靶法,對助推-滑翔導彈的最大射程和可達區(qū)域進行了分析。直接打靶法只對控制量進行離散,其優(yōu)化參數(shù)少,便于理解,但是直接打靶法需要使用高階積分算法求解狀態(tài)變量,計算量大,效率不高。根據(jù)FRABIEN的研究,直接打靶法對初始猜測敏感度很高。偽譜法是另一種最常用的數(shù)值解法。陳琦等針對滑翔制導炮彈,采用Gauss偽譜法離散彈道,以最短飛行時間為性能指標,利用序列二次規(guī)劃算法(SQP)規(guī)劃了滑翔段彈道。黃詰等通過Radau偽譜法離散導彈軌跡,分別取最大終端速度和最大落角為性能指標,采用SQP方法規(guī)劃出最優(yōu)攻擊軌跡。張科南等利用SQP方法,為高超聲速飛行器規(guī)劃了多約束突防彈道。偽譜離散與序列二次規(guī)劃法結(jié)合的偽譜法在求解非線性最優(yōu)化問題時十分有效,但其求解速度較慢。造成這一問題的主要原因是序列二次規(guī)劃方法在每一次迭代中不僅要更新Jacobian矩陣,還需要更新Hessian矩陣。

凸問題在數(shù)學上已經(jīng)被證明,使用內(nèi)點法可以保證在多項式時間內(nèi)獲得全局最優(yōu)解。由于這一固有的優(yōu)點,近年來,利用凸優(yōu)化求解軌跡優(yōu)化問題成為研究熱點。陳嘉澍等采用Chebyshev偽譜法離散高超聲速飛行器動態(tài)方程,結(jié)合無損凸化方法,將原問題轉(zhuǎn)化為標準二階錐規(guī)劃問題求解。ACIKMESE等以最小燃料消耗為指標,利用序列凸優(yōu)化方法規(guī)劃最優(yōu)火星著陸軌跡。同樣以火星著陸為研究背景,文獻[13-14]以最小著陸誤差為性能指標,利用凸優(yōu)化算法規(guī)劃著陸軌跡。王金波等利用翻轉(zhuǎn)Radau偽譜法離散一級火箭回收軌跡,將空心非凸推力約束轉(zhuǎn)換為二階錐約束,以燃料消耗量為性能指標,通過序列凸優(yōu)化方法得到最優(yōu)回收軌跡。文獻[16]采用梯形離散方法離散連續(xù)變量,利用序列凸優(yōu)化方法規(guī)劃了空地導彈的最大存速軌跡,但精度較差。ROH等采用梯形離散方法,離散導彈簡化三自由度模型,將動態(tài)方程線性凸化,采用帶L1懲罰的序列凸優(yōu)化算法(LPSCP)實現(xiàn)實時彈道規(guī)劃。

文獻[9]提出的(LPSCP)算法信賴域?qū)挾容^大,且恒定;算法中L1懲罰系數(shù)初值較大,且呈指數(shù)增長。該算法在解決簡單軌跡規(guī)劃問題時表現(xiàn)良好,但面對較復雜的彈道規(guī)劃問題時表現(xiàn)較差。較大的信賴域?qū)挾仍斐删€性近似精度低;指數(shù)增長的懲罰系數(shù)會在迭代過程中引起浮點數(shù)計算誤差,尤其制導炮彈軌跡優(yōu)化問題中包含大量小數(shù)量級決策變量,嚴重時會導致無法求解。

基于上述考慮,本文提出一種改進的L1懲罰序列凸規(guī)劃算法(ILPSCP),算法中采用指數(shù)衰減策略對相對信賴域?qū)挾冗M行更新,并采用指數(shù)增長帶上界的L1懲罰系數(shù)。仿真結(jié)果顯示:與傳統(tǒng)L1懲罰序列凸規(guī)劃算法(LPSCP)相比,本文提出的ILPSCP算法收斂速度更快,穩(wěn)定性更強。

1 建模

本文以控制能量最優(yōu)為性能指標,規(guī)劃制導炮彈打擊固定目標的滑翔彈道。為實現(xiàn)一定的打擊效能,要求落角不大于給定上界值,末速度應不小于給定下界值。

1.1 制導炮彈縱向平面內(nèi)滑翔段彈道模型

本文對制導炮彈縱向二維平面內(nèi)方案彈道進行規(guī)劃,對應的動態(tài)方程為

(1)

式中:為速度,為彈道傾角,為阻力,為升力。阻力和升力的計算公式為

式中:為特征面積,0為零升阻力系數(shù),為誘導阻力系數(shù),為升力系數(shù)導數(shù),為攻角。

1.2 模型約束

假設(shè)制導炮彈在滑翔起控點速度為,起控點彈道傾角為,起控點坐標為(,),目標點坐標為(,)。為實現(xiàn)較好的打擊效能,末速度約束下界為,落角約束上界為。

于是建立初始狀態(tài)約束為

(2)

末狀態(tài)約束:

(3)

考慮彈體穩(wěn)定性,攻角不能過大,建立控制量約束:

≤≤

(4)

1.3 最優(yōu)化模型

參考文獻[5],以滑翔能量消耗最低為性能指標,對應二維方案彈道的性能指標函數(shù)為

需要注意的是,初始時間和終止時間是可變的,也是問題的決策變量。

綜上,制導炮彈滑翔軌跡最優(yōu)化問題模型為

2 偽譜離散和凸化

2.1 模型一般化

為了方便后續(xù)公式推導,本文先建立制導炮彈滑翔段控制能量最優(yōu)軌跡優(yōu)化模型對應的一般化模型。一般化模型中狀態(tài)變量、控制變量以及狀態(tài)函數(shù)都以向量形式表示。

狀態(tài)向量∈×1,控制向量=×1,本文中=4,=1。狀態(tài)向量:

=()

控制向量為

=

于是動態(tài)方程可以寫為

(5)

式中:∈×1為狀態(tài)函數(shù)向量。

目標函數(shù)可以寫為

(6)

式中:()為控制向量的Euclidean范數(shù)。狀態(tài)向量、控制向量和時間變量的上下界約束可以寫為

≤≤

(7)

≤≤

(8)

≤≤

(9)

式中:L表示下邊界約束,U表示上邊界約束。

綜上,制導炮彈滑翔段控制能量最優(yōu)軌跡規(guī)劃模型對應的一般化最優(yōu)化模型為

2.2 Radau偽譜離散

偽譜法采用全局多項式插值,相對于傳統(tǒng)的等距梯形離散方法,偽譜離散可以獲得更高的精度。Radau偽譜法使用階Radau勒讓德多項式的根為節(jié)點,其節(jié)點是非對稱、非等間距的,采用它的根作為插值節(jié)點可以避免等距高階多項式插值經(jīng)常出現(xiàn)的龍格現(xiàn)象。相比Gauss偽譜離散,Radau偽譜離散沒有附加的積分約束,更為簡潔。Radau多項式是階勒讓德多項式和-1階勒讓德多項式-1的差:

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

階Radau多項式有個根,這個根包含=-1,<1?？梢钥闯?配置點包含初始點,但不包含終止點。

定義∈[-1,1] 為偽譜時間,∈[,]為物理時間。偽譜時間和物理時間之間存在如下映射關(guān)系:

(10)

Radau偽譜法的節(jié)點包含階Radau多項式的個根和一個附加終止點,共計+1個節(jié)點。

定義為第個節(jié)點對應的拉格朗日多項式函數(shù)的基函數(shù),寫作:

于是,狀態(tài)向量的第個分量()可以用經(jīng)過這+1個節(jié)點的拉格朗日多項式函數(shù)來近似。

(11)

對式(11)兩邊求導,可得:

于是第個狀態(tài)變量在第個配置點處微分的值為

(12)

式中:∈(+1)×為狀態(tài)矩陣,為第個狀態(tài)變量在第個離散節(jié)點的狀態(tài)。在這里定義控制矩陣為∈×,為第個控制變量在第個離散節(jié)點的控制量。∈×(+1)被稱為微分矩陣,其每一行對應著相應配置點,每一列對應著相應拉格朗日基函數(shù)。

對式(10)兩邊求導,并代入式(5)可得:

定義偽譜域函數(shù):

即:

離散形式為

(13)

式中:∈×。

=(,,,)

式中矩陣下標對應第個離散節(jié)點。結(jié)合式(12)和式(13)可得:

=(,,,)

(14)

由式(14)可知,本文研究的最優(yōu)化問題的決策變量包含:狀態(tài)變量、控制量、起始時間和結(jié)束時間。為了方便數(shù)值求解,需要根據(jù)標準數(shù)值規(guī)劃模型,將所有決策變量轉(zhuǎn)換成列向量。

定義決策變量列向量為

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

式中:vec表示矩陣列變換操作,∈×1,為決策變量總數(shù)。

=(+1)++2

式中:(+1)對應狀態(tài)決策變量總數(shù),對應控制決策變量總數(shù),常數(shù)2對應初末時間,為控制變量數(shù)目。

定義狀態(tài)函數(shù)列向量為∈×1:

()=vec(())

于是式(14)可以改寫為

=()

(15)

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

Radau偽譜法的離散求積公式為

式中:為被積函數(shù)在第個配置點的值;為對應配置點的積分權(quán)系數(shù)。于是式(6)的離散形式為

式中:表示第個控制變量。引入輔助變量建立二階錐約束:

(16)

(17)

2.3 線性凸化

(18)

=

(19)

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

彈道規(guī)劃問題中,決策變量分為狀態(tài)變量、控制變量和時間變量3個大類。每一個大類之間數(shù)量級相差較大,而大類之中決策變量的數(shù)量級也相差懸殊。本文采用比例縮放的思想,指定不同決策變量的縮放比例系數(shù)。比例縮放系數(shù)定義為決策變量上邊界值與下邊界值的差,即:

設(shè)定統(tǒng)一的相對信賴域?qū)挾?從而實現(xiàn)簡單地調(diào)整一個相對信賴域?qū)挾?完成對所有決策變量信賴域?qū)挾鹊恼{(diào)整。

相對信賴域約束定義為

(20)

式中:相對信賴域?qū)挾葹橐粋€足夠小的數(shù)。

同理,定義相對誤差為

式(18)、式(20)的物理意義為:新軌跡點在原軌跡點附近,由于軌跡由軌跡點連接而成,即新軌跡在原軌跡附近,從變分學角度看,新軌跡在原軌跡的鄰域內(nèi)。

式(17)為線性目標函數(shù),式(19)為線性等式約束方程,式(16)和式(20)為錐約束,于是式(16)～式(17)、式(19)～式(20)組成了標準凸優(yōu)化模型,稱該凸優(yōu)化模型為原子問題:

3 改進L1懲罰序列凸規(guī)劃算法

由文獻[9]可知,采用普通序列凸規(guī)劃算法(SCP)求解彈道規(guī)劃問題時,其對初始預測十分敏感。如果不能提供合適的初始預測,子問題的解空間將為空,無法求解。合適的初始預測依靠豐富的經(jīng)驗,當問題較為復雜時,經(jīng)驗預測也難以適用。

L1懲罰序列凸規(guī)劃算法(LPSCP)成功解決了序列凸規(guī)劃算法(SCP)對初始預測敏感這個問題。LPSCP算法通過引入松弛輔助變量,,,,在原問題的基礎(chǔ)上建立L1懲罰子問題。子問題定義為

式中:∈(·)×1和∈(·)×1為一組用于松弛線性等式約束方程的輔助變量;∈×1和∈×1為一組用于松弛信賴域不等式約束的輔助變量;為違反等式約束的懲罰系數(shù),為違反信賴域不等式約束的懲罰系數(shù);=1,2,…,。

2個懲罰系數(shù)的獲得方式為

式中:為違反原子問題約束時,隨著迭代次數(shù)呈指數(shù)增長的總懲罰系數(shù),增長速率由系數(shù)確定?！蔥0,1]為一個權(quán)系數(shù),代表對2種約束懲罰權(quán)重。越大,對違反動態(tài)方程等式約束的懲罰越大,對違反信賴域不等式約束懲罰越少,反之同理。當=1時,代表只對動態(tài)約束采取L1懲罰;當=0時,代表只對信賴域約束采取L1懲罰。的具體取值根據(jù)實際情況而定,本文研究的制導炮彈軌跡優(yōu)化問題動態(tài)約束直接影響初始預測的敏感性,相比之下更為重要,所以懲罰權(quán)系數(shù)取值最好不小于0.5。傳統(tǒng)L1懲罰序列凸規(guī)劃(LPSCP)算法流程如圖1所示。

圖1 傳統(tǒng)LPSCP算法流程圖

從式(18)、式(20)可以看出,信賴域?qū)挾炔粌H影響線性近似的精度,還將決定每一次迭代優(yōu)化的最大收斂步長。相對信賴域?qū)挾仍叫?則最大收斂步長越小,但線性近似精度越高;相對信賴域?qū)挾仍酱?則迭代最大收斂步長越大,但線性近似精度越低。傳統(tǒng)序列凸規(guī)劃算法使用固定信賴域?qū)挾?如果相對信賴域?qū)挾忍?則無法滿足線性近似條件,會導致算法失效。為了兼具收斂速度和近似精度,改進的L1懲罰序列凸規(guī)劃(ILPSCP)算法采用指數(shù)衰減策略更新相對信賴域?qū)挾?

=max(×,)

式中:為初始相對信賴域?qū)挾?由于采用了指數(shù)衰減策略,可以取得適當大一些,有利于前期迭代快速收斂;為衰減系數(shù),越大相對信賴域?qū)挾人p越慢,的取值一般在08～098之間;為迭代次數(shù);為相對信賴域?qū)挾认陆?。算法的終止條件為前后兩次迭代決策變量之間的相對差值不大于相對誤差容錯系數(shù)。信賴域?qū)挾仍O(shè)置下界的原因在于相對信賴域?qū)挾缺仨毚笥诘扔谙鄬φ`差容錯系數(shù),即:

≥

否則當<,會導致解的精度并未達到預期,迭代過程便提前結(jié)束。

理論上,當總懲罰系數(shù)足夠大時,子問題和子問題求得的解相同。然而在研究一個數(shù)值算法時,有必要考慮計算機浮點數(shù)計算精度。計算機運算中,數(shù)量級差距太大的數(shù)進行運算時可能會出現(xiàn)誤差,這對于決策變量本身數(shù)量級差距較大的彈道規(guī)劃問題是不利的,懲罰系數(shù)太大時,將無法求解。如圖2所示,當懲罰系數(shù)達到10時,求得的所有的決策變量都在0附近,顯然不是期望的彈道,是無效解。

圖2 懲罰系數(shù)cex=107 時的仿真結(jié)果

綜合考慮,為懲罰系數(shù)設(shè)定一個上界,使懲罰系數(shù)滿足:

=min(×,)

改進的L1懲罰序列凸規(guī)劃(ILPSCP)算法流程如圖3所示。

圖3 ILPSCP算法流程圖

4 仿真和分析

本節(jié)以制導炮彈縱向平面內(nèi)滑翔段彈道規(guī)劃問題為例,分別采用傳統(tǒng)L1懲罰序列凸規(guī)劃算法(LPSCP)和本文提出的改進L1懲罰序列凸規(guī)劃算法(ILPSCP)進行仿真對比。為驗證本文提出的ILPSCP算法的準確性,使用非線性問題最優(yōu)化通用工具箱GPOPS2進行仿真,作為參照組。

4.1 仿真參數(shù)設(shè)置

制導炮彈參數(shù)如表1所示。表2為各決策變量的上下邊界值。

表1 制導炮彈參數(shù)

表2 決策變量的上下邊界值

直接打靶法對初始預測敏感,需要遴選較好的初始預測才能求解,然而本文提出的ILPSCP算法具備強大的冷啟動能力。正常情況下,彈道的始末端點坐標預測值會選擇制導炮彈滑翔起點坐標和目標的坐標,這種初始預測經(jīng)過拉格朗日插值得到的初始預測彈道軌跡比較接近期望的最優(yōu)彈道軌跡,有利于算法的求解。表3為本文仿真使用的始末變量預測值,下標“0”對應變量起始值,下標“f”對應變量末值,下標“g”表示預測值。為了驗證ILPSCP算法的冷啟動能力,表3提供的始末端點坐標預測值極其惡劣,是2個重合的點,并且這個重合點即不是制導炮彈滑翔起點坐標,也不是目標坐標。

表3 始末變量初始預測值

表4中為彈道規(guī)劃實例參數(shù),主要包含各狀態(tài)變量的初始點約束、終點約束以及規(guī)劃算法中需要的各參數(shù)。為保證制導炮彈在飛行末段具備一定的機動能力,設(shè)置飛行速度下界為200 m/s;考慮到制導炮彈飛行穩(wěn)定性;攻角絕對值最大值為π/12;其他參數(shù)見表4。

表4 實例仿真參數(shù)表

4.2 仿真結(jié)果與分析

仿真結(jié)果如圖4～圖9所示。圖4為飛行軌跡對比圖;圖5為飛行速度對比圖;圖6為飛行過程中彈道傾角對比圖;圖7為攻角曲線對比圖,也是該模型的控制曲線;圖8為本文提出的ILPSCP算法和傳統(tǒng)LPSCP算法的目標函數(shù)收斂曲線對比圖;圖9為GPOPS2和ILPSCP算法解的節(jié)點分布圖。

圖4 飛行軌跡對比圖

圖5 速度曲線對比圖

圖6 彈道傾角曲線對比圖

圖7 攻角曲線對比圖

圖8 目標函數(shù)收斂曲線圖

圖9 離散節(jié)點分布對比圖

圖4～圖7中,給出了傳統(tǒng)LPSCP算法第199次和第200次迭代的結(jié)果,可以看出,這兩組曲線都沒能和參照組曲線重合,即沒有收斂到最優(yōu)解。結(jié)合圖8可以看出,LPSCP算法對應的目標函數(shù),僅需20次迭代便趨于一個值附近,并在該值附近震蕩。但在經(jīng)過50次、100次甚至150次迭代之后,該目標函數(shù)曲線仍然持續(xù)震蕩,難以收斂,最后在0.607和0.690之間持續(xù)跳變?？梢灶A計,即使進行更多迭代仍然難以收斂。形成鮮明對比的是,圖8中ILPSCP算法對應的目標函數(shù)曲線在經(jīng)過30次迭代后便基本穩(wěn)定,之后變化量極小。

由圖4～圖7和圖9可以看出,ILPSCP算法的節(jié)點和GPOPS2的節(jié)點雖然不同,但ILPSCP算法規(guī)劃出的各變量曲線與參照組GPOPS2高度重合,這證明本文提出的ILPSCP算法是有效的。

本文的彈道規(guī)劃模型的性能指標為攻角平方積分,可以發(fā)現(xiàn),圖7中LPSCP算法的2條攻角曲線的平方積分比其他2種算法的更小,從圖8中也可以看出,LPSCP算法的性能指標小于ILPSCP算法,這與本文描述的ILPSCP算法優(yōu)于LPSCP算法這一觀點不符。出現(xiàn)這種現(xiàn)象的原因是:LPSCP算法使用固定相對信賴域=03,較大的信賴域?qū)挾葘τ诤唵螁栴}影響不大,但對于本文較復雜的彈道規(guī)劃問題,會造成線性近似假設(shè)條件不成立,近似誤差大。過大的線性近似誤差導致近似模型無法反映制導炮彈物理模型,所以,對于制導炮彈模型,圖7中LPSCP算法的攻角曲線和圖8中LPSCP算法的目標函數(shù)曲線都是無效的。過大的線性近似誤差,也是導致LPSCP算法無法收斂的原因。采用=03的LPSCP算法在2個無效解之間跳變,而=03的ILPSCP算法收斂到最優(yōu)解,顯然ILPSCP算法比LPSCP算法更適合解復雜的彈道規(guī)劃問題。

LPSCP算法在=03時便因為信賴域?qū)挾冗^大無法求解,而ILPSCP算法在更大的初始相對信賴域?qū)挾惹闆r下,仍然能穩(wěn)定求解。增大后,ILPSCP算法收斂到指定精度(最優(yōu)解與=03時相同)所需迭代次數(shù),如表5所示。隨著初始相對信賴域的增大,所需迭代次數(shù)增加,但相比信賴域?qū)挾鹊脑鲩L倍數(shù),其增長幅度較小。相比于傳統(tǒng)的LPSCP算法,即使設(shè)置了更大的初始相對信賴域?qū)挾?ILPSCP算法仍然能快速穩(wěn)定收斂,這進一步證明了ILPSCP算法的優(yōu)越性。

表5 不同δ0下ILPSCP算法迭代次數(shù)

5 結(jié)束語

本文推導了時間自由且以控制能量最優(yōu)為指標的彈道規(guī)劃一般化模型,針對該模型采用偽譜離散和線性凸化完成凸優(yōu)化建模。針對傳統(tǒng)L1序列凸規(guī)劃算法(LPSCP)在解決復雜的初末時間自由、非線性強的彈道規(guī)劃問題時存在穩(wěn)定性較差等問題,本文提出了一種改進的L1偽譜凸規(guī)劃算法(ILPSCP)。ILPSCP算法中引入指數(shù)衰減的相對信賴域?qū)挾群蛶辖绲腖1懲罰系數(shù)。仿真結(jié)果證明,ILPSCP算法收斂速度快,精度高,冷啟動能力強,很好地解決了傳統(tǒng)L1序列凸規(guī)劃算法(LPSCP)存在的問題。