沙建科，施雨陽，萬自明，徐敏

(1.西北工業(yè)大學(xué) 航天學(xué)院，陜西 西安 710072； 2.北京電子工程總體研究所，北京 100854)

0 引言

沖壓發(fā)動機以質(zhì)量輕、比沖高、具有良好的超聲速特性而成為中遠程戰(zhàn)術(shù)導(dǎo)彈首選裝置，其可為導(dǎo)彈超聲速飛行提供加速爬升所需要的富裕推力及續(xù)航推力[1]。以沖壓發(fā)動機為動力的導(dǎo)彈武器系統(tǒng)研制工作是一項系統(tǒng)工程，總體設(shè)計決定導(dǎo)彈的性能和成本。優(yōu)化設(shè)計是從多種方案中選擇最佳方案的設(shè)計方案，是導(dǎo)彈總體設(shè)計的一個重要手段。導(dǎo)彈總體優(yōu)化設(shè)計的重要任務(wù)之一是對導(dǎo)彈設(shè)計參數(shù)進行尋優(yōu)，使導(dǎo)彈以最小代價滿足戰(zhàn)術(shù)技術(shù)指標的要求[2]。

隨著優(yōu)化理論的不斷完善，分系統(tǒng)優(yōu)化設(shè)計方法如導(dǎo)彈總體參數(shù)和彈道優(yōu)化設(shè)計已經(jīng)很成熟[3-7]。沖壓發(fā)動機的性能參數(shù)(如比沖、推力)與導(dǎo)彈飛行參數(shù)緊密相關(guān)，而導(dǎo)彈的飛行性能與給定的彈道形狀有關(guān)，各分系統(tǒng)單獨開展的優(yōu)化設(shè)計并不能滿足總體參數(shù)最優(yōu)的要求，要想獲得總體參數(shù)優(yōu)化最優(yōu)必須結(jié)合飛行軌跡進行一體化優(yōu)化設(shè)計?？傮w參數(shù)/飛行軌跡一體化優(yōu)化才能使系統(tǒng)參數(shù)最優(yōu)，從而最大程度地挖掘?qū)椀脑O(shè)計潛能。

本文從一體化優(yōu)化設(shè)計思想出發(fā)，對以沖壓發(fā)動機為動力導(dǎo)彈軌跡/總體參數(shù)一體化設(shè)計中的相關(guān)數(shù)學(xué)模型、優(yōu)化目標、一體化優(yōu)化方法等幾個問題進行了深入研究。提出了采用一種精度高的參數(shù)化方法-高斯偽譜法將軌跡最優(yōu)控制問題轉(zhuǎn)換為非線性規(guī)劃問題，然后用序列二次規(guī)劃方法(sequential quadratic programming, SQP)求解一體化優(yōu)化問題。

1 導(dǎo)彈總體優(yōu)化數(shù)學(xué)模型

數(shù)學(xué)模型包括質(zhì)量計算模型、彈道及氣動力計算模型和發(fā)動機性能計算模型。

1.1 導(dǎo)彈質(zhì)量計算模型

導(dǎo)彈的質(zhì)量計算是軌跡/總體參數(shù)一體化優(yōu)化設(shè)計基礎(chǔ)，導(dǎo)彈總質(zhì)量由各部分質(zhì)量組成。在導(dǎo)彈總體方案論證階段宜采用導(dǎo)出型質(zhì)量計算模型。本文的質(zhì)量計算按照部件分析方法進行估算，全彈質(zhì)量M0由助推器質(zhì)量M1和二級導(dǎo)彈質(zhì)量M2組成：

M0=M1+M2.

(1)

助推器采用固體火箭發(fā)動機，為了簡化模型，認為助推器質(zhì)量由推進劑Mp和結(jié)構(gòu)質(zhì)量Ms組成。依據(jù)有關(guān)統(tǒng)計數(shù)據(jù)，發(fā)動機推進劑質(zhì)量為

Ms=0.85M1.

(2)

二級導(dǎo)彈的質(zhì)量可分為有效載荷質(zhì)量mp、發(fā)動機推進劑質(zhì)量mpp、發(fā)動機結(jié)構(gòu)質(zhì)量mss1、能源電池質(zhì)量mcs、導(dǎo)彈結(jié)構(gòu)質(zhì)量mss和彈上設(shè)備質(zhì)量mps。即

M2=mp+mpp+mss1+mcs+mss+mps.

(3)

各部分質(zhì)量具體展開表達式可參考文獻[2].

1.2 發(fā)動機性能計算模型

總體參數(shù)優(yōu)化建模中必須根據(jù)實際要求建立發(fā)動機性能的計算模型。發(fā)動機為導(dǎo)彈提供飛行動力，以保證導(dǎo)彈獲得所需的速度和射程。沖壓發(fā)動機最大的缺點是在靜止狀態(tài)下不能自行啟動，必須借助助推器將其加速到一定的速度后才能工作。因此，本文的動力系統(tǒng)由助推器(固體火箭發(fā)動機)和沖壓發(fā)動機組成。

為了簡化計算，在總體方案設(shè)計階段，認為固體火箭發(fā)動機的推力為定值，燃料流率為一確定值。采用等熵流動，一維計算模型對沖壓發(fā)動機進行熱力計算，具體可見文獻[8]。通過仿真得到?jīng)_發(fā)動機的推力、耗油量隨飛行高度、速度、攻角、沖壓發(fā)動機余氣系數(shù)、喉道收縮比的變化關(guān)系，即

F=f1(h,v,α,β,η)，

(4)

mc=f2(h,v,α,β,η)，

(5)

式中：F為推力；mc為耗油量；h為飛行高度；v為飛行速度；α為飛行攻角；β為余氣系數(shù)；η為喉道收縮比。

1.3 彈道計算模型

導(dǎo)彈的彈道由起飛助推段、爬升段和巡航段組成。助推器在地面點火加速到?jīng)_壓發(fā)動機啟動馬赫數(shù)以后脫落，沖壓發(fā)動機開始接力工作使導(dǎo)彈加速爬升到巡航高度、巡航速度后進行水平巡航。為了研究方便，采用了以下假設(shè)：①將導(dǎo)彈看作可控質(zhì)點；②為簡化問題，僅研究導(dǎo)彈在垂直平面內(nèi)的運動；③略去飛行過程中的隨機干擾，即控制系統(tǒng)是理想的，導(dǎo)彈運動完全按照程序飛行。

基于以上假設(shè)，導(dǎo)彈在鉛垂平面內(nèi)的質(zhì)心運動方程組為

(6)

大氣模型采用1976美國標準大氣。導(dǎo)彈的氣動力計算采用部件組合的工程估算方法，基本能滿足導(dǎo)彈總體方案設(shè)計的精度要求。阻力D和升力L分別為

(7)

(8)

式中：S為參考面積；ρ為大氣密度；CD和CL分別為阻力系數(shù)和升力系數(shù)，可表示為攻角和馬赫數(shù)的插值函數(shù)，也可近似表示為攻角的函數(shù)。

2 導(dǎo)彈軌跡/總體一體化優(yōu)化問題描述

2.1 目標函數(shù)

優(yōu)化設(shè)計的任務(wù)就是要對各個設(shè)計方案進行尋優(yōu)比較，得出最佳方案，而對設(shè)計方案進行優(yōu)劣比較的準則就是目標函數(shù)。目標函數(shù)應(yīng)能充分反映主要設(shè)計性能指標。本文以起飛質(zhì)量最小為目標函數(shù)，該目標函數(shù)既能反映了導(dǎo)彈性能的好壞，也在一定程度上反映了成本。

2.2 優(yōu)化變量

優(yōu)化設(shè)計變量是指能夠用來描述設(shè)計方案特征的獨立變量。優(yōu)化變量應(yīng)選取與目標函數(shù)有關(guān)系，對目標有較大影響的變量。本文的優(yōu)化變量由總體參數(shù)優(yōu)化和軌跡優(yōu)化兩者的優(yōu)化變量組成。根據(jù)問題性質(zhì)和設(shè)計經(jīng)驗，選取下述的優(yōu)化變量。

(1) 助推器工作時間t1

在裝藥量一定時，助推器工作時間越短，助推段末速就越大，但相應(yīng)的導(dǎo)彈結(jié)構(gòu)質(zhì)量也增加，且受到發(fā)動機燃燒室壓強的限制，發(fā)動機工作時間也不可能無限制的縮短，因此，需通過優(yōu)化助推工作時間以提高導(dǎo)彈總體性能。

(2) 初始彈道傾角θ0

增大初始彈道傾角可使導(dǎo)彈快速爬升到預(yù)定的巡航高度，減少阻力的消耗，增加巡航射程，但是也增加了爬升段重力的消耗，同時減少了水平飛行距離，為了提高總體性能應(yīng)優(yōu)化初始彈道傾角。

(3) 接力馬赫數(shù)Mat1f

接力馬赫數(shù)為助推段的末速，二級導(dǎo)彈的初速，是導(dǎo)彈總體設(shè)計指標。接力馬赫數(shù)值增大時，則助推段質(zhì)量增加，反之，二級導(dǎo)彈的質(zhì)量增加，導(dǎo)彈的起飛質(zhì)量M0會因Mat1f的不同而發(fā)生變化。優(yōu)化接力馬赫數(shù)可以使得導(dǎo)彈總體性能得到進一步的優(yōu)化。

(4) 沖壓發(fā)動機余氣系數(shù)β及喉道收縮比η

發(fā)動機的性能參數(shù)取決于余氣系數(shù)及喉道收縮比，因此β和η影響著導(dǎo)彈的飛行性能，為達到給定的射程，總可能找到合適的β和η使得導(dǎo)彈的起飛質(zhì)量最小。

(5) 軌跡最優(yōu)攻角變化規(guī)律α(t)

飛行攻角是軌跡的控制變量，優(yōu)化飛行攻角使得軌跡最優(yōu)減少燃油消耗，從而達到起飛質(zhì)量最小的目的。

2.3 約束條件

導(dǎo)彈在飛行過程中經(jīng)歷助推段、爬升段和巡航段，跨越的時間和空間范圍比較大。飛行過程中受到過程不等式約束、終端約束。此外，相鄰的2段軌跡(助推段與爬升段、爬升段與巡航段)之間還需要滿足結(jié)點連接條件。本文的約束條件如下：

不等式約束：

αmin≤αmax，Mat1f>Ma0，

(9)

等式約束：

hti3=hcruise，θti3=θcruise，vti3=vcruise，L=Lmax.

(10)

相鄰兩段連接約束：

(11)

式中：Mat1f為沖壓發(fā)動機機轉(zhuǎn)級馬赫數(shù)；tf1為助推段關(guān)機時刻，為固定值；ti2和tf2分別為爬升段初始時刻和終端時刻；ti3為巡航初始時間。

爬升段轉(zhuǎn)到巡航段的標志是爬升段彈道傾角為0。

2.4 一體化優(yōu)化方法

導(dǎo)彈軌跡/總體參數(shù)一體化優(yōu)化變量分為2類，一類是導(dǎo)彈總體參數(shù)靜態(tài)優(yōu)化；另一類是動態(tài)軌跡優(yōu)化。由以上分析可知，總體參數(shù)優(yōu)化問題屬于非線性規(guī)劃問題，而軌跡優(yōu)化設(shè)計本質(zhì)上屬于最優(yōu)控制問題。采用基于pontryagin極小值原理的間接求解軌跡優(yōu)化問題比較繁瑣，且只能與總體參數(shù)優(yōu)化分開單獨優(yōu)化，有可能失去整體的最優(yōu)解，本文利用基于非線性規(guī)劃的直接法來求解。軌跡最優(yōu)控制問題的搜索空間是泛函空間，不能直接利用非線性規(guī)劃算法進行求解。為了統(tǒng)一用非線性規(guī)劃算法求解軌跡/總體一體化優(yōu)化問題，應(yīng)該用參數(shù)化方法將軌跡最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題。參數(shù)化方法是將最優(yōu)軌跡問題轉(zhuǎn)化為非線性規(guī)劃問題求解的橋梁，其精度直接關(guān)系到最優(yōu)軌跡的好壞程度。本文提出采用高斯偽譜法將軌跡最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題。

高斯偽譜法是近年來發(fā)展迅速、應(yīng)用較多的一種同時離散控制變量和狀態(tài)變量的函數(shù)逼近參數(shù)化方法[9-11]。高斯偽譜法屬于配點法的范疇，其基本原理是：以Legendre多項式的零點(也稱為高斯點)作為節(jié)點同時將狀態(tài)方程中的控制變量和狀態(tài)變量進行離散，然后以這些節(jié)點作為插值點構(gòu)造全局的拉格朗日多項式來逼近狀態(tài)變量及控制變量，通過對全局插值多項求導(dǎo)來近似動力學(xué)方程中的狀態(tài)變量對時間的導(dǎo)數(shù)，從而將動力學(xué)方程轉(zhuǎn)化為一組代數(shù)方程。目標函數(shù)中的積分項由高斯積分計算得到。經(jīng)過上述變換，可將軌跡最優(yōu)控制問題轉(zhuǎn)化為受一系列代數(shù)約束的參數(shù)優(yōu)化問題，即非線性規(guī)劃問題，基本流程圖如圖1所示。

高斯偽譜法的優(yōu)點主要表現(xiàn)在：①可估算出原最優(yōu)控制問題的協(xié)態(tài)變量信息，因此能保證所獲得的非線性規(guī)劃的最優(yōu)解是原最優(yōu)控制問題的解。②因同時離散控制變量及狀態(tài)變量，其優(yōu)化變量的維數(shù)遠高于其他方法，但以此為代價，降低了目標函數(shù)的病態(tài)方程，提高了收斂性和精度，且對初值猜測值不敏感。

高斯偽譜方法將軌跡最優(yōu)控制問題參數(shù)化，轉(zhuǎn)化為非線性規(guī)劃問題后， 理論上各種非線性規(guī)劃算法均適用該問題求解。本文提出采用SQP算法對該問題進行求解。SQP優(yōu)化算法[12]是在每一迭代步中，通過求解一個二次規(guī)劃子問題，來確定一個下降方向，以減少價值函數(shù)來取得步長，重復(fù)這些步驟，直到求得原問題的解，該方法在局部整體收斂性的同時，保持局部超1次收斂性，是求解約束優(yōu)化問題最有效的算法之一，也是應(yīng)用最為成功和廣泛的一種算法。

3 仿真結(jié)果與討論

根據(jù)本文所建立的軌跡/總體一體化優(yōu)化方法，以某遠程沖壓發(fā)動機導(dǎo)彈為例，分別對不考慮軌跡優(yōu)化的總體參數(shù)優(yōu)化(總體優(yōu)化方案)和考慮軌跡/總體參數(shù)一體化優(yōu)化(一體化優(yōu)化方案)進行了仿真分析。本算例的初值為：(h,x,v)=(0,0,0)；約束為Mat1f>1.8，α(t)∈[-3°,8°]，t1∈[3,7]。計算中巡航高度為16 km，巡航馬赫數(shù)為3，射程為300 km。優(yōu)化結(jié)果如表1所示。

圖1 高斯偽譜法基本流程圖Fig.1 Basic block diagram of Gauss pseudo-spectral method

表1 總體參數(shù)優(yōu)化結(jié)果Table 1 Overall parameter optimization results

參數(shù)總體優(yōu)化一體化優(yōu)化初始彈道傾角/(°)68.4262.457余氣系數(shù)β1.8791.45～2.0喉道收縮比η0.610.57助推時間/s4.0483.868爬升時間/s145.35354.144巡航時間/s235.655300.06總飛行時間/s385.056358.072轉(zhuǎn)級馬赫數(shù)1.8431.837轉(zhuǎn)級高度/km1.0730.923平均速度/(m·s-1)780.325837.82射程/km300.469300.02助推器裝藥量/kg283.36270.76助推器總質(zhì)量/kg333.364318.54二級導(dǎo)彈裝藥量/kg164.936135.18二級導(dǎo)彈總質(zhì)量/kg764.936735.18起飛質(zhì)量/kg1 098.301 053.72

從上述優(yōu)化結(jié)果可以看出，在沖壓發(fā)動機接力工作時，2種方案的導(dǎo)彈飛行馬赫數(shù)分別為1.834和1.837，均滿足接力馬赫數(shù)的條件。接力馬赫數(shù)的大小，主要取決于助推器和二級發(fā)動機的比沖，沖壓發(fā)動機的比沖遠遠大于固體火箭發(fā)動機的比沖，因此為了使起飛質(zhì)量最小，接力馬赫數(shù)取值在沖壓發(fā)動機能正常工作附近。導(dǎo)彈起飛質(zhì)量減少了4.06%，主要是助推段和爬升段燃油減少的結(jié)果，在總體參數(shù)優(yōu)化中，增加軌跡的優(yōu)化可以大大提高導(dǎo)彈的設(shè)計潛力。彈道曲線和速度隨時間變化曲線如圖2和圖3所示，一體化優(yōu)化方案中，為了減少助推段重力的損失，助推器工作時間減小到3.818 s；導(dǎo)彈以較大的推力迅速爬升到巡航高度，在爬升轉(zhuǎn)平飛時，馬赫數(shù)為2.95，迅速爬升到預(yù)定高度，可減少阻力消耗，增加巡航段射程，充分發(fā)揮沖壓發(fā)動優(yōu)勢。隨著飛行高度的增加，空氣密度減少，進入進氣道的空氣流量減小，使得燃料秒流量減少，從而達到了減輕起飛質(zhì)量的目的。飛行攻角隨時間的變化如圖4所示，滿足約束要求，變化比較平滑，容易控制。

圖2 導(dǎo)彈彈道曲線Fig.2 Change curve of trajectory vs. time

圖3 導(dǎo)彈飛行速度時間曲線Fig.3 Change curve of velocity vs. time

圖4 導(dǎo)彈飛行攻角時間曲線Fig.4 Attack angle curve of vs. time

4 結(jié)束語

本文針對沖壓發(fā)動機導(dǎo)彈的特點，建立了導(dǎo)彈軌跡/總體參數(shù)一體化優(yōu)化設(shè)計模型及方法。為了用非線性規(guī)劃方法求解一體化優(yōu)化問題，提出了采用高斯偽譜法將一體化中的軌跡最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題。分別對軌跡/總體參數(shù)一體化優(yōu)化和未考慮軌跡優(yōu)化2種方案進行了仿真分析。仿真結(jié)果表明，起飛質(zhì)量減少了4.06%，采用一體化優(yōu)化設(shè)計思想，可以最大程度地挖掘?qū)椏傮w設(shè)計的潛力，其結(jié)果可為沖壓發(fā)動機導(dǎo)彈總體優(yōu)化設(shè)計研究提供理論參考。

參考文獻：

[1] 鮑福廷，黃熙君, 張振鵬. 固體火箭沖壓組合發(fā)動機[M]. 北京：中國宇航出版社，2006.

BAO Fu-ting, HUANG Xi-jun, ZHANG Zhen-peng. Integral Solid Propellant Ramjet Rocket Motor [M].Beijing: China Astronautic Publishing House, 2004.

[2] 谷良賢，溫炳恒.導(dǎo)彈總體設(shè)計原理[M],西安：西北工業(yè)大學(xué)出版社，2004.

GU Liang-xian, WEN Bing-heng. Missile Overall Design Principle [M].Xi’an: NWPU Publishing House, 2004.

[3] 張弫，鄭時鏡，于本水. 遺傳算法在遠程防空導(dǎo)彈總體優(yōu)化設(shè)計中的應(yīng)用[J].系統(tǒng)工程與電子技術(shù),2003,25(1):34-37.

ZHANG Zhen, ZHENG Shi-jing, YU Ben-shui. Application of the Genetic Algorithms in the Optimization Design of the Long-Range Air-Defense Missile System [J].Systems Engineering and Electronics,2003, 25(1):34-37.

[4] 李向林，于本水.GA與ST在防空導(dǎo)彈總體參數(shù)優(yōu)化設(shè)計中的應(yīng)用[J].系統(tǒng)工程與電子技術(shù)，2000,22(10):14-16.

LI Xiang-lin, YU Ben-shui. Application of Genetic Algorithm and Stepwise Tolerance Method in the Optimal Design of Systemic Parameter of Antiaircraft Missile [J]. Systems Engineering and Electronics, 2000, 22(10):14-16.

[5] 張弫，鄭時鏡，于本水.遠程防空導(dǎo)彈彈道設(shè)計技術(shù)研究[J]. 系統(tǒng)工程與電子技術(shù)，2003,25(3):304-307.

ZHANG Zhen, ZHENG Shi-jing,YU Ben-shui.Study on the Trajectory Design Technology for the Long-Range Air-Defense Missile [J].Systems Engineering and Electronics，2003,25(3):304-307.

[6] 王志剛,嚴輝，陳士櫓.軌跡/飛行器總體參數(shù)的一體化優(yōu)化方法研究[J].飛行力學(xué)，1997,15(2):19-26.

WANG Zhi-gang, YAN Hui, CHEN Shi-lu. Study of Integrated Optimization Method of Trajectory/Vehicle Overall Parameter [J]. Flight Dynamics, 1997, 15(2):19-26.

[7] 胡凡，楊希祥，江振宇，等，固體運載火箭軌跡/總體參數(shù)一體化優(yōu)化設(shè)計研究[J].固體火箭技術(shù)，2010,33(6):599-602.

HU Fan,YANG Xi-xiang, JIANG Zhen-yu，et al. Integrated Optimization Design of Trajectory/System Parameters for Solid Launch Vehicles[J].Journal of Solid Rocket Technology, 2010,33(6):599-602.

[8] 劉興洲.飛航導(dǎo)彈動力裝置(上)[M].北京：中國宇航出版社，1992.

LIU Xing-zhou. Cruise Missile Power Device (Volume 1) [M]. Beijing: China Astronautic Publishing House, 1992.

[9] BENSON D A, HUNTINGTON G T, THORVALDSEN T P,et al. Direct Trajectory Optimization and Costate Estimation via an Orthogonal Collocation Method[J]. Journal of Guidance, Control and Dynamics, 2006, 29(6):1435-1440.

[10] BENSON D A. A Gauss Pseudospectral Transcription for Optimal Control [D]. Boston: Massachusetts Institute of Technology, 2004.

[11] 宗群，田柏苓，竇立謙. 基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化[J].宇航學(xué)報,2010,31(7):1775-1781.

ZONG Qun,TIAN Bai-ling,DOU Li-qian.Ascent Phase Trajectory Optimization for Near Space Vehicle Based on Gauss Pseudospectral Method[J].Journal of Astronautics, 2010,31(7):1775-1781.

[12] BOGGS P T, TOLLE J W. Sequential Quadratic Programming [J].Journal of Acta Numerica, 1995(4):1-51.