李怡昕，李旭，劉少波，康鵬

(北京電子工程總體研究所，北京 100854)

0 引言

推力可調(diào)固體火箭發(fā)動(dòng)機(jī)可根據(jù)作戰(zhàn)任務(wù)、作戰(zhàn)平臺(tái)實(shí)時(shí)調(diào)節(jié)推力大小，實(shí)現(xiàn)推進(jìn)劑能量的合理分配，從而實(shí)現(xiàn)武器裝備作戰(zhàn)能力多樣化，因而成為固體火箭發(fā)動(dòng)機(jī)發(fā)展的趨勢(shì)[1]。

國外已在推力可控固體火箭發(fā)動(dòng)機(jī)的應(yīng)用方面取得一些進(jìn)展。美國Aerojet公司的EIHAWK導(dǎo)彈集成了喉栓式推力可調(diào)噴管，在裝藥量一定的條件下成功實(shí)現(xiàn)導(dǎo)彈增程。美陸軍精確攻擊導(dǎo)彈PAM采用了Aerojet公司的喉栓式變推力發(fā)動(dòng)機(jī)，可根據(jù)任務(wù)目標(biāo)實(shí)時(shí)調(diào)節(jié)推力，從而擴(kuò)大作戰(zhàn)范圍，且在2007～2009年間，成功進(jìn)行多次飛行試驗(yàn)，展示出固體發(fā)動(dòng)機(jī)推力可調(diào)技術(shù)的可行性[2]。

近年來，國內(nèi)外學(xué)者對(duì)于偽譜法用于最優(yōu)控制問題的應(yīng)用研究愈加廣泛，發(fā)展愈加成熟，而飛行器軌跡優(yōu)化中的偽譜法應(yīng)用研究是偽譜法應(yīng)用最成功的領(lǐng)域之一。

目前，偽譜法在飛行器軌跡優(yōu)化中的應(yīng)用已取得相當(dāng)可觀的成果。其中，Legendre偽譜法已被應(yīng)用于解決運(yùn)載火箭上升段軌跡快速優(yōu)化問題[3]、通用航空飛行器性能優(yōu)化問題[4]、航天飛機(jī)應(yīng)急下降軌跡優(yōu)化問題[5]、航天器最優(yōu)轉(zhuǎn)移軌道問題[6]、月球軟著陸軌跡快速優(yōu)化問題[7]；Gauss偽譜法已被應(yīng)用于解決小推力航天器軌跡優(yōu)化問題和Delta運(yùn)載火箭上升段軌跡優(yōu)化問題[8]、日-火Halo轉(zhuǎn)移軌道快速優(yōu)化設(shè)計(jì)問題[9]、多級(jí)固體運(yùn)載火箭上升段軌跡快速優(yōu)化問題[10]、空空導(dǎo)彈最優(yōu)中制導(dǎo)律設(shè)計(jì)中的應(yīng)用問題[11]；月球軟著陸軌道優(yōu)化中的應(yīng)用問題[12]、高超聲速飛行器的軌跡優(yōu)化問題[13]?？偟膩碚f，國內(nèi)外對(duì)于偽譜法在飛行器軌跡優(yōu)化設(shè)計(jì)中的應(yīng)用正處于快速發(fā)展向逐步成熟的邁進(jìn)階段[14]。

針對(duì)攔截近程目標(biāo)，本文采用變推力控制，選用固體火箭發(fā)動(dòng)機(jī)實(shí)現(xiàn)推力可調(diào)，提出考慮末速最大的目標(biāo)函數(shù)，給出相應(yīng)的約束條件，并使用高斯偽譜法進(jìn)行彈道優(yōu)化設(shè)計(jì)，最后給出相應(yīng)的仿真校驗(yàn)結(jié)果。其中，選用的變推力固體火箭發(fā)動(dòng)機(jī)推進(jìn)劑質(zhì)量給定，可根據(jù)實(shí)際需要調(diào)節(jié)推力大小。本文將高斯偽譜法應(yīng)用于處理可變推力近程導(dǎo)彈彈道優(yōu)化設(shè)計(jì)，在推進(jìn)劑一定的條件下考慮合理分配能量，提出了一種設(shè)計(jì)變推力彈道的新思路。

1 建立高斯偽譜優(yōu)化模型

1.1 導(dǎo)彈運(yùn)動(dòng)學(xué)模型

本文作如下假設(shè)：導(dǎo)彈在整個(gè)飛行期間都處于瞬時(shí)平衡狀態(tài)；導(dǎo)彈繞彈體軸的轉(zhuǎn)動(dòng)是無慣性的；導(dǎo)彈的控制系統(tǒng)理想工作，無誤差，無時(shí)間延遲。

導(dǎo)彈在鉛垂平面內(nèi)的運(yùn)動(dòng)方程組可以簡(jiǎn)化為

(1)

式中：m為導(dǎo)彈質(zhì)量；v為速度；FY為升力；FX為阻力；G為重力；P為推力；α為攻角；θ為彈道傾角;重力加速度g=9.81 m/s2。

改寫成狀態(tài)方程的形式

(2)

式中：狀態(tài)變量x=(v,θ,x,y,m)T；控制變量u=(P,α)T；t0為初始時(shí)刻；tf為末端時(shí)刻。

1.2 基于高斯偽譜法的離散表示

Gauss偽譜法實(shí)質(zhì)上是一種離散變量表示法，其核心在于將狀態(tài)變量和控制變量在高斯點(diǎn)上離散，并以離散點(diǎn)為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式逼近原連續(xù)函數(shù)。下面將利用高斯偽譜法完成系統(tǒng)微分方程式(2)到代數(shù)方程式(11)的轉(zhuǎn)化[15]：

首先，將原系統(tǒng)的狀態(tài)方程式(2)做時(shí)域變換

(3)

將原時(shí)間區(qū)間t∈[t0,tf]映射到τ∈[-1,1]，則變換后的最優(yōu)控制問題狀態(tài)方程為

(4)

性能指標(biāo)、邊界條件、不等式約束為

(5)

Φ(x(τ0),t0,x(τf),tf)=0,

(6)

C(x(τ),u(τ),τ;t0,tf)≤0,τ∈[-1,1],

(7)

式中：容許控制u(τ)∈U;τ∈[-1,1]。

其次，在區(qū)間τ∈[-1,1)內(nèi)選取Gauss點(diǎn)τk(k=1,2,…,N)作為節(jié)點(diǎn)構(gòu)造Lagrange多項(xiàng)式逼近原連續(xù)狀態(tài)變量x(τ)和控制變量u(τ)，即離散化為

(8)

(9)

記Xk=X(τk)=x(τk)，Uk=U(τk)=u(τk)。式(8)對(duì)時(shí)間求導(dǎo)則有

(10)

式中：D∈RN∧(N+1)為微分矩陣，從而系統(tǒng)微分方程式(2)轉(zhuǎn)化為代數(shù)約束：

(11)

式中：k=1,2,…,N。

原性能指標(biāo)、邊界條件和不等式約束式(5)到式(7)可離散表示為

(12)

Φ(X0,t0,Xf,tf)=0,

(13)

C(Xk,Uk,τk;t0,tf)≤0,

(14)

式中：ωk為高斯權(quán)重系數(shù);k=1,2,…,N。

至此，原最優(yōu)控制問題最終轉(zhuǎn)化為非線性規(guī)劃問題，式(11)到式(14)定義的非線性規(guī)劃問題的解就是原最優(yōu)控制問題的近似解[16]。轉(zhuǎn)換所得NLP問題采用序列二次規(guī)劃(SQP)方法求解。

2 彈道優(yōu)化方案設(shè)計(jì)

本文以近距導(dǎo)彈攔截目標(biāo)為特定研究背景，在推進(jìn)劑質(zhì)量給定的情況下采用Gauss偽譜法進(jìn)行彈道優(yōu)化設(shè)計(jì)，考慮選取J=vf末速最大為優(yōu)化目標(biāo)函數(shù)，旨在提高其末端過載能力。

2.1 發(fā)動(dòng)機(jī)方案

為了進(jìn)一步驗(yàn)證推力可調(diào)的彈道性能，本節(jié)給出變推力發(fā)動(dòng)機(jī)的主要參數(shù)，且將在此基礎(chǔ)上用高斯偽譜法進(jìn)行彈道優(yōu)化設(shè)計(jì)并與傳統(tǒng)推力形式的彈道作比較和分析。

給出某變推力發(fā)動(dòng)機(jī)參數(shù)見表1，推力調(diào)節(jié)范圍500 N～15 kN。

可控推力固體火箭發(fā)動(dòng)機(jī)的技術(shù)途徑包括調(diào)節(jié)噴管喉部面積、實(shí)時(shí)改變推進(jìn)劑燃速、加質(zhì)發(fā)動(dòng)機(jī)、膠狀推進(jìn)劑發(fā)動(dòng)機(jī)等，其中喉部面積可調(diào)的方案理論和試驗(yàn)基礎(chǔ)好，能夠?qū)崿F(xiàn)推力大小無級(jí)調(diào)節(jié)[1]。

2.2 優(yōu)化模型

狀態(tài)方程(2)中，狀態(tài)變量x=(v,θ,x,y,m)T、控制變量u=(P,α)T為未知量，即需按一定的制導(dǎo)規(guī)律生成制導(dǎo)指令。即原運(yùn)動(dòng)方程聯(lián)立控制方程則方程組封閉，給出相應(yīng)的初始條件即可求解。

本節(jié)提出一種高斯偽譜法的優(yōu)化模型，以J=vf最大為目標(biāo)函數(shù)，約束條件為攻角約束-40°≤α≤40°和法向過載約束-60≤ny≤60，選取u=(P)作為待優(yōu)化控制變量，由比例導(dǎo)引給出α值。

比例導(dǎo)引給出α值的原理如下：

第1步：由比例導(dǎo)引

(15)

給出N個(gè)高斯點(diǎn)τk(k=1,2,…,N)上的過載指令；

第2步：根據(jù)預(yù)先設(shè)定的狀態(tài)變量、控制變量猜測(cè)區(qū)間將狀態(tài)變量、控制變量在各高斯點(diǎn)τk(k=1,2,…,N)上初始化，即此時(shí)各高斯點(diǎn)上x=(v,θ,x,y,m)T，u=(P)為一系列確定值；

第3步：由動(dòng)力學(xué)方程

ny=

(16)

得出各高斯點(diǎn)上的一系列α值。

隨著高斯偽譜法優(yōu)化的逐步進(jìn)行，α值同步更新，該最優(yōu)控制問題可以求解。

3 仿真校驗(yàn)

導(dǎo)彈初始狀態(tài)x0=0，y0=10 km，v0=400 m/s，滿載質(zhì)量mfull=82.3 kg，空載質(zhì)量memp=55 kg。

仿真圖中*號(hào)各點(diǎn)為高斯偽譜法優(yōu)化得到的離散狀態(tài)量/控制量，另將高斯偽譜法得到的P，α作為輸入進(jìn)行數(shù)值積分，其結(jié)果用細(xì)實(shí)線表示，以驗(yàn)證高斯偽譜法優(yōu)化結(jié)果的正確性。

3.1 發(fā)動(dòng)機(jī)方案對(duì)比

目標(biāo)設(shè)定為(xt,yt)=(5,12) km定點(diǎn)目標(biāo)，導(dǎo)彈發(fā)射角θ=90°。

首先在表1給出的發(fā)動(dòng)機(jī)參數(shù)的基礎(chǔ)上，給定幾組工程常用推力比，通過仿真計(jì)算其在比例導(dǎo)引下的彈道，并比較末速如表2所示，仿真結(jié)果見圖1。

表2 傳統(tǒng)發(fā)動(dòng)機(jī)參數(shù)Table 2 Parameters of traditional engines

圖1 傳統(tǒng)發(fā)動(dòng)機(jī)方案的彈道仿真Fig.1 Trajectory simulation of traditional engine schemes

由表2可知，近距目標(biāo)考察末速最大時(shí)推力比1∶1性能最優(yōu)(末速最大)，將其與高斯偽譜法優(yōu)化所得的變推力方案作對(duì)比，給出仿真結(jié)果如圖2所示。

分析仿真結(jié)果可知，對(duì)比最佳的傳統(tǒng)發(fā)動(dòng)機(jī)方案比例導(dǎo)引彈道，推力可調(diào)的發(fā)動(dòng)機(jī)方案下，以末速最大為優(yōu)化目標(biāo)，高斯偽譜法優(yōu)化得到的彈道末速提升9.97%，表明采用變推力控制可顯著提升導(dǎo)彈的彈道性能，對(duì)比結(jié)果見表3。

圖2 兩種發(fā)動(dòng)機(jī)方案的彈道仿真Fig.2 Trajectory simulation of two engine schemes

表3 2種發(fā)動(dòng)機(jī)方案的性能Table 3 Trajectory performance comparison of two engine schemes

由圖2可看出，高斯偽譜法優(yōu)化結(jié)果與數(shù)值積分結(jié)果有一定誤差，這是由于高斯偽譜法給出的最優(yōu)解為離散值，由式(9)可知，最優(yōu)控制變量應(yīng)以這些離散點(diǎn)作為節(jié)點(diǎn)構(gòu)造Lagrange插值多項(xiàng)式來逼近。而細(xì)實(shí)線表示的曲線則是在離散推力值的基礎(chǔ)上線性插值作為輸入量進(jìn)行數(shù)值積分的結(jié)果，因此在推力值波動(dòng)較大時(shí)，二者出現(xiàn)誤差。

3.2 有效性驗(yàn)證

為進(jìn)一步驗(yàn)證高斯偽譜法應(yīng)用于可變推力導(dǎo)彈彈道優(yōu)化的有效性，本節(jié)考慮近程導(dǎo)彈全向發(fā)射情況，選用比例導(dǎo)引控制變量u=(P)的高斯偽譜法優(yōu)化方案，在表1給出的推力可調(diào)發(fā)動(dòng)機(jī)參數(shù)的基礎(chǔ)上進(jìn)行彈道仿真。

目標(biāo)設(shè)定為(xt,yt)=(5,12) km定點(diǎn)目標(biāo)，選取導(dǎo)彈發(fā)射角θ=0°,30°,…,150°，同3.1節(jié)給出推力可調(diào)與最優(yōu)傳統(tǒng)推力的彈道性能對(duì)比結(jié)果見表4。

表4 全向發(fā)射兩種發(fā)動(dòng)機(jī)方案的性能Table 4 Trajectory performance comparison of two omnidirectional attack engine schemes

高斯偽譜法給出導(dǎo)彈在不同的發(fā)射角θ=0°,30°,…,150°下的變推力彈道優(yōu)化結(jié)果如圖3所示，高斯偽譜法優(yōu)化得到的彈道末速分別提升1.67%，1.25%，2.92%，9.97%，30.44%，35.91%。區(qū)別于傳統(tǒng)推力方案，圖3d)推力曲線圖體現(xiàn)出可變推力控制下的推力形式多樣性，說明了高斯偽譜法應(yīng)用于可變推力導(dǎo)彈彈道優(yōu)化設(shè)計(jì)的可行性和有效性。

圖3 全向攻擊彈道仿真Fig.3 Trajectory simulation of omnidirectional attack

4 結(jié)束語

對(duì)于可變推力導(dǎo)彈彈道優(yōu)化問題，本文提出基于比例導(dǎo)引的變推力導(dǎo)彈彈道優(yōu)化模型，研究在推進(jìn)劑質(zhì)量給定的條件下進(jìn)行能量的合理分配。仿真結(jié)果表明了高斯偽譜法應(yīng)用于可變推力導(dǎo)彈彈道優(yōu)化設(shè)計(jì)的有效性。本文采用比例導(dǎo)引給出制導(dǎo)指令，總體計(jì)算量小，可以實(shí)時(shí)在線運(yùn)行。后續(xù)將對(duì)打擊中遠(yuǎn)程目標(biāo)及機(jī)動(dòng)目標(biāo)的情況進(jìn)行研究，同時(shí)考慮采用hp-自適應(yīng)方法對(duì)其進(jìn)行改進(jìn)，逐步細(xì)化網(wǎng)格并提高基函數(shù)階次以獲得更高的精度。