基于改進hp自適應偽譜法的高超聲速飛行器上升段軌跡規劃

劉瑞帆, 于云峰, 閆斌斌

(西北工業大學 航天學院, 陜西 西安 710072)

?

基于改進hp自適應偽譜法的高超聲速飛行器上升段軌跡規劃

劉瑞帆, 于云峰, 閆斌斌

(西北工業大學 航天學院, 陜西 西安 710072)

針對高超聲速飛行器的上升段軌跡優化問題,提出了一種改進的hp自適應Radau偽譜法。該算法將軌跡劃分為多個子區間,采用雙層優化策略迭代調整子區間個數和區間內的配點數。首先,以微分-代數約束在采樣點處的相對誤差作為解的誤差評估準則;在需要提高精度的區間,將相對誤差作為迭代判據,若通過相對誤差求得的多項式階數小于允許的最大階數,增加區間內的配點數(p法);否則,細化時間區間(h法)來提高求解精度。仿真結果表明,相較于h方法和p方法,文中的算法在求解最優問題時配點總數更少,求解效率更高。

上升段;軌跡優化;hp自適應;偽譜法；高超聲速飛行器

飛行器軌跡優化是指在特定的約束條件下,尋找飛行器從初始點到目標點滿足某種性能指標最優的運動軌跡。在數學本質上,飛行器軌跡優化可以抽象為在包含微分方程、代數方程和不等式等約束下求解泛函極值的開環最優控制問題。一般將軌跡優化方法分為直接法和間接法,間接法基于Pontryagin極大值原理將最優控制問題轉換為Hamiltonian邊值問題(HBVP),其解的精度高,滿足一階最優必要條件,但HBVP問題的收斂域小,需要估計協態變量的初值,且難以求解存在路徑約束的最優控制問題[1-3];直接法無需求解最優必要條件,而是將連續的最優控制問題離散并參數化,用數值方法對性能指標尋優,其解的收斂半徑大,近年來得到廣泛應用[4-7]。作為直接法中的一個分支,全局偽譜法發展極為迅速,宗群等[7]將Gauss偽譜法與序列二次規劃算法相結合,對存在邊值及加速度約束的軌跡優化問題進行了求解。相比于傳統的間接法,全局偽譜法收斂速度快,然而在求解非光滑問題時,其收斂速度變慢且無法捕捉狀態突變點導致求解精度下降,運用效果不理想[8]。Jiang Zhao等[9]通過分段優化以改善軌跡的非光滑特性,然而在解出最優軌跡之前難以確定要分多少段以及在哪分段。趙吉松等[10-11]基于局部配點法提出自適應網格細化算法,在非光滑區域加密網格以提高求解精度,能夠有效控制問題規模,但失去了全局快速收斂的優勢。Darby等[12-13]首次提出hp自適應偽譜法并將其應用到最優控制問題上,該方法汲取了全局偽譜法和網格細分法的優點,采用雙層優化策略,具有更合理的配點分布,從而提高了算法的優化效率。文獻[14-15]將hp自適應偽譜法分別應用于滑翔彈道快速優化問題和再入軌跡優化問題。

本文基于Radau偽譜法提出了一種改進的hp自適應偽譜法,并將該算法應用于高超聲速飛行器的最省燃料爬升問題。結果表明，與目前已有的偽譜法相比，改進的hp自適應偽譜法能夠以較少的時間代價獲得較高精度的解。

1 上升段軌跡優化數學模型

高超聲速飛行器上升段軌跡優化問題的數學模型包括動力學模型、推力模型和氣動力模型,以及飛行器爬升過程中必須滿足的路徑約束條件。

1.1 動力學模型

僅考慮飛行器在縱向平面內的質心運動,忽略地球曲率和自轉的影響,得到高超聲速飛行器上升段的動力學方程如下

(1)

式中,[V,θ,x,y,m]T為系統的狀態變量,分別表示飛行器的速度、航跡傾角、航程、飛行高度和質量;P為發動機推力,D為阻力,L為升力,α為攻角,Isp為燃料比沖,μ為地球引力參數,Re為地球半徑。

1.2 推力模型

本文采用渦輪沖壓發動機與超燃沖壓發動機相結合的推力模型。當馬赫數小于4時,渦輪沖壓發動機提供推力;當馬赫數大于4時,超燃沖壓發動機提供推力。推力系數和特征沖量的計算公式如下[16]:

1) 推力系數Cp為

(2)

推力P為

(3)

2) 特征沖量Isp為

(4)

式中,h為飛行高度,q為動壓,Ae為進氣口面積,M為飛行馬赫數。

1.3 氣動力模型

飛行器飛行過程中受到升力L和阻力D的作用,計算公式如下

(5)

式中,Sref為飛行器參考面積,ρ為大氣密度,計算公式為ρ=ρ0e-βh,其中ρ0為海平面大氣密度,β為密度系數常數,h為飛行器高度。升力系數CL和阻力系數CD由攻角α和馬赫數Ma表示

(6)

式中的系數CLi和CDi基于通用航空飛行器數據[17]。

1.4 約束條件

1.4.1 熱流率約束

爬升過程中,為了滿足嚴苛的熱防護條件,必須對熱流率Qd添加如下約束

(7)

式中,KQ為熱流率標準化常數。

1.4.2 動壓約束

動壓的上限由飛行器執行機構的鉸鏈力矩確定,下限值則是由超燃沖壓發動機啟動條件所限制。動壓q的路徑約束可表示為

(8)

1.4.3 過載約束

為保證飛行器的結構不被破壞,需對過載進行約束,即

(9)

式中,g0為海平面重力加速度。

1.4.4 控制量和狀態量約束

受機動性能等因素的影響,飛行器在上升段飛行過程中,要求滿足一定的控制約束和狀態約束。因此建立如下約束條件

(10)

1.4.5 終端約束

飛行器巡航段所需要的初始高度、速度和彈道傾角確定了其爬升段的終端約束,即

(11)

1.5 目標函數

飛行器爬升過程中節省的燃料作為巡航段的燃油增加,可以大大增加巡航段的飛行航程。據此,選取上升段燃料最省為優化目標,即

(12)

1.6 多區間軌跡優化問題描述

基于高超聲速飛行器的動力學模型,最省燃料軌跡優化問題的本質為尋找最優控制量α使得目標函數J最小,同時滿足多種約束條件((7)～(11)式)。為表述方便,將該問題描述為一般形式的Bolza型最優控制問題:定義狀態變量x=[h,V,θ,m]T、控制變量u=α,使代價函數最小化

(13)

約束條件為

(14)

多區間最優控制問題是指將上述Bolza型問題的時間區間[t0,tf]劃分為K個子區間,用t0

(15)

將時間區間t∈[tk-1,tk]變換為τ=[-1,1]。

2 Radau偽譜法

在RPM中,連續Bolza問題的狀態變量近似為

(16)

(17)

(18)

(19)

(20)

3 改進的hp自適應偽譜法

3.1 誤差評估準則

在數值求解過程中,若使數值方法的解能夠精確地近似原問題的確切解,則需要在任意點處滿足微分-代數約束方程。因此,以微分-代數約束方程在配點間的滿足程度作為解的近似誤差評估準則。

(21)

(22)

(23)

(24)

3.2 多項式階數的估計

(25)

為保證求得的配點數為整數,對(25)式右邊項向上取整,即

(26)

3.3 hp自適應迭代策略

(27)

每個區間內的配點數均設為配點數初始值Nmin。

可以看出,本文所提出的hp迭代策略優先使用p法來提高精度,只有當配點數超過所設定的上限值時才進行區間劃分,劃分后每個區間內繼續采用p法,以此規律不斷迭代直到滿足解的精度要求。

基于以上的迭代策略,hp自適應算法流程總結如下:

①初始化。令迭代次數M=0,設定區間的初始劃分以及區間內配點數上下限Nmin,Nmax;

②將最優控制問題離散化,轉化為NLP問題進行求解,并按照3.1節內容計算解的相對最大誤差emax;

③若所有區間的emax均小于ε或M>Mmax,停止迭代;否則,繼續④;

④在emax>ε的區間內采用3.3節的自適應策略更新網格,直到所有的網格完成更新;

⑤令M=M+1,轉向步驟②進行下一次迭代。

4 仿真分析

本文以某高超聲速飛行器為例,飛行器初始質量為3 600 kg,氣動參考面積為0.6 m2,給定初始高度y0=17 km,初始速度v0=1 100 m/s,初始彈道傾角θ0=4°,終端高度yf=27 km,終端速度vf=1 870 m/s,終端航跡傾角θf=0°。爬升過程中狀態變量、控制變量和路徑約束的邊界限制如下

表1 約束邊界設置

優化計算過程中軟件環境為:Windows7 32位操作系統和MATLAB R2013a;硬件環境為:Intel Core i5-2450處理器,4G內存。在精度要求為10-5時,求得的爬升段最小燃料消耗為165.56 kg,整個計算過程用時7.069 8 s,優化結果如圖1～圖4中a)圖的圓點所示。

為驗證解的最優性,將最優控制量通過三次樣條插值代入到動力學方程中,求得飛行器的實際飛行狀態,計算結果圖1～圖4中a)圖的實線所示。離散點處的狀態誤差如圖中星型線所示。圖中可以看出兩者的相對誤差在10-4～10-6之間,解的最優性得到驗證。

圖1 高度及高度誤差曲線 圖2 速度及速度誤差曲線

圖3 航跡傾角及航跡角誤差曲線 圖4 質量及質量誤差曲線

在精度要求為ε的情況下,采用3種不同的偽譜方法對同一航跡優化問題進行求解,優化結果如表2所示。

表2 不同偽譜法的優化比較

表中,p代表全局Radau偽譜法,h代表分段定階Radau偽譜法,hp代表改進的hp自適應Radau偽譜法。分別表示初始子區間個數和區間內的初始配點數。K0、N0分別表示優化完成時的最大微分-代數約束的誤差、計算時間區間個數、整個時間區間的配點總數、迭代次數以及目標函數。需要注意的是在采用全局偽譜法優化時,由于求解問題的不光滑性,精度要求較高時,過多的配點導致優化時間過長,無法得到優化結果。

分析表中數據,我們可以得出如下結論:①精度要求較低時,3種方法均可以優化出結果,且在優化速度和配點總數上無明顯差異;②的精度要求下,hp偽譜法逐漸地顯現出其優勢,不論是優化速度還是配點總數均優于其他2種方法;③隨著精度要求的提高,p方法僅通過增加配點數無法得到最優解,而h方法和hp方法成功地克服了不光滑帶來的問題,得到最優解。兩者相對而言,hp偽譜法在配點總數,子區間劃分和優化速度上均優于h方法。

圖5 不同偽譜法的配點分布

圖5給出了10-3精度下不同偽譜法的配點分布(圖中縱坐標無實際物理意義,僅為圖示方便,直線上點的橫坐標分別表示p偽譜法、h偽譜法和hp偽譜法的配點位置),從圖中可以看出:①p方法為全局偽譜法,配點分布具有兩端密、中間疏的特點;②h方法通過細化區間來提高精度,每個區間內配點相同,圖中可以看出區間的劃分;③hp方法不具備統一的配點分布標準,相比于h方法,區間劃分較少。

本文中所采用的hp自適應策略為:設定區間配點數上限,當配點數未達到上限時,通過增加區間內配點數(p法)提高精度,當配點數達到上限仍未滿足求解精度時,再利用細化時間區間(h法)來提高精度。文獻[12]提出另一種自適應策略:通過比較相對曲率rk與曲率閾值rmax的大小關系來決定采取細分時間區間(h法)還是增加區間內的配點數(p法)來提高算法精度。表3為這兩種細化策略在不同的精度要求下的優化結果比較,其中hp1表示本文細化策略,hp2為文獻[12]中細化策略。

從表3中可以看出,精度要求較低時,hp2方法所需的配點總數小于hp1方法,計算效率也高于hp1;精度要求提高到10-5時,hp2方法逐漸失去了優勢;精度為10-6時,hp1方法在計算效率上明顯優于hp2方法。分析其可能的原因,文獻[12]中hp偽譜法在初始迭代中,NLP問題的解與原問題的解存在較大的誤差,此時根據曲率進行分段可能導致分段過多且分段點不準確,且這種情況隨著精度要求的提高表現得更為突出,而本文中所采用的hp偽譜法優先采用p法來提高精度,在一定程度上可以消除該問題對求解效率的影響。

表3 不同hp自適應策略的優化對比

5 結 論

本文提出了一種基于hp自適應Radau偽譜法的hp自適應網格算法,并將其應用于高超聲速飛行器上升段的軌跡優化問題。仿真結果表明:①與全局偽譜法和固定階偽譜法相比,hp偽譜法充分結合了p法的快速收斂性和p法的計算稀疏性,具有更合理的配點分布,能夠以較少的時間代價獲得更高精度的解;②與以曲率為迭代判據的hp偽譜法相比,本文所提出的偽譜法優先采用p方法,在一定程度上能夠避免分段過多且不準確的問題,更適合求解高精度的最優控制問題,具有一定的工程應用價值。

[1] Betts J T. Survey of Numerical Methods for Trajectory Optimization[J]. Journal of Guidance Control And Dynamics, 1998, 21: 193-207

[2] 雍恩米,陳磊,唐國金.飛行器軌跡優化數值方法綜述[J]. 宇航學報,2008(2):397-406

Yong Enmi, Chen Lei, Tang Guojin. A Survey of Numerical Methods for Trajectory Optimization of Spacecraft[J]. Journal of Astronautics, 2008(2): 397-406 (in Chinese)

[3] 黃長強,國海峰,丁達理. 高超聲速滑翔飛行器軌跡優化與制導綜述[J]. 宇航學報, 2014(4):369-379

Huang Changqiang, Guo Haifeng, Ding Dali. A Survey of Trajectory Optimization and Guidance for Hypersonic Gliding Vehicle[J]. Journal of Astronautics, 2014(4): 369-379 (in Chinese)

[4] Ma L, Shao Z J, Chen W F, et al. Three-Dimensional Trajectory Optimization for Lunar Ascent Using Gauss Pseudospectral Method[R]. AIAA-2016-1372

[5] 白瑞光. 基于Gauss偽譜法的多UAV協同航跡規劃[J]. 宇航學報,2014(9):1022-1029

Bai Ruiguang. Multiple UAV Cooperative Trajectory Planning Based on Gauss Pseudospectral Method[J]. Journal of Astronautics, 2014(9):1022-1029 (in Chinese)

[6] 彭祺擘,李海陽,沈紅新. 基于高斯-偽譜法的月球定點著陸軌道快速優化設計[J]. 宇航學報, 2010, 31(4): 1012-1016

Peng Qibo, Li Haiyang, Shen Hongxin. Rapid Lunar Exact-Landing Trajectory Optimization via Guass Pseudospectral Method[J]. Journal of Astronautics, 2010, 31(4): 1012-1016 (in Chinese)

[7] 宗群,田栢苓,竇立謙. 基于Gauss偽譜法的臨近空間飛行器上升段軌跡優化[J]. 宇航學報,2010(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(7): 1775-1781 (in Chinese)

[8] 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: 432-436

[9] Zhao J, Zhou R. Reentry Trajectory Optimization for Hypersonic Vehicle Satisfying Complex Constraints[J]. Chinese Journal of Aeronautics, 2013, 26: 1544-1553

[10] 趙吉松,谷良賢,佘文學. 配點法和網格細化技術用于非光滑軌跡優化[J]. 宇航學報, 2013(11): 1442-1450

Zhao Jisong, Gu Liangxian, She Wenxue. Application of Local Collocation Method and Mesh Refinement to Non-Smooth Trajectory Optimization[J]. Journal of Astronautics, 2013(11): 1442-1450

[11] 陳琦. 求解非光滑最優控制問題的自適應網格優化[J]. 系統工程與電子技術, 2015,37(6): 1377-1383

Chen Qi. Adaptive Mesh Refinement for Solving Non-Smooth Optimal Control Problem[J]. Systems Engineering and Electronics, 2015, 37(6): 1377-1383 (in Chinese)

[12] Arby C L, Hager W W, Rao A V. An hp-Adaptive Pseudospectral Method for Solving Optimal Control Problems[J]. Optimal Control Applications & Methods, 2011, 32: 476-502

[13] Darby C L, Hager W W, Rao A V. Direct Trajectory Optimization Using a Variable Low-Order Adaptive Pseudospectral Method[J]. Journal of Spacecraft and Rockets, 2011, 48: 433-445

[14] 洪蓓,辛萬青. hp自適應偽譜法在滑翔彈道快速優化中的應用[J]. 計算機測量與控制, 2012(05): 1283-1286.

Hong Bei, Xin Wanqing. Application of hp-Adaptive Pseudospectral Method in Rapid Gliding Trajectory Optimization[J]. Computer Measurement & Control, 2012(05): 1283-1286 (in Chinese)

[15] 王海濤. 基于hp自適應Radau偽譜法的再入飛行器軌跡優化[J]. 科學技術與工程, 2015(02): 165-171

Wang Haitao. Track Optimizing for Reentry Vehicle Based on hp-Adaptive Radau Pseudospectral Method[J]. Science Technology and Engineering, 2015(02): 165-171 (in Chinese)

[16] Morimoto H. Trajectory Optimization for a Hypersonic Vehicle with Constraint[D]. Atlanta, Georgia Institute of Technology,

1997

[17] Phillips T H. A Common Aero Vehicle(CAV) Model, Description, and Employment Guide[R]. Technical Report, Schafer Corporation for AFRL and AF-SPC, 2003

[18] Bengson D. A Gauss Pseudospectral Transcription for Optimal Control[D]. Colorado, Massachusetts Institute of Technology,

2005

[19] Hou H, Hager W W, Rao A V. Convergence of a Gauss Pseudospectral Transcription for Optimal Control[C]∥2012 AIAA Guidance, Navigation, and Control Conference, Minneapolis, MN, 2012

Ascent Phase Trajectory Optimization for Hypersonic Vehicle Based on hp-Adaptive Pseudospectral Method

Liu Ruifan, Yu Yunfeng,Yan Binbin

(School of Astronautics, Northwestern Polytechnical University, Xi′an 710072, China)

An improved hp-adaptive pseudospectral method is developed in the present paper for solving hypersonic vehicle trajectory optimization in ascent phase. The method adopts double-layer optimization strategy to adjust both the mesh spacing and the degree of the polynomial on each interval. First, a relative error estimate is derived based on the residual of the differential-algebraic constraints at sample points. The derived relative error estimate is then used to predict the polynomial degree within a mesh interval. If the polynomial degree estimate remains below a maximum allowable degree, the degree of polynomial is increased to degree estimate. Otherwise, the mesh interval is divided into subinterval. Simulation results show the proposed method is more efficient than either h method or p method.

non-smooth; ascent phase; hp-adaptive; pseudospectral method; double-layer optimization

2016-04-02

國家自然科學基金(61104195、61503302)與西北工業大學基礎研究基金(GKEY1003)資助

劉瑞帆(1995—)，女，西北工業大學碩士研究生，主要從事高超聲速飛行器制導研究。

V412.4

A

1000-2758(2016)05-0790-08