趙吉松

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

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

趙吉松

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

直接配點法通過對控制變量和狀態(tài)變量都進(jìn)行離散將軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃(NLP)問題。為了提高NLP的求解效率，需要利用其偏導(dǎo)數(shù)的稀疏特性并建立偏導(dǎo)數(shù)的高效計算方法。本文研究了局部配點法離散得到的NLP的一階偏導(dǎo)數(shù)的稀疏特性，建立了一階偏導(dǎo)數(shù)的高效計算方法。推導(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ù)的計算量。含有離散氣動力和推力數(shù)據(jù)的仿真算例驗證了本文方法的有效性。仿真結(jié)果表明，與有限差分法直接計算NLP的偏導(dǎo)數(shù)相比，本文方法能夠?qū)?yōu)化耗時減小至4%以內(nèi)，隨著離散節(jié)點數(shù)目的增加，計算效率的提升更為顯著。

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

0 引 言

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

雖然配點法具有多方面優(yōu)勢，但是其具體實施方法對于提高優(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ù)的計算量通常比較大，甚至超過優(yōu)化算法本身。因此，提高NLP的一階/二階偏導(dǎo)數(shù)的計算效率對于提高軌跡優(yōu)化效率有重要意義。

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

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

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

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)

端點條件為

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格式的配點法離散

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

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

(5)

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

記xi=x(τi),ui=u(τi), 對于狀態(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)時，離散格式為顯式格式，否則為隱式格式。采用類似的方法，可將目標(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，簡記HS格式)，以及經(jīng)典四階Runge-Kutta格式(q= 4，簡記RK格式)。

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

(13)

(14)

約束條件為

(15)

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

(16)

(17)

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

(18)

其中

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

Δt=tf-t0>0

(19)

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

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

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

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

(20)

其中G1的每一項仍為矩陣，以?f/?xT為例，

(21)

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

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

(22)

記

S1=struct(G1)

(23)

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

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

(24)

S2=struct(G2)

(25)

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

3.2 變量記法

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

以狀態(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ù)對優(yōu)化變量的偏導(dǎo)數(shù)，具體定義如下

(30)

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

(31)

(32)

對角陣h=diag(h0,h1,…,hN-1)，其中hi(i=0, 1, …,N-1)為積分步長，矩陣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的約束對優(yōu)化變量的偏導(dǎo)數(shù)矩陣，對于HS格式，形式如下

(34)

式中向量F和z的定義參見方程(29)。雅克比矩陣GF的展開形式遵循向量求偏導(dǎo)數(shù)運算規(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端點約束的偏導(dǎo)數(shù)

(39)

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

(40)

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

(41)

4 仿真算例

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

(42)

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

式中L為升力，D為阻力，CL為升力系數(shù)，CD為阻力系數(shù)，S為參考面積，ρ為大氣密度。氣動力的相關(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 氣動力數(shù)據(jù)Table 2 Aerodynamic data

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

J=tf

(43)

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

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

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

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

5 結(jié) 論

本文以HS格式為例，研究了局部配點離散得到的NLP問題的稀疏特性，推導(dǎo)出NLP一階偏導(dǎo)數(shù)的高效計算方法。首先引入一種新的變量記法將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ù)的計算量。采用帶有離散推力數(shù)據(jù)和氣動力數(shù)據(jù)的仿真算例驗證了本文方法的有效性。仿真結(jié)果表明，與采用有限差分法直接計算NLP的偏導(dǎo)數(shù)相比，采用本文方法能夠大幅度減小優(yōu)化耗時(對于算例，減小至4%以下)，并且隨著離散節(jié)點數(shù)量的增加，本文方法計算效率的提升更為顯著。雖然本文的研究工作基于HS格式，但是所給出的方法可容易推廣至局部配點法的其它離散格式，比如梯形格式，經(jīng)典四階RK格式等。

[1] 唐國金, 羅亞中, 雍恩米. 航天器軌跡優(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] 雍恩米, 唐國金, 陳磊. 基于Gauss偽譜方法的高超聲速飛行器再入軌跡快速優(yōu)化 [J]. 宇航學(xué)報, 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é)報, 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] 鐘睿, 徐世杰. 基于直接配點法的繩系衛(wèi)星系統(tǒng)變軌控制 [J]. 航空學(xué)報, 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é)報, 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é). 配點法和網(wǎng)格細(xì)化技術(shù)用于非光滑軌跡優(yōu)化 [J]. 宇航學(xué)報, 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

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

V412

A

1000-1328(2017)12- 1263- 10

10.3873/j.issn.1000- 1328.2017.12.002

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

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

電話：18260412336

E-mail: zhaojisongjinling@163.com