張蕊 孫國慶 楊澤山

（北京航空航天大學(xué)宇航學(xué)院，北京 100191）

1 引言

臨近空間飛行器再入時，由于需要較長一段時間在大氣層中飛行，將受到熱流密度、動壓、過載的嚴(yán)格約束，以及需要進(jìn)行主動姿態(tài)控制，給飛行器材料和內(nèi)部設(shè)施熱絕緣帶來巨大的挑戰(zhàn)，所以安全可靠的再入成為重要問題，再入軌跡的優(yōu)化也成為一項關(guān)鍵技術(shù)。

現(xiàn)有的軌跡設(shè)計方法有工程規(guī)劃法，依據(jù)工程經(jīng)驗將軌跡分為數(shù)段，然后分段求解［1］。軌跡優(yōu)化設(shè)計主要有基于極大值原理的間接法［2-3］和基于非線性規(guī)劃理論的直接法，直接法中的偽譜法［4-5］應(yīng)用最為廣泛，但無論是直接法還是間接法，均采用傳統(tǒng)的尋優(yōu)方法，根據(jù)工程經(jīng)驗設(shè)置參數(shù)初值，然后通過試湊法對這些參數(shù)進(jìn)行調(diào)整。由于目標(biāo)函數(shù)對待優(yōu)化參數(shù)的微小變化十分敏感，這樣在全局尋優(yōu)時給軌跡優(yōu)化帶來很大問題，很難對整個的解空間進(jìn)行搜索并找到全局最優(yōu)解，即使能找到全局最優(yōu)解，也需花費大量的時間尋優(yōu)，而且優(yōu)化結(jié)果的好壞，在很大程度上取決于某些參數(shù)的初始給定值，因而不便于工程應(yīng)用。除了傳統(tǒng)的直接法和間接法之外，還有許多其他優(yōu)化方法應(yīng)用也較廣泛，如動態(tài)規(guī)劃方法［6］、快速搜索 隨機樹 法［7］、遺傳 算法［8］等，這些算法在OPTIMUS軟件平臺中都有體現(xiàn)。

OPTIMUS軟件平臺是一種過程集成與多學(xué)科優(yōu)化的軟件平臺，主要功能包括過程自動化、設(shè)計空間開發(fā)、先進(jìn)的概率和統(tǒng)計方法、數(shù)值優(yōu)化算法等。OPTIMUS已經(jīng)被廣泛應(yīng)用于航空航天領(lǐng)域的結(jié)構(gòu)優(yōu)化、渦輪氣動優(yōu)化設(shè)計［9］、彈道仿真控制系統(tǒng)參數(shù)優(yōu)化［10］等方面。在再入軌跡優(yōu)化問題中，基于OPTIMUS平臺，可以充分應(yīng)用其集成的數(shù)值優(yōu)化算法，快速確定初始猜想值，并選取最佳優(yōu)化方案進(jìn)行求解。

本文針對再入軌跡優(yōu)化中初值選取和計算效率的問題，利用OPTIMUS 尋優(yōu)平臺，找出軌跡優(yōu)化的初值點，并確定最佳的優(yōu)化算法，再依據(jù)算法編寫優(yōu)化程序，這種計算機集成的輔助優(yōu)化設(shè)計可以取代大量的低效手動設(shè)計。仿真結(jié)果表明，基于OPTIMUS的軌跡優(yōu)化方案計算快速、尋優(yōu)準(zhǔn)確，且計算結(jié)果精度高。

2 再入軌跡優(yōu)化問題

飛行器再入過程中，由于各種狀態(tài)變量的絕對值相差較大，這對于再入運動方程的數(shù)值計算非常不利，為提高計算精度、增強算法收斂性與快速性，需將實際再入運動方程進(jìn)行歸一化處理。因此，在考慮地球旋轉(zhuǎn)的情況下，飛行器再入無量綱運動方程組為

上述再入運動方程組是關(guān)于無量綱時間τ＝t／的微分方程組，其中，R0為地球半徑，g0是海平面重力加速度。式（1）～（6）方程組中：r為飛行器質(zhì)心-地心的距離與地球半徑R0的比值；θ和φ分別為再入過程中的經(jīng)度緯度；V為飛行器相對地球的速度與歸一化變量（圓周速度）的比值；γ為航跡傾角；ψ為航向角；Ω為地球自轉(zhuǎn)角速度，σ為滾轉(zhuǎn)角；L為飛行器的總升力加速度；D為阻力加速度。L、D計算公式為

式中：ρ為大氣密度；m為飛行器質(zhì)量；Sref為飛行器的參考面積；CL和CD分別為飛行器的升力系數(shù)和阻力系數(shù)，其與飛行器的馬赫數(shù)Ma和攻角α有關(guān)，也即CL＝CL（Ma，α），CD＝CD（Ma，α）。

再入軌跡優(yōu)化問題，即尋找最優(yōu)控制律σ（t），使飛行器在再入過程中，滿足要求的初始狀態(tài)、目標(biāo)狀態(tài)和運動學(xué)方程，并使給定的性能泛函具有極值。本文中設(shè)計的性能泛函為在給定縱程下使再入橫程最大。

為簡化問題，再入軌跡優(yōu)化時采用以下假設(shè)條件：重力加速度g為常值g0；滿足準(zhǔn)平衡滑翔條件γ≈0，γ≈0；由于高度相對地球半徑很小，令r為常值且r＝1；不考慮地球自轉(zhuǎn)，將哥氏加速度與牽連加速度影響因素忽略。簡化后的再入運動方程為

過程約束考慮氣動加熱、過載及動壓3種不等式約束。熱流率約束為

動壓約束q主要取決于機體結(jié)構(gòu)和運載器表面防熱材料的強度，qmax為允許的最大動壓值。

過載約束n主要取決于飛行器結(jié)構(gòu)強度和器載設(shè)備的過載范圍，nmax為允許的最大過載值。

給定縱程下最大橫程問題的描述如下（見圖1）：設(shè)飛行器的初始再入點為O點，其對應(yīng)的經(jīng)緯度為θO和φO；要求到達(dá)的目標(biāo)點為T，其對應(yīng)的經(jīng)緯度為θT和φT；P點為經(jīng)過O點和T點的大圓弧上的點，且滿足O點到P點的距離為給定的縱程值SP，其對應(yīng)的經(jīng)緯度為θP和φP；F點即為給定縱程值SP所對應(yīng)的最大橫程點，其對應(yīng)的經(jīng)緯度為θF和φF。再入軌跡優(yōu)化問題即為尋求最優(yōu)滾轉(zhuǎn)角控制量σ（t），使飛行器在滿足初始狀態(tài)（O點）和給定的縱程SP的情況下，到達(dá)最大橫程點F處，且滿足再入過程約束和再入運動方程。

圖1 最大橫程示意圖Fig．1 Maximum crossrange diagram

本文主要是基于最優(yōu)控制理論，設(shè)計滾轉(zhuǎn)角控制量，通過選擇合適的求解方法求解最優(yōu)控制問題，使再入飛行器在給定縱程SP下，獲得最大橫程SF。優(yōu)化目標(biāo)是橫程SF最大，優(yōu)化變量為滾轉(zhuǎn)角σ（｜σ｜＜π／2），性能泛函即為

終點約束條件主要有2個：①達(dá)到預(yù)設(shè)的末端能量；②保證連接∠OPF為直角。

式中：rF為終端高度；VF為終端速度；eF為終端能量，哈密爾頓函數(shù)H寫為

式中：pθ，pφ，pV，pψ為共軛狀態(tài)變量，共軛狀態(tài)方程為

根據(jù)最優(yōu)條件，最優(yōu)滾轉(zhuǎn)角控制律要滿足

根據(jù)運動的積分以及整理，具體參見文獻(xiàn)［11］，可得到滾轉(zhuǎn)角σ最優(yōu)控制律為

在滾轉(zhuǎn)角表達(dá)式（22）中，可以看出分子分母是常數(shù)c1，c2和c3的線性表達(dá)式，由此可以得出，滾轉(zhuǎn)角的值只與其中2個獨立變量有關(guān)，最終優(yōu)化目標(biāo)也即通過搜索找到2個獨立變量的值。假設(shè)：c2＝sinλ，c3＝cosλ。終點約束條件2可以改寫為

至此，給定縱程下最大橫程優(yōu)化問題可表述為：在滿足終點約束條件的前提下，利用尋優(yōu)算法求解最優(yōu)控制參數(shù)c1和λ的值。c1作用域為［-1，1］，λ作用域為［-π，π］。對實際飛行器軌跡優(yōu)化而言，可以取

其中K1，K2為權(quán)值，可取為K1＝1，K2＝1。需要尋找兩個參數(shù)c1和λ，使Δ（c1，λ）最小。利用OPTIMUS集成的優(yōu)化平臺，可以方便快速地搜索出最優(yōu)參數(shù)c1和λ。

3 基于OPTIMUS軟件平臺的再入軌跡優(yōu)化

在OPTIMUS軟件平臺中，集成了多種優(yōu)化算法，分為兩大類：全局優(yōu)化算法和局部優(yōu)化算法。全局優(yōu)化算法主要包括：差分進(jìn)化、自適應(yīng)進(jìn)化、模擬退火、遺傳算法［12］等；局部優(yōu)化算法主要包括：模式搜索法［13］、序列二次規(guī)劃、非線性二次規(guī)劃、廣義簡約梯度等。全局優(yōu)化算法往往收斂速度較慢，而局部優(yōu)化算法收斂速度較快，但如果初始點選取不夠準(zhǔn)確，局部優(yōu)化算法將有可能使搜索進(jìn)入局部極值，導(dǎo)致結(jié)果不準(zhǔn)確。本文結(jié)合全局優(yōu)化算法和局部優(yōu)化算法各自的優(yōu)缺點，提出了基于全局優(yōu)化算法和局部優(yōu)化算法相結(jié)合的混合優(yōu)化的概念，全局優(yōu)化算法采用遺傳算法，局部優(yōu)化算法采用模式搜索法，將遺傳算法得到的結(jié)果作為模式搜索法的初始點，并進(jìn)一步精確搜索得到最優(yōu)解。采用混合優(yōu)化算法，不僅縮短了仿真時間、降低了錯誤率，而且提高了優(yōu)化精度，滿足了全局最優(yōu)性。下面主要介紹在OPTIMUS軟件平臺下優(yōu)化設(shè)計的詳細(xì)步驟。

3．1 建立工作流

為了應(yīng)用OPTIMUS軟件平臺，通過可視化的建模方式，將仿真程序和工作流集成在OPTIMUS平臺的框架中，首先建立M 文件global．variable．dat，其中包含有待優(yōu)化的控制參數(shù)c1和λ的值。接著將事先編好的數(shù)學(xué)模型仿真程序?qū)隣PTIMUS平臺，搭建系統(tǒng)工作流。在OPTIMUS 中選擇控件，并將所有控件用連接線進(jìn)行連接，搭建OPTIMUS工作流。工作流控件中包含設(shè)計變量、輸入輸出文本文件、仿真程序模塊、輸出向量、輸出變量以及連接線。搭建工作流的目的，在于以圖像化的方式描述仿真中計算的邏輯和數(shù)據(jù)流程，仿真工作流在輸入文件中更新輸入?yún)?shù)、進(jìn)行仿真計算，從結(jié)果文件中讀入計算結(jié)果。

在手動輸入完成后，OPTIMUS 就開始自動化的計算。先對c1和λ取值，修改global．variable．dat中相應(yīng)的數(shù)值。之后運行M 文件搭建的動力學(xué)仿真模塊，待程序運行結(jié)束后讀取Output．dat中的數(shù)據(jù)，根據(jù)事先設(shè)置的算法，自動計算出c1和λ更新值，回到第一步進(jìn)行數(shù)據(jù)刷新，然后重復(fù)這個工作流程，直至得到符合條件的解。整個工作不需要人工干預(yù)，由計算機獨立完成長時間的運算，提高了工作效率。

3．2 試驗設(shè)計與響應(yīng)面模型的搭建

完成建模后要進(jìn)行試驗設(shè)計，通過較少的試驗次數(shù)，比較準(zhǔn)確地反映輸入?yún)?shù)和輸出參數(shù)的關(guān)系，并且通過響應(yīng)面模型使結(jié)果可視化。OPTIMUS提供了多種試驗設(shè)計（Design of Experiment，DOE）方法，包括全參數(shù)、隨機試驗設(shè)計、部分參數(shù)等。試驗設(shè)計的目的在于分析參數(shù)相關(guān)性，為響應(yīng)面的建立提供樣本點，保證在設(shè)計變量變化時（特別是邊界上）不發(fā)生干涉、參數(shù)不匹配及命令不執(zhí)行等情況。此外還應(yīng)對網(wǎng)格尺度、計算步長、差分格式及數(shù)等進(jìn)行匹配調(diào)試，以便在計算時間、精度及穩(wěn)定性等方面得到平衡。

優(yōu)化變量c1可行域為［-1，1］，λ可行域為［-π，π］，圖2為拉丁超立方樣本點分布，圖3為參數(shù)相關(guān)性。

圖2 拉丁超立方樣本點分布Fig．2 Latin super cubic sample point distribution

圖3 拉丁超立方參數(shù)相關(guān)性Fig．3 Latin super cubic parameters correlation

通過實驗設(shè)計分析得出：該飛行器軌跡優(yōu)化的目標(biāo)泛函對選擇的參數(shù)敏感性非常高，參數(shù)響應(yīng)面存在很多局部的急劇變化，因此選擇一種合適的優(yōu)化算法顯得尤其重要。試驗點的選擇應(yīng)針對不同模型有側(cè)重地進(jìn)行考慮，雖然盡可能地將試驗點個數(shù)減少，但前提是體現(xiàn)整個參數(shù)可行域內(nèi)的分布情況。

選擇相應(yīng)的設(shè)計實驗方法可以得到響應(yīng)面模型，來模擬一系列輸入與輸出參數(shù)之間的響應(yīng)關(guān)系。由于該飛行器優(yōu)化模型非線性極強，對于參數(shù)的響應(yīng)面波動十分劇烈，所以需要相當(dāng)精確的擬合方法對其進(jìn)行擬合，以減少后續(xù)分析處理中由于模型誤差引入的不確定性。以插值法中的RBF 神經(jīng)網(wǎng)絡(luò)方法為例，形成的插值響應(yīng)面如圖4所示。

圖4 RBF神經(jīng)網(wǎng)絡(luò)方法插值響應(yīng)面Fig．4 RBF neural network method interpolation response surface

響應(yīng)面模型搭建的出發(fā)點在于利用一系列的擬合算法，在若干試驗點的基礎(chǔ)上，對待優(yōu)化的精確響應(yīng)面模型進(jìn)行離散樣本點擬合，用數(shù)理分析統(tǒng)計的方法代替實際飛行器模型。雖然舍棄了小部分精度，但換來的是后續(xù)的高效運算，極大地減少了數(shù)據(jù)存儲和運算量，基于該擬合模型，就可以將輸入?yún)?shù)和輸出參數(shù)關(guān)聯(lián)起來，從而具有較快的響應(yīng)速度和較好的響應(yīng)精度，降低了運算成本和響應(yīng)時間。

完成了上述的工作流的建立、以及試驗設(shè)計和響應(yīng)面模型的搭建之后，就可以選擇OPTIMUS軟件平臺內(nèi)部自帶的遺傳算法和模式搜索法對參數(shù)c1和λ進(jìn)行尋優(yōu)計算，并最終得到參數(shù)的最優(yōu)值為c1＝0．784 7，λ＝0．524 6。

4 編程實現(xiàn)軌跡優(yōu)化與仿真結(jié)果

將OPTIMUS軟件平臺優(yōu)化出的參數(shù)c1和λ的值，代入到公式（22）的滾轉(zhuǎn)角最優(yōu)控制律中，并編寫相應(yīng)的MATLAB 再入軌跡仿真程序，即得到滿足性能泛函和過程約束的再入軌跡。

MATLAB 仿真用到的再入飛行器參數(shù)為：參考面積149．388m2，飛行器質(zhì)量37 362．9kg，最大熱流率794 005 W／m2，最大動壓15 425N／m2，最大過載5gn。攻角α變化范圍為［0°，45°］；馬赫數(shù)Ma變化范圍為［3，25］；高度變化范圍為［0km，120km］。初始條件為：高度123．104km，經(jīng)度和緯度分別為-122．498°和-33．262 9°，速度為7 625．15m／s，航跡傾角為-1．249 4°，航跡偏角為46．061 7°。終端條件為：高度30．427km，經(jīng)度和緯度分別-62．346 5°，35．270 6°，速度為908．151 6m／s，縱程為9 871．4km。圖5～8是仿真結(jié)果，對任意的給定縱程，都有左右兩邊的兩個最大橫程。圖5是高度-速度空間中的再入軌跡，從圖中可以看到滿足熱流率、動壓和過載約束；圖6是對應(yīng)的最優(yōu)滾轉(zhuǎn)角；圖7是航跡傾角；圖8是星下點軌跡（給定縱程下對應(yīng)的最大橫程有左右2條）。

圖5 再入軌跡圖Fig．5 Entry trajectory

圖6 滾轉(zhuǎn)角控制曲線Fig．6 Bank angle profile

圖7 航跡傾角變化曲線Fig．7 Flight path angle profile

圖8 星下點軌跡Fig．8 Trajectory of ground track

5 結(jié)束語

本文基于OPTIMUS優(yōu)化軟件平臺，提出了臨近空間飛行器的再入軌跡優(yōu)化方案。與傳統(tǒng)的優(yōu)化方法不同，先利用OPTIMUS平臺優(yōu)化求出初值，并選出最佳優(yōu)化算法，再編寫優(yōu)化程序進(jìn)行求解。OPTIMUS平臺在軌跡優(yōu)化中的應(yīng)用，可以快速地搭建起整個數(shù)據(jù)流程，通過應(yīng)用多種不同的優(yōu)化算法，對同一個問題進(jìn)行求解分析，可以更加全面地了解軌跡優(yōu)化問題各個環(huán)節(jié)的細(xì)節(jié)，并獲得優(yōu)化問題的初始猜想，然后針對具體問題設(shè)計優(yōu)化方案，所開發(fā)的算法程序具有計算高效、便捷準(zhǔn)確的特點。

（References）

［1］趙漢元，飛行器再入動力學(xué)和制導(dǎo)［M］．長沙：國防科學(xué)技術(shù)大學(xué)出版社，1997

Zhao Hanyuan．Flight vehicle entry dynamics and guidance［M］．Changsha：University of National Defense Science and Technology Press，1997（in Chinese）

［2］Vinh N X．Optimal trajectories in atmospheric flight［M］．NewYork：Elsevier Scientific Publishing Company，1981

［3］Lu P，Xue S B．Rapid generation of accurate entry landing footprints［J］．Journal of Guidance，Control and Dynamics，2010，33（3）：756-767

［4］雍恩米，唐國金，陳磊．基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化［J］．宇航學(xué)報，2008，29（6）：1766-1772

Yong Eenmi，Tang Guojin，Chen Lei．Rapid trajectory optimization for hypersonic reentry vehicle via Guass pseudospectral method［J］．Journal of Astronautics，2008，29（6）：1766-1772（in Chinese）

［5］宗群，田栢苓，竇立謙．基于Gauss偽譜法的臨近空間飛行器上升段軌跡優(yōu)化［J］．宇航學(xué)報，2010，31（7）：1775-1781

Zong Qun，Tian Bailing，Dou Liqian．Ascent phase trajectory optimization for near space vehicle based on Gauss pseudospectral method［J］．Journal of Astronautics，2010，31（7）：1775-1781（in Chinese）

［6］John T B．Survey of numerical methods for trajectory optimization［J］．Journal of Guidance，Control and Dynamics，1998，21（3）：193-207

［7］周克強，高曉光，白奕．改進(jìn)的快速擴(kuò)展隨機樹在航跡規(guī)劃中的應(yīng)用［J］．系統(tǒng)工程與電子技術(shù)，2006，28（10）：1538-1540

Zhou Keqiang，Gao Xiaoguang，Bai Yi．Trajectory planning using improved rapidly exploring random tree［J］．System Engineering and Electronics，2006，28（10）：1538-1540（in Chinese）

［8］周浩，陳萬春．高超聲速飛行器滑行航跡優(yōu)化［J］．北京航空航天大學(xué)學(xué)報，2006，32（5）：513-517

Zhou Hao，Chen Wanchun．Optimization of glide trajectory for a hypersonic vehicle［J］．Journal of Beijing University of Aeronautics and Astronautics，2006，32（5）：513-517（in Chinese）

［9］嚴(yán)俊峰，吳寶元，逯婉若．基于OPTIMUS的渦輪氣動優(yōu)化設(shè)計［J］．火箭推進(jìn)，2008，34（2）：13-17

Yan Junfeng，Wu Baoyuan，Lu Wanruo．Turbo aerodynamics optimization design based on OPTIMUS flat［J］．Journal of Rocket Propulsion，2008，34（2）：13-17（in Chinese）

［10］周健，姜曉菡，周喻虹．基于OPTIMUS的控制系統(tǒng)參數(shù)優(yōu)化計算［J］．彈箭與制導(dǎo)學(xué)報，2008，28（3）：13-18

Zhou Jian，Jiang Xiaohan，Zhou Yuhong．Optimization of parameters for control system based on OPTIMUS［J］．Journal of Projectiles，Rockets，Missiles and Guidance，2008，28（3）：13-18（in Chinese）

［11］Lu P．Entry trajectory optimization with analytical feedback bank angle law，AIAA-2008-7268［R］．Washington：AIAA，2008

［12］袁麟博，章衛(wèi)國，李廣文．一種基于遺傳算法-模式搜索法的無人機路徑規(guī)劃［J］．彈箭與制導(dǎo)學(xué)報，2009，29（3）：279-282

Yuan Linbo，Zhang Weiguo，Li Guangwen．A path planning for UAV based on genetic-pattern searching algorithm［J］．Journal of Projectiles，Rockets，Missiles and Guidance，2009，29（3）：279-282（in Chinese）

［13］張煜東，吳樂南，王水花．基于遺傳算法與模式搜索的混合優(yōu)化算法［J］．南京信息工程大學(xué)學(xué)報，2012，4（1）：34-39

Zhang Yudong，Wu Lenan，Wang Shuihua．A hybrid optimization method based on genetic algorithm and pattern search［J］．Journal of Nanjing University of Information Science ＆ Technology，2012，4（1）：34-39（in Chinese）