基于子集模擬優(yōu)化的空天飛機(jī)再入軌跡混合優(yōu)化方法

朱俊杰 余雄慶

1.飛行器先進(jìn)設(shè)計(jì)技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南京210016

2.南京航空航天大學(xué)，南京210016

空天飛機(jī)軌跡優(yōu)化是在飛行動(dòng)力學(xué)和物理學(xué)約束的條件下，尋求某種意義下最優(yōu)解的過程，它貫穿于整個(gè)飛行器設(shè)計(jì)過程中，影響著總體、氣動(dòng)布局、制導(dǎo)控制、動(dòng)力和結(jié)構(gòu)等多個(gè)分系統(tǒng)的設(shè)計(jì)。

軌跡優(yōu)化是一類最優(yōu)控制問題，屬于復(fù)雜的大規(guī)模非線性優(yōu)化范疇。軌跡優(yōu)化設(shè)計(jì)方法主要有基于極大值原理的間接法和基于非線性規(guī)劃理論的直接法。間接法的優(yōu)點(diǎn)是計(jì)算精度高，且能滿足最優(yōu)必要性條件，但由于繁雜的數(shù)學(xué)推導(dǎo)和難于求解的兩點(diǎn)邊值問題，目前已很少使用間接法求解軌跡優(yōu)化問題。隨著計(jì)算機(jī)技術(shù)的發(fā)展，目前的研究更多傾向于使用直接法求解。采用直接法求解軌跡優(yōu)化問題時(shí)，先將最優(yōu)控制問題轉(zhuǎn)化為非線性規(guī)劃問題(nonlinear programming，NLP)，對(duì)控制變量和狀態(tài)變量進(jìn)行離散，對(duì)狀態(tài)方程、端點(diǎn)約束以及路徑約束進(jìn)行節(jié)點(diǎn)轉(zhuǎn)化，再選擇優(yōu)化算法求解決策向量和目標(biāo)函數(shù)［1］。

優(yōu)化算法一般可分為2類:1)經(jīng)典優(yōu)化算法，例如序列二次規(guī)劃法(簡(jiǎn)稱SQP)，梯度下降法和罰函數(shù)方法［2］等;2)智能優(yōu)化算法，例如遺傳算法(Genetic Algorithm，簡(jiǎn)稱GA)，進(jìn)化策略(Evolution Strategy，簡(jiǎn)稱ES)和粒子群優(yōu)化(Particle Swarm Optimization，簡(jiǎn)稱PSO)［3］等。基于梯度的優(yōu)化算法的優(yōu)點(diǎn)是收斂速度快，但對(duì)初值較為敏感，且容易陷入局部最優(yōu)點(diǎn)。對(duì)于空天飛機(jī)再入軌跡問題，有數(shù)百甚至數(shù)千個(gè)設(shè)計(jì)變量，面對(duì)這么多變量，以人工方式給出一個(gè)合理的初始值，幾乎不可能。智能優(yōu)化算法的優(yōu)點(diǎn)是無需給出初值，穩(wěn)健性好，但收斂速度慢，優(yōu)化結(jié)果不夠精確。為了克服經(jīng)典優(yōu)化算法和智能優(yōu)化算法的缺點(diǎn)，近年來，有學(xué)者將GA與SQP相結(jié)合來求解軌跡優(yōu)化問題，如Vavrina采用GA與梯度迭代算法相結(jié)合的方法求解低推力軌跡優(yōu)化問題［4］。但是，用GA-SQP混合算法求解再入軌跡優(yōu)化問題時(shí)，通常選取的配點(diǎn)數(shù)為30～50［1］，配點(diǎn)數(shù)繼續(xù)增加的話，GA為SQP提供的初始點(diǎn)的質(zhì)量會(huì)越來越差，導(dǎo)致SQP最終不易收斂。而在某些軌跡優(yōu)化問題中，取30～50個(gè)配點(diǎn)可能精度上不能滿足設(shè)計(jì)要求。

本文嘗試采用一種新的隨機(jī)優(yōu)化算法-子集模擬優(yōu)化算法［5］(Subset Simulation Optimization，簡(jiǎn)稱SSO)，來改進(jìn)現(xiàn)有的混合優(yōu)化方法，這種方法簡(jiǎn)稱為SSO-SQP。有關(guān)研究表明:SSO與GA相比，具有更快的收斂速度［6］。本文以典型空天飛機(jī)再入軌跡混合優(yōu)化問題為算例，測(cè)試SSO-SQP的有效性。

1 再入軌跡優(yōu)化問題的描述

1.1 運(yùn)動(dòng)方程及設(shè)計(jì)變量

飛行器再入大氣層時(shí)的三自由度運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)方程組如下［7］:

式中:m為空天飛機(jī)質(zhì)量(kg)，R為飛行器質(zhì)心到地心的距離(m)，Φ為經(jīng)度(rad)，θ為緯度(rad)，v為速度(m/s)，γ為航跡角(rad)，ψ為航跡方向角(rad)，α為迎角(rad)，β為傾側(cè)角(rad)，D為阻力，L為升力，g為不同高度的重力加速度。定義x=［R，Φ ，θ，v，γ，ψ］為狀態(tài)變量，u=［α，β］為控制變量，設(shè)計(jì)變量包括狀態(tài)變量和控制變量。D和L由下式確定:

式中:g0為地球表面的重力加速度，大小為9.8m/s2;H為空天飛機(jī)離地面高度;R0為地球半徑，大小為6357km;S為空天飛機(jī)參考面積;CD和CL分別為阻力系數(shù)與升力系數(shù)。

在大氣飛行過程的熱流密度為:

1.2 歸一化處理

對(duì)于再入軌跡優(yōu)化問題，不同變量之間的數(shù)值量級(jí)相差較大，例如R是106級(jí)別，而經(jīng)度和緯度化為弧度后是1～10級(jí)別，兩者數(shù)值上相差太大，會(huì)導(dǎo)致數(shù)值奇異，優(yōu)化計(jì)算失敗。為了增加優(yōu)化計(jì)算的穩(wěn)健性，需要對(duì)原始模型歸一化處理。令:

1.3 優(yōu)化目標(biāo)

優(yōu)化目標(biāo)的選擇主要由空天飛機(jī)的設(shè)計(jì)及任務(wù)要求確定，通常有如下4種優(yōu)化目標(biāo):1)以飛行航程最大作為性能指標(biāo);2)以再入過程中總吸熱量最小作為性能指標(biāo)［9］;3)以最小能量消耗作為性能指標(biāo)［10］;4)以到達(dá)指定目標(biāo)點(diǎn)飛行時(shí)間最短作為性能指標(biāo)。

1.4 約束條件

再入軌跡優(yōu)化中的約束包括終端約束、路徑約束、動(dòng)壓及熱流峰值約束、過載約束和控制量約束［9］等。

2 SSO-SQP混合優(yōu)化方法

SSO-SQP混合優(yōu)化方法的思路是:首先在整個(gè)設(shè)計(jì)空間使用SSO算法進(jìn)行全局搜索，迭代適當(dāng)步數(shù)。待SSO迭代終止后，將SSO計(jì)算結(jié)果作為SQP優(yōu)化算法的初始點(diǎn)，使用SNOPT軟件包進(jìn)行再次優(yōu)化，得出收斂后的結(jié)果。以下對(duì)SSO和SQP算法進(jìn)行簡(jiǎn)要介紹。

2.1 SSO算法

子集模擬優(yōu)化算法的基本思想是極值問題(優(yōu)化問題)，可以看作小失效概率問題(可靠性問題)［5］。在進(jìn)行轉(zhuǎn)換時(shí)，首先考慮如下無約束單目標(biāo)全局優(yōu)化問題:

其中f:Rn→Ω是一個(gè)實(shí)值函數(shù)，x是設(shè)計(jì)變量并且Ω是一個(gè)封閉有邊界的集合。全局最小值fopt為滿足下式的解xopt:

為了在可靠性分析框架中研究?jī)?yōu)化問題，人為設(shè)定設(shè)計(jì)變量為隨機(jī)變量，這個(gè)“隨機(jī)化”處理使得目標(biāo)函數(shù)f也變?yōu)橐粋€(gè)隨機(jī)變量，且擁有自己的概率密度函數(shù)和累積分布函數(shù)。根據(jù)累積分布函數(shù)的數(shù)學(xué)定義可知，累積分布函數(shù)曲線是單調(diào)非減，并且右連續(xù)，有且僅有唯一的最小值fopt。構(gòu)造一個(gè)可靠性問題如下:

其中，失效事件為F={f(x)≤fopt}。顯然，與方程(14)有關(guān)的失效概率是0，因?yàn)閒opt是全局最小。然而所關(guān)心的不是這個(gè)失效概率，而是目標(biāo)函數(shù)最小值點(diǎn)或者包含目標(biāo)函數(shù)最小值的區(qū)域。失效事件的變量區(qū)域與函數(shù)取得最小值的變量區(qū)域是對(duì)應(yīng)的，從而實(shí)現(xiàn)了最優(yōu)化問題向可靠性問題的轉(zhuǎn)化，進(jìn)而可用子集模擬方法解決優(yōu)化問題［6］。

SSO算法的計(jì)算步驟［5］如下:

首先，定義一個(gè)約束違反函數(shù)用于處理約束條件并將樣本排序。對(duì)于一個(gè)不等式約束條件，gi(x)≤0對(duì)應(yīng)約束違反函數(shù)定義為

2)根據(jù)人工分布類型，利用直接的Monte Carlo模擬來生成 N個(gè)獨(dú)立同分布的樣本{x1，x2，…，xN}。計(jì)算所有隨機(jī)樣本的目標(biāo)函數(shù)值，并將其用雙標(biāo)準(zhǔn)排序方法升序排列，即有{W1，l:l=1，…，N}。此處，下標(biāo)“1”表示這些目標(biāo)函數(shù)的值對(duì)應(yīng)于第1層模擬。設(shè)定xN為當(dāng)前序列中的最佳求解方案，而x1為最差求解方案。合理選取p1和N，使得N(1-p1)為整數(shù)，從唯一序列中得到第N(1-p1)個(gè)樣本xN(1-p1)以及相應(yīng)的W1，N(1-p1)和Fcon(xN(1-p1))。選取第1組中間事件的臨界值為:W1=W1，N(1-p1)以及 Fc1=F(xN(1-p1))。其中，下標(biāo)“1”表示第1層模擬，則第1個(gè)中間事件F1定義為:

中間事件F1的定義也是滿足雙標(biāo)準(zhǔn)排序準(zhǔn)則的，以約束條件優(yōu)先。分開考慮{Fc1=0}和{Fc1＜0}兩種情況，容易得到P(F1)的估計(jì)值為p1。無論何種情況，總有p1N個(gè)樣本屬于F1，為下一層模擬生成新樣本提供“種子”樣本。

3)采用改進(jìn)的Metropolis-Hasting算法從中間事件中生成條件樣本［5］。在第k層模擬中，從Fk-1中每個(gè)樣本開始，可以生成相同分布的Markov鏈。再一次計(jì)算全局約束函數(shù)值和目標(biāo)函數(shù)值，對(duì)該層內(nèi)的所有隨機(jī)樣本進(jìn)行雙標(biāo)準(zhǔn)排序，確定序列中第N(1-pk)個(gè)樣本xN(1-pk)。然后根據(jù)Fk選擇樣本，F(xiàn)k定義如下:

2.2 SQP算法

SQP法的思想是通過一系列二次規(guī)劃(QP)子問題來求解非線性問題。SQP的基本結(jié)構(gòu)包括主迭代和次迭代，主迭代逐漸收斂于問題的最優(yōu)解，主迭代中的QP子問題產(chǎn)生下次主迭代的搜索方向。QP子問題本身是一個(gè)迭代過程，即SQP的次迭代是QP子問題迭代。本文的SQP算法來源于SNOPT程序［11］。SNOPT程序采用了一種改進(jìn)的序列二次規(guī)劃算法(SQP)，適用于求解大規(guī)模非線性稀疏的優(yōu)化問題。

3 算例

空天飛機(jī)再入軌跡優(yōu)化算例取自文獻(xiàn)［12］，并在該算例基礎(chǔ)上增加了熱流密度約束、動(dòng)壓約束和過載約束。

飛行器的氣動(dòng)模型如下:

空天飛機(jī)參考面積S=250.237m2，質(zhì)量m=92162kg，前緣半徑RN=1m。不同高度下大氣參數(shù)參見文獻(xiàn)［13］。

本例的設(shè)計(jì)變量包括狀態(tài)變量和控制變量。x=［R，Φ ，θ，v，γ，ψ］為狀態(tài)變量，u=［α，β］為控制變量。該優(yōu)化問題的初始狀態(tài)是:

優(yōu)化目標(biāo)是通過控制迎角α和傾側(cè)角β，使航程最大，對(duì)于本例來說，與緯度最大等價(jià)。本例以緯度最大為最優(yōu)目標(biāo):J=maxθ(tf)。在本例中SSO算法的參數(shù)設(shè)定為:單層樣本量取100，條件概率取0.51，最大層數(shù)取100，終止精度取10-6。離散方法采用Hermite-Simpson法。

為了驗(yàn)證SSO-SQP算法的魯棒性，計(jì)算了不同配點(diǎn)數(shù)情況下的軌跡優(yōu)化問題。計(jì)算機(jī)配置為3.10GHz，內(nèi)存2GB。SSO迭代100步，7s內(nèi)就能計(jì)算完畢。最終結(jié)果見表1。

表1 不同配點(diǎn)數(shù)的SSO-SQP優(yōu)化結(jié)果

由表1可見:1)隨著配點(diǎn)數(shù)的增加，優(yōu)化結(jié)果少量增加，但是當(dāng)配點(diǎn)數(shù)增加到160時(shí)，最優(yōu)值基本不變了;2)隨著配點(diǎn)數(shù)的增加，SQP迭代步數(shù)和優(yōu)化時(shí)間也逐漸增加，當(dāng)配點(diǎn)數(shù)增加到160時(shí)，SQP迭代步數(shù)和優(yōu)化時(shí)間開始大幅度增加，所以在優(yōu)化計(jì)算時(shí)，配點(diǎn)數(shù)可以選在100～140之間，這樣精度上既能滿足要求，又相對(duì)節(jié)省計(jì)算時(shí)間。

當(dāng)配點(diǎn)數(shù)為140時(shí)，迎角和傾側(cè)角隨時(shí)間變化曲線如圖1和2所示。熱流密度、動(dòng)壓和過載曲線如圖3～5所示，可以看出優(yōu)化結(jié)果能很好地滿足約束要求。

圖1 迎角隨時(shí)間的變化

圖2 傾側(cè)角隨時(shí)間的變化

圖3 熱流密度隨時(shí)間的變化

圖4 動(dòng)壓隨時(shí)間的變化

圖5 過載隨時(shí)間的變化

4 結(jié)論

提出了一種SSO-SQP混合優(yōu)化方法來求解空天飛機(jī)再入軌跡優(yōu)化問題。SSO算法屬于隨機(jī)優(yōu)化算法，不需要初始值，魯棒性強(qiáng)。算例證明，僅需迭代100步(7s內(nèi)就能計(jì)算完畢)就能獲得一個(gè)較好的初始值，因此用SSO作為獲得SQP算法初始值的手段是有效的。在SSO給出了合理初始值后，利用SQP算法收斂速度快、精度高的優(yōu)點(diǎn)，能高效地尋到最優(yōu)軌跡。算例表明:

1)該算法能有效地求解空天飛機(jī)軌跡優(yōu)化問題，能根據(jù)精度和計(jì)算時(shí)間要求，設(shè)定配點(diǎn)數(shù);

2)SSO-SQP算法魯棒性較強(qiáng)，當(dāng)配點(diǎn)從20增加到180時(shí)，均能優(yōu)化出正確結(jié)果;

3)通常情況下，選擇配點(diǎn)數(shù)為100～140，算例的設(shè)計(jì)變量數(shù)為809～1129，約束方程有903～1263個(gè)。經(jīng)SSO優(yōu)化得到初始值后，SQP僅需迭代300步左右就能收斂，可見該混合算法優(yōu)化效率較高。

［1］ 張鼎逆，劉毅.基于改進(jìn)遺傳算法和序列二次規(guī)劃的再入軌跡優(yōu)化［J］.浙江大學(xué)學(xué)報(bào)(工學(xué)版)，2014，48(1):161-167.(Zhang Dingni，Liu Yi.Reentry trajectory optimization based on improved genetic algorithm and sequential quadratic programming［J］.Journal of Zhejiang University(Engineering Science)，2014，48(1):161-167.)

［2］ 陳功，傅瑜，郭繼峰.飛行器軌跡優(yōu)化方法綜述［J］.飛行力學(xué)，2011，29(4):1-5.(Chen Gong，F(xiàn)u Yu，Guo JiFeng.Survey of aircraft trajectory optimization methods［J］.Flight Dynamics，2011，29(4):1-5.)

［3］ 楊希祥，李曉斌，肖飛，等.智能優(yōu)化算法及其在飛行器優(yōu)化設(shè)計(jì)領(lǐng)域的應(yīng)用綜述［J］.宇航學(xué)報(bào)，2009，30(6):2051-2061.(Yang Xixiang，Li Xiaobin，Xiao Fei，et al.Overview of intelligent optimization algorithm and its application in flight vehicles optimization design［J］.Journal of Astronautics，2009，30(6):2051-2061.)

［4］ Vavrina M A，Howell K C.Global low-thrust trajectory optimization through hybridization of a genetic algorithm and a direct method［C］//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Honolulu:［s.n.］，2008:1-27.

［5］ 李紅雙，馬遠(yuǎn)卓.結(jié)構(gòu)可靠性分析與隨機(jī)優(yōu)化設(shè)計(jì)的統(tǒng)一方法［M］.北京:國(guó)防工業(yè)出版社，2015:223-260.(Li Hongshuang，Ma Zhuoyuan.Unified Methods for Structual Reliability Analysis and Stochastic Optimization Design［M］.Beijing:National Defense industry press，2015:223-260.)

［6］ Li H S，Au S K.Design optimization using subset simulation algorithm［J］.Structural Safety，2010，32(6):384-392.

［7］ Rahimi A，Kumar K D，Alighanbari H.Particle swarm optimization applied to spacecraft reentry trajectory［J］.Journal of Guidance，Control，and Dynamics，2013，36(1):307-310.

［8］ 雍恩米，唐國(guó)金，陳磊.基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化［J］.宇航學(xué)報(bào)，2008，29(6):1766-1772.(Yong Enmi，Tang Guojin，Chen Lei.Rapid trajectory optimization for hypersonic reentry vehicle［J］.Journal of Astronautics，2008，29(6):1766-1772.)

［9］ 湯亮，楊建民，陳風(fēng)雨.多約束條件下的升力滑翔式再入軌跡優(yōu)化［J］.導(dǎo)彈與航天運(yùn)載技術(shù)，2013(1):1-5.(Tang Liang，Yang Jianmin，Chen Fengyu.Trajectory optimization with multi-constraints for a lifting reentry vehicle［J］.Missiles and Space Vehicles，2013(1):1-5.)

［10］ 陳剛，萬自明，徐敏，等.飛行器軌跡優(yōu)化應(yīng)用遺傳算法的參數(shù)化與約束處理方法研究［J］.系統(tǒng)仿真學(xué)報(bào)，2005，17(11):2737-2740.(Chen Gang，Wan Ziming，Xu min，et al.Parameterization method and constraints transformation method for RLV reentry trajectory optimization using genetic algorithm［J］.Journal of System Simulation，2005，17(11):2737-2740.)

［11］ Gill P E，Murray W，Saunders M A.User's Guide for SNOPT Version 7:Software for Large-Scale Nonlinear Programming［R］.Department of Mathematics University of California，San Diego.2010.06.16.

［12］ Zondervan K P，Bauer T P，Betts J T，Huff-man W P.Solving the Optimal Control Problem Using a Nonlinear Programming Technique Part 3:Optimal Shuttle Reentry Trajectories［C］.AIAA-84-2039.

［13］ 楊炳尉.標(biāo)準(zhǔn)大氣參數(shù)的公式表示［J］.宇航學(xué)報(bào)，1983，(1):83-86.(Yang Bingwei.Formulization of standard atmospheric parameters［J］.Journal of Astronautics，1983，(1):83-86.)