郭尚生，楊榮軍，王良明

（1.遼沈工業(yè)集團(tuán)有限公司，沈陽110045；2.中國工程物理研究院 總體工程研究所，四川 綿陽621900；3.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，南京210094）

1 制導(dǎo)炮彈增程方案彈道優(yōu)化問題

制導(dǎo)炮彈在進(jìn)入末制導(dǎo)之前將沿預(yù)先設(shè)定的彈道飛行，為了使炮彈具有更遠(yuǎn)的打擊能力，需要對方案彈道進(jìn)行規(guī)劃。預(yù)先設(shè)計(jì)的方案彈道需要考慮對炮彈運(yùn)動(dòng)有重要影響的重力、空氣動(dòng)力等因素，并且為了使炮彈在各種干擾情況下依然能夠跟蹤該方案彈道，需要在彈道優(yōu)化時(shí)限定攻角控制量的幅值，即在規(guī)劃彈道時(shí)使制導(dǎo)炮彈有一定的控制余量。

1.1 制導(dǎo)炮彈增程彈道優(yōu)化性能指標(biāo)函數(shù)

考慮制導(dǎo)炮彈的無動(dòng)力增程彈道優(yōu)化問題，狀態(tài)變量主要包括飛行速度、彈道傾角、飛行高度和射程，控制變量取飛行中的攻角。針對增程優(yōu)化的設(shè)計(jì)要求，為使射程最大化，并考慮到終端約束條件，選擇的性能指標(biāo)函數(shù)為

式中：t0，tf分別為優(yōu)化段起始時(shí)刻、終端時(shí)刻；xg（t0），xg（tf）分別為起始時(shí)刻、終端時(shí)刻的位置射程分量。

1.2 制導(dǎo)炮彈的運(yùn)動(dòng)模型

制導(dǎo)炮彈的飛行彈道主要由初始條件和作用在彈上的力決定，并且描述制導(dǎo)炮彈運(yùn)動(dòng)的各狀態(tài)變量是相互關(guān)聯(lián)的，受到飛行力學(xué)建立的運(yùn)動(dòng)模型的約束。在控制系統(tǒng)理想的情況下，忽略彈體繞質(zhì)心的短周期姿態(tài)運(yùn)動(dòng)過程，制導(dǎo)炮彈縱向平面質(zhì)點(diǎn)彈道模型為

式中：m為制導(dǎo)炮彈質(zhì)量；v為制導(dǎo)炮彈速度，S為彈體特征面積；θ為彈道傾角；q＝0.5ρv2，為動(dòng)壓；Cx，Cy分別為制導(dǎo)炮彈氣動(dòng)阻力系數(shù)和升力系數(shù)，主要與Ma和攻角α相關(guān)。在某一Ma下Cx，Cy可分別用關(guān)于攻角的二次多項(xiàng)式模型和線性模型近似：Cx其中：Cx0，Cx2分別為零升阻力系數(shù)和誘導(dǎo)阻力系數(shù)為升力系數(shù)導(dǎo)數(shù)。

1.3 最優(yōu)控制問題描述

本質(zhì)上求解滿足式（2）和各種飛行狀態(tài)約束的制導(dǎo)炮彈彈道尋優(yōu)是一種動(dòng)態(tài)優(yōu)化問題，以攻角為控制量u（t），使下述目標(biāo)函數(shù)最?。?/p>

式中：x（t）是系統(tǒng)的狀態(tài)變量；t是時(shí)間變量；標(biāo)量目標(biāo)函數(shù)J（x，u，t）由終端懲罰項(xiàng)Φ（x（tf））和積分項(xiàng)L（x，u，t）構(gòu)成；t0，tf分別為積分項(xiàng)的積分下限和上限＝f（x，u，t）是系統(tǒng) 狀 態(tài) 方 程；c（x，u，t）是與狀態(tài)、控制輸入相關(guān)的等式約束函數(shù)；d（x，u，t）是與狀態(tài)、控制輸入相關(guān)的不等式約束函數(shù)；x（t0）＝x0是系統(tǒng)初始狀態(tài)；Y（x（tf）），G（x（tf））分別是與狀態(tài)相關(guān)的終端等式約束函數(shù)和終端不等式約束函數(shù)。

由表2可得出：在陰極極化條件下，該材料——介質(zhì)體系在自腐蝕電壓（-658mV）的電位下斷裂壽命最長，其應(yīng)力腐蝕敏感性最小、抗拉強(qiáng)度最大、應(yīng)變量也最大。在比自腐蝕電壓更正的-1 300、-1 000mV的電位下，也能激發(fā)10＃鋼的應(yīng)力腐蝕開裂。

2 變時(shí)間尺度的控制向量參數(shù)化

2.1 非均勻控制向量參數(shù)化

最優(yōu)控制問題（OCP）式（3）所描述的是一個(gè)兩點(diǎn)邊值問題，直接求解十分困難。目前，該問題的數(shù)值求解常采用參數(shù)化方法：先將控制時(shí)域離散化為等長的時(shí)間段，然后在優(yōu)化時(shí)域內(nèi)對系統(tǒng)模型積分來消除OCP中的系統(tǒng)模型動(dòng)態(tài)約束，轉(zhuǎn)化為一般的非線性規(guī)劃問題（NLP），進(jìn)而采用常用的NLP方法求解[8－11]?？刂葡蛄繀?shù)化方法是利用分段函數(shù)近似連續(xù)控制作用，通常采用分段常函數(shù)的形式，即每個(gè)時(shí)間段的控制量為常值。

將固定時(shí)間區(qū)間[t0，tf]分成N個(gè)時(shí)間段，該分段方法對應(yīng)有N＋1個(gè)時(shí)間節(jié)點(diǎn)，每段時(shí)間內(nèi)的控制量為ui，i＝0，1，…，N。傳統(tǒng)的方法是將各時(shí)間分段均勻化，即每個(gè)節(jié)點(diǎn)間等距，各時(shí)間段的長度為（tf－t0）/N。若考慮到分段的影響，那么應(yīng)該采用非均勻的分段方法，除了初始和終端時(shí)間節(jié)點(diǎn)外，每個(gè)時(shí)間節(jié)點(diǎn)ti都需要進(jìn)行優(yōu)化。控制變量經(jīng)參數(shù)化處理后，式（3）的最優(yōu)控制問題即可轉(zhuǎn)化為含有N－1個(gè)時(shí)間節(jié)點(diǎn)和N個(gè)控制參數(shù)的NLP問題。

引入標(biāo)準(zhǔn)化的時(shí)間變量τ，在時(shí)間段[ti，ti＋1]上定義dτ，滿足條件：

即，pi為時(shí)域[ti，ti－1]到標(biāo)準(zhǔn)化時(shí)域[τi，τi－1]的尺度因子，原控制時(shí)域[t0，tf]就轉(zhuǎn)化為標(biāo)準(zhǔn)化時(shí)域[0，1]，稱式（5）的變換為時(shí)間尺度變換。對式（5）積分可得到：

因此，每個(gè)標(biāo)準(zhǔn)化的時(shí)間段長度均為

這樣就可以用傳統(tǒng)的均勻控制向量參數(shù)法來解決非均勻參數(shù)法的問題。

定義函數(shù)：

且滿足約束條件則原時(shí)域[t0，tf]上的連續(xù)的無限維 OCP式（3）即轉(zhuǎn)化為標(biāo)準(zhǔn)化時(shí)域[0，1]內(nèi)的有限維優(yōu)化問題：

式（9）中的各標(biāo)準(zhǔn)化時(shí)間節(jié)點(diǎn)τi是等間距的，這樣每個(gè)時(shí)間子區(qū)間的狀態(tài)轉(zhuǎn)變?yōu)槌踔祮栴}求解，可以采用均勻參數(shù)化的方法來求解。

為了有效地計(jì)算目標(biāo)函數(shù)值，引入一個(gè)新變量l（τ），并 使 其 滿 足＝p（tf－t0）L（x，u，τ），初值l（0）＝0。式（9）的系統(tǒng)微分方程、參數(shù)和初始狀態(tài)已知，系統(tǒng)的后續(xù)狀態(tài)實(shí)際上只由控制輸入u（τ）確定，即x（τ）＝（ui，τ）。對 方 程＝p（tf－t0）×L（x，u，τ）和＝f（x，u，τ）在時(shí)間區(qū)間[0，1]數(shù)值積分可獲得系統(tǒng)狀態(tài)（ui，τ）和目標(biāo)函數(shù)J。將描述標(biāo)準(zhǔn)化時(shí)間分段的參數(shù)pi寫成向量形式p＝（p1p2…pN）T，則可取k維的控制向量ui（i＝1，2，…，N）和p作為待優(yōu)化的參數(shù)，用一個(gè)kN＋N維的參數(shù)向量表示：

通過上述變換，用離散的控制量在優(yōu)化時(shí)域內(nèi)對系統(tǒng)狀態(tài)方程積分，從而將動(dòng)態(tài)系統(tǒng)模型的約束完全消去，式（9）轉(zhuǎn)化成通用、易求解的有限維參數(shù)NLP問題：

2.2 參數(shù)化的制導(dǎo)炮彈最大射程彈道優(yōu)化問題

制導(dǎo)炮彈滑翔彈道優(yōu)化是一個(gè)初始狀態(tài)已知，但終端時(shí)刻和終端狀態(tài)未知或部分未知的過程優(yōu)化問題。由于制導(dǎo)炮彈飛行的終端時(shí)刻tf一般是未知的，這樣給微分形式彈道模型的積分計(jì)算帶來不便。設(shè)飛行時(shí)間ts＝tf－t0，將ts也作為待優(yōu)化的參數(shù)，可將彈道優(yōu)化問題轉(zhuǎn)化成如下形式：

用離散化的攻角在彈道優(yōu)化時(shí)域內(nèi)對彈道方程積分，轉(zhuǎn)化為形如式（11）的NLP問題，待優(yōu)化參數(shù)為＝（α1…αNp1…pNts）T。針對轉(zhuǎn)化后的NLP問題，可利用適于約束優(yōu)化問題的序列二次規(guī)劃方法[12]高效地求解。

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

設(shè)某制導(dǎo)炮彈的基本參數(shù)、初始發(fā)射條件和約束條件：火炮炮口動(dòng)能E＝1.34×107J，火箭發(fā)動(dòng)機(jī)點(diǎn)火前的炮彈質(zhì)量m0＝40 kg；火箭助推劑質(zhì)量mp＝4kg，比沖I＝2 200 Ns/kg，點(diǎn)火時(shí)刻tp＝5s，工作時(shí)間tpw＝3s，被動(dòng)段制導(dǎo)炮彈的質(zhì)量m＝m0－mp；制導(dǎo)炮彈的初始發(fā)射角θ0＝55°；增程彈道優(yōu)化中飛行攻角約束αmin＝－8°，αmax＝8°；終端時(shí)刻高度約束yg，tf＝2km；終端時(shí)刻的傾角約束θtf＝－20°；終端時(shí)刻的最小速度vtf＝150m/s。

按初始條件對炮彈無控段彈道解算，得到彈道頂點(diǎn)的時(shí)刻t0＝42.91s，飛行速度vt0＝234.35m/s，射程xg，t0＝13.140km，高度yg，t0＝11.524km。利用非均勻參數(shù)化方法與序列二次規(guī)劃相結(jié)合對降弧段彈道進(jìn)行優(yōu)化。為便于對比觀察，將增程段的彈道按時(shí)間分割為N＝25，即有25個(gè)攻角控制量參數(shù)，加上標(biāo)準(zhǔn)化的時(shí)間段參數(shù)和終端時(shí)間參數(shù)，共有2N＋1，即51個(gè)優(yōu)化參數(shù)。

圖1為對應(yīng)于每個(gè)時(shí)間段的無量綱攻角參數(shù)，該無量綱參數(shù)乘以最大攻角幅值即為攻角。圖中i的含義和式（4）相同。圖2為無量綱的時(shí)間分段參數(shù)，該無量綱參數(shù)通過式（5）的時(shí)間尺度變換，即乘以參數(shù)ts/N，可得到各段的時(shí)間長度（優(yōu)化的時(shí)間總長度ts＝180.18s）。通過觀察時(shí)間分段參數(shù)可以看出，非均勻分段后的各時(shí)間長度是不同的，這也是與時(shí)間均勻分段的本質(zhì)區(qū)別。圖3（a）為無控彈道與增程彈道的對比曲線，圖3（b）為對比曲線的局部放大圖，圖4、圖5分別為無控彈道與增程彈道的速度對比曲線和彈道傾角對比曲線，圖6為時(shí)間均勻分段與時(shí)間非均勻分段的攻角對比曲線。

圖1 無量綱的攻角參數(shù)

圖2 無量綱的時(shí)間分段參數(shù)

圖3 無控段彈道與增程彈道對比曲線

圖3顯示以增程彈道優(yōu)化的終端條件yg，tf＝2km為終點(diǎn)，無控彈道的射程為22 466.2m，時(shí)間均勻分段增程彈道優(yōu)化的射程為50 721.2m，非均勻分段增程彈道優(yōu)化的射程為51 450.8m。相比無控彈道，優(yōu)化后的彈道可顯著提高制導(dǎo)炮彈的射程，其中采用非均勻分段的射程比采用均勻分段的射程增加了729.6m。由圖4所示的速度對比曲線看出，受到重力的作用無控彈在彈道頂點(diǎn)后速度會(huì)有所增大，增程優(yōu)化彈道通過控制攻角產(chǎn)生向上的升力對重力進(jìn)行一定的補(bǔ)償，速度變化緩慢，并在增程彈道的末段具有一定的速度存量。采用時(shí)間均勻分段增程彈道優(yōu)化的終端速度為151.2 m/s，非均勻分段增程彈道優(yōu)化的終端速度為157.5 m/s，后者更加滿足對終端速度的約束。圖5所示增程彈道的彈道傾角幅值較小、變化較慢，相比無控彈的飛行時(shí)間也較長。時(shí)間均勻分段優(yōu)化彈道的飛行時(shí)間為220.9s，非均勻分段優(yōu)化彈道的飛行總時(shí)間為223.1s，終端傾角均為θtf＝－20.0°，而相比非均勻分段的增程彈道，均勻分段增程彈道的傾角波動(dòng)較大，前者更利于彈道控制。圖6表明時(shí)間非均勻分段優(yōu)化彈道的攻角控制量在相應(yīng)時(shí)刻迅速調(diào)整，可以較好地反映出攻角的變化規(guī)律。

圖4 無控段彈道與增程彈道速度對比曲線

圖5 無控段彈道與增程彈道傾角對比曲線

圖6 均勻分段與非均勻分段攻角對比曲線

4 結(jié)束語

本文基于縱向質(zhì)點(diǎn)彈道模型，引入了基于時(shí)間尺度變換的時(shí)間非均勻分段思想，對攻角變量參數(shù)化，將動(dòng)態(tài)規(guī)劃問題轉(zhuǎn)化為通用的有限維NLP問題，然后利用序列二次規(guī)劃方法進(jìn)行求解，獲得了滿足各種約束的增程優(yōu)化彈道。通過仿真驗(yàn)證了該方法的有效性，相比傳統(tǒng)的時(shí)間均勻分段方法，非均勻分段優(yōu)化的彈道射程更遠(yuǎn)、終端速度更大，且更為平滑，有利于進(jìn)行彈道跟蹤控制。

[1]JEPPS G.Linearised optimal control and application to a gliding projectile[C]//AIAA Atmospheric Flight Mechanics Conference.Canberrn：AIAA，1985.

[2]FLECK V.Increase of range for an artillery projectile by using the lifting force[C]//19th International Symposium on Ballistics.San Francisco：IBC，1996.

[3]BETTS J T.Survey of numerical methods for trajectory optimization[J].Journal of Guidance，Control and Dynamics，1998，21（2）：193－207.

[4]MICHAEL R I，F(xiàn)ARIBA F.A perspective on methods for trajectory optimization [C ]//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Monterey，California：AIAA，2002.

[5]李瑜，楊志紅，崔乃剛.助推－滑翔導(dǎo)彈彈道優(yōu)化研究[J].宇航學(xué)報(bào)，2008，29（1）：66－71.LI Yu，YANG Zhi－h(huán)ong，CUI Nai－gang.A study of optimal trajectory for boost－glide missile[J].Journal of Astronautics，2008，29（1）：66－67.（in Chinese）

[6]CHEN Gang，HU Ying.Optimization design on RLV reentry trajectory based on genetic algorithm [J].Journal of Solid Rocket Technology，2006，29（4）：235－238.

[7]CHEN Gang，XU Min，WANG Zi－ming，et al.RLV reentry trajectory multi－objective optimization design based on NSGA－Ⅱalgorithm.AIAA－2005－6131[R].2005.

[8]LOXTON R C，TEOA K L，REHBOCKA V.Optimal control problems with a continuous inequality constraint on the state and the control[J].Automatica，2009，45（10）：2 250－2 257.

[9]LOXTON R C，TEOA K L.Optimal control problems with multiple characteristic time points in the objective and constraints[J].Automatica，2008，44（11）：2 923－2 929.

[10]雷陽.一種求解最優(yōu)控制問題的非均勻控制向量參數(shù)化方法[J].中國石油大學(xué)學(xué)報(bào)，2011，35（5）：180－184.LEI Yang.A non－uniform control vector parameterization approach for optimal control problems[J].Journal of China University of Petroleum，2011，35（5）：180－184.（in Chinese）

[11]BRIAN C，F(xiàn)ABIEN.Some tools for the direct solution of optimal control problems[J].Advances in Engineering Software，1998，29（1）：45－61.

[12]GOH C J，TEO L K.Control parameterization：a unified approach to optimal control problems with general constraints[J].Automatica，1988，24（1）：3－18.