陳 琦，王中原，常思江

（南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，南京210094）

隨著戰(zhàn)爭(zhēng)理念的進(jìn)步和發(fā)展，特別是現(xiàn)代非對(duì)稱作戰(zhàn)中“遠(yuǎn)程壓制、精確打擊、高效毀傷”等概念的提出，傳統(tǒng)火炮射程不遠(yuǎn)、精度不高等缺點(diǎn)逐漸暴露，已不能適應(yīng)現(xiàn)代戰(zhàn)爭(zhēng)的需求。在這種背景下，滑翔增程制導(dǎo)炮彈以其反應(yīng)速度快、效費(fèi)比高、使用靈活并且可打擊靜止和運(yùn)動(dòng)目標(biāo)等優(yōu)勢(shì)逐漸受到各國(guó)的重視。滑翔増程彈由火炮發(fā)射后，尾翼張開(kāi)保持穩(wěn)定飛行，之后火箭發(fā)動(dòng)機(jī)啟動(dòng)并開(kāi)始助推，發(fā)動(dòng)機(jī)工作結(jié)束后，彈丸像普通尾翼彈一樣繼續(xù)在升弧段上飛行，到達(dá)彈道頂點(diǎn)附近，鴨舵展開(kāi)，開(kāi)始對(duì)彈丸實(shí)施控制，最終命中目標(biāo)。文獻(xiàn)［1－4］闡述了滑翔増程彈的飛行原理；文獻(xiàn)［5］研究了制導(dǎo)炮彈最優(yōu)滑翔彈道與控制問(wèn)題，分別采用龐特里亞金極大值原理和共軛梯度法求解了控制變量有約束的最優(yōu)滑翔彈道；文獻(xiàn)［6］對(duì)滑翔彈道分3段進(jìn)行了分析，研究了滑翔増程炮彈的彈道模型，并對(duì)滑翔彈道進(jìn)行了設(shè)計(jì)。文獻(xiàn)［7］以最大升阻比的方法對(duì)滑翔彈道進(jìn)行了優(yōu)化設(shè)計(jì)；文獻(xiàn)［8］分別采用了序列二次規(guī)劃法和最大升阻比法優(yōu)化了俯仰舵偏角規(guī)律。以上文獻(xiàn)主要以射程最遠(yuǎn)為性能指標(biāo)優(yōu)化方案彈道，這對(duì)于考察彈丸的滑翔能力、打擊距離均具有十分重要的意義。但是這種方法也往往使得優(yōu)化出的彈道飛行時(shí)間過(guò)長(zhǎng)，易被敵方防御系統(tǒng)發(fā)現(xiàn)和攔截；存速較低，削弱了末端攻擊的機(jī)動(dòng)性。這在很大程度上限制了制導(dǎo)炮彈快速精確打擊、及時(shí)火力支援等能力的發(fā)揮。此外，上述研究均假設(shè)攻角能瞬時(shí)響應(yīng)指令信號(hào)，導(dǎo)致設(shè)計(jì)出的方案彈道過(guò)于理想，給實(shí)際應(yīng)用帶來(lái)了較大的難度。針對(duì)以上問(wèn)題，本文以飛行時(shí)間為性能指標(biāo)對(duì)滑翔彈道進(jìn)行了優(yōu)化設(shè)計(jì)，考慮了動(dòng)力學(xué)滯后，將其簡(jiǎn)化為一階慣性環(huán)節(jié)，引入了虛擬控制量，并將其作為優(yōu)化變量，同時(shí)為了保證攻擊效果，對(duì)末端存速、落角進(jìn)行一定的限制，在縱向平面內(nèi)建立了彈道優(yōu)化模型；引入了偽譜理論對(duì)模型進(jìn)行離散，這在很大程度上降低了彈道優(yōu)化耗時(shí)，縮短了滑翔彈的作戰(zhàn)反應(yīng)時(shí)間，同時(shí)也為今后實(shí)現(xiàn)算法的實(shí)時(shí)性提供了理論依據(jù)。

1 彈道優(yōu)化問(wèn)題描述

1.1 滑翔増程彈運(yùn)動(dòng)方程

為研究問(wèn)題的本質(zhì)，本文只考慮炮彈在縱向平面內(nèi)的運(yùn)動(dòng)，并將動(dòng)力學(xué)滯后簡(jiǎn)化為一階慣性環(huán)節(jié)，大氣為標(biāo)準(zhǔn)氣象條件。在這種情況下，滑翔増程制導(dǎo)炮彈的滑翔段運(yùn)動(dòng)方程可簡(jiǎn)化為

Fx，F(xiàn)y分別為阻力和升力，有：

式中：θ為彈道傾角，q為動(dòng)壓，Sref為參考面積，k和k′分別為彈翼組合體和舵面的誘導(dǎo)阻力系數(shù)，Cx0BW和Cx0δ分別為彈翼組合體和舵面的零升阻力系數(shù)，和Cδz

y分別為彈翼組合體和舵面的升力系數(shù)導(dǎo)數(shù)，以上各參數(shù)的值由吹風(fēng)實(shí)驗(yàn)獲得。αB和δzB分別為平衡攻角和平衡舵偏角，Td為一階慣性環(huán)節(jié)的時(shí)間常數(shù)，uα為虛擬控制量，ε（uα）＝0為虛擬量uα的控制方程。方程中其它變量的意義參見(jiàn)文獻(xiàn)［9］。

1.2 彈道優(yōu)化模型

為了實(shí)現(xiàn)快速打擊的要求，滑翔増程制導(dǎo)炮彈在彈道最高點(diǎn)處張開(kāi)鴨舵，對(duì)彈體實(shí)施控制，經(jīng)過(guò)滑翔飛行，在最短的時(shí)間內(nèi)命中目標(biāo)，同時(shí)，為了保證對(duì)目標(biāo)的毀傷效果，還要對(duì)彈道末端速度和落角進(jìn)行一定的限制。因此，滑翔彈道優(yōu)化問(wèn)題可以具體描述為：在［0，tf］時(shí)間內(nèi)（tf為未知的彈道末端時(shí)刻），確定出控制量uα，使得由式（1）確定的動(dòng)態(tài)系統(tǒng)從初始狀態(tài)轉(zhuǎn)移到終端狀態(tài)，在滿足規(guī)定約束的條件下，使得參數(shù)tf最小。由以上分析可知，彈道優(yōu)化問(wèn)題實(shí)際上是一個(gè)動(dòng)態(tài)優(yōu)化問(wèn)題，其優(yōu)化模型可表示如下：

①性能指標(biāo)。為了實(shí)現(xiàn)快速打擊的要求，性能指標(biāo)取為攻擊時(shí)間最短，即：

②初始狀態(tài)為無(wú)控段末端狀態(tài)，本文取為彈道最高點(diǎn)處狀態(tài)。

③終端約束條件。為了保證最終的攻擊效果，滑翔段終端狀態(tài)應(yīng)滿足一定的約束，包括終端位置約束、速度約束和落角約束，即：

④攻角約束。考慮到彈體的控制性能，同時(shí)避免彈體失速，需要對(duì)攻角的幅值進(jìn)行約束，即：

⑤系統(tǒng)狀態(tài)方程為式（1）描述的彈體縱向平面運(yùn)動(dòng)模型。

可見(jiàn)，滑翔彈道優(yōu)化問(wèn)題是一個(gè)起始狀態(tài)固定，終端狀態(tài)受約束，終端時(shí)刻自由的非線性動(dòng)態(tài)優(yōu)化問(wèn)題。

2 Gauss偽譜法

作為一種求解最優(yōu)控制問(wèn)題的直接方法，Gauss偽譜法的主要思想是將連續(xù)的無(wú)限維最優(yōu)控制問(wèn)題轉(zhuǎn)換為非線性規(guī)劃問(wèn)題進(jìn)行求解，在一系列的Legendre－Gauss（LG）節(jié)點(diǎn)上將狀態(tài)變量和控制變量進(jìn)行離散，并且利用這些離散點(diǎn)構(gòu)造全局Lagrange插值多項(xiàng)式來(lái)近似系統(tǒng)的動(dòng)力學(xué)方程。該方法具有收斂速度快、對(duì)初值不敏感且無(wú)需猜測(cè)協(xié)態(tài)變量等優(yōu)點(diǎn)，2005年經(jīng)Benson D在其博士論文中介紹并完善后［10］，便受到了廣泛的關(guān)注和應(yīng)用。

2.1 系統(tǒng)動(dòng)力學(xué)方程離散

考慮積分形式的系統(tǒng)動(dòng)力學(xué)方程：

式中：X（t）∈Rn為狀態(tài)向量，U（t）∈Rm為控制向量，f：Rn×Rm×R→Rn為連續(xù)向量函數(shù)。采用Gauss偽譜法需要將時(shí)間區(qū)域［t0，tf］轉(zhuǎn)換到［－1，1］上，為此，引入變量τ對(duì)時(shí)間t進(jìn)行變換［6］，即：

因此，系統(tǒng)動(dòng)力學(xué)方程變?yōu)?/p>

由數(shù)值分析的知識(shí)可知，選取N個(gè)Lagrange插值基函數(shù)Lk（ζ）近似式（4）中的積分項(xiàng)，其代數(shù)精度可達(dá)2N－1次［11］，有：

這樣，結(jié)合式（4）和式（5），第i個(gè)節(jié)點(diǎn)處的狀態(tài)向量便可表示為

若記X0＝X（－1），Ui＝U（τi），Xi＝X（τi），1≤i≤N，根據(jù)式（6），便可將積分形式的系統(tǒng)動(dòng)力學(xué)方程在N個(gè)LG節(jié)點(diǎn)上轉(zhuǎn)換為代數(shù)方程，即：

式中：Wi為第i個(gè)節(jié)點(diǎn)處的積分權(quán)重，Wi＝

2.2 邊界條件和性能指標(biāo)離散

式（8）僅在區(qū)間的內(nèi)點(diǎn)計(jì)算狀態(tài)量，并未包含終端時(shí)刻節(jié)點(diǎn)，終端狀態(tài)應(yīng)滿足動(dòng)力學(xué)方程約束，即：

將終端條件離散并利用Gauss積分來(lái)近似，可得：

由于本文所選取的性能指標(biāo)為Mayer型指標(biāo)，不包含積分項(xiàng)。因此無(wú)需進(jìn)行特殊處理，即：

根據(jù)以上的數(shù)學(xué)變換，本文的問(wèn)題可以描述為：在［t0，tf］時(shí)間內(nèi)（tf為未知的彈道末端時(shí)刻），確定離散點(diǎn)上的狀態(tài)向量Xi（i＝1，2，…，N）、控制向量Ui（i＝1，2，…，N）和終端時(shí)刻tf，使得性能指標(biāo)（11）最小，并滿足系統(tǒng)動(dòng)力學(xué)方程約束（8）、終端狀態(tài)約束（10），以及原問(wèn)題的邊界條件和過(guò)程約束：

從而將連續(xù)的無(wú)限維最優(yōu)控制問(wèn)題轉(zhuǎn)換為一般的非線性規(guī)劃問(wèn)題，寫(xiě)成標(biāo)準(zhǔn)形式為

3 數(shù)值仿真算例

以某滑翔增程制導(dǎo)炮彈為對(duì)象進(jìn)行數(shù)值仿真。經(jīng)火箭發(fā)動(dòng)機(jī)助推后彈丸質(zhì)量m＝37.27kg，炮彈達(dá)到彈道最高點(diǎn)時(shí)的高度為Hp0＝14 746.6m，速度為vp0＝299.3m／s，彈道傾角為θp0＝0；滑翔段舵偏角約束為｜αB（t）｜≤15°，滑翔末端高度取為yf＝0，彈道傾角取為｜θ（tf）｜≥θf(wàn)＝89°，同時(shí)，為了保證攻擊效果，也要對(duì)末端存速進(jìn)行限制，v（tf）≥vf＝200m／s，一階慣性環(huán)節(jié)的時(shí)間常數(shù)Td取為1.5s。取彈道最高點(diǎn)處x0＝0，對(duì)目標(biāo)xT＝40km，50km，60km分別進(jìn)行仿真計(jì)算，LG節(jié)點(diǎn)個(gè)數(shù)取為35，優(yōu)化算法工具包為SNOPT，仿真結(jié)果如圖1～圖5所示。

從圖1可見(jiàn)，制導(dǎo)炮彈在經(jīng)過(guò)彈道最高點(diǎn)后，近乎以直線滑翔，當(dāng)快達(dá)到末端時(shí)，彈道開(kāi)始俯沖，這是由于為了保證最后的攻擊效果，對(duì)末端條件進(jìn)行約束所致。圖2描述了滑翔過(guò)程中速度的變化規(guī)律，可以發(fā)現(xiàn)，速度在中間階段變化非常緩慢，但當(dāng)接近彈道末端時(shí)開(kāi)始衰減，這是因?yàn)閺椡柙谀┒我暂^大的負(fù)攻角進(jìn)行俯沖，增大了飛行阻力，因此速度衰減較快。

圖1 不同射程情況下的彈道曲線

圖2 不同射程情況下的速度曲線

圖3 也很明顯地反映出彈道傾角在經(jīng)過(guò)短暫的減小后，在較長(zhǎng)的時(shí)間內(nèi)近乎保持不變，驗(yàn)證了彈丸近似直線飛行，當(dāng)接近彈道末端時(shí)，彈道傾角先增加后急劇減小，說(shuō)明彈丸開(kāi)始俯沖，且最終落角均大于89°，滿足設(shè)計(jì)指標(biāo)。

圖3 不同射程情況下的彈道傾角曲線

圖4 為不同射程條件下虛擬控制量uα和實(shí)際控制量αB的對(duì)比曲線，可以看出，因?yàn)榭紤]了動(dòng)力學(xué)滯后，為了使攻角迅速增加，uα在初始時(shí)刻的值較大，之后慢慢減小，逐漸與αB重合。由于在構(gòu)建優(yōu)化模型時(shí)對(duì)末端彈道傾角進(jìn)行了約束，｜θ（tf）｜≥θf(wàn)，為了滿足這種約束，升力需在彈道末端變?yōu)樨?fù)值，根據(jù)升力計(jì)算公式，此時(shí)攻角也應(yīng)該為負(fù)值，且θf(wàn)越大，攻角的絕對(duì)值也越大，從而攻角曲線在末端出現(xiàn)了圖4所示的陡峭現(xiàn)象。此外，由于動(dòng)力學(xué)滯后的存在，uα在末端的下降過(guò)程中要超前于αB。

圖4 不同射程情況下的實(shí)際控制量和虛擬控制量曲線

圖5 展示了考慮動(dòng)力學(xué)滯后的攻角改善效果?？梢悦黠@地看出，考慮動(dòng)力學(xué)滯后的攻角變化較為平緩，以xT＝60km為例，不考慮動(dòng)力學(xué)滯后時(shí)，攻角在初始時(shí)刻從0瞬間增加到接近6°，這種突變往往會(huì)給今后控制系統(tǒng)的準(zhǔn)確跟蹤帶來(lái)較大的困難。而在考慮動(dòng)力學(xué)滯后時(shí)，這種現(xiàn)象便不再出現(xiàn)，優(yōu)化得到的攻角曲線變化非常平緩。

圖5 考慮動(dòng)力學(xué)滯后攻角改善效果

為了進(jìn)一步驗(yàn)證Gauss偽譜法的高效性，利用直接打靶法［14］在同樣的計(jì)算平臺(tái)上分別對(duì)xT＝40km，50km，60km的彈道進(jìn)行了優(yōu)化計(jì)算，離散控制量所需的節(jié)點(diǎn)設(shè)置為均勻分布，且個(gè)數(shù)取為35，積分步長(zhǎng)取為0.5s，對(duì)于優(yōu)化模型中的各種約束，采用罰函數(shù)法進(jìn)行處理，優(yōu)化結(jié)果對(duì)比見(jiàn)表1。從中可看出，直接打靶法優(yōu)化出的彈道飛行時(shí)間大于Gauss偽譜法，這說(shuō)明其尋優(yōu)效果要低于Gauss偽譜法。此外，直接打靶法的優(yōu)化結(jié)果并沒(méi)有完全滿足落角約束，可見(jiàn)，在終端位置、速度和落角等諸多約束的條件下，相比于Gauss偽譜法，直接打靶法的精度有所欠缺。最后，直接打靶法的優(yōu)化耗時(shí)tcost要遠(yuǎn)大于Gauss偽譜法，這是因?yàn)橹苯哟虬蟹ㄔ诿看蔚鷷r(shí)都會(huì)利用控制量積分彈道，雖然選取了較大的積分步長(zhǎng)（0.5s），但是大量的迭代還是產(chǎn)生了巨大的計(jì)算耗時(shí)，明顯地降低了算法的效率，而Gauss偽譜法無(wú)需在每次迭代中積分彈道，因此在計(jì)算效率上具有很大的優(yōu)勢(shì)。圖6為兩種方法優(yōu)化得到的攻角的對(duì)比圖，從中可看出二者基本吻合，但Gauss偽譜法的結(jié)果要更為光滑。從以上的分析可以看出，Gauss偽譜法不論在優(yōu)化效率還是精度上均強(qiáng)于直接打靶法。

表1 直接打靶法和Gauss偽譜法優(yōu)化結(jié)果對(duì)比

圖6 直接打靶法和Gauss偽譜法攻角對(duì)比曲線

4 結(jié)束語(yǔ)

本文基于Gauss偽譜法求解了遠(yuǎn)程制導(dǎo)炮彈滑翔彈道快速優(yōu)化問(wèn)題。以快速精確打擊為指標(biāo)對(duì)滑翔彈道進(jìn)行了優(yōu)化設(shè)計(jì)，考慮動(dòng)力學(xué)滯后，引入虛擬控制量，并將其作為優(yōu)化變量，在縱向平面內(nèi)建立了彈道優(yōu)化模型；利用Gauss偽譜法將優(yōu)化模型在一系列LG節(jié)點(diǎn)上離散，將連續(xù)的無(wú)限維最優(yōu)控制問(wèn)題轉(zhuǎn)換為一般的非線性規(guī)劃問(wèn)題。仿真結(jié)果表明，彈丸能以一定的速度和大著角攻擊目標(biāo)，所有約束條件均能較好地滿足，引入虛擬控制量后，優(yōu)化得到的攻角曲線更為平滑。與傳統(tǒng)的直接打靶法相比，Gauss偽譜法計(jì)算時(shí)間短，優(yōu)化效率高，具有在線優(yōu)化的潛力，促進(jìn)了制導(dǎo)炮彈快速精確打擊、及時(shí)火力支援等能力的有效發(fā)揮，可為遠(yuǎn)程制導(dǎo)炮彈的方案彈道優(yōu)化設(shè)計(jì)提供一定的參考。

［1］COSTELLO M.Extended range of a gun launched smart projectile using controllable canards［J］.Journal of Guidance，Control，and Dynamics，2006，29（6）：1 404－1 414.

［2］COSTELLO M.Extended range of a gun launched smart projectile using controllable canards［J］.Journal of Shock and Vibration，2001，8：203－213.

［3］史金光，王中原，易文俊.滑翔増程彈方案彈道特性研究［J］.彈道學(xué)報(bào)，2003，15（1）：51－54.SHI Jin－guang，WANG Zhong－yuan，YI Wen－jun.A study on the projectile trajectory characteristics of gliding extended range projectile［J］.Journal of Ballistics，2003，15（1）：51－54.（in Chinese）

［4］史金光，王中原，易文俊.滑翔増程彈飛行彈道［J］.火力與指揮控制，2007，32（11）：88－90.SHI Jin－guang，WANG Zhong－yuan，YI Wen－jun.Study of fight trajectory for gliding extended range projectile［J］.Fire Control and Command Control，2007，32（11）：88－90.（in Chinese）

［5］卜奎晨，劉莉.復(fù)合制導(dǎo)炮彈最優(yōu)滑翔彈道與控制［J］.彈道學(xué)報(bào)，2007，19（4）：23－25.BU Kui－chen，LIU Li.Optimal glide trajectory and control for combined guidance munition［J］.Journal of Ballistics，2007，19（4）：23－25.（in Chinese）

［6］汪小娜，王樹(shù)宗，朱華兵.火箭助推滑翔増程炮彈彈道模型研究［J］.艦船科學(xué)技術(shù)，2005，27（1）：77－79.WANG Xiao－na，WANG Shu－zong，ZHU Hua－bing.Study on ballistic models of rocket－aid gliding extended projectiles［J］.Ship Science and Technology，2005，27 （1）：77 － 79.（in Chinese）

［7］易文俊，王中原，史金光，等.帶鴨舵滑翔増程炮彈方案彈道研究［J］.南京理工大學(xué)學(xué)報(bào)，2008，32（3）：322－326.YI Wen－jun，WANG Zhong－yuan，SHI Jin－guang.Trajectory of glide extended－range projectile with canards configuration［J］.Journal of Nanjing University of Science and Technology，2008，32（3）：322－326.（in Chinese）

［8］史金光，王中原，涂泗華.求解滑翔増程彈較優(yōu)舵偏角方法［J］.彈箭與制導(dǎo)學(xué)報(bào)，2005，25（4）：717－719.SHI Jin－guang，WANG Zhong－yuan，TU Si－h(huán)ua.Calculating method on a suitable angle of lift－canards of gliding extended range projectile［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2005，25（4）：717－719.（in Chinese）

［9］錢(qián)杏芳，林瑞雄，趙亞男.導(dǎo)彈飛行力學(xué)［M］.北京：北京理工大學(xué)出版社，2008.QIAN Xing－fang，LIN Rui－xiong，ZHAO Ya－nan.Missile flight dynamics［M］.Beijing：Beijing Institute of Technology Press，2008.（in Chinese）

［10］BENSON D.A Gauss pseudospectral transcription for optimal control［D］.Cambridge：Massachusetts Institute of Technology，2005.

［11］關(guān)冶，陸金甫.數(shù)值分析基礎(chǔ)［M］.北京：高等教育出版社，2005.GUAN Ye，LU Jin－fu.Fundamentals of numerical analysis［M］.Beijing：High Education Press，2005.（in Chinese）

［12］GILL P E，MURRAY W，SAUNDERS M A.SNOPT：an SQP algorithm for large－scale constrained optimization［J］.SIAM Review，2005，47（1）：99－131.

［13］ACHTER A W，BIEGLER L T.On the implementation of a primal－dual interior point filter line search algorithm for largescale nonlinear programming［J］.Mathematical Programming，2006，106（1）：25－57.

［14］GERDTS M.Direct shooting method for the numerical solution of higher－index DAE optimal control problems［J］.Journal of Optimization Theory and Applications，2003，17（2）：267－294.