權(quán)申明，左 坤，晁 濤，楊 明

（哈爾濱工業(yè)大學(xué) 控制與仿真中心，哈爾濱 150080）

單個(gè)飛行器在對(duì)目標(biāo)探測(cè)與機(jī)動(dòng)能力上均存在一定的局限性，傳統(tǒng)單飛行器對(duì)目標(biāo)制導(dǎo)效果有限。相比之下，多飛行器協(xié)同探測(cè)目標(biāo)存在更大的優(yōu)勢(shì)，多飛行器可同時(shí)進(jìn)行探測(cè)識(shí)別與信息交流，顯著提高對(duì)真實(shí)目標(biāo)的識(shí)別能力。多飛行器協(xié)同還能降低飛行器的彈載計(jì)算和探測(cè)設(shè)備的要求，從而降低成本。

在協(xié)同末制導(dǎo)階段，相關(guān)學(xué)者針對(duì)不同約束條件進(jìn)行制導(dǎo)律設(shè)計(jì)。趙久奮等人[1]針對(duì)推力可控的導(dǎo)彈，在視線和視線法線方向上，結(jié)合超螺旋算法和有限時(shí)間滑模設(shè)計(jì)分布式協(xié)同制導(dǎo)律，并針對(duì)機(jī)動(dòng)目標(biāo)進(jìn)行仿真驗(yàn)證。肖惟等人[2]研究了考慮過載約束的多枚弱機(jī)動(dòng)能力導(dǎo)彈強(qiáng)目標(biāo)協(xié)同問題，結(jié)合飛行器可達(dá)域、可行域以及目標(biāo)逃逸域，設(shè)計(jì)了協(xié)同攔截策略。司玉潔等人[3]研究了考慮通訊拓?fù)錀l件下的協(xié)同末制導(dǎo)律，并基于有限時(shí)間理論設(shè)計(jì)快收斂的滑模制導(dǎo)律。此外，強(qiáng)化學(xué)習(xí)等方法也用于協(xié)同末制導(dǎo)律設(shè)計(jì)中，并取得較好的制導(dǎo)效果。

模型預(yù)測(cè)靜態(tài)規(guī)劃（Model Predictive Static Programming,MPSP）自Padhi 等人[4]提出以來，在飛行器制導(dǎo)控制[5-10]方面有較多的應(yīng)用，該方法在計(jì)算效率和制導(dǎo)精度上具有一定優(yōu)勢(shì)。MPSP 根據(jù)靈敏度矩陣不斷更新控制量，能夠滿足終端約束[11]，后續(xù)研究中，輸入與狀態(tài)約束逐漸引起關(guān)注[12,13]。文獻(xiàn)[14]考慮中制導(dǎo)過程中，零控?cái)r截和交班視窗角約束，推導(dǎo)了時(shí)間固定的廣義模型預(yù)測(cè)靜態(tài)規(guī)劃算法。文獻(xiàn)[15]基于優(yōu)化模型預(yù)測(cè)靜態(tài)規(guī)劃算法，針對(duì)吸氣式高超聲速飛行器設(shè)計(jì)了滿足終端碰撞角約束的制導(dǎo)律。文獻(xiàn)[16]提出了一種初、中制導(dǎo)聯(lián)合規(guī)劃制導(dǎo)方法，用于解決多階段目標(biāo)攔截問題。在MPSP 的基礎(chǔ)上，廣義MPSP 算法被提出[17,18]，該方法無需對(duì)動(dòng)力學(xué)模型離散化，提高算法求解速度與精度。

MPSP 算法同模型預(yù)測(cè)控制有相似之處，離散步長通常為定值，因此難以處理終端時(shí)間自由的問題。然而，在協(xié)同末制導(dǎo)問題中，終端時(shí)間難以事先獲得；同時(shí)，初始軌跡的選取影響了優(yōu)化結(jié)果與優(yōu)化效率。

針對(duì)現(xiàn)有序列二階錐規(guī)劃（Second-Order Cone Programming,SOCP）求解協(xié)同制導(dǎo)問題速度較慢難以滿足實(shí)時(shí)性要求的問題，本文結(jié)合MPSP 與SOCP，提出一種基于模型預(yù)測(cè)靜態(tài)二階錐規(guī)劃的協(xié)同末制導(dǎo)算法。首先，建立考慮過載約束的飛行器協(xié)同末制導(dǎo)模型；然后，在MPSP 算法基礎(chǔ)上，提出終端時(shí)間自由的MPSP 算法；進(jìn)一步將有狀態(tài)約束的MPSP 問題轉(zhuǎn)化為模型預(yù)測(cè)靜態(tài)二階錐規(guī)劃問題，可以有效處理過程約束。仿真結(jié)果表明，所提出的模型預(yù)測(cè)靜態(tài)二階錐規(guī)劃能兼?zhèn)銶PSP 算法的快速性與SOCP 的處理約束與最優(yōu)性特點(diǎn)，通過減少求解變量個(gè)數(shù)，提升求解速度，易于實(shí)現(xiàn)在線快速計(jì)算。通過不同仿真場(chǎng)景、不同優(yōu)化方法的仿真與分析，驗(yàn)證了本文所提算法的計(jì)算效率與精度。

1 問題描述

建立發(fā)射系下飛行器運(yùn)動(dòng)模型為：

其中，i為每個(gè)飛行器的編號(hào)，xi、yi、zi、Vi分別表示第i個(gè)飛行器的在空間坐標(biāo)系下的位置分量、速度大小，θi、σi、γi分別是第i個(gè)飛行器的彈道傾角、彈道偏角、傾側(cè)角。Di、Li分別是第i個(gè)飛行器在飛行過程中受到的阻力與升力，速度系下氣動(dòng)力分量為：

其中ρ=ρ0e-y/hs是大氣密度，ρ0為標(biāo)準(zhǔn)大氣密度，hs為標(biāo)準(zhǔn)高度。Sref為飛行器的最大橫截面積。CL和CD分別是升力與阻力系數(shù)，由飛行攻角α與飛行馬赫數(shù)M計(jì)算得到，馬赫數(shù)的數(shù)值大小由飛行器的飛行速度進(jìn)行計(jì)算。

在飛行器協(xié)同飛行過程中，需滿足初始狀態(tài)約束、終端狀態(tài)約束以及過載約束。其中，過載約束為：

2 算法設(shè)計(jì)

本節(jié)首先進(jìn)行終端時(shí)間固定、無過程約束MPSP算法介紹。在此基礎(chǔ)上，提出考慮過載約束的終端時(shí)間自由協(xié)同打擊末制導(dǎo)算法。

2.1 終端時(shí)間固定的MPSP 算法

記離散時(shí)間的非線性動(dòng)力學(xué)系統(tǒng)為：

其中，k=1,2,......N-1為一系列的離散點(diǎn)，離散步長為tN為設(shè)定的終端時(shí)間，x∈Rn是系統(tǒng)的狀態(tài)變量，u∈Rm是系統(tǒng)的控制變量，F(xiàn)關(guān)于x與u連續(xù)可導(dǎo)，與飛行器的動(dòng)力學(xué)方程有關(guān)，MPSP的目的是找到合適的du，對(duì)于先前的輸入up進(jìn)行校正（上角標(biāo)p表示上一次迭代結(jié)果），得到新的控制量u=up+du，狀態(tài)變量隨之得到更新。即：x=xp+dx。其中du、dx分別是控制變量與狀態(tài)變量的增量。

假設(shè)狀態(tài)變量的變化dx與控制變量的變化du較小，那么在更新k+1時(shí)刻的狀態(tài)時(shí)，可以對(duì)離散的狀態(tài)方程進(jìn)行泰勒展開，考慮一階展開形式：

同理可以得到：

將式(7)代入式(6)，得：

不斷迭代可得：

其中，

其中j=1,2…k-1。

當(dāng)j=k時(shí)，有：

初始狀態(tài)通常為給定值，那么dx1=0，則式(9)可寫為：

式(13)表明，第k+1個(gè)時(shí)間點(diǎn)的狀態(tài)是由前k個(gè)控制量進(jìn)行校正的，為簡(jiǎn)化B陣的計(jì)算，引入中間變量：

至此B計(jì)算完成，稱B為靈敏度矩陣。

2.2 基于二階錐規(guī)劃終端時(shí)間自由的MPSP

考慮終端時(shí)間自由的離散時(shí)間的動(dòng)力學(xué)系統(tǒng)方程：

其中，k表示第k個(gè)離散點(diǎn)，系統(tǒng)的非線性動(dòng)力學(xué)方程為f(xk,uk)，本文離散點(diǎn)總數(shù)N保持不變，終端時(shí)刻變化后，每一個(gè)離散點(diǎn)表示的時(shí)間將不斷改變，離散時(shí)間間隔tk為：

將模型(17)在前一次結(jié)果(xkp,ukp,tkp)的基礎(chǔ)上進(jìn)行線性化，可以得到：

其中，終端時(shí)間的變化量為dtN=tN-tNp，每個(gè)離散時(shí)間間隔的變化量為表示前一次迭代得到的終端時(shí)間。

滿足上述式子的du很多，難以獲得解析解，因此使用優(yōu)化的思想來進(jìn)行求解。

構(gòu)造目標(biāo)函數(shù)為：

其中，L為過程中的優(yōu)化目標(biāo)，φ為終端時(shí)刻的優(yōu)化目標(biāo)。本文不僅考慮終端約束的影響，還要考慮過程約束。過程約束與控制約束線性化后可以寫成如下形式：

其中，Lk、Ck與具體的約束有關(guān)，已知代入式(22)可以得到：

由上述不等式可知，上述約束為一個(gè)凸約束，可以使用二階錐規(guī)劃方法解決。

綜上，初始狀態(tài)為x1*，期望終端狀態(tài)為xf*，可以將求控制變量du與終端時(shí)間tN的問題轉(zhuǎn)化為一個(gè)模型預(yù)測(cè)靜態(tài)二階錐規(guī)劃（MPSP&SOCP）問題：

由于在MPSP&SOCP 中，優(yōu)化變量只有終端時(shí)間和輸入量的偏差，因此無法直接計(jì)算迭代過程中的過載大小，采用一階泰勒展開近似線性化后，通過逐級(jí)求導(dǎo)可得到攻角對(duì)速度、高度的影響，進(jìn)而通過序列迭代求出過載值，具體如式(25)：

設(shè)計(jì)MPSP&SOCP 算法求解流程如圖1所示。

圖1 MPSP&SOCP 算法流程圖Fig.1 MPSP&SOCP algorithm flow chart

MPSP&SOCP 算法需要設(shè)置初始軌跡，傳統(tǒng)MPSP 求解末制導(dǎo)問題時(shí)，往往通過比例導(dǎo)引進(jìn)行初始軌跡的計(jì)算。本文中，為了降低算法對(duì)初值的依賴性，提高算法的適應(yīng)性，采用常值控制量，開環(huán)積分后得到初值。在求解帶有過程約束最優(yōu)問題時(shí)，第一次迭代時(shí)，僅加入終端位置約束，在求得結(jié)果基礎(chǔ)上，增加終端角度約束。通過將約束條件逐步加入，可以有效提升求解速度與算法穩(wěn)定性。計(jì)算靈敏度矩陣后，將過程約束轉(zhuǎn)化為SOCP 問題，僅需求解控制量的更新值即可，無需進(jìn)行狀態(tài)量的優(yōu)化，大大減少了求解時(shí)間。

3 數(shù)值仿真及分析

為了驗(yàn)證本文所提算法的有效性與適應(yīng)性，本節(jié)分別進(jìn)行無過程約束MPSP 末制導(dǎo)仿真、考慮側(cè)窗約束MPSP 末制導(dǎo)仿真、對(duì)比算法仿真以及適應(yīng)性仿真實(shí)驗(yàn)。其中，選取文獻(xiàn)[19]中SOCP 算法作為對(duì)比算法。

飛行器的升力系數(shù)與阻力系數(shù)是關(guān)于攻角和馬赫數(shù)的多項(xiàng)式，如式(26)所示：

協(xié)同飛行器數(shù)量為3（M1、M2、M3），飛行器質(zhì)量為907.8 kg，面積為0.355 m2，最大過載約束為1.5g，仿真中初始條件如表1所示。

表1 仿真初始條件Tab.1 Initial conditions of simulation

3.1 二階錐規(guī)劃方法仿真結(jié)果

在相同仿真條件下，采用二階錐規(guī)劃方法進(jìn)行對(duì)比仿真，仿真結(jié)果如圖2-7所示（圖中Initial 表示第一次求解得到的軌跡，Opt 表示最終優(yōu)化軌跡）。

圖4 速度傾角變化曲線Fig.4 Flight path angles

圖5 攻角變化曲線Fig.5 Angles of attack

圖 傾 化曲Fig.6 Curves of bank angles

由仿真結(jié)果可以看出，基于二階錐規(guī)劃方法的協(xié)同末制導(dǎo)可以實(shí)現(xiàn)對(duì)目標(biāo)的精確打擊。該算法首先加入終端位置約束，因此圖7中第一次求解結(jié)果中峰值過載接近2.4g，然后考慮終端速度傾角約束，在上述結(jié)果基礎(chǔ)上，考慮過載約束能夠?qū)崿F(xiàn)收斂，同時(shí)過程約束滿足要求。

圖7 過載變化曲線Fig.7 Curves of overload

3.2 MPSP&SOCP 算法仿真結(jié)果

基于MPSP&SOCP 算法仿真結(jié)果如圖8-14所示（圖中Guess 表示初始猜測(cè)軌跡，Opt 表示計(jì)算得到的最優(yōu)軌跡）。

圖8 三維軌跡Fig.8 Curves of 3D trajectories

圖9 速度變化曲線Fig.9 Curves of speed

圖10 速度傾角變化曲線Fig.10 Flight path angles

圖11 攻角變化曲線Fig.11 Angles of attack

圖12 傾側(cè)角變化曲線Fig.12 Curves of bank angles

圖13 過載變化曲線Fig.13 Curves of overload

圖14 相對(duì)距離變化曲線Fig.14 Curves of relative distances

由圖8可以看出，3 個(gè)飛行器由猜測(cè)控制量開環(huán)積分得到的軌跡不滿足終端位置約束，飛行器1 和飛行器3 相距目標(biāo)點(diǎn)較遠(yuǎn)，3 個(gè)飛行器終端高度為負(fù)，圖11和圖12中開環(huán)控制量相同，因此各飛行器初始軌跡相似。在第一次迭代中，僅加入終端位置約束，從第二次迭代開始，加入終端彈道傾角約束如圖10。由于末制導(dǎo)階段飛行器高度下降較快，阻力減速效果沒有重力勢(shì)能轉(zhuǎn)化為動(dòng)能的效果明顯，因此圖9中速度略有增加。經(jīng)過MPSP&SOCP 迭代計(jì)算后，能夠?qū)崿F(xiàn)對(duì)目標(biāo)的精確打擊。

本仿真算例中，過載約束為1.5g，初始軌跡中由于控制量變化率為常值，因此過載較小，最大過載為1.1g，在優(yōu)化中考慮過載約束，求得圖13中最優(yōu)軌跡滿足過載約束。

3.3 對(duì)比分析

將兩種方法得到的最優(yōu)軌跡進(jìn)行對(duì)比，如圖15所示?？梢钥闯龌赟OCP 與MPSP&SOCP 的求解三維軌跡相吻合。

圖15 不同方法求解結(jié)果Fig.15 Results of different methods

對(duì)終端位置約束進(jìn)行小范圍改變，變化范圍為原點(diǎn)附近水平面內(nèi)1000 m 的半徑圓區(qū)域，進(jìn)行50 次蒙特卡洛仿真。統(tǒng)計(jì)終端時(shí)刻脫靶量，如圖16所示?？梢钥闯鰞煞N方法求解精度均在2 m 以內(nèi)。

圖16 求解精度對(duì)比Fig.16 Solving accuracy comparison

將上述MPSP&SOCP 算法與二階錐規(guī)劃的求解效率進(jìn)行對(duì)比仿真分析，結(jié)果如圖17所示，可以看出在同等迭代次數(shù)、求解精度情況下，MPSP&SOCP算法能夠節(jié)省40%的時(shí)間，有望實(shí)現(xiàn)在線實(shí)時(shí)規(guī)劃求解。

圖17 求解效率對(duì)比Fig.17 Solving efficiency comparison

綜上，基于MPSP&SOCP 算法的協(xié)同末制導(dǎo)可以在保障求解精度的前提下，耗費(fèi)更短的時(shí)間，易于實(shí)現(xiàn)在線計(jì)算，工程前景廣闊。

4 結(jié) 論

本文針對(duì)現(xiàn)有序列二階錐規(guī)劃求解協(xié)同打擊問題速度較慢難以滿足實(shí)時(shí)性要求的問題，結(jié)合MPSP 與SOCP，提出一種基于模型預(yù)測(cè)靜態(tài)二階錐規(guī)劃的協(xié)同末制導(dǎo)算法。主要結(jié)論如下：

(1) 在MPSP 算法基礎(chǔ)上，提出終端時(shí)間自由的MPSP 算法，能夠解決傳統(tǒng)MPSP 問題終端時(shí)間固定的缺陷，求解時(shí)間最優(yōu)問題。

(2) 為了解決傳統(tǒng)MPSP 算法過于依賴初始軌跡的局限性，提出逐步增加約束條件的迭代求解流程，從而使得小偏差線性化更具合理性，算法穩(wěn)定性更強(qiáng)。

(3) 所提出的模型預(yù)測(cè)靜態(tài)二階錐規(guī)劃能兼?zhèn)銶PSP算法的快速性與SOCP 的處理約束與最優(yōu)性特點(diǎn)，通過減少求解變量個(gè)數(shù)，提升求解速度，易于實(shí)現(xiàn)在線快速計(jì)算。

本文目前工作從軌跡優(yōu)化的角度進(jìn)行，后續(xù)將開展在線實(shí)時(shí)優(yōu)化-跟蹤的閉環(huán)制導(dǎo)算法研究，同時(shí)考慮更多的約束條件，以增加本文算法的適用性。