楊良, 周浩, 陳萬(wàn)春

(北京航空航天大學(xué) 宇航學(xué)院, 北京 100191)

高超聲速飛行器末段軌跡快速優(yōu)化

楊良, 周浩, 陳萬(wàn)春

(北京航空航天大學(xué) 宇航學(xué)院, 北京 100191)

利用高斯偽譜法求解具有路徑約束和落角、落速等終端約束條件下的高超聲速飛行器末段軌跡優(yōu)化問題。分析了不同LG節(jié)點(diǎn)數(shù)目對(duì)該問題求解精度和計(jì)算時(shí)間的影響,以及不同LG節(jié)點(diǎn)的擬合多項(xiàng)式與實(shí)際積分結(jié)果之間的相對(duì)誤差關(guān)系。根據(jù)精度需要,選擇合適的LG節(jié)點(diǎn)進(jìn)行優(yōu)化,分析了狀態(tài)變量與積分結(jié)果的相對(duì)誤差傳播情況。仿真結(jié)果表明,采用節(jié)點(diǎn)分析后的高斯偽譜法在求解該問題上具有較高的計(jì)算精度和計(jì)算效率。

高超聲速飛行器; 軌跡優(yōu)化; 高斯偽譜法

0 引言

高超聲速飛行器對(duì)控制變量比較敏感,并且需要考慮最大動(dòng)壓約束和最大過載約束等情況,因此,可行解就在一個(gè)非常狹窄的范圍內(nèi),許多學(xué)者對(duì)這一問題進(jìn)行了研究。近十幾年來(lái),隨著偽譜法在優(yōu)化領(lǐng)域的應(yīng)用,能夠較快地尋到最優(yōu)解[1]。文獻(xiàn)[2]使用間接法對(duì)高超聲速飛行器的滑行軌跡進(jìn)行了優(yōu)化。文獻(xiàn)[3-4]分別利用高斯偽譜法對(duì)高超聲速飛行器的再入軌跡進(jìn)行了優(yōu)化,文獻(xiàn)[4]同時(shí)提出了分段串行優(yōu)化策略。文獻(xiàn)[5]利用Legendre偽譜法對(duì)高超聲速飛行器的再入軌跡進(jìn)行了優(yōu)化,同時(shí)利用Hamilton函數(shù)分析了解的最優(yōu)性。文獻(xiàn)[6]利用高斯偽譜法對(duì)路徑點(diǎn)與禁飛區(qū)等多過程約束的高超聲速飛行器進(jìn)行了軌跡優(yōu)化。這些研究都集中在高超聲速飛行器的滑翔段,很少有研究對(duì)不同LG節(jié)點(diǎn)與求解精度進(jìn)行定量的分析。本文針對(duì)高超聲速飛行器末段軌跡優(yōu)化這一特定問題進(jìn)行研究,分析了不同節(jié)點(diǎn)對(duì)計(jì)算時(shí)間和計(jì)算精度的影響,利用Hamilton函數(shù)分析了求解結(jié)果的最優(yōu)性,計(jì)算得到不同LG節(jié)點(diǎn)的擬合多項(xiàng)式與實(shí)際積分結(jié)果之間的相對(duì)誤差關(guān)系。通過相對(duì)誤差的關(guān)系,根據(jù)精度要求,選擇合適的LG節(jié)點(diǎn),仿真結(jié)果表明,選取合適LG節(jié)點(diǎn)的高斯偽譜法[7]在求解該問題上具有較高的計(jì)算精度和計(jì)算效率,最后分析了不同狀態(tài)變量之間的誤差傳播關(guān)系。對(duì)于高超聲速飛行器的末段軌跡優(yōu)化問題,準(zhǔn)確的選點(diǎn)能夠保證較高的計(jì)算精度和計(jì)算效率,這為高斯偽譜法的在線制導(dǎo)提供了較好的仿真數(shù)據(jù)。

1 高斯偽譜法介紹

考慮下面的最優(yōu)控制問題,狀態(tài)量x(τ)∈Rn,控制量u(τ)∈Rm,初始時(shí)刻和終端時(shí)刻為t0和tf。

J[u,x,τf]=Φ[x(τ0),x(τf),τ0,τf]+

φ[x(-1),x(1),t0,tf]=0∈Rq

C[x(-1),x(1),t0,tf]≤0∈Rc

式中,φ為過程等式約束;C為過程不等式約束。一般的最優(yōu)控制問題可以通過下面的公式轉(zhuǎn)化為時(shí)間區(qū)域?yàn)棣印蔥-1,1]的問題。

t=(tf-t0)τ/2+(tf-t0)/2

高斯偽譜法與文獻(xiàn)[8-9]中介紹的Legendre偽譜法和Chebyshev偽譜法一樣,都是采用插值多項(xiàng)式來(lái)逼近狀態(tài)變量和控制變量。

在LG節(jié)點(diǎn)的每個(gè)Lagrange多項(xiàng)式的微分方程可以用一個(gè)微分矩陣來(lái)逼近,這個(gè)矩陣就是D∈RN×N+1,微分矩陣的元素可以表示為:

式中,k=1,…,N;i=0,…,N。這樣狀態(tài)微分方程就可以通過微分矩陣轉(zhuǎn)換為代數(shù)約束,可表示為:

式中,ωk為高斯權(quán)函數(shù)。連續(xù)的性能指標(biāo)函數(shù)也可以通過積分公式表示為:

J=Φ[x(-1),x(1),t0,tf]+

φ(X0,Xf,t0,tf)=0,C(Xk,Uk,τk;t0,tf)≤0

不規(guī)范處方主要包括：①處方超量，如普通處方超過7天量，急診處方超過3天量，省特約處方超過7天量，無(wú)特殊治療需要處方超過1個(gè)月用量等；②臨床診斷書寫不全、病人信息不全；③醫(yī)生輸入錯(cuò)誤；④醫(yī)師處方未簽名（蓋章）或與留樣不一致；⑤處方上藥品品種數(shù)量超過規(guī)定（5種）；⑥處方修改未簽名蓋章，也未注明修改日期；⑦處方用箋不規(guī)范，如門診手寫處方用的是住院處方或自費(fèi)藥房處方。

這樣就把連續(xù)的Bolza問題轉(zhuǎn)換為非線性規(guī)劃(NLP)問題,通過成熟的非線性規(guī)劃求解器就可以找到該問題的最優(yōu)解。

2 高超聲速飛行器動(dòng)力學(xué)模型

2.1 動(dòng)力學(xué)方程

大氣相對(duì)于地球靜止的無(wú)動(dòng)力三自由度質(zhì)點(diǎn)彈道運(yùn)動(dòng)學(xué)方程[10]為:

式中,x,y,z為大地坐標(biāo);V為地球相對(duì)速度;γ為航跡角;φ為航向角;α為迎角;σ為側(cè)傾角;g=μ/r2為重力加速度;μ為地球重力加速度常數(shù)。升力L和阻力D如下:

L=qSrefCL/m,D=qSrefCD/m

式中,q=ρV2/2為飛行過程中的動(dòng)壓;CL為升力系數(shù);CD為阻力系數(shù);Sref為飛行器的參考面積;m為飛行器質(zhì)量;ρ為大氣密度,ρ=ρ0exp(-βH),其中ρ0為海平面大氣密度,H=r-Re為高度,Re為地球半徑,β為密度常數(shù)。

2.2 約束條件

從結(jié)構(gòu)和熱防護(hù)的角度出發(fā),要求高超聲速飛行器在飛行過程中滿足動(dòng)壓約束以及過載約束,其約束的表達(dá)式為:

q=(1/2)ρV2≤qmax,nz=L/(mg)≤nzmax

式中,qmax,nzmax根據(jù)具體飛行器的設(shè)計(jì)進(jìn)行選取。

|u1|≤u1max, |u2|≤u2max

2.3 優(yōu)化性能指標(biāo)

為了保證優(yōu)化得到的控制變量能夠足夠光滑,選取各個(gè)控制變量的平方和最小這一性能指標(biāo)作為優(yōu)化性能指標(biāo)。各個(gè)控制變量的平方和最小也意味著飛行器在飛行過程中以消耗能量最小的方式達(dá)到所滿足的條件,其具體表達(dá)形式為:

3 數(shù)值計(jì)算與分析

3.1LG節(jié)點(diǎn)數(shù)目的影響

選取高斯偽譜法的LG節(jié)點(diǎn)數(shù)為60,30,15,10,6分別對(duì)該問題進(jìn)行優(yōu)化。所有計(jì)算均在Window 7的系統(tǒng)環(huán)境下運(yùn)行MATLAB 2008B進(jìn)行計(jì)算,非線性規(guī)劃問題采用著名的snopt軟件包進(jìn)行求解。

采用不同LG節(jié)點(diǎn)進(jìn)行優(yōu)化的仿真結(jié)果如圖1所示。由圖可知,這些結(jié)果非常相似,多次優(yōu)化的結(jié)果都同時(shí)收斂到最優(yōu)解。從各狀態(tài)變量的曲線來(lái)看,其最優(yōu)結(jié)果都非常光滑。

通過對(duì)不同擬合多項(xiàng)式的相對(duì)節(jié)點(diǎn)誤差比較,采用階數(shù)相對(duì)較少的多項(xiàng)式直接逼近各狀態(tài)變量就可以達(dá)到較高的逼近精度,這也意味著能夠選取很少的LG節(jié)點(diǎn)達(dá)到較高的計(jì)算精度。LG節(jié)點(diǎn)的減少同時(shí)也能夠帶來(lái)通過高斯偽譜法離散后的非線性規(guī)劃問題的規(guī)模減少。但隨著節(jié)點(diǎn)數(shù)目的減少,計(jì)算精度也會(huì)下降,從圖中可以看出,采用6個(gè)LG節(jié)點(diǎn)的優(yōu)化結(jié)果明顯出現(xiàn)較大的偏差。

圖1 優(yōu)化結(jié)果仿真曲線Fig.1 Curve of optimization simulation

表1分別對(duì)不同LG節(jié)點(diǎn)數(shù)進(jìn)行優(yōu)化的計(jì)算時(shí)間進(jìn)行對(duì)比。隨著節(jié)點(diǎn)的減少,優(yōu)化變量和計(jì)算時(shí)間大大減少,優(yōu)化問題的規(guī)模得到了很大的改善。從表中可以看出,采用10個(gè)LG節(jié)點(diǎn)的計(jì)算時(shí)間相對(duì)于60個(gè)LG節(jié)點(diǎn)也從520.71 s減少到4.975 1 s。

表1 不同LG節(jié)點(diǎn)的優(yōu)化結(jié)果Table 1 Optimization of various LG nodes

3.2 最優(yōu)性驗(yàn)證

根據(jù)文獻(xiàn)[7],可以通過約束的Lagrange乘子估算出最優(yōu)控制問題的協(xié)態(tài)變量。這樣就提供了一種對(duì)所得結(jié)果進(jìn)行最優(yōu)性驗(yàn)證的方法。

根據(jù)最優(yōu)控制理論,在約束不受限的情況下,Hamilton方程恒等于零;在受限情況下,其Hamilton函數(shù)應(yīng)該盡可能地接近于不受限的Hamilton函數(shù),這樣就為解的最優(yōu)性提供了一個(gè)驗(yàn)證(見圖2)。在Hamilton函數(shù)曲線上,Hamilton函數(shù)的數(shù)值結(jié)果基本接近于零,通過Hamilton函數(shù)的近似性,保證了解的最優(yōu)性。

圖2 Hamilton函數(shù)Fig.2 Hamilton function

3.3 優(yōu)化結(jié)果比較與誤差分析

圖3為不同LG節(jié)點(diǎn)與積分彈道的相對(duì)誤差之間的關(guān)系。從圖中可以看出,隨著LG節(jié)點(diǎn)的增加,誤差精度越來(lái)越小,但變化趨勢(shì)趨于平緩,在較少的LG節(jié)點(diǎn)時(shí),精度出現(xiàn)比較大的變化,因此,根據(jù)精度要求,合理選擇LG節(jié)點(diǎn)個(gè)數(shù)非常重要。圖4為采用10個(gè)LG節(jié)點(diǎn)的高斯偽譜法優(yōu)化彈道與積分彈道進(jìn)行對(duì)比的結(jié)果,通過MATLAB 2008B中的Simulink搭建動(dòng)力學(xué)模型實(shí)現(xiàn)彈道積分,積分方法選用四階R-K法。從圖中可以看出,優(yōu)化結(jié)果與積分結(jié)果基本相同。

圖3 相對(duì)誤差精度與LG節(jié)點(diǎn)的關(guān)系Fig.3 Relationship between relative error and LG nodes

圖4 優(yōu)化結(jié)果對(duì)比曲線Fig.4 Comparison of optimization results

圖5為各個(gè)狀態(tài)變量的相對(duì)誤差傳播變化情況,是通過各個(gè)狀態(tài)變量的積分結(jié)果與優(yōu)化結(jié)果的Lagrange插值多項(xiàng)式相減,然后除以相對(duì)系數(shù)。從圖中可以看出,各個(gè)狀態(tài)量之間的相對(duì)誤差都很小,且呈振蕩分布,在LG節(jié)點(diǎn)上達(dá)到優(yōu)化算法的可行性精度。x和y方向的相對(duì)誤差基本上保持在10-5數(shù)量級(jí)上,而z方向、速度和航跡角的相對(duì)誤差基本保持在10-4數(shù)量級(jí)上,最小的為航向角相對(duì)誤差,保持在10-8數(shù)量級(jí)上。因此,采用10個(gè)LG節(jié)點(diǎn)的高斯偽譜法在處理高超聲速飛行器末段軌跡快速優(yōu)化問題上,能夠很快地收斂到全局最優(yōu)解,并且具有很高的計(jì)算效率與計(jì)算精度。

圖5 誤差的仿真結(jié)果Fig.5 Results of error simulation

4 結(jié)束語(yǔ)

本文針對(duì)高超聲速飛行器的末段軌跡快速優(yōu)化問題進(jìn)行了研究,采用高斯偽譜法進(jìn)行優(yōu)化。分析了不同LG節(jié)點(diǎn)對(duì)計(jì)算時(shí)間和計(jì)算效率的影響,不同LG節(jié)點(diǎn)的擬合多項(xiàng)式與實(shí)際積分結(jié)果之間的相對(duì)誤差關(guān)系以及狀態(tài)的誤差傳播情況。結(jié)果表明,合理選擇LG節(jié)點(diǎn)的高斯偽譜法具有很高的計(jì)算精度和計(jì)算效率。高精度、高效率的優(yōu)化結(jié)果也為高超聲速飛行器的末段軌跡在線優(yōu)化實(shí)現(xiàn)提供了很好的計(jì)算依據(jù),下一步將以本文工作為基礎(chǔ),重點(diǎn)研究高斯偽譜法的在線優(yōu)化問題。

[1] Betts J T.Survey of numerical methods for trajectory optimization[J].Journal of Guidance,Control,and Dynamics,1998,21(2):193-207.

[2] 周浩,陳萬(wàn)春,殷興良.高超聲速飛行器滑行軌跡優(yōu)化[J].北京航空航天大學(xué)學(xué)報(bào),2006,32(5):513-517.

[3] 周文雅,楊滌,李順利.利用高斯偽譜法求解升力航天器最優(yōu)再入軌跡[J].南京理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,34(1):85-90.

[4] 雍恩米,唐國(guó)金,陳磊.基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化[J].宇航學(xué)報(bào),2008,29(6):1765-1772.

[5] Rao A V,Clarke K A.Performance optimization of a maneuvering re-entering vehicle suing a legendre pseudospectral method [R].AIAA-2002-4885,2002.

[6] Jorris T R,Rao A V.Three-dimensional trajectory optimization satisfying waypoint and no-fly zone constraints[J].Journal of Guidance,Control,and Dynamics,2009,32(2):551-572.

[7] Benson D A,Huntington G T. Direct trajectory optimization and costate estimation via an orthogonal collocation method[J].Journal of Guidance,Control,and Dynamics,2006,29(6):1435-1440.

[8] Fahroo F,Ross I M.Costate estimation by a Legendre pseudospectral method[J].Journal of Guidance,Control,and Dynamics,2001,24(2):270-277.

[9] Fahroo F,Ross I M.Direct trajectory optimization by a Chebyshev pseudospectral method [J].Journal of Guidance,Control,and Dynamics,2002,25(1):160-166.

[10] Vinh N X.Flight mechanics of high-performance aircraft[M].Cambridge,England:Cambridge University Press,1993:22-25.

Fastterminaltrajectoryoptimizationforhypersonicspacecraft

YANG Liang, ZHOU Hao, CHEN Wan-chun

(School of Astronautics, Beijing University of Aeronautics and Astronautics, Beijing 100191, China)

In this paper, the terminal trajectory optimization with the constraints on process, impact angles and impact speed for hypersonic spacecraft is solved by Gauss pesudospctral method. By various Lambert-Gauss nodes operated to optimize this problem, the effect of the number of nodes on calculation time and calculation efficiency is studied. By comparing optimal results of various nodes with integral results, the relative error is obtained. According to specified accuracy, appropriate LG node is selected to optimize this problem and the extensive analysis of propagation relationship of the relative errors between Lagrange interpolation polynomial and integrated states is carried out. The results show that Gauss pseudospectral method applying nodes analysis is of high calculation accuracy and high calculation time to solve this problem.

hypersonic spacecraft; trajectory optimization; Gauss pseudospectral method

V412.4

A

1002-0853(2013)04-0341-04

2012-11-14;

2013-03-24; < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間:2013-06-06 12:25

楊良(1985-)，男，湖南常德人，博士研究生，研究方向?yàn)轱w行動(dòng)力學(xué)及彈道優(yōu)化。

(編輯：方春玲)