張鼎逆，劉 毅

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

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

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

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

1.1 飛行運(yùn)動(dòng)學(xué)方程

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

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

1.2 發(fā)動(dòng)機(jī)模型

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

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

圖1 火箭發(fā)動(dòng)機(jī)比沖和推力

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

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

圖2 沖壓發(fā)動(dòng)機(jī)空氣流量與比沖

1.3 氣動(dòng)力模型

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

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

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

1.4 約束條件

式中:RN為駐點(diǎn)曲率半徑.

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

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

2.1 配點(diǎn)離散

考慮滿足運(yùn)動(dòng)方程、端點(diǎn)和過(guò)程條件等約束的RLV軌跡優(yōu)化問(wèn)題為

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

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

2.2 多鄰域改進(jìn)PSO

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

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

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

1)初始化

2)鄰域劃分

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

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

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

4)速度更新

按照如下公式計(jì)算速度:

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

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

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

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

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

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

5)位置更新

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

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

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

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

2.3 HPSO混合策略

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

圖4 HPSO流程圖

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

3.1 最小燃料消耗問(wèn)題

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

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

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

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

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

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

3.2 RLV亞軌道上升問(wèn)題

上升段的高度、速度和迎角等參數(shù)對(duì)巡航以及再入有很大影響，速度越高表示飛行動(dòng)能越大，入軌所需動(dòng)力越少;高度越高，則其飛行勢(shì)能越大，入軌所需動(dòng)力也越少;而迎角則是需要優(yōu)化的參變量.如何最大程度增加上升段RLV的有效載荷，是總體設(shè)計(jì)和性能分析的重要研究點(diǎn).除了通過(guò)改進(jìn)動(dòng)力裝置和采用輕型材料2種途徑以外，在起飛重量恒定的條件下還可通過(guò)減少燃料消耗來(lái)實(shí)現(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è)熱流密度、動(dòng)壓和過(guò)載約束分別為420 kW·m－2，29 kPa 和2.5g.初始總重 m0=5000 kg，機(jī)翼參考面積S=12 m2，駐點(diǎn)曲率半徑RN=0.25 m.海平面大氣 密度 ρ0=1.225 kg·m－3，大 氣常量 β =0.000137.

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

圖6 上升軌跡曲線

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

4 結(jié)論

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

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

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

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ìn)粒子群算法［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.