基于子集模擬優化的空天飛機再入軌跡混合優化方法

朱俊杰 余雄慶

1.飛行器先進設計技術國防重點學科實驗室，南京210016

2.南京航空航天大學，南京210016

空天飛機軌跡優化是在飛行動力學和物理學約束的條件下，尋求某種意義下最優解的過程，它貫穿于整個飛行器設計過程中，影響著總體、氣動布局、制導控制、動力和結構等多個分系統的設計。

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

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

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

1 再入軌跡優化問題的描述

1.1 運動方程及設計變量

飛行器再入大氣層時的三自由度運動學和動力學方程組如下［7］:

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

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

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

1.2 歸一化處理

對于再入軌跡優化問題，不同變量之間的數值量級相差較大，例如R是106級別，而經度和緯度化為弧度后是1～10級別，兩者數值上相差太大，會導致數值奇異，優化計算失敗。為了增加優化計算的穩健性，需要對原始模型歸一化處理。令:

1.3 優化目標

優化目標的選擇主要由空天飛機的設計及任務要求確定，通常有如下4種優化目標:1)以飛行航程最大作為性能指標;2)以再入過程中總吸熱量最小作為性能指標［9］;3)以最小能量消耗作為性能指標［10］;4)以到達指定目標點飛行時間最短作為性能指標。

1.4 約束條件

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

2 SSO-SQP混合優化方法

SSO-SQP混合優化方法的思路是:首先在整個設計空間使用SSO算法進行全局搜索，迭代適當步數。待SSO迭代終止后，將SSO計算結果作為SQP優化算法的初始點，使用SNOPT軟件包進行再次優化，得出收斂后的結果。以下對SSO和SQP算法進行簡要介紹。

2.1 SSO算法

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

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

為了在可靠性分析框架中研究優化問題，人為設定設計變量為隨機變量，這個“隨機化”處理使得目標函數f也變為一個隨機變量，且擁有自己的概率密度函數和累積分布函數。根據累積分布函數的數學定義可知，累積分布函數曲線是單調非減，并且右連續，有且僅有唯一的最小值fopt。構造一個可靠性問題如下:

其中，失效事件為F={f(x)≤fopt}。顯然，與方程(14)有關的失效概率是0，因為fopt是全局最小。然而所關心的不是這個失效概率，而是目標函數最小值點或者包含目標函數最小值的區域。失效事件的變量區域與函數取得最小值的變量區域是對應的，從而實現了最優化問題向可靠性問題的轉化，進而可用子集模擬方法解決優化問題［6］。

SSO算法的計算步驟［5］如下:

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

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

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

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

2.2 SQP算法

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

3 算例

空天飛機再入軌跡優化算例取自文獻［12］，并在該算例基礎上增加了熱流密度約束、動壓約束和過載約束。

飛行器的氣動模型如下:

空天飛機參考面積S=250.237m2，質量m=92162kg，前緣半徑RN=1m。不同高度下大氣參數參見文獻［13］。

本例的設計變量包括狀態變量和控制變量。x=［R，Φ ，θ，v，γ，ψ］為狀態變量，u=［α，β］為控制變量。該優化問題的初始狀態是:

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

為了驗證SSO-SQP算法的魯棒性，計算了不同配點數情況下的軌跡優化問題。計算機配置為3.10GHz，內存2GB。SSO迭代100步，7s內就能計算完畢。最終結果見表1。

表1 不同配點數的SSO-SQP優化結果

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

當配點數為140時，迎角和傾側角隨時間變化曲線如圖1和2所示。熱流密度、動壓和過載曲線如圖3～5所示，可以看出優化結果能很好地滿足約束要求。

圖1 迎角隨時間的變化

圖2 傾側角隨時間的變化

圖3 熱流密度隨時間的變化

圖4 動壓隨時間的變化

圖5 過載隨時間的變化

4 結論

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

1)該算法能有效地求解空天飛機軌跡優化問題，能根據精度和計算時間要求，設定配點數;

2)SSO-SQP算法魯棒性較強，當配點從20增加到180時，均能優化出正確結果;

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

［1］ 張鼎逆，劉毅.基于改進遺傳算法和序列二次規劃的再入軌跡優化［J］.浙江大學學報(工學版)，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］ 陳功，傅瑜，郭繼峰.飛行器軌跡優化方法綜述［J］.飛行力學，2011，29(4):1-5.(Chen Gong，Fu Yu，Guo JiFeng.Survey of aircraft trajectory optimization methods［J］.Flight Dynamics，2011，29(4):1-5.)

［3］ 楊希祥，李曉斌，肖飛，等.智能優化算法及其在飛行器優化設計領域的應用綜述［J］.宇航學報，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］ 李紅雙，馬遠卓.結構可靠性分析與隨機優化設計的統一方法［M］.北京:國防工業出版社，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］ 雍恩米，唐國金，陳磊.基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優化［J］.宇航學報，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］ 湯亮，楊建民，陳風雨.多約束條件下的升力滑翔式再入軌跡優化［J］.導彈與航天運載技術，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］ 陳剛，萬自明，徐敏，等.飛行器軌跡優化應用遺傳算法的參數化與約束處理方法研究［J］.系統仿真學報，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］ 楊炳尉.標準大氣參數的公式表示［J］.宇航學報，1983，(1):83-86.(Yang Bingwei.Formulization of standard atmospheric parameters［J］.Journal of Astronautics，1983，(1):83-86.)