劉瑞帆, 于云峰, 閆斌斌

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

?

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

劉瑞帆, 于云峰, 閆斌斌

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

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

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

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

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

1 上升段軌跡優(yōu)化數(shù)學模型

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

1.1 動力學模型

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

(1)

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

1.2 推力模型

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

1) 推力系數(shù)Cp為

(2)

推力P為

(3)

2) 特征沖量Isp為

(4)

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

1.3 氣動力模型

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

(5)

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

(6)

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

1.4 約束條件

1.4.1 熱流率約束

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

(7)

式中,KQ為熱流率標準化常數(shù)。

1.4.2 動壓約束

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

(8)

1.4.3 過載約束

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

(9)

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

1.4.4 控制量和狀態(tài)量約束

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

(10)

1.4.5 終端約束

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

(11)

1.5 目標函數(shù)

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

(12)

1.6 多區(qū)間軌跡優(yōu)化問題描述

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

(13)

約束條件為

(14)

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

(15)

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

2 Radau偽譜法

在RPM中,連續(xù)Bolza問題的狀態(tài)變量近似為

(16)

(17)

(18)

(19)

(20)

3 改進的hp自適應偽譜法

3.1 誤差評估準則

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

(21)

(22)

(23)

(24)

3.2 多項式階數(shù)的估計

(25)

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

(26)

3.3 hp自適應迭代策略

(27)

每個區(qū)間內(nèi)的配點數(shù)均設為配點數(shù)初始值Nmin。

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

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

①初始化。令迭代次數(shù)M=0,設定區(qū)間的初始劃分以及區(qū)間內(nèi)配點數(shù)上下限Nmin,Nmax;

②將最優(yōu)控制問題離散化,轉化為NLP問題進行求解,并按照3.1節(jié)內(nèi)容計算解的相對最大誤差emax;

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

④在emax>ε的區(qū)間內(nèi)采用3.3節(jié)的自適應策略更新網(wǎng)格,直到所有的網(wǎng)格完成更新;

⑤令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°。爬升過程中狀態(tài)變量、控制變量和路徑約束的邊界限制如下

表1 約束邊界設置

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

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

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

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

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

表2 不同偽譜法的優(yōu)化比較

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

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

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

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

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

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

表3 不同hp自適應策略的優(yōu)化對比

5 結 論

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

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

[2] 雍恩米,陳磊,唐國金.飛行器軌跡優(yōu)化數(shù)值方法綜述[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] 黃長強,國海峰,丁達理. 高超聲速滑翔飛行器軌跡優(yōu)化與制導綜述[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協(xié)同航跡規(guī)劃[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] 彭祺擘,李海陽,沈紅新. 基于高斯-偽譜法的月球定點著陸軌道快速優(yōu)化設計[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偽譜法的臨近空間飛行器上升段軌跡優(yōu)化[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] 趙吉松,谷良賢,佘文學. 配點法和網(wǎng)格細化技術用于非光滑軌跡優(yōu)化[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] 陳琦. 求解非光滑最優(yōu)控制問題的自適應網(wǎng)格優(yōu)化[J]. 系統(tǒng)工程與電子技術, 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自適應偽譜法在滑翔彈道快速優(yōu)化中的應用[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偽譜法的再入飛行器軌跡優(yōu)化[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)與西北工業(yè)大學基礎研究基金(GKEY1003)資助

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

V412.4

A

1000-2758(2016)05-0790-08