陳永信

（上海機電工程研究所，上海 201109）

0 引言

滑翔飛行器通常需要通過動力裝置達到一定飛行高度和速度，而后依靠自身氣動升力進行無動力飛行［1］。該類飛行器的飛行速度高、機動范圍大且突防生存能力強，可應用于遠程物資的快速運送或精確打擊任務，具有重要的經濟和軍事價值［2］。

與傳統(tǒng)飛行器相比，滑翔飛行器的飛行環(huán)境更為嚴酷，約束條件更為嚴格，因此在該類飛行器的總體設計中，各學科更加緊密耦合［3］。傳統(tǒng)總體設計方法采用解析公式、經驗數據或簡化的學科模型進行分析，難以體現(xiàn)多學科耦合因素的影響，在多約束條件下獲得可行優(yōu)化解的難度大［4］。多學科設計優(yōu)化（multidisciplinary design optimization，MDO）方法可充分考慮多學科間的耦合效應，近年來成為航空航天領域的研究熱點［4-7］。但是MDO 方法中往往涉及耗時長的學科分析和大量的設計變量，直接進行總體設計優(yōu)化的難度極高且計算時間成本巨大，優(yōu)化效果通常并不理想，甚至難以獲得可行解［8］。代理優(yōu)化方法能夠有效針對優(yōu)化難度高且耗時長的問題，成為MDO研究中的關鍵技術［8-10］。目前代理優(yōu)化方法多應用于涉及學科較少的設計優(yōu)化問題中，例如飛行器氣動外形設計優(yōu)化［8，10-12］、結構設計優(yōu)化［13-14］等，對于滑翔飛行器與軌跡一體化設計優(yōu)化的相關研究尚較少。

本文針對滑翔飛行器與軌跡一體化設計優(yōu)化問題，建立飛行器無動力縱向運動方程和主要的飛行約束模型；為獲得飛行器的氣動特性參數，建立滑翔飛行器的氣動外形參數化模型，并采用非結構三角形面元法和流線追蹤技術實現(xiàn)考慮黏性影響的氣動特性計算；采用自適應偽譜法獲得多約束條件下飛行器的最大航程及最優(yōu)軌跡；基于Kriging代理模型和多采樣點高效全局優(yōu)化算法提升一體化設計優(yōu)化問題的求解效率，實現(xiàn)滑翔飛行器與軌跡一體化設計優(yōu)化，并對計算結果進行分析。

1 問題描述

1.1 優(yōu)化問題

針對尾部采用全動空氣舵進行姿態(tài)控制的面對稱升力體式滑翔飛行器，其初始飛行條件確定，通過一體化設計優(yōu)化滑翔飛行器的氣動外形參數和飛行軌跡，在滿足飛行器外形尺寸、容積率和飛行過程中主要約束的情況下，提升飛行器滑翔飛行段的最大航程能力。

1.2 飛行器縱向運動方程

關注飛行器的航程能力，為了方便描述滑翔飛行器的運動特性，基于瞬時平衡假設，考慮地球為無旋轉圓球，建立極坐標形式下無動力無側滑的飛行器縱向運動方程為

式中：h、v、λ、s為狀態(tài)變量，分別代表飛行器的飛行高度、速度、航跡角和航程；m為飛行器的質量，飛行過程中保持不變；引力常數μE=3.986×1014m3/s2；地球半徑Re=6 378 137 m；FLb和FDb分別為飛行器所受的配平升力和配平阻力，可表示為

式中：CLb和CDb分別為飛行器的配平升力系數和配平阻力系數，由氣動外形和飛行狀態(tài)共同決定；α、Ma和δ為分別為飛行器的迎角、馬赫數和俯仰舵偏角；SR為參考面積；q為飛行器的動壓；ρ為大氣密度，由飛行高度h決定，本文采用美國1976 標準大氣模型［15］進行計算。

1.3 飛行約束條件

1.3.1 俯仰力矩配平約束

初步設計階段，僅考慮飛行器的靜態(tài)俯仰力矩配平能力，對配平俯仰舵偏角進行限制，并留有一定余量，即

式中：Cmz為飛行器的俯仰力矩系數。對于給定的迎角α和馬赫數Ma，需要通過數值迭代求解式（3）獲得δ，則俯仰力矩配平約束可表示為

另一方面為減輕姿態(tài)控制系統(tǒng)的壓力，還需要對迎角變化速率進行限制，即

在軌跡設計優(yōu)化中，為了方便地描述控制約束，將式（1）中的狀態(tài)變量擴展為x=［h，v，λ，s，α］T，并將作為軌跡設計中待優(yōu)化的控制變量u，即

1.3.2 能量管理約束

為保證飛行器在滑翔段末端的能量管理要求，需要對末端的飛行高度和速度進行限制

式中：tf為滑翔段的末端時刻；hf和vf為滑翔段的末端高度和速度。

1.3.3 飛行過程約束

飛行器在大氣層內高速飛行，其載荷環(huán)境嚴酷，需要考慮飛行器的防熱、結構強度等約束，因此對飛行過程中飛行器的駐點熱流密度、法向過載ny和動壓q進行限制，即

式中：RN為駐點曲率半徑；完全氣體比熱比γ=1.4，高溫氣體比熱比γh=1.2；海平面引力加速度g0=9.806 6 m/s2。

2 滑翔飛行器氣動特性計算

2.1 滑翔飛行器氣動外形參數化建模

基于類型函數/形狀函數的幾何參數化建模方法（class function/shape function transformation technique，CST）具有控制參數少且外形表達豐富的優(yōu)點［16］，近年來在飛行器氣動外形參數化建模中得到了廣泛的應用［17-19］。通常參數化曲線方程可以統(tǒng)一表示為

式中：x、y和z為三維外形表面的物理坐標分量；ξ，η∈［0，1］為軸向歸一化變量；L0為三維外形的軸向特征長度，H0為曲線特征高度，W0為曲線特征寬度；C為類型函數；S為形狀函數，通常采用Bernstein 多項式或B樣條函數表示［20］；H為厚度函數。由式（9）可知，對于三維外形沿軸線方向的任意截面，其形狀類型和大小由軸向和橫向變量共同控制決定。因此，該方法可通過選擇不同的類型、形狀和厚度函數來獲得豐富的三維外形。

本文以文獻［19］中的面對稱升力體構型為基礎進行設計優(yōu)化，該類升力體外形示意圖見圖1。

采用式（9）形式對圖1 中的升力體外形進行參數化表達，其參數化方程為

圖1 升力體外形示意圖Fig.1 Schematic diagram of lifting body shape

式中：Lb為升力體的長度；Wb為底部截面的寬度；nb為升力體俯視面輪廓冪指數；yu和yd分別為升力體上和下表面的高度坐標；Hbu和Hbd分別為底部截面上和下輪廓曲線的最大高度；Nu和Nd分別為底部截面上和下輪廓曲線的形狀函數指數。

升力體尾部配置全動空氣舵以進行飛行器姿態(tài)控制，其幾何外形見圖2。

圖2 空氣舵外形示意圖Fig.2 Schematic diagram of air rudder shape

圖2中：Lw為舵面前緣未鈍化時的根弦長度，比例參數k1=Lw/Lb，用來控制根弦長度；hw為舵面后緣厚度，由結構載荷學科確定；rw為舵面前緣鈍化半徑，由氣動防熱學科確定；k2、k3和k4為比例參數，其中k2用來控制根梢比和后掠角，k3用來控制展弦比，k4用來控制舵面前緣的傾角。

為了安裝如圖2所示的全動空氣舵，需要對升力體尾部進行削平，見圖1。圖1 中WG為升力體削平后的寬度；Ls和hs為全動空氣舵的舵軸安裝位置參數。另一方面，考慮到滑翔飛行器頭部駐點的防熱要求，需要對頭部進行鈍化，定義鈍化半徑為rG，由氣動防熱要求確定。為了保證頭部區(qū)域的外形光滑且連續(xù)，需要對球頭與升力體連接處進行平滑過渡。此外，考慮到滑翔飛行器的內部需要裝載一定儀器部件和有效載荷，因此總體設計對飛行器的容積率RV提出要求，其定義為

式中：VL為削平后升力體的體積；SL為削平后升力體在xz平面的投影面積。

2.2 滑翔飛行器氣動特性計算

本文采用非結構三角形面元法和流線追蹤技術實現(xiàn)考慮黏性影響的氣動特性計算，具體計算方法見文獻［21］，本文不再贅述。對于傾斜轉彎（back-toturn，BTT）控制的面對稱滑翔飛行器，無側滑角，即側向力系數為0，因此，當獲得考慮黏性修正的升力系數CLb和阻力系數CDb，并給定參考坐標系下的飛行器質心后，飛行器繞質心的氣動力矩系數可表達為

式中：Cmx、Cmy、Cmz分別為飛行器繞質心的滾轉、偏航和俯仰力矩系數；Opx、Opy、Opz分別為參考坐標系下各三角形面元中心的坐標分量；Ogx、Ogy、Ogz分別為參考坐標系下飛行器質心的坐標分量；LR為飛行器的參考長度。

3 滑翔飛行器最大航程計算

當滑翔飛行器的質量和氣動特性確定后，其最大航程計算問題可描述為一個具有多種約束條件的連續(xù)時間最優(yōu)控制問題（continuous-time optimal control problem，COCP），即

式中：Eq.（1）（6）為式（1）和（6）組成的微分約束；Eq.（4）（5）（7）（8）為式（4）、（5）、（7）和（8）組成的不等式約束。

對于形如式（13）這類的多約束COCP，可采用自適應Legendre-Gauss-Radau（LGR）偽譜法進行有效的離散求解［22］。本文采用文獻［22］中基于Legendre 近似理論的自適應算法，并引入配點區(qū)間縮減技術［22］，該算法相比于傳統(tǒng)算法具有更高的求解效率和穩(wěn)健性。

4 基于代理優(yōu)化的滑翔飛行器設計優(yōu)化

4.1 Kriging代理模型

Kriging 代理模型的預測精度高，且能夠提供目標預測值的平均標準差信息，因此被廣泛應用于代理優(yōu)化中［8］。其假設目標函數為高斯靜態(tài)隨機過程的具體實現(xiàn)，可表示為

式中：f（X）為目標函數，X為f（X）的輸入變量；B（X）為回歸基函數矩陣，b為回歸系數矩陣，BT（X）b代表f（X）的數學期望；G（X）為高斯靜態(tài)隨機過程，其均值為0，方差為σ2，且對于不同的輸入變量Xi和Xj，隨機變量間存在相關性，由相關函數R描述。R隨Xi和Xj間距離的增大而減小，且滿足：|Xi-Xj|=0 時，R=1；|Xi-Xj|→∞時，R=0。R的表達式為

式中：nx為輸入變量的維數；下標k代表輸入變量第k維所對應的參數；Rk為相關基函數。

Kriging 代理模型利用已計算的ns個輸入變量（樣本點）處的目標函數值進行線性加權，根據無偏估計要求，采用拉格朗日乘子法推導獲得均方差MSE最小的權重系數，從而利用插值來預測未知輸入變量處的目標函數估計值，則最優(yōu)目標函數估計值的表達式為

式中：ω（X）為權重系數矩陣；fs為樣本點的目標函數值列向量；r（X）為預測點的相關函數列向量；Bs和Rs分別為樣本點的回歸基函數和相關函數矩陣。

式（16）可以化簡為式（14）的形式，即為

根據式（17）可以發(fā)現(xiàn)，b*和Gs僅取決于樣本點，因此對預測點的目標函數值進行預測時，僅需計算預測點的B（X）和r（X）。上述插值預測的計算耗時相對于原分析模型的計算耗時可以忽略不計。此外，預測點處目標函數值的平均標準差估計值為

回歸基函數和相關基函數的選擇對Kriging 代理模型的預測精度和算法魯棒性有較大影響。本文中回歸基函數選擇為一階回歸函數，相關基函數選擇為各向異性的指數高斯函數，其控制參數較多，能夠有效改善Kriging 模型的魯棒性［8］。則相應模型的表達式為

式中：θk和pk（1≤pk≤2）為模型訓練超參數。

根據“極大似然估計”準則，可建立無約束優(yōu)化問題來確定2nx個模型訓練超參數，則優(yōu)化問題為

4.2 代理優(yōu)化算法

序列代理優(yōu)化（sequential surrogate-based optimization，SSBO）算法是一類求解高計算成本優(yōu)化問題的算法，該算法首先采用較少的樣本點建立初始代理模型代替計算成本很高的目標函數或約束函數，并求解代理優(yōu)化問題，當求解結果不滿足求解要求時，按一定準則加入新樣本點來重構代理模型，并重復進行優(yōu)化直至結果滿足求解要求。

加點準則是序列代理優(yōu)化算法的核心［11］，對于基于Kriging代理模型的代理優(yōu)化算法，目前廣泛采用的加點準則主要包括：最小目標函數預測值（minimum predicted-value，MP）加點準則、最大目標函數值改善期望（expected improvement，EI）加點準則等。其中高效全局優(yōu)化（efficient global optimization，EGO）算法即采用EI加點準則［23］。

對于約束優(yōu)化問題，MP 準則通過求解約束優(yōu)化問題（21）獲得新樣本點［11］，EI準則通過求解約束優(yōu)化問題（22）獲得新樣本點［23］，即

式中：nk為采用代理模型的約束數量；ng為約束的總數量為第i個約束gi（X）的代理模型；E（X）為目標函數值改善的期望，其表達式為

式中：fmin為當前樣本點中目標函數真實值的最小值；Φ為標準正態(tài)累積分布函數；φ為標準正態(tài)概率密度函數。

MP 準則具有較好的局部收斂性，但容易收斂至局部最優(yōu)解［11］；EI 準則具有全局收斂性，但局部收斂速度較慢［23］。因此，為使得序列代理優(yōu)化算法兼具一定全局和局部搜索能力，本文在每次迭代過程中，在原樣本點集中同時加入MP 準則和EI 準則獲得的兩種樣本點來重構Kriging代理模型，并稱算法為多采樣點高效全局優(yōu)化算法（multi-sampling efficient global optimization，MEGO）。MEGO 算法的迭代收斂條件為

式中：ε1、ε2和ε3分別為相對誤差精度、代理模型精度和約束滿足精度；為第i次序列代理優(yōu)化的最優(yōu)解。

4.3 滑翔飛行器氣動外形與軌跡一體化設計優(yōu)化

在MEGO算法中，為了降低初始采樣時樣本點間的相關性，以提高Kriging 代理模型的近似精度，本文采用優(yōu)化拉丁超立方體試驗設計（optimal latin hypercube design，OLHD）算法［24］在設計空間內生成10nx個初始樣本點。對于式（21）和（22）的約束優(yōu)化問題，本文采用遺傳算法（genetic algorithm，GA）和序列二次規(guī)劃（seqential quadratic programming，SQP）算法的組合算法進行求解，以避免陷入局部最優(yōu)解。

對于本文滑翔飛行器氣動外形與軌跡一體化設計優(yōu)化問題，采用MEGO算法的求解流程見圖3。

圖3 MEGO算法流程Fig.3 Algorithm flow of MEGO method

5 仿真校驗

5.1 仿真環(huán)境

本文基于PYTHON 語言驅動ABAQUS 軟件，實現(xiàn)滑翔飛行器氣動外形參數化建模和非結構三角形表面網格的自動化劃分；基于MATLAB R2019b 平臺編寫氣動特性計算模塊、自適應LGR偽譜法軌跡優(yōu)化模塊以及代理優(yōu)化模塊程序，采用SNOPT 軟件包求解非線性規(guī)劃（nonlinear programming，NLP）問題；計算環(huán)境為Windows 10 2.60GHz，16.0 GB內存。

5.2 仿真結果

本文中滑翔飛行器的質量m為固定值，假設其質心在縱軸上，并定義質心到頭部頂點的距離與飛行器長度Lb的比值為質心位置系數k0。本文優(yōu)化變量X=［Lb，nb，Hbu，Hbd，Nu，Nd，k0，k1，k2，k3，k4］，包括10 個滑翔飛行器氣動外形參數以及位置質心系數，優(yōu)化設計變量的搜索范圍設置見表1。

表1 優(yōu)化變量搜索空間Tab.1 Searching space of optimization variables

其他氣動外形參數則根據發(fā)射平臺幾何尺寸、結構載荷和熱防護等要求預先確定為固定值，其中包括：Wb=1.1 m，WG=0.9 m，rG=0.03 m，rw=0.003 m和hw=0.02 m?？諝舛娑孑S安裝位置參數設置為：Ls=0.51Lw和hs=Hbu。在氣動特性計算時設置：SR=0.151 8 m2和LR=Lb。滑翔飛行器飛行初始條件和約束邊界參數為固定值，且，相應參數的設置見表2。

表2 飛行初始條件和約束定義Tab.2 Definition of flight initial conditions and constraints

為了分析滑翔飛行器容積率與最大航程的關系，本文分別對10 種不同容積率約束RV≥RVmin的情況進行了求解，其中RVmin分別取0.24、0.26、0.28、0.30、0.32、0.34、0.36、0.38、0.40 和0.42。MEGO 算法的收斂精度設置為：ε1=ε2=ε3=0.001?；栾w行器氣動外形與軌跡一體化設計優(yōu)化結果見圖4～13。為了驗證MEGO 算法的有效性，采用EGO 算法作為對比算法，對相同問題進行求解，算法性能對比結果見表3，表中：Rf為優(yōu)化結果相對基準方案提升的百分比；Nf為圖3中學科分析模塊的調用次數。其中基準方案選擇為初始樣本點中滿足全部約束條件的最優(yōu)點。

表3 不同算法的優(yōu)化結果Tab.3 Optimization results of different methods

由圖4 和5 可知，優(yōu)化結果中飛行器末端飛行高度和速度均滿足端點約束要求。由圖7～10 可知，整個飛行過程中，最大迎角變化率不超過2（°）/s，即迎角變化平緩，滿足約束要求；熱流密度、法向過載和動壓均滿足過程約束要求，其中最大熱流密度均達到5 500 kW/m2的約束邊界值，即熱流密度為本文滑翔飛行器最大航程的限制因素。由圖11可知，俯仰力矩系數均不大于4×10-4，滿足配平精度要求。俯仰舵偏角配平約束限制了飛行器的可用配平迎角，由圖12可知，在考慮常值干擾俯仰力矩系數后，整個飛行過程中的配平俯仰舵偏角仍能滿足約束要求，即設計優(yōu)化結果具有一定魯棒性。

圖4 高度與航程Fig.4 Altitude versus range

圖5 速度與時間Fig.5 Speed versus time

圖6 迎角與時間Fig.6 Angle of attack versus time

圖7 迎角變化率與時間Fig.7 Change rate of angle of attack versus time

圖8 熱流密度與時間Fig.8 Heat flux versus time

圖9 法向過載與時間Fig.9 Normal overload versus time

圖11 俯仰力矩系數與時間Fig.11 Pitching moment coefficient versus time

由圖13可知，滑翔飛行器的最大航程與容積率之間存在指標沖突，即最大航程隨容積率的增大而減小，且對于本文飛行器，最大航程與容積率呈現(xiàn)近似線性關系。對于不同容積率約束下的優(yōu)化結果，容積率均達到約束下限值，因此容積率亦為飛行器最大航程的限制因素。此外，對于優(yōu)化結果，隨著容積率的增大，Hbu隨之明顯增大以滿足容積率要求；Hbd均為下限值，即增大Hbd不利于提升最大射程；Lb均為上限值，即較大長細比有利于提升最大射程；nb均為較小值，表明飛行器的頭部寬度較??；空氣舵的外形尺寸由軌跡優(yōu)化確定，即在滿足配平能力約束下降低阻力，優(yōu)化結果中，空氣舵均具有較小的根梢比，且舵面前緣具有一定后掠角。質心位置系數k0分布在上限值附近，結合圖12 可知，飛行器飛行過程中先處于臨界穩(wěn)定狀態(tài)，而后處于靜穩(wěn)定狀態(tài)，且采用較小的舵俯仰偏角即可實現(xiàn)配平。

圖10 動壓與時間Fig.10 Dynamic pressure versus time

圖12 配平俯仰舵偏角與時間Fig.12 Equilibrium pitching rudder angle versus time

圖13 射程與容積率Fig.13 Range versus volume efficiency

由表3 可知，經過MEGO 算法優(yōu)化后，飛行器的最大航程相對基準方案的最大航程提升了近10%左右，當容積率約束邊界為RVmin=0.40 和RVmin=0.42時，初始采樣無法得到滿足全部約束條件的基準方案，但采用MEGO算法即可獲得最優(yōu)可行解。根據表3 可以發(fā)現(xiàn)，MEGO 和EGO 算法得到的最優(yōu)結果一致，且均滿足全部約束要求，但MEGO 算法對學科分析模塊的調用次數均少于EGO 算法對學科分析模塊的調用次數，因此MEGO 算法的優(yōu)化計算時間較EGO 算法的優(yōu)化計算時間也有顯著降低。由于EI 準則兼顧指標函數和代理模型的預測誤差，經過數次加點后，最優(yōu)解附近的EI 函數值可能較小，使得后續(xù)加點多處于樣本點較為稀疏的區(qū)域，導致最優(yōu)解附近的代理模型精度改善較慢，而MEGO 算法加入MP 點則能有效改善EGO算法局部收斂速度，驗證了本文算法的有效性。

6 結束語

本文針對滑翔飛行器氣動外形與軌跡一體化設計優(yōu)化問題，建立滑翔飛行器氣動特性計算和最大航程評估模型，并基于MEGO算法求解了相應一體化設計優(yōu)化問題，計算結果表明：

1）MEGO算法的局部收斂速度相比于EGO算法的局部收斂速度有明顯提升，能夠有效降低學科分析模塊的調用次數，從而顯著提升優(yōu)化效率。

2）MEGO 算法能夠在多學科約束條件下，有效地獲得可行優(yōu)化方案，且性能指標提升明顯，該方法具有在初始方案設計階段縮短設計周期的潛力。