面向航程能力的固體火箭發(fā)動(dòng)機(jī)方案設(shè)計(jì)優(yōu)化

宋少倩，陳永信，任鵬飛，周文勇，李偉喆

(1.中國(guó)運(yùn)載火箭技術(shù)研究院,北京 100076；2.上海機(jī)電工程研究所,上海 201109；3.北京宇航系統(tǒng)工程研究所,北京 100076)

為適應(yīng)未來高威脅作戰(zhàn)環(huán)境的迫切需要，對(duì)飛行器的航程提出了越來越高的要求[1]。由于發(fā)射平臺(tái)的約束，飛行器的質(zhì)量與尺寸往往有著嚴(yán)苛的限制，要實(shí)現(xiàn)高性能的戰(zhàn)術(shù)技術(shù)指標(biāo)，除了先進(jìn)總體布局及新理論、新材料、新技術(shù)的應(yīng)用外，必須有高性能的推進(jìn)系統(tǒng)與之相匹配[2]。

固體火箭發(fā)動(dòng)機(jī)(Solid rocket motor,SRM)具有結(jié)構(gòu)簡(jiǎn)單、快速反應(yīng)能力強(qiáng)、工作可靠性高等優(yōu)點(diǎn)，廣泛應(yīng)用于運(yùn)載火箭、導(dǎo)彈等飛行器系統(tǒng)中[3]。傳統(tǒng)上在飛行器設(shè)計(jì)初期，發(fā)動(dòng)機(jī)設(shè)計(jì)單位按給定質(zhì)量與尺寸要求給出初步的性能參數(shù)，供總體設(shè)計(jì)評(píng)估，在不滿足航程要求時(shí)再迭代設(shè)計(jì)，由于傳統(tǒng)的設(shè)計(jì)方法沒有直接引入航程能力等總體性能的影響[4]，故難以得到系統(tǒng)最優(yōu)需求的發(fā)動(dòng)機(jī)設(shè)計(jì)方案。因此在飛行器設(shè)計(jì)初期，除了不斷提升SRM自身性能外，基于總體性能驅(qū)動(dòng)的一體化設(shè)計(jì)優(yōu)化是提高飛行器航程的有效途徑之一。后者需要結(jié)合飛行器總體布局合理設(shè)計(jì)發(fā)動(dòng)機(jī)的內(nèi)彈道特性，并通過飛行器軌跡合理分配與使用發(fā)動(dòng)機(jī)的能量。

早期的一體化設(shè)計(jì)優(yōu)化方法受限于計(jì)算資源和時(shí)間成本，往往采用大量近似解析公式、經(jīng)驗(yàn)數(shù)據(jù)，適用范圍受限，對(duì)于新型高性能飛行器，其設(shè)計(jì)結(jié)果往往與后續(xù)詳細(xì)設(shè)計(jì)的性能差別較大，容易導(dǎo)致設(shè)計(jì)的反復(fù)[4-5]。為保證設(shè)計(jì)優(yōu)化結(jié)果的可信度，需要在分析模型中采用大量的高精度數(shù)值計(jì)算，這導(dǎo)致進(jìn)行一次方案評(píng)估的計(jì)算時(shí)間成本很高。若直接采用傳統(tǒng)參數(shù)優(yōu)化方法調(diào)用精確分析模型進(jìn)行設(shè)計(jì)優(yōu)化，則優(yōu)化總時(shí)間難以承受。因此，有必要采用近似優(yōu)化(Approximate optimization,AO)方法[6-11]進(jìn)行求解，即采用計(jì)算代價(jià)較小的代理模型來代替復(fù)雜且耗時(shí)的數(shù)值分析，從而以較小的時(shí)間成本獲得較優(yōu)的設(shè)計(jì)方案。

AO方法的策略主要分為：直接法(One-shot approach,OA)[6-8]和序列迭代法(Updating approach, UA)[7-11]。OA策略僅構(gòu)建一次代理模型，并通過一定準(zhǔn)則校驗(yàn)代理模型的全局精度，而后在該代理模型上進(jìn)行優(yōu)化。UA策略則采用較少的初始樣本點(diǎn)構(gòu)建低精度代理模型，并根據(jù)一定加點(diǎn)準(zhǔn)則，在序列優(yōu)化過程中不斷更新代理模型的信息，直至滿足優(yōu)化停止條件。相比于OA策略，UA策略無(wú)需校驗(yàn)代理模型的全局精度，特別是對(duì)于高維優(yōu)化問題，UA策略調(diào)用耗時(shí)模型的次數(shù)更少。近年來，基于UA策略的序列近似優(yōu)化(Sequential approximate optimization,SAO)方法得到了廣泛發(fā)展[7-11]。例如，基于Kriging代理模型的高效全局優(yōu)化(Efficient global optimization,EGO)算法[8]，其定義了目標(biāo)函數(shù)值改善期望(Expected improvement,EI)，并選擇EI最大值點(diǎn)作為新的樣本點(diǎn)，該算法具有較好的全局搜索性，但在最優(yōu)解附近收斂速度較慢[9]；基于目標(biāo)函數(shù)預(yù)測(cè)值最小(Minimizing the predicted objective function,MP)準(zhǔn)則的優(yōu)化算法，其屬于貪婪算法，收斂速度較快，但容易陷入局部最優(yōu)解；結(jié)合多種加點(diǎn)準(zhǔn)則的并行加點(diǎn)算法[7,9,11]，其能夠彌補(bǔ)單一加點(diǎn)準(zhǔn)則的不足，目前仍是研究的熱點(diǎn)。

本文以固體動(dòng)力助推滑翔式飛行器為背景，針對(duì)其規(guī)模約束嚴(yán)格、飛行環(huán)境嚴(yán)酷及高性能需求的特點(diǎn)，提出面向航程能力的SRM方案設(shè)計(jì)優(yōu)化方法，建立復(fù)雜三維裝藥發(fā)動(dòng)機(jī)性能分析的參數(shù)化模型，以及基于自適應(yīng)偽譜法的多約束飛行器航程能力評(píng)估模型，提出基于Kriging代理模型的改進(jìn)混合整數(shù)序列優(yōu)化算法，開展發(fā)動(dòng)機(jī)與飛行器軌跡一體化設(shè)計(jì)優(yōu)化。

1 模型建立

1.1 參數(shù)優(yōu)化問題

不失一般性，參數(shù)優(yōu)化問題可表示為

minf(z)

(1)

式中：z為設(shè)計(jì)變量，f(z)為目標(biāo)函數(shù)，gi(z)為第i個(gè)不等式約束，ng為不等式約束數(shù)量，hj(z)為第j個(gè)等式約束，nh為等式約束數(shù)量。

對(duì)于SRM總體方案設(shè)計(jì)優(yōu)化問題，本文以單級(jí)助推滑翔式飛行器為例，在不改變飛行器氣動(dòng)外形的條件下，通過設(shè)計(jì)優(yōu)化發(fā)動(dòng)機(jī)的裝藥參數(shù)和噴管參數(shù)，即設(shè)計(jì)變量z，來提升飛行器的航程能力，即目標(biāo)函數(shù)f(z)，其中要求發(fā)動(dòng)機(jī)總質(zhì)量不大于要求值，即不等式約束gi(z)，發(fā)動(dòng)機(jī)的外徑、總長(zhǎng)度和裝藥質(zhì)量為固定值，即等式約束hj(z)。

(2)

式中：s為飛行器的航程；dM、dMd分別為發(fā)動(dòng)機(jī)外徑和要求值；lM、lMd分別為發(fā)動(dòng)機(jī)總長(zhǎng)度和要求值；mp、mpd分別為裝藥質(zhì)量和要求值；mM為發(fā)動(dòng)機(jī)總質(zhì)量；下標(biāo)max表示最大值或上限值；|z表示設(shè)計(jì)方案z所對(duì)應(yīng)的結(jié)果。

值得注意的是，設(shè)計(jì)變量z中除了實(shí)數(shù)型變量外，往往還包括整數(shù)型變量，例如翼柱型裝藥的翼槽數(shù)量等。因此，式(1)所示的參數(shù)優(yōu)化問題為混合整數(shù)優(yōu)化問題。

1.2 發(fā)動(dòng)機(jī)性能分析模型

1.2.1 裝藥模型

考慮到翼柱形裝藥結(jié)構(gòu)完整性好，體積裝填分?jǐn)?shù)高，燃面可調(diào)范圍大，目前廣泛應(yīng)用于各類SRM中[12-13]。本文發(fā)動(dòng)機(jī)裝藥采用前、后三角翼柱形構(gòu)型，如圖1所示。

1.2.2 殼體模型

發(fā)動(dòng)機(jī)前、后封頭均采用橢球比相同的橢球封頭，殼體圓筒段的厚度δc1采用最大應(yīng)變能強(qiáng)度理論確定，橢球封頭厚度沿橢球表面從邊緣向中心逐漸加厚，頂點(diǎn)處的厚度δc2采用最大應(yīng)力強(qiáng)度理論確定，其計(jì)算式[13]為

(3)

式中：pc為燃燒室壓強(qiáng)；kc、φc、ξc分別為安全系數(shù)、壓力波動(dòng)系數(shù)和焊縫強(qiáng)度系數(shù)，本文分別取為1.25、1.10和0.90；[σc]為殼體材料許用應(yīng)力；d1為圓筒的內(nèi)徑；km為前、后封頭的橢球比，取為2。

圖1 翼柱形裝藥幾何示意

1.2.3 噴管模型

噴管采用特型噴管，入口段采用橢圓形型面，擴(kuò)張段采用優(yōu)化3次多項(xiàng)式型面，入口段與擴(kuò)張段之間采用圓弧過渡。其中擴(kuò)張段型面的表達(dá)式為[13]

(4)

式中：xn為擴(kuò)張段型面上的點(diǎn)到噴管入口段前端面的距離；yn為擴(kuò)張段型面上的點(diǎn)到旋成軸線的距離；lnb、lnt、lne分別為噴管收斂段、過渡段和擴(kuò)張段的長(zhǎng)度；rnt、rne分別為喉部半徑和擴(kuò)張段出口半徑；rnm為過渡段圓弧半徑，本文取rnm=rnt；αnm、αne分別為擴(kuò)張段入口和出口的擴(kuò)張半角，本文分別取值為26°和22°。

圖2 噴管構(gòu)型

1.2.4 質(zhì)量模型

發(fā)動(dòng)機(jī)總質(zhì)量mM可以表示為

(5)

式中：ρp、Vp分別為裝藥的密度和體積；ρc、Vc分別為燃燒室殼體的材料密度和體積；ρq、Vq分別為燃燒室絕熱層的材料密度和體積；ρnb、Vnb分別為喉襯的材料密度和體積；ρne、Vne分別為擴(kuò)張段絕熱層的材料密度和體積；ma為其他裝置質(zhì)量，包括：點(diǎn)火裝置、發(fā)火裝置、連接裝置等，按經(jīng)驗(yàn)估算。其中Vnb和Vne采用文獻(xiàn)[14]方法計(jì)算，其余部件體積由幾何模型計(jì)算獲得。不同時(shí)刻發(fā)動(dòng)機(jī)相對(duì)飛行器頭部頂點(diǎn)處的質(zhì)心位置xGM亦由幾何模型計(jì)算獲得。

1.2.5 推力模型

對(duì)于長(zhǎng)徑比較小的SRM，采用零維內(nèi)彈道模型進(jìn)行性能預(yù)示可滿足工程設(shè)計(jì)的精度要求[12]。零維內(nèi)彈道模型假設(shè)燃燒室內(nèi)的燃速、壓強(qiáng)、溫度等參數(shù)隨空間分布均勻，根據(jù)質(zhì)量守恒關(guān)系，則可建立燃燒室內(nèi)的流動(dòng)控制方程為

(6)

式中：Vce為燃燒室空腔自由容積；Rg、Tg分別為燃?xì)獾臍怏w常數(shù)和燃燒溫度；ep為燃燒厚度；Ap為裝藥在燃燒厚度為ep時(shí)的燃面面積；a0、n0分別為裝藥在標(biāo)準(zhǔn)設(shè)計(jì)狀態(tài)下的燃速系數(shù)和燃速壓強(qiáng)指數(shù)；C*為裝藥的特征速度；ant為喉部燒蝕速率，取為0.24 mm/s。

則發(fā)動(dòng)機(jī)推力大小FT為

(7)

式中：FT0為動(dòng)推力，γg為燃?xì)獗葻岜?，pe為噴管出口壓力，pa為環(huán)境壓力，εe為噴管膨脹比，其中εe與pe存在關(guān)系為

(8)

可采用牛頓迭代法求解式(8)以獲得pe。

對(duì)于復(fù)雜三維裝藥的燃面計(jì)算問題，最小距離函數(shù)(Minimum distance function,MDF)法能夠自動(dòng)處理燃燒過程中裝藥幾何特征消失/增加等拓?fù)渥兓?，具有較好的穩(wěn)定性、通用性和計(jì)算精度[13,15]。該方法首先建立裝藥體網(wǎng)格到初始燃面面網(wǎng)格的最小距離函數(shù)，根據(jù)平行層燃燒定律，距離值等于燃燒厚度ep的等值面即為當(dāng)前燃面Ap，從而為式(6)提供輸入數(shù)據(jù)。詳細(xì)的計(jì)算步驟見文獻(xiàn)[15]，本文不再贅述。

1.3 航程能力評(píng)估模型

1.3.1 飛行器運(yùn)動(dòng)模型

為有效進(jìn)行航程能力分析，在方案設(shè)計(jì)優(yōu)化階段可忽略地球扁率和自轉(zhuǎn)影響，建立縱向平面有動(dòng)力無(wú)側(cè)滑的質(zhì)點(diǎn)運(yùn)動(dòng)方程。同時(shí)，在方程中引入迎角變化率，以便考慮控制系統(tǒng)限制的影響。則運(yùn)動(dòng)方程可表示為

(9)

式中：v、θ、h、s、α分別為飛行器的速度、航跡角、高度、航程和迎角；ud為迎角變化率控制指令；mT為飛行器有效載荷質(zhì)量，其值為500 kg；μ、RE分別為地球引力常數(shù)和半徑，其值分別為3.986 005×1014m3/s2和6 371 004 m；FD、FL分別為飛行器的阻力和升力，其可表示為

(10)

式中：CA、CN分別為軸向力系數(shù)和法向力系數(shù)；CD、CL分別為阻力系數(shù)和升力系數(shù)；SR為飛行器的參考面積；Ma為馬赫數(shù)；q為動(dòng)壓；ρ、a分別為當(dāng)?shù)卮髿饷芏群彤?dāng)?shù)芈曀伲沪臑楦┭龆嫫牵鶕?jù)“瞬時(shí)平衡假設(shè)”有

(11)

式中：CMZ、xGT分別為相對(duì)飛行器頭部頂點(diǎn)處的俯仰力矩系數(shù)和有效載荷質(zhì)心位置，xGT= 1.8 m；LR為飛行器的參考長(zhǎng)度，值為6 m。由于式(10)難以解析求解，故本文采用牛頓迭代法來獲得當(dāng)前Ma和α所對(duì)應(yīng)的δ。

1.3.2 飛行約束條件

助推滑翔式飛行器在大氣層內(nèi)高速飛行，其力熱環(huán)境嚴(yán)酷，因此，在航程評(píng)估時(shí)需要考慮防熱、載荷、控制的限制約束，即

(12)

式中：RN為駐點(diǎn)曲率半徑，值為0.03；γ為完全氣體比熱比，其值為1.4；γh為高溫氣體比熱比，其值為1.2；VQ為駐點(diǎn)熱流密度；ny為法向過載；g0為海平面引力加速度，其值為9.806 6 m/s2；下標(biāo)min為最小值或下限值。

此外，飛行末端還需要進(jìn)行能量管理，即

(13)

式中：tf、vf、θf(wàn)分別為末段的時(shí)刻、速度和航跡角。

1.3.3 基于最優(yōu)控制方法的航程能力評(píng)估

當(dāng)給定一組設(shè)計(jì)方案z時(shí)，根據(jù)發(fā)動(dòng)機(jī)性能分析模型方法可獲得所對(duì)應(yīng)的發(fā)動(dòng)機(jī)推力FT、發(fā)動(dòng)機(jī)質(zhì)量mM、發(fā)動(dòng)機(jī)質(zhì)心位置xGM隨時(shí)間t的變化關(guān)系。至此，對(duì)于式(2)中目標(biāo)函數(shù)的計(jì)算，則可轉(zhuǎn)化為對(duì)連續(xù)時(shí)間最優(yōu)控制問題的求解，即

f(z)=smaxz=-s(tf)

(14)

式中：狀態(tài)變量x=[vθhsα]T；控制變量u=ud；C為路徑約束，見式(12)；φ為端點(diǎn)約束，見式(13)。

上述航程能力評(píng)估的實(shí)質(zhì)為求解最優(yōu)控制問題，即為內(nèi)層優(yōu)化。為了保證式(1)所示的參數(shù)優(yōu)化問題能夠穩(wěn)健有效進(jìn)行，在保證式(14)的求解精度前提下，需要對(duì)求解方法的快速性和穩(wěn)定性提出要求。自適應(yīng)偽譜法將連續(xù)時(shí)間最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題，根據(jù)一定準(zhǔn)則來保證問題轉(zhuǎn)化的精度，這類方法對(duì)于初值不敏感，近年來廣泛應(yīng)用于飛行器軌跡優(yōu)化中[16-19]。

本文采用改進(jìn)自適應(yīng)Legendre-Gauss-Radau (LGR)偽譜法[19]作為內(nèi)層優(yōu)化求解器。文獻(xiàn)[19]驗(yàn)證了該方法對(duì)于多約束飛行器軌跡優(yōu)化問題的求解快速性和可靠性，該方法根據(jù)當(dāng)前解的相對(duì)誤差來調(diào)整配點(diǎn)區(qū)間數(shù)量，在精度不滿足的區(qū)間內(nèi)根據(jù)Legendre近似理論增加配點(diǎn)或進(jìn)行分割，在精度滿足的區(qū)間內(nèi)根據(jù)多項(xiàng)式系數(shù)分析減少配點(diǎn)或進(jìn)行合并[18-19]，直至計(jì)算精度滿足要求。

1.4 發(fā)動(dòng)機(jī)方案評(píng)估模型

綜上所述，進(jìn)行單次發(fā)動(dòng)機(jī)總體方案評(píng)估的流程如圖3所示。

圖3 單次SRM設(shè)計(jì)評(píng)估流程

2 改進(jìn)混合整數(shù)序列近似優(yōu)化算法

2.1 Kriging代理模型

Kriging代理模型能夠提供任意位置處函數(shù)的預(yù)測(cè)值以及該預(yù)測(cè)值的統(tǒng)計(jì)學(xué)誤差估計(jì)。該誤差估計(jì)信息可用于指導(dǎo)新樣本點(diǎn)的選擇。

(15)

(16)

本文中回歸基函數(shù)采用一次多項(xiàng)式模型[9,20]，相關(guān)基函數(shù)采用3次樣條模型[9,20]，即

(17)

式中：θ(k)為模型控制參數(shù)，可根據(jù)極大似然原則[20]，求解參數(shù)優(yōu)化問題式(18)獲得，即

(18)

2.2 混合整數(shù)LHD初始采樣

提高投影均勻性和空間均布性，能夠有效豐富設(shè)計(jì)空間的信息獲取[7]。拉丁超立方設(shè)計(jì)(Latin hypercube design,LHD)采樣自然滿足投影均勻性，已廣泛應(yīng)用于各類近似優(yōu)化算法的初始采樣中。LHD采樣點(diǎn)可表示為

(19)

由于LHD采樣具有隨機(jī)性，難以保證空間均布性，故本文將πi和ρi作為優(yōu)化參數(shù)以改善空間均布性，其中優(yōu)化準(zhǔn)則采用Maximin準(zhǔn)則[7]，即

(20)

對(duì)于整型變量，可將其視作實(shí)數(shù)型變量處理，并通過四舍五入取整以確保變量始終為整數(shù)[21]。

2.3 局部增強(qiáng)的EI加點(diǎn)準(zhǔn)則

(21)

則I(z)的數(shù)學(xué)期望E[I(z)]為

(22)

式中：Φ[·]為標(biāo)準(zhǔn)正態(tài)累積分布函數(shù)，φ[·]為標(biāo)準(zhǔn)正態(tài)概率密度函數(shù)。

EI加點(diǎn)準(zhǔn)則即選擇E[I(z)]值(稱為EI值)最大的樣本點(diǎn)。對(duì)于混合整數(shù)參數(shù)優(yōu)化問題，可在滿足等式約束、不等式約束以及整數(shù)約束[21]的可行集內(nèi)搜索EI值最大點(diǎn)作為新的采樣點(diǎn)，即

min(-E[I(z)])

(23)

本文針對(duì)上述問題提出一種局部增強(qiáng)的EI加點(diǎn)準(zhǔn)則(簡(jiǎn)稱為L(zhǎng)EEI加點(diǎn)準(zhǔn)則)。其基本思想是：在當(dāng)前最優(yōu)解附近限定一個(gè)小的局部搜索區(qū)域，在該區(qū)域內(nèi)采用EI加點(diǎn)準(zhǔn)則，以提升局部搜索性能；同時(shí)在全局搜索區(qū)域內(nèi)采用EI加點(diǎn)準(zhǔn)則，以保留原算法的全局收斂性；引入比例參數(shù)對(duì)局部搜索范圍進(jìn)行動(dòng)態(tài)調(diào)整，以控制算法收斂。

則引入局部搜索范圍的EI加點(diǎn)優(yōu)化問題(23)可改寫為

min(-E[I(z)])

(24)

式中：round(·)為四舍五入取整算子，rD為搜索空間縮放因子，z*(n)為當(dāng)前最優(yōu)解z*為的第n維值。

2.4 算法流程

基于LEEI加點(diǎn)準(zhǔn)則的改進(jìn)混合整數(shù)SAO算法的具體計(jì)算步驟為：

1)確定初始采樣點(diǎn)數(shù)量(可取ns=10nz)，求解參數(shù)優(yōu)化問題式(20)獲得初始樣本點(diǎn)集{z1,z2,…,zns}，計(jì)算各樣本點(diǎn)真實(shí)的目標(biāo)函數(shù)值和約束函數(shù)值；初始化迭代次數(shù)k=0，并設(shè)置最大迭代次數(shù)kmax(可取kmax=200)和rD的初始值(可取rD=0.2)。

(25)

3)令k=k+1，基于當(dāng)前樣本點(diǎn)集構(gòu)建目標(biāo)函數(shù)和約束函數(shù)的Kriging代理模型，分別求解參數(shù)優(yōu)化問題式(23)和式(24)獲得全局EI點(diǎn)zG和局部EI點(diǎn)zL。

4)若zG在局部搜索區(qū)域內(nèi)(即zG滿足式(24)中的約束3、4)，則將zL加入樣本點(diǎn)集中,計(jì)算zL真實(shí)的目標(biāo)函數(shù)值和約束函數(shù)值；否則將zG和zL同時(shí)加入樣本點(diǎn)集中，分別計(jì)算zG和zL真實(shí)的目標(biāo)函數(shù)值和約束函數(shù)值。

6)判斷收斂條件式(26)，其中εi(i=1,2,3,4)為收斂精度(可根據(jù)目標(biāo)函數(shù)和約束函數(shù)的值域?qū)?yīng)選擇εi的大小)。若滿足則停止；否則返回步驟3)。

(26)

3 仿真驗(yàn)證

本文燃燒室殼體材料采用40SiMnCrMoV，絕熱層材料采用EPDM，推進(jìn)劑采用HTPB/AP/Al，噴管喉襯和擴(kuò)張段材料采用碳/碳復(fù)合材料。共優(yōu)化17個(gè)設(shè)計(jì)變量，包括前、后翼槽數(shù)量共2個(gè)整型設(shè)計(jì)變量，以及裝藥參數(shù)、噴管參數(shù)、燃速系數(shù)共15個(gè)實(shí)數(shù)型設(shè)計(jì)變量，為保證幾何參數(shù)化模型的拓?fù)湔_，定義設(shè)計(jì)變量與幾何參數(shù)的關(guān)系為

(27)

設(shè)計(jì)變量搜索范圍見表1。設(shè)計(jì)約束、飛行條件定義見表2，其中上標(biāo)P和H分別為發(fā)動(dòng)機(jī)工作和不工作時(shí)的參數(shù)，其余均為全程參數(shù)。

表1 設(shè)計(jì)變量搜索范圍

表2 設(shè)計(jì)約束、飛行條件定義

本文基于PYTHON腳本語(yǔ)言驅(qū)動(dòng)ABAQUS 6.14實(shí)現(xiàn)發(fā)動(dòng)機(jī)殼體、絕熱層、裝藥和噴管的幾何參數(shù)化建模和四面體/三角形網(wǎng)格自動(dòng)劃分；采用文獻(xiàn)[15]方法實(shí)現(xiàn)燃面計(jì)算；采用文獻(xiàn)[19]方法實(shí)現(xiàn)改進(jìn)自適應(yīng)LGR偽譜法；采用經(jīng)典的DACE工具箱實(shí)現(xiàn)Kriging代理模型的構(gòu)建；基于MATLAB R2019b平臺(tái)編寫混合整數(shù)LHD初始采樣程序，以及基于LEEI加點(diǎn)準(zhǔn)則的改進(jìn)混合整數(shù)SAO算法(簡(jiǎn)稱為L(zhǎng)EEI算法)程序，其中混合整數(shù)參數(shù)優(yōu)化問題均采用MATLAB自帶的混合整數(shù)遺傳算法(Genetic algorithm, GA)進(jìn)行求解，LEEI算法中的收斂精度分別設(shè)置為：ε1=ε2=10-3，ε3=0，ε4=10-1)；仿真計(jì)算環(huán)境為Windows 7 2.10 GHz，32.0 GB內(nèi)存，優(yōu)化結(jié)果如圖4所示。

紅色為發(fā)動(dòng)機(jī)工作段，黑色為發(fā)動(dòng)機(jī)不工作段，實(shí)線為參數(shù)曲線，虛線為約束邊界

為說明設(shè)計(jì)優(yōu)化的有效性，以初始采樣中的最優(yōu)方案作為優(yōu)化前的方案(由算法流程步驟2)獲得)進(jìn)行對(duì)比，其中優(yōu)化前的最大航程為554.90 km，裝藥前、后翼槽數(shù)nw1、nw2分別為8和6，裝藥質(zhì)量mp為593.82 kg，不滿足等式約束要求；采用LEEI算法優(yōu)化后，最大航程為626.42 km(如圖4(a)所示)，相比優(yōu)化前提升了12.89%，裝藥的前、后翼槽數(shù)nw1、nw2均為8個(gè)，發(fā)動(dòng)機(jī)外徑dM、總長(zhǎng)度lM和裝藥質(zhì)量mp分別為0.44，2.90，589.97 kg，發(fā)動(dòng)機(jī)總質(zhì)量mM為737.40 kg，其中等式約束函數(shù)違反程度為0.03(≤ε4=10-1)，不等式約束函數(shù)違反程度為0(≤ε3=0)，均滿足式(2)中的等式和不等式約束要求，驗(yàn)證了LEEI算法對(duì)于等式和不等式約束處理的有效性。

由圖4(b)～(e)可知，飛行器的迎角、迎角變化率、駐點(diǎn)熱流密度、動(dòng)壓、法向過載和配平舵偏角均全程滿足相應(yīng)的約束要求。其中迎角在初始爬升時(shí)達(dá)到了約束上界，即最大爬升迎角是航程的限制因素之一；法向過載在發(fā)動(dòng)機(jī)工作時(shí)達(dá)到了約束上界，而駐點(diǎn)熱密度全程均未達(dá)到約束邊界，可知對(duì)于本文飛行器，法向過載為主要限制因素，其制約了可用迎角和發(fā)動(dòng)機(jī)推力大小；最高點(diǎn)動(dòng)壓達(dá)到了約束下限，由于控制系統(tǒng)對(duì)動(dòng)壓的限制，導(dǎo)致最大飛行高度受到制約，進(jìn)而影響了航程；配平舵偏角全程變化平穩(wěn)，僅在滑翔末段達(dá)到了約束下界，即限制了可用迎角。由圖4(d)可知，為提高航程，飛行末端的飛行速度達(dá)到了約束下界，而航跡角達(dá)到了上界，兩者均滿足相應(yīng)約束要求。由此可知，優(yōu)化結(jié)果滿足全部飛行約束要求，驗(yàn)證了LEEI算法結(jié)果的正確性。

下面為對(duì)比分析優(yōu)化結(jié)果，設(shè)置2個(gè)整型變量為若干固定值，并分別對(duì)剩余的15個(gè)實(shí)數(shù)型變量進(jìn)行優(yōu)化，優(yōu)化對(duì)比結(jié)果如圖5(a)～(d)所示。

圖5 不同翼槽數(shù)對(duì)應(yīng)的優(yōu)化結(jié)果

由圖5(a)可知，優(yōu)化方案的最大航程優(yōu)于其他情況；由圖5(b)～(d)可知，前8后8翼槽方案的推力接近于雙推力的形式，即采用短時(shí)間大推力輔助飛行器爬升，而后采用長(zhǎng)時(shí)間小推力持續(xù)加速，由于發(fā)動(dòng)機(jī)長(zhǎng)度、外徑以及裝藥質(zhì)量被嚴(yán)格限制，因此當(dāng)翼槽數(shù)較少時(shí)，則需要翼面和翼厚更大，以滿足尺寸質(zhì)量的約束條件，從而導(dǎo)致發(fā)動(dòng)機(jī)的初始推力相對(duì)更大，工作時(shí)間相對(duì)更短；對(duì)于翼槽數(shù)較多時(shí)，結(jié)果相反；此外，可以發(fā)現(xiàn)優(yōu)化方案對(duì)應(yīng)的關(guān)機(jī)點(diǎn)高度和速度也相對(duì)更大，從而其最大航程也更大。

為進(jìn)一步驗(yàn)證LEEI算法的性能，分別采用基于目標(biāo)函數(shù)預(yù)測(cè)值最小加點(diǎn)準(zhǔn)則的SAO算法(簡(jiǎn)稱為MP算法)、基于傳統(tǒng)EI加點(diǎn)準(zhǔn)則的SAO算法(簡(jiǎn)稱為EI算法)和GA算法對(duì)本文問題進(jìn)行求解與對(duì)比，其中MP算法來自文獻(xiàn)[14]，EI算法來自MATLAB軟件包Surrogate Model Toolbox。為剝離初始采樣對(duì)SAO算法的性能影響，LEEI、MP、EI算法均采用本文混合整數(shù)LHD初始采樣算法進(jìn)行初始采樣，采樣數(shù)量為170個(gè)，算法性能對(duì)比結(jié)果如圖6、7和表3所示。

圖6 目標(biāo)函數(shù)的收斂過程

圖7 約束函數(shù)的收斂過程

表3 算法性能對(duì)比

由圖6、7可知，隨著序列加點(diǎn)與迭代優(yōu)化，3種SAO算法均得到收斂。MP算法收斂速度最快，滿足全部約束要求，但收斂時(shí)的目標(biāo)函數(shù)值較大，即優(yōu)化的最大航程較小。LEEI和EI算法的結(jié)果與GA算法結(jié)果一致，均滿足全部約束要求，且優(yōu)于MP算法的結(jié)果，其中EI算法初始收斂速度較快，但后期收斂速度很慢，耗時(shí)函數(shù)調(diào)用次數(shù)約為MP算法的2倍；GA算法雖然優(yōu)化結(jié)果較好，但耗時(shí)函數(shù)調(diào)用次數(shù)過大，約為SAO算法的30～60倍；本文LEEI算法的耗時(shí)函數(shù)調(diào)用次數(shù)多于MP算法，但相比于EI算法顯著減少，由于增加了局部搜索控制策略，從而有效改善了收斂速度，同時(shí)保留了原算法的全局搜索性，驗(yàn)證了LEEI算法的有效性。

4 結(jié) 論

1)LEEI算法能夠有效求解SRM方案設(shè)計(jì)優(yōu)化的混合整數(shù)參數(shù)優(yōu)化問題，優(yōu)化后飛行器的最大航程提升了12.89%，且設(shè)計(jì)結(jié)果滿足全部約束條件。

2)LEEI算法的優(yōu)化結(jié)果與GA算法的優(yōu)化結(jié)果相一致，但耗時(shí)函數(shù)的調(diào)用次數(shù)僅為后者的2%。

3)LEEI算法具有一定全局搜索能力，且相比于EI算法有效提高了局部收斂速度。