求解軌跡優(yōu)化問題的局部配點(diǎn)法的稀疏性研究

趙吉松

(南京航空航天大學(xué)航天學(xué)院，南京 210016)

求解軌跡優(yōu)化問題的局部配點(diǎn)法的稀疏性研究

趙吉松

(南京航空航天大學(xué)航天學(xué)院，南京 210016)

直接配點(diǎn)法通過對(duì)控制變量和狀態(tài)變量都進(jìn)行離散將軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃(NLP)問題。為了提高NLP的求解效率，需要利用其偏導(dǎo)數(shù)的稀疏特性并建立偏導(dǎo)數(shù)的高效計(jì)算方法。本文研究了局部配點(diǎn)法離散得到的NLP的一階偏導(dǎo)數(shù)的稀疏特性，建立了一階偏導(dǎo)數(shù)的高效計(jì)算方法。推導(dǎo)了NLP的目標(biāo)函數(shù)梯度和約束雅克比矩陣的數(shù)學(xué)表達(dá)式，得到了NLP偏導(dǎo)數(shù)的稀疏型，并且將NLP的偏導(dǎo)數(shù)分解為原始軌跡優(yōu)化問題的偏導(dǎo)數(shù)。由于原始軌跡優(yōu)化問題的約束和變量的數(shù)量遠(yuǎn)少于NLP的約束和變量的數(shù)量，從而顯著減小了NLP的一階偏導(dǎo)數(shù)的計(jì)算量。含有離散氣動(dòng)力和推力數(shù)據(jù)的仿真算例驗(yàn)證了本文方法的有效性。仿真結(jié)果表明，與有限差分法直接計(jì)算NLP的偏導(dǎo)數(shù)相比，本文方法能夠?qū)?yōu)化耗時(shí)減小至4%以內(nèi)，隨著離散節(jié)點(diǎn)數(shù)目的增加，計(jì)算效率的提升更為顯著。

軌跡優(yōu)化；局部配點(diǎn)法；非線性規(guī)劃；一階偏導(dǎo)數(shù)；稀疏特性

0 引 言

軌跡優(yōu)化對(duì)于飛行器設(shè)計(jì)有著十分重要的意義和工程實(shí)際價(jià)值[1]。軌跡優(yōu)化本質(zhì)上屬于最優(yōu)控制問題，其求解方法主要分為間接法和直接法。其中，直接法中的配點(diǎn)法通過對(duì)控制變量和狀態(tài)變量進(jìn)行離散將軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃(NLP)問題，降低了對(duì)初值的敏感性并且具有很好的收斂性，近年來得到廣泛的研究和應(yīng)用[2-7]。

雖然配點(diǎn)法具有多方面優(yōu)勢(shì)，但是其具體實(shí)施方法對(duì)于提高優(yōu)化效率具有重要影響。目前，以序列二次規(guī)劃(SQP)為代表的基于梯度算法的NLP求解器需要提供NLP的目標(biāo)函數(shù)和約束的一階偏導(dǎo)數(shù)，甚至二階偏導(dǎo)數(shù)。其中，一階偏導(dǎo)數(shù)NLP求解器需要提供一階偏導(dǎo)數(shù)，然后采用擬牛頓法(DFP法或者BFGS法等)構(gòu)造近似的二階偏導(dǎo)數(shù)，比如SNOPT[8]；二階偏導(dǎo)數(shù)NLP求解器除了需要一階偏導(dǎo)數(shù)，還需要準(zhǔn)確的二階偏導(dǎo)數(shù)，比如IPOPT[9]。但是，偏導(dǎo)數(shù)的計(jì)算量通常比較大，甚至超過優(yōu)化算法本身。因此，提高NLP的一階/二階偏導(dǎo)數(shù)的計(jì)算效率對(duì)于提高軌跡優(yōu)化效率有重要意義。

研究發(fā)現(xiàn)，由軌跡優(yōu)化離散得到的NLP是非常稀疏的，即NLP的一階偏導(dǎo)數(shù)和二階偏導(dǎo)數(shù)含有大量零元素[2,10-12]。Betts等[10]較早地研究了局部配點(diǎn)法的稀疏特性，得到了梯形格式、Hermite-Simpson格式和Runge-Kutta格式的狀態(tài)方程離散殘差約束的偏導(dǎo)數(shù)的稀疏型，并將其中的非零元素轉(zhuǎn)化為最優(yōu)控制問題的偏導(dǎo)數(shù)，有效地減小了計(jì)算量。但是，Betts等[10]的研究存在一些不足之處：一是只給出了狀態(tài)方程離散殘差約束的偏導(dǎo)數(shù)稀疏特性和非零元素的計(jì)算方法，沒有給出目標(biāo)函數(shù)、路徑約束和端點(diǎn)約束的偏導(dǎo)數(shù)的稀疏型和計(jì)算方法；二是對(duì)于最常用的Hermite-Simpson格式，沒有完全探索出緊湊形式的狀態(tài)方程離散殘差的偏導(dǎo)數(shù)的稀疏型，而是通過在離散節(jié)點(diǎn)中間位置添加離散格式約束和狀態(tài)變量將其轉(zhuǎn)換為分離形式才得到完全的稀疏型。在Betts等的研究工作的基礎(chǔ)上，Patterson等[11]研究了Radau偽譜法的一階和二階偏導(dǎo)數(shù)的稀疏性，得到了完整的稀疏型，并且推導(dǎo)出非零元素的高效計(jì)算方法(將NLP偏導(dǎo)數(shù)分解為最優(yōu)控制問題的偏導(dǎo)數(shù))。因此，如果能夠基于上述研究，完全探索出局部配點(diǎn)法的NLP偏導(dǎo)數(shù)的稀疏型并建立偏導(dǎo)數(shù)的高效計(jì)算方法，對(duì)于提高局部配點(diǎn)法的優(yōu)化效率具有重要意義。此外，與全局配點(diǎn)法(又稱偽譜法，離散節(jié)點(diǎn)是正交多項(xiàng)式的根)相比，局部配點(diǎn)法的離散節(jié)點(diǎn)可以根據(jù)需要任意布置，在網(wǎng)格細(xì)化方面具有更好的靈活性[13-19]，適合求解非光滑軌跡優(yōu)化問題。因而，提高局部配點(diǎn)法的優(yōu)化效率還有能夠促進(jìn)局部配點(diǎn)法在非光滑軌跡優(yōu)化領(lǐng)域的應(yīng)用。

在實(shí)際應(yīng)用中，一階偏導(dǎo)數(shù)NLP求解器比二階偏導(dǎo)數(shù)NLP求解器更為常用，因?yàn)槎A偏導(dǎo)數(shù)的計(jì)算通常比較繁瑣并且計(jì)算量較大。以NLP的一階偏導(dǎo)數(shù)為例，其計(jì)算方法可分為兩大類。一類方法是直接計(jì)算NLP的偏導(dǎo)數(shù)，包括自動(dòng)微分法[20]、復(fù)變量微分法[21]、有限差分法等。其中自動(dòng)微分法計(jì)算量小，精度高，得到了廣泛應(yīng)用，其局限性在于要求優(yōu)化模型解析可導(dǎo)，不適用于帶有離散數(shù)據(jù)的軌跡優(yōu)化問題。對(duì)于帶有離散數(shù)據(jù)的優(yōu)化模型，只能采用有限差分法。但是，采用有限差分直接計(jì)算NLP偏導(dǎo)數(shù)的效率較低，需要耗費(fèi)大量的機(jī)時(shí)。另一類方法是將NLP的偏導(dǎo)數(shù)分解為最優(yōu)控制問題的偏導(dǎo)數(shù)，然后采用各種微分算法計(jì)算最優(yōu)控制問題的偏導(dǎo)數(shù)，并組裝得到NLP的偏導(dǎo)數(shù)。因?yàn)榕cNLP相比，最優(yōu)控制的約束和變量的數(shù)量大幅減少，因而這樣處理可以顯著提高偏導(dǎo)數(shù)的計(jì)算效率。

本文在Betts等[10]和Patterson等[11]的研究工作的基礎(chǔ)上，以Hermite-Simpson格式為例，研究了局部配點(diǎn)法離散得到的NLP的一階偏導(dǎo)數(shù)(目標(biāo)函數(shù)梯度和約束雅克比矩陣)的稀疏性，建立非零元素的高效計(jì)算方法。采用帶有離散參數(shù)模型的優(yōu)化算例驗(yàn)證了所述方法的有效性。仿真結(jié)果表明，與采用有限差分法直接計(jì)算NLP的偏導(dǎo)數(shù)相比，本文方法能夠?qū)?yōu)化耗時(shí)減小至4%以下，并且隨著離散節(jié)點(diǎn)數(shù)量的增加，計(jì)算效率的提升更為顯著。

1 軌跡優(yōu)化問題數(shù)學(xué)描述

軌跡優(yōu)化問題本質(zhì)上屬于最優(yōu)控制問題，以Bolza型最優(yōu)控制問題為例，可描述為：求解控制變量u(t)∈Rm，使得如下目標(biāo)函數(shù)最小化

(1)

式中:M:Rn×R×Rn×R→R,L:Rn×Rm×R→R,x∈Rn,u∈Rm,t∈[t0,tf]?R。

狀態(tài)方程為

(2)

端點(diǎn)條件為

E(x(t0),t0,x(tf),tf)=0

(3)

路徑約束為

C(x(t),u(t),t)≤0,t∈[t0,tf]

(4)

式中f:Rn×Rm×R→Rn,E:Rn×R×Rm×R→Re,1:Rn×Rm×R→Rc。方程(1)-(4)所描述的問題稱為連續(xù)Bolza型最優(yōu)控制問題。

2 基于Runge-Kutta格式的配點(diǎn)法離散

首先利用積分變換τ=(t-t0)/(tf-t0)將軌跡優(yōu)化問題(方程(1)-(4))變換至?xí)r間區(qū)間τ([0,1]。假設(shè)單位區(qū)間[0,1]上的N個(gè)離散區(qū)間的節(jié)點(diǎn)為

τN=τf=1;τi<τi+1,i=0,1,…,N-1}

(5)

式中τi稱為節(jié)點(diǎn)或網(wǎng)格點(diǎn)，τi在[0, 1]上可以均勻分布，也可以非均勻分布。

記xi=x(τi),ui=u(τi), 對(duì)于狀態(tài)方程，基于q階Runge-Kutta(RK)方法的離散格式為

(6)

式中：Δt=tf-t0,hi=τi+1-τi,fij=f(xij,uij,τij;t0,tf),xij,uij和τij為中間變量，xij由下式給出

(7)

式中：τij=τi+hiρj,uij=u(τij).ρj,βj,αjl均為已知常數(shù)并且滿足0≤ρ1≤ρ2≤…≤ρq≤1。當(dāng)αjl=0(l≥j)時(shí)，離散格式為顯式格式，否則為隱式格式。采用類似的方法，可將目標(biāo)函數(shù)可離散化。常用的離散格式包括梯形格式(q=2)，Hermite-Simpson格式(q=3)，以及經(jīng)典四階Runge-Kutta格式(q=4)。

(8)

并且滿足如下約束

(9)

Ci=C(xi,ui,τi;t0,tf)≤0

(10)

(11)

E(x0,t0,xf,tf)=0

(12)

式中

常用的離散格式有梯形格式(q= 2)，Hermite- Simpson格式(q= 3，簡(jiǎn)記HS格式)，以及經(jīng)典四階Runge-Kutta格式(q= 4，簡(jiǎn)記RK格式)。

以HS格式為例，該格式需要用到區(qū)間中點(diǎn)的變量和函數(shù)值，為此將區(qū)間中點(diǎn)的控制變量作為優(yōu)化變量，并且在區(qū)間中點(diǎn)添加路徑約束，即

(13)

(14)

約束條件為

(15)

Ci=C(xi,ui,τi;t0,tf)≤0

(16)

(17)

E(x0,t0,xf,tf)=0

(18)

其中

在數(shù)值優(yōu)化時(shí)，為了使問題具有實(shí)際物理意義，還需要添加時(shí)間差約束

Δt=tf-t0>0

(19)

3 NLP偏導(dǎo)數(shù)計(jì)算方法

3.1 依賴關(guān)系矩陣

在推導(dǎo)NLP一階偏導(dǎo)數(shù)的稀疏特性時(shí)，需要用到原始軌跡優(yōu)化問題對(duì)自變量的依賴關(guān)系。

由于狀態(tài)方程、路徑約束和目標(biāo)函數(shù)Lagrange積分項(xiàng)都定義在整個(gè)時(shí)域區(qū)間，因而本文將這三項(xiàng)對(duì)自變量的偏導(dǎo)數(shù)定義在一起，

(20)

其中G1的每一項(xiàng)仍為矩陣，以?f/?xT為例，

(21)

易知，G1為(n+c+1)×(n+m+1)維矩陣。

通常情況下，G1是稀疏矩陣。為了描述G1的稀疏型，定義如下struct函數(shù)

(22)

記

S1=struct(G1)

(23)

式中struct (G1)表示對(duì)G1的每個(gè)元素進(jìn)行struct運(yùn)算。S1表示G1的稀疏型。為了得到S1，不需要計(jì)算G1的每個(gè)元素的具體值，只需要判斷是否為0。

類似可以定義端點(diǎn)約束和目標(biāo)函數(shù)的Mayer項(xiàng)對(duì)自變量的依賴關(guān)系矩陣和稀疏型

(24)

S2=struct(G2)

(25)

易知，G2為(e+1)×2(n+1)維矩陣。

3.2 變量記法

前述離散格式將同一個(gè)節(jié)點(diǎn)處的變量或約束記為一個(gè)列向量，這種記法與數(shù)值積分格式的形式一致，但是不利于推導(dǎo)偏導(dǎo)數(shù)矩陣的稀疏特性。為此，本文定義一種新的變量記法，將變量或約束的同一個(gè)分量在不同節(jié)點(diǎn)的值記為一個(gè)新的向量。

以狀態(tài)方程離散殘差為例，定義

ξ:,j=(ξ0,j,ξ1,j,…,ξN-1,j)T, (j=1,2,…,n)

(26)

J(z)

(27)

并且滿足約束條件

Fmin≤F(z)≤Fmax

(28)

式中：目標(biāo)函數(shù)J(z)的表達(dá)式參見方程，優(yōu)化變量z和約束函數(shù)F(z)的定義如下

(29)

3.3 目標(biāo)函數(shù)梯度

目標(biāo)函數(shù)梯度是指目標(biāo)函數(shù)對(duì)優(yōu)化變量的偏導(dǎo)數(shù)，具體定義如下

(30)

將目標(biāo)函數(shù)寫成矩陣乘積形式可得到

(31)

(32)

對(duì)角陣h=diag(h0,h1,…,hN-1)，其中hi(i=0, 1, …,N-1)為積分步長(zhǎng)，矩陣D1和D2定義如下

D1和D2均為N×(N+1)矩陣，其中空白元素為0。

(33)

式中

可見，NLP的目標(biāo)函數(shù)梯度可以分解為軌跡優(yōu)化問題的目標(biāo)函數(shù)和狀態(tài)方程的偏導(dǎo)數(shù)。

3.4 雅克比矩陣

NLP的雅克比矩陣定義為NLP的約束對(duì)優(yōu)化變量的偏導(dǎo)數(shù)矩陣，對(duì)于HS格式，形式如下

(34)

式中向量F和z的定義參見方程(29)。雅克比矩陣GF的展開形式遵循向量求偏導(dǎo)數(shù)運(yùn)算規(guī)則(參見方程(21))。下文推導(dǎo)GF的數(shù)學(xué)表達(dá)式。

3.4.1狀態(tài)方程離散殘差約束的偏導(dǎo)數(shù)

結(jié)合前述定義的變量記法，將狀態(tài)方程離散殘差約束即方程寫成如下形式

(35)

(36)

式中

3.4.2路徑約束的偏導(dǎo)數(shù)

(37)

(38)

3.4.3端點(diǎn)約束的偏導(dǎo)數(shù)

(39)

3.4.4時(shí)間差約束的偏導(dǎo)數(shù)

(40)

可見，NLP的雅克比矩陣可分解為軌跡優(yōu)化問題的狀態(tài)方程、路徑約束、端點(diǎn)約束和時(shí)間約束的偏導(dǎo)數(shù)。計(jì)算出這些約束在離散節(jié)點(diǎn)和區(qū)間中點(diǎn)的偏導(dǎo)數(shù)(對(duì)于f和C)以及端點(diǎn)處的偏導(dǎo)數(shù)(對(duì)于E和Δt)之后，采用本節(jié)的方法組裝得到雅克比矩陣。對(duì)于HS格式，以N=4為例，其雅克比矩陣的稀疏型如圖1所示。其中，空白元素表示恒為零，“·”表示非零元素，“×”表示可能不為零的元素。

(41)

4 仿真算例

飛行器最短時(shí)間爬升問題最初由Bryson等[22]提出，此后得到廣泛研究[2]。該算例的特色之處是推力和氣動(dòng)力數(shù)據(jù)以離散表格形式給出，與實(shí)際工程問題比較接近?；締栴}是求解最優(yōu)攻角α(t)，使得飛行器從跑道起飛爬升至指定高度消耗的時(shí)間最短。在縱向平面內(nèi)，飛行器的運(yùn)動(dòng)方程組為[22]

(42)

式中h為高度，v為速度，γ為航跡角，m為質(zhì)量，μ為引力常數(shù)，Re為地球半徑。發(fā)動(dòng)機(jī)推力T(Ma,h)由表 1給出(缺少的數(shù)據(jù)實(shí)際飛行不會(huì)用到)，Isp=1600 s，g0=9.80665 m/s2。氣動(dòng)力由下式給出[22]

式中L為升力，D為阻力，CL為升力系數(shù)，CD為阻力系數(shù)，S為參考面積，ρ為大氣密度。氣動(dòng)力的相關(guān)參數(shù)CLα，CD0和η由表 2給出。

初始條件是h(t0)=0 m,v(t0)=129.314 m/s,γ(t0)=0°,m(t0)=19050.864 kg；終端約束是h(tf)=19994.88 m,v(tf)= 295.092 m/s,γ(tf)= 0°。

Ma高度h/km01．5243．0484．5726．0967．629．14412．19215．2421．3360．024．20．228．024．621．118．115．212．810．70．428．325．221．918．715．913．411．27．34．40．630．827．223．820．517．314．712．38．14．90．834．530．326．623．219．816．814．19．45．61．11．037．934．330．426．823．319．816．811．26．81．41．236．138．034．931．327．323．620．113．48．31．71．436．638．536．131．628．124．216．210．02．21．638．735．732．028．119．311．92．91．834．631．121．713．33．1

表2 氣動(dòng)力數(shù)據(jù)Table 2 Aerodynamic data

目標(biāo)函數(shù)為飛行器爬升消耗的時(shí)間最短，即

J=tf

(43)

本算例的特點(diǎn)是推力數(shù)據(jù)和氣動(dòng)力數(shù)據(jù)是離散形式，并且推力數(shù)據(jù)不完整。文獻(xiàn)[2]通過對(duì)推力數(shù)據(jù)進(jìn)行最小曲率樣條擬合，得到了完整、光滑的曲面，但是處理過程比較復(fù)雜，難以通用。本文采用線性外插法將表1中的推力數(shù)據(jù)補(bǔ)充完整，采用三次樣條插值計(jì)算飛行過程的推力和氣動(dòng)力參數(shù)。大氣模型采用美國(guó)1976 年標(biāo)準(zhǔn)大氣進(jìn)行插值。

采用局部配點(diǎn)法將該軌跡優(yōu)化問題轉(zhuǎn)化為NLP問題，采用SNOPT 7[8]求解NLP問題。對(duì)NLP一階偏導(dǎo)數(shù)的兩種計(jì)算方案進(jìn)行了研究。一種方案是不提供NLP的一階偏導(dǎo)數(shù)，此時(shí)SNOPT內(nèi)部采用有限差分法直接計(jì)算NLP的偏導(dǎo)數(shù)。另一種方案即本文方法，采用有限差分法計(jì)算原始軌跡優(yōu)化問題的一階偏導(dǎo)數(shù)，采用第3節(jié)給出的方法組裝得到NLP的一階偏導(dǎo)數(shù)矩陣并提供給SNOPT。采用Matlab 2009a語言編程，表 3給出采用32、64和128個(gè)均勻離散節(jié)點(diǎn)時(shí)兩種方案的優(yōu)化效率對(duì)比。采用的計(jì)算機(jī)為MacBook Air(處理器Intel Core i5-5250U 1.6 GHz，操作系統(tǒng)Windows 10企業(yè)版，內(nèi)存4 GB 1600 MHz DDR3)，計(jì)算耗時(shí)為10次運(yùn)行的平均耗時(shí)。從表 3可以看出，與采用有限差分法直接計(jì)算NLP的偏導(dǎo)數(shù)相比，本文方法能夠?qū)⒂?jì)算耗時(shí)減小到4%以下，并且隨著離散節(jié)點(diǎn)數(shù)目的增加，本文方法的優(yōu)化效率提升更為顯著。這是因?yàn)殡S著離散節(jié)點(diǎn)數(shù)量增加，NLP一階偏導(dǎo)數(shù)矩陣具有更大的稀疏度(稀疏度定義為NLP的一階偏導(dǎo)數(shù)矩陣中的零元素?cái)?shù)量所占的比例)，參見表 3。

本算例的最優(yōu)飛行攻角在部分區(qū)域變化比較劇烈。為了高精度捕捉最優(yōu)攻角，本文將前述建立的偏導(dǎo)數(shù)高效計(jì)算方法與作者最近發(fā)展的網(wǎng)格細(xì)化技術(shù)[16]相結(jié)合求解該問題，優(yōu)化結(jié)果如圖 2和圖 3所示。圖中的圓圈為離散最優(yōu)解，實(shí)線為基本無振蕩插值(ENO)結(jié)果(對(duì)于控制變量)或者數(shù)值積分結(jié)果(對(duì)于狀態(tài)變量)。網(wǎng)格細(xì)化算法迭代6次，從641個(gè)均勻離散節(jié)點(diǎn)中選取82個(gè)節(jié)點(diǎn)求解該問題，總耗時(shí)29.2 s，優(yōu)化的最短爬升時(shí)間J* = 320.25 s。文獻(xiàn)[2]的優(yōu)化結(jié)果為324.98 s，二者差異主要是由采用的大氣模型不同和對(duì)推力數(shù)據(jù)的處理方式不同造成的。文獻(xiàn)[18]同樣采用局部配點(diǎn)法和網(wǎng)格細(xì)化技術(shù)求解該問題，但是采用了有限差分法直接計(jì)算NLP的偏導(dǎo)數(shù)，優(yōu)化耗時(shí)長(zhǎng)達(dá)12.1分鐘(CPU頻率2.59 GHZ)?？梢姡疚姆椒ㄅc網(wǎng)格細(xì)化技術(shù)相結(jié)合，能夠快速、高精度地求解軌跡優(yōu)化問題。

表3 兩種方法計(jì)算NLP偏導(dǎo)數(shù)的優(yōu)化效率對(duì)比Table 3 Efficiency comparison of two derivative methods

5 結(jié) 論

本文以HS格式為例，研究了局部配點(diǎn)離散得到的NLP問題的稀疏特性，推導(dǎo)出NLP一階偏導(dǎo)數(shù)的高效計(jì)算方法。首先引入一種新的變量記法將NLP問題寫成向量形式，然后應(yīng)用向量鏈?zhǔn)角髮?dǎo)規(guī)則將NLP的偏導(dǎo)數(shù)轉(zhuǎn)化為原始軌跡優(yōu)化問題的偏導(dǎo)數(shù)。由于與NLP相比，軌跡優(yōu)化問題的約束和自變量的數(shù)量大幅度減少，因而這樣處理可以顯著減小NLP的一階偏導(dǎo)數(shù)的計(jì)算量。采用帶有離散推力數(shù)據(jù)和氣動(dòng)力數(shù)據(jù)的仿真算例驗(yàn)證了本文方法的有效性。仿真結(jié)果表明，與采用有限差分法直接計(jì)算NLP的偏導(dǎo)數(shù)相比，采用本文方法能夠大幅度減小優(yōu)化耗時(shí)(對(duì)于算例，減小至4%以下)，并且隨著離散節(jié)點(diǎn)數(shù)量的增加，本文方法計(jì)算效率的提升更為顯著。雖然本文的研究工作基于HS格式，但是所給出的方法可容易推廣至局部配點(diǎn)法的其它離散格式，比如梯形格式，經(jīng)典四階RK格式等。

[1] 唐國(guó)金, 羅亞中, 雍恩米. 航天器軌跡優(yōu)化理論、方法及應(yīng)用 [M]. 北京: 科學(xué)出版社, 2012.

[2] Betts J T. Practical methods for optimal control and estimation using nonlinear programming, advances in design and control series [M]. Philadelphia: Soc. for Industrial and Applied Mathematics, 2009.

[3] 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.

[4] Benson, David. A Gauss pseudospectral transcription for optimal control [D]. Cambridge: Massachusetts Institute of Technology, 2005.

[5] 雍恩米, 唐國(guó)金, 陳磊. 基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化 [J]. 宇航學(xué)報(bào), 2008, 29(6): 1766-1772. [Yong En-mi, Tang Guo-jin , Chen Lei. Rapid trajectory optimization for hypersonic reentry vehicle via Gauss pseudospectral method [J]. Journal of Astronautics, 2008, 29(6): 1766-1772.]

[6] 丁洪波, 曹淵, 佟衛(wèi)平, 等. 亞軌道飛行器上升段軌跡優(yōu)化與快速重規(guī)劃 [J]. 宇航學(xué)報(bào), 2009, 30(3): 877-883. [Ding Hong-bo, Cao Yuan, Tong Wei-ping, et al. Ascent trajectory optimization and fast-reconstruction for suborbital launch vehicle [J]. Journal of Astronautics, 2009, 30(3): 877-883.]

[7] 鐘睿, 徐世杰. 基于直接配點(diǎn)法的繩系衛(wèi)星系統(tǒng)變軌控制 [J]. 航空學(xué)報(bào), 2010, 31(3): 572-578. [Zhong Rui, Xu Shi-jie. Orbit-transfer control for TSS using direct collocation method [J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(3): 572-578.]

[8] Gill P E, Murray W, Saunders M A. SNOPT: An SQP algorithm for large-scale constrained optimization [J]. Siam Journal on Optimization, 2002, 12(4): 979-1006.

[9] Biegler L T, Zavala V M. Large-scale nonlinear programming using IPOPT: An integrating framework for enterprise-wide dynamic optimization [J]. Computers & Chemical Engineering, 2009, 33(3): 575-582.

[10] Betts J T, Huffman W P. Exploiting sparsity in the direct transcription method for optimal control [J]. Computational Optimization & Applications, 1999, 14(2): 179-201.

[11] Patterson M A, Rao A V. Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems [J]. Journal of Spacecraft and Rockets, 2012, 49(2): 364-377.

[12] 童科偉, 周建平, 何麟書. 稀疏擬譜最優(yōu)控制法求解Goddard火箭問題 [J]. 固體火箭技術(shù), 2009, 32(4): 360-364. [Tong Ke-wei, Zhou Jian-ping, He Lin-shu. Sparse pseudospectral optimal control method for solving Goddard rocket problem. Journal of Solid Rocket Technology, 2009, 32(4): 360-364.]

[13] Betts J T, Huffman W P. Mesh refinement in direct transcription methods for optimal control [J]. Optimal Control Applications and Methods, 1998, 19(1): 1-21.

[14] Jain S, Tsiotras P. Trajectory optimization using multiresolution techniques [J]. Journal of Guidance Control and Dynamics, 2008, 31(5): 1424-1436.

[15] Zhao Y, Tsiotras P. Density functions for mesh refinement in numerical optimal control [J]. Journal of Guidance Control and Dynamics, 2011, 34(1): 271-277.

[16] Zhao J, Li S. Modified multiresolution technique for mesh refinement in numerical optimal control [J]. Journal of Guidance, Control, and Dynamics, 2017, Online: https://arc.aiaa.org/doi/10.2514/1.G002796.

[17] 陳小慶, 侯中喜, 劉建霞. 基于多分辨率技術(shù)的滑翔飛行器軌跡優(yōu)化算法 [J]. 宇航學(xué)報(bào), 2010, 08): 1944-1950. [Chen Xiao-qing, Hou Zhong-xi, Liu Jian-xia. A multiresolution technique-based gliding vehicle trajectory optimization algorithm [J]. Journal of Astronautics, 2010, 31(8): 1944-195.]

[18] 趙吉松, 谷良賢, 佘文學(xué). 配點(diǎn)法和網(wǎng)格細(xì)化技術(shù)用于非光滑軌跡優(yōu)化 [J]. 宇航學(xué)報(bào), 2013, 34(11): 1442-1450. [Zhao Ji-song, Gu Liang-xian, She Wen-xue. Application of local collocation method and mesh refinement to nonsmooth trajectory optimization [J]. Journal of Astronautics, 2013, 34(11): 1442-1450. ]

[19] 張松, 侯明善. 軌跡優(yōu)化的LASSO網(wǎng)格自適應(yīng)加密方法 [J]. 系統(tǒng)工程與電子技術(shù), 2016, 38(5): 1195-1200. [Zhang Song, Hou Ming-shan. LASSO-based node adaptive refinement in trajectory optimization [J]. System Engineering and Electronics, 2016, 38(5): 1195-1200.]

[20] Csendes T. Developments in Reliable Computing [M]. Dordrecht: Kluwer Academic Publishers, 1999, 77-104.

[21] Squire W, Trapp G. Using complex variables to estimate derivatives of real functions [J]. SIAM Review, 1998, 40(1): 110-112.

[22] Bryson J A E, Desai M N, Hoffman W C. Energy-state approximation in performance optimization of supersonic aircraft [J]. Journal of Aircraft, 1969, 6(6): 481-488.

ExploitingSparsityinLocalCollocationMethodsforSolvingTrajectoryOptimizationProblems

ZHAO Ji-song

(College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

In a direct local collocation method, a trajectory optimization problem is transcribed into a nonlinear programming (NLP) problem. Solving this NLP as efficiently as possible requires that the sparsity of the NLP derivatives should be exploited and the derivatives should be efficiently calculated. In this paper, a computational efficient method is developed for computing the first derivatives of the NLP functions arising from a local discretization of a trajectory optimization problem. Specifically, the expressions are derived for the NLP objective function gradient and constraint Jacobian. It is shown that the NLP derivatives can be reduced to the first derivatives of the functions in the trajectory optimization problem. As a result, the method derived in this paper reduces significantly the amount of computation required to obtain the first-derivatives required by a NLP solver. The approach derived in this paper is demonstrated by an example with discrete aerodynamic data and thrust data where it is forund that the time required to solve the NLP is reduced to less than 4% compared with the direct differentiation of the NLP functions using a finite difference method, and the efficiency improvement is more significant as the number of the grid points increases.

Trajectory optimization; Local collocation method; NLP; First derivatives; Sparsity

2017- 07- 06;

2017- 10- 10

國(guó)家自然科學(xué)基金(11602107)；中國(guó)博士后科學(xué)基金(一等資助，168884)

V412

A

1000-1328(2017)12- 1263- 10

10.3873/j.issn.1000- 1328.2017.12.002

趙吉松(1984-)，男，博士，南京航空航天大學(xué)航天學(xué)院講師，主要從事飛行器總體設(shè)計(jì)與軌跡優(yōu)化等方面的研究。

通信地址：南京市秦淮區(qū)御道街29號(hào)(210016)

電話：18260412336

E-mail: zhaojisongjinling@163.com