張鼎逆，劉 毅

(同濟大學(xué)航空航天與力學(xué)學(xué)院，上海 200092)

可重復(fù)使用運載器(reusable launch vehicle，RLV)軌跡優(yōu)化問題屬于非線性、帶有狀態(tài)約束和控制約束的最優(yōu)控制問題.與巡航段和再入段不同，大氣環(huán)境、運動狀態(tài)以及發(fā)動機推力的劇烈變化，使得RLV的上升軌跡呈高度非線性，且受到控制約束、路徑約束和端點約束等條件的限制，增加了軌跡優(yōu)化設(shè)計的難度［1－2］.優(yōu)化求解方法的設(shè)計是上升段軌跡優(yōu)化的難點之一，求解的基本思路是用一個或多個有限維子問題來代替無限維的最優(yōu)控制問題.傳統(tǒng)的最速下降法、擬牛頓法和乘子罰函數(shù)法等方法被廣泛用于最小熱載，最大荷載等軌跡優(yōu)化問題，但普遍存在對初始值敏感等局限性.近年來，序列二次規(guī)劃法(sequential quadratic programming，SQP)和一些智能優(yōu)化方法如遺傳算法(genetic algorithm，GA)、粒子群算法(particle swarm optimization，PSO)等也越來越多被用來求解軌跡優(yōu)化問題［3－4］.智能算法具有很好的全局搜索性能，但是在局部最優(yōu)值附近收斂緩慢［5－6］.基于梯度的算法在局部搜索方面效率較高，可以獲得高精度的解，然而此類方法很大程度上依賴于初始點的選取，可行域的結(jié)構(gòu)和目標(biāo)函數(shù)的類型，不合適的初始點往往導(dǎo)致迭代發(fā)散或收斂至局部最優(yōu).為了提升算法性能，可以將多種優(yōu)化策略結(jié)合使用來彌補單獨使用的不足.如為了避免智能算法結(jié)果的抖動現(xiàn)象，往往采用梯度算法做后續(xù)計算.這種組合技術(shù)通常被稱為混合優(yōu)化方法，并已成功應(yīng)用于各類工程問題和研究中［7－9］.美國德克薩斯大學(xué)的 K.Subbarao 和 B.M.Shippey［10］采用混合遺傳算法求解了最短時間空間軌道轉(zhuǎn)移問題，通過打靶法和遺傳算法能夠獲得合適的解初值，但是遺傳算法在處理空間軌跡問題時明顯增加了問題的復(fù)雜度.英國南安普頓大學(xué)的 J.S.Alsumait等人［11］將GA、模式搜索和SQP的混合算法用于經(jīng)濟學(xué)優(yōu)化問題，對多種方法的混合進行了一系列嘗試.

文中根據(jù)軌跡優(yōu)化問題的分析方法、動態(tài)系統(tǒng)的離散技術(shù)、智能算法與基于梯度算法的混合優(yōu)化方案等理論，發(fā)展一種結(jié)合粒子群算法和序列二次規(guī)劃法的混合優(yōu)化方法，即HPSO法，并將其應(yīng)用于RLV上升段的軌跡優(yōu)化問題中.

1 上升段軌跡優(yōu)化問題

1.1 飛行運動學(xué)方程

RLV上升段屬于有動力的爬升飛行，上升過程的狀態(tài)量和控制量直接決定了巡航和再入的高度和速度等軌跡狀態(tài)，也間接影響了末端區(qū)域能量管理段和自動著陸段的能量控制和軌跡控制.因此，上升段在RLV任務(wù)飛行中有著重要的作用，有必要對其進行軌跡優(yōu)化設(shè)計.針對機載空中發(fā)射起飛的RLV，假設(shè)地球坐標(biāo)系靜止，考慮縱向平面內(nèi)的三自由度運動，上升段的軌跡模型如下［1］:

式中:狀態(tài)向量 X=［h，r，v，γ，m］包括高度 h，航程r，速度 v，航跡角 γ 和質(zhì)量 m;控制向量 U=［α，σ］包括迎角α，滾轉(zhuǎn)角σ;Re為地球半徑;L為升力;D為阻力;T為推力;g為重力加速度;Isp為發(fā)動機比沖.

1.2 發(fā)動機模型

采用火箭發(fā)動機和沖壓發(fā)動機的組合作為推進系統(tǒng)，以便在不同的馬赫數(shù)下達到最優(yōu)推力狀態(tài).組合式發(fā)動機包括如下2種工作模態(tài):①0≤M≤3，火箭發(fā)動機工作;②M＞3，沖壓發(fā)動機工作.

火箭發(fā)動機的比沖和推力如圖1所示.

圖1 火箭發(fā)動機比沖和推力

沖壓模態(tài)下發(fā)動機的推力由空氣流量˙mair和比沖 Isp決定，如下式［12］:

式中:下標(biāo)32.5表示該物理量在高度32.5 km處的值;ρ為大氣密度;vsound32.5為聲速;v32.5為 RLV 速度.沖壓發(fā)動機在32.5 km的空氣流量和比沖如圖2所示.

圖2 沖壓發(fā)動機空氣流量與比沖

1.3 氣動力模型

升力、阻力公式表示為如下標(biāo)準形式:

式中:CL為升力系數(shù);CD為阻力系數(shù);S為機翼參考面積.其中CL，CD與M和α的關(guān)系如圖3所示.

圖3 CL，CD與M和α的關(guān)系

1.4 約束條件

式中:RN為駐點曲率半徑.

2 優(yōu)化設(shè)計方法

本節(jié)對最優(yōu)控制問題的離散預(yù)處理，多鄰域改進PSO的求解步驟以及HPSO混合策略做簡要概述.SQP是一種應(yīng)用相對成熟的優(yōu)化方法，其詳細方法參考文獻［13］.

2.1 配點離散

考慮滿足運動方程、端點和過程條件等約束的RLV軌跡優(yōu)化問題為

上式可看作是對無限維狀態(tài)量X(t)與控制量U(t)的非線性優(yōu)化，屬于泛函優(yōu)化問題.由于粒子群算法和序列二次規(guī)劃法的搜索空間是歐幾里得空間，所以需要對X(t)和U(t)進行參數(shù)離散.文中采用直接配點法將軌跡優(yōu)化問題轉(zhuǎn)化為多維變量x的非線性NLP，即為

直接配點法是基于連續(xù)最優(yōu)控制系統(tǒng)與離散最優(yōu)控制系統(tǒng)的變換關(guān)系發(fā)展而來的，具有嚴格的數(shù)學(xué)理論推導(dǎo)［4］.

2.2 多鄰域改進PSO

傳統(tǒng)PSO在處理高復(fù)雜度優(yōu)化問題時存在著粒子信息共享緩慢，速度和位置更新容易超出界限等不足，因此，采用多鄰域共享方案，時變慣性因子和粒子限速等策略對傳統(tǒng)方法進行改進，可有效改善PSO的優(yōu)化性能.

在n維搜索空間中，m個粒子組成一個群體，每個粒子pi具有位置xi和速度vi等屬性，其中:xi=(xi1，xi2，…，xij，…，xin)，vi=(vi1，vi2，…，vij，…，vin)，1≤i≤m，1≤j≤n.

運用多鄰域改進PSO計算非線性規(guī)劃問題的初始值，算法步驟如下［14］:

1)初始化

2)鄰域劃分

將群體劃分成多個鄰域，鄰域內(nèi)的粒子共享此鄰域信息，特定粒子傳遞全局信息至其他鄰域，從而既保證了粒子搜索的獨立性，也容易得到全局最優(yōu)解.

3)適應(yīng)度計算與處理

令pbest，i為第i個粒子的歷史最優(yōu)值;gbest為群體的歷史最優(yōu)值;plbest，l為第l個鄰域的歷史最優(yōu)值.通過適應(yīng)度函數(shù)計算各粒子的適應(yīng)度值并進行最優(yōu)值的統(tǒng)計與更新.

4)速度更新

按照如下公式計算速度:

式中:w為速度慣性因子;c1，c2，c3為學(xué)習(xí)因子，通常取c1=c2=c3=2;εrand為0和1之間的隨機數(shù).采用時變慣性因子，即迭代初期w取值較大，隨著各粒子適應(yīng)度的接近而逐漸減小:

式中:kmax為最大迭代次數(shù);k為當(dāng)前迭代次數(shù).如果在最優(yōu)點附近速度較大，粒子將飛離最優(yōu)點，所以根據(jù)優(yōu)化狀態(tài)對w做進一步處理:

如果粒子的適應(yīng)度值變優(yōu)，速度將按照原有的慣性進行更新;若適應(yīng)度值變差，降低粒子速度以提高局部搜索能力.

Ti為學(xué)習(xí)能力函數(shù)，用于判別粒子獲取全局或鄰域知識的能力，Ti的表示為

試驗的荷載(p)-沉降(s)曲線如圖3、圖4所示。在較低的荷載范圍內(nèi)(200 kPa)，地基基本處于彈性變形狀態(tài)。

式中:H為可以獲取全局知識的粒子集合.

5)位置更新

按照如下公式計算位置可以起到限速作用，并且不會喪失粒子的搜索能量:

對于超出位置限制的粒子，通過如下公式可以滿足位置區(qū)間約束:

當(dāng)粒子位于局部最優(yōu)點附近時，速度的大幅衰減可導(dǎo)致其搜索能量的不足，采用變異率σ進行位置的隨機更新，可使粒子重新獲得搜索能力.速度判別函數(shù)與位置更新函數(shù)如下:

式中:σ∈［0，1］;常數(shù)τ是最小速度限制，大小和問題的求解精度有關(guān).

2.3 HPSO混合策略

HPSO混合算法的計算步驟如下:首先利用PSO的突變性搜索全局初始最優(yōu)解，當(dāng)不滿足問題精度時，跳出局部極值點繼續(xù)迭代計算;當(dāng)達到PSO終止條件后，利用SQP的數(shù)值解析計算得到更優(yōu)的最終值.PSO和SQP分別進行優(yōu)化計算，不僅迭代結(jié)果優(yōu)于PSO算法得到的初始解，而且提高了收斂速度和計算精度.HPSO的算法流程如圖4所示.

圖4 HPSO流程圖

3 上升軌跡優(yōu)化問題

3.1 最小燃料消耗問題

考慮DeltaIII運載火箭的上升段軌跡優(yōu)化問題，目標(biāo)函數(shù)為 max J=mf，狀態(tài)向量 X=［r，v，m］包括位置r，速度v和質(zhì)量m.控制向量 U=［ux，uy，uz］為推力矢量.運動方程、設(shè)計參數(shù)、動力參數(shù)和氣動參數(shù)等詳見文獻［15］.

經(jīng)過多次計算，取離散配點數(shù)20，PSO群體規(guī)模20，鄰域數(shù)3，最大最小慣性因子分別取0.7和0.1，位置變異率 0.01，迭代終止的誤差閥值設(shè)為10－6.表1為各算法的計算時間和最優(yōu)值比較.

表1 各次主要計算結(jié)果比較

HPSO法與文獻［15］以及單獨采用SQP的計算結(jié)果比較見圖5.

圖5 狀態(tài)與控制變量曲線

由圖5可見，HPSO的優(yōu)化結(jié)果與文獻［15］近乎重合，可以看出該算法的優(yōu)化精度很高.而SQP在沒有合適初始估計的情況下，結(jié)果相對較差.仿真結(jié)果表明，HPSO法優(yōu)化的軌跡滿足約束條件和目標(biāo)函數(shù)要求，并且在計算精度和效率上較SQP有了較大提高.

3.2 RLV亞軌道上升問題

上升段的高度、速度和迎角等參數(shù)對巡航以及再入有很大影響，速度越高表示飛行動能越大，入軌所需動力越少;高度越高，則其飛行勢能越大，入軌所需動力也越少;而迎角則是需要優(yōu)化的參變量.如何最大程度增加上升段RLV的有效載荷，是總體設(shè)計和性能分析的重要研究點.除了通過改進動力裝置和采用輕型材料2種途徑以外，在起飛重量恒定的條件下還可通過減少燃料消耗來實現(xiàn).

RLV上升段的初始與末端條件如下:h(t0)=9 km，v(t0)=238 m·s－1，γ(t0)=4°，h(tf)=32 km，v(tf)=2240 m·s－1，γ(tf)= －1°.

給出狀態(tài)量與控制量的有效范圍如下:－45°≤γ≤45°，－45°≤α≤45°，－45°≤σ≤45°.

設(shè)熱流密度、動壓和過載約束分別為420 kW·m－2，29 kPa 和2.5g.初始總重 m0=5000 kg，機翼參考面積S=12 m2，駐點曲率半徑RN=0.25 m.海平面大氣 密度 ρ0=1.225 kg·m－3，大 氣常量 β =0.000137.

取離散配點數(shù)n=30，HPSO的群體規(guī)模100，迭代次數(shù)20000.計算的末端重量mf為3272.6 kg.軌跡上升時間為218.74 s，部分狀態(tài)和控制曲線如圖6所示.

圖6 上升軌跡曲線

由圖6a高度曲線可見，火箭發(fā)動機工作位于9～16 km的高度范圍內(nèi)，隨后進入沖壓發(fā)動機的大推力狀態(tài)，高度為16～32 km.上升末端高度空氣稀薄，沖壓推力處于較小值，此時不能產(chǎn)生足夠的加速度，因此高度變化趨于平緩.由圖6a速度曲線可見，由于忽略了動力轉(zhuǎn)換的過渡段，火箭動力結(jié)束后速度曲線存在著一個拐點，拐點前后具有不同的加速度斜率.從高度和速度的曲線看出，到達末端高度時，速度為2240 m·s－1，高度為32 km，滿足上升段末端條件.由圖6b迎角控制曲線可見，發(fā)動機點火時迎角較大，隨著飛行過程中高度增加，迎角逐漸減小以有利于RLV的姿態(tài)控制，從而可以平滑過渡到下一階段.同時，控制迎角滿足約束以防止上升率過大，否則容易影響RLV的操縱性和穩(wěn)定性.由圖6b過載曲線可見，升力L和阻力D隨著迎角的增大而增加，法向過載也隨之增大.兩次點火位置對應(yīng)的局部迎角相對最大，因此局部最大過載也發(fā)生在這兩點.

4 結(jié)論

1)HPSO法不需要指定初始迭代點，SQP所需的合適初始估計由PSO優(yōu)化階段自動獲得，因此克服了SQP對初始估計過分依賴的缺點，避免了對初值敏感而導(dǎo)致陷入局部極值.隨著優(yōu)化系統(tǒng)規(guī)模和復(fù)雜度的增加，HPSO算法的這個特點更是顯而易見的.

2)HPSO法可以得到更優(yōu)的數(shù)值解，并且能有效縮短計算時間，減少大尺度復(fù)雜優(yōu)化問題的計算量.

3)對RLV上升段軌跡優(yōu)化問題的計算表明文中的混合算法可以很好地解決工程優(yōu)化問題.

References)

［1］Hirschel E H，Weiland C.Design of hypersonic flight vehicles:some lessons from the past and future challenges［J］.CEAS Space Journal，2010，1(1/2/3/4):3－22.

［2］Williams P.Hermite-Legendre-Gauss-Lobatto direct transcription in trajectory optimization［J］.Journal of Guidance， Control， andDynamics， 2009， 32(4):1392－1395.

［3］Huntington G T，Rao A V.Comparison of global and local collocation methods for optimal control［J］.Journal of Guidance，Control，and Dynamics，2008，31(2):432－436.

［4］Jain S，Tsiotras P.Trajectory optimization using multiresolution techniques［J］.Journal of Guidance，Control，and Dynamics，2008，31(5):1424 －1436.

［5］Wuerl A，Crain T，Braden E.Genetic algorithm and calculus of variations-based trajectory optimization technique［J］.Journal of Spacecraft and Rockets，2003，40(6):882－888.

［6］Pontani M，Conway B A.Particle swarm optimization applied to space trajectories［J］.Journal of Guidance，Control，and Dynamics，2010，33(5):1429 －1441.

［7］Fesanghary M，Mahdavi M，Minary-Jolandan M，et al.Hybridizing harmony search algorithm with sequential quadratic programming for engineering optimization problems［J］.Computer Methods in Applied Mechanics and Engineering，2008，197(33/…/40):3080－3091.

［8］Sentinella M R，Casalino L.Cooperative evolutionary algorithm for space trajectory optimization［J］.Celestial Mechanics＆ Dynamical Astronomy，2009，105(1/2/3):211－227.

［9］Vavrina M A，Howell K C.Global low-thrust trajectory optimization through hybridization of a genetic algorithm and a direct method［C］∥Proceedings of AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Reston:American Institute ofAeronauticsand Astronautics Inc.，2008，article number:2008 －6614.

［10］Subbarao K，Shippey B M.Hybrid genetic algorithm collocation method for trajectory optimization［J］.Journal of Guidance，Control，and Dynamics，2009，32(4):1396－1403.

［11］Alsumait J S，Sykulski J K，Al-Othman A K.A hybrid GA-PS-SQP method to solve power system valve-point economic dispatch problems［J］.Applied Energy，2010，87(5):1773－1781.

［12］Prasanna H M，Ghose D，Bhat M S，et al.Interpolationaware trajectory optimization for a hypersonic vehicle using nonlinear programming［C］∥Proceedings of Collection of Technical Papers-AIAA Guidance，Navigation，and Control Conference.Reston:American Institute of Aeronautics and Astronautics Inc.，2005:2160 －2171.

［13］謝 政，李建平，陳 摯.非線性最優(yōu)化理論與方法［M］.北京:高等教育出版社，2010.

［14］楊雪榕，梁加紅，陳 凌，等.多鄰域改進粒子群算法［J］.系統(tǒng)工程與電子技術(shù)，2010，32(11):2453－2458.Yang Xuerong，Liang Jiahong，Chen Ling，et al.Multineighborhood improved particle swarm optimization algorithm［J］.Systems Engineering and Electronics，2010，32(11):2453－2458.(in Chinese)

［15］Benson D.A Gauss pseudospectral transcription for optimal control［D］.Cambridge:Massachusetts Institute of Technology，2005.