一種求解軌跡優(yōu)化問(wèn)題的改進(jìn)多分辨率技術(shù)

趙吉松，李 爽

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

0 引 言

軌跡優(yōu)化對(duì)于飛行器設(shè)計(jì)有著重要意義和工程實(shí)際價(jià)值。軌跡優(yōu)化屬于最優(yōu)控制問(wèn)題，其求解方法主要分為間接法和直接法[1]。間接法借助變分法或龐特里亞金最大值原理，把泛函優(yōu)化轉(zhuǎn)化為邊值問(wèn)題求解；直接法通過(guò)對(duì)控制變量和/或狀態(tài)變量進(jìn)行離散把泛函優(yōu)化轉(zhuǎn)化為非線性規(guī)劃(Nonlinear programming，NLP)求解。直接法中的配點(diǎn)法由于不需要推導(dǎo)最優(yōu)性必要條件，并且對(duì)初值的敏感性較低，容易收斂，近年來(lái)得到廣泛研究。

如果需要高精度求解軌跡優(yōu)化問(wèn)題，均勻加密不是一種較好的方案。因?yàn)檫@樣處理會(huì)導(dǎo)致計(jì)算量顯著增加，需要更多的計(jì)算資源，并且隨著節(jié)點(diǎn)數(shù)目的增加，算法的收斂性通常變差。這種情況下有必要引入網(wǎng)格細(xì)化技術(shù)對(duì)離散節(jié)點(diǎn)進(jìn)行局部加密或者調(diào)整，在軌跡變化平緩區(qū)域采用稀疏網(wǎng)格，在軌跡變化劇烈的區(qū)域采用稠密網(wǎng)格。目前國(guó)內(nèi)外研究者在網(wǎng)格細(xì)化方面已做了不少工作。對(duì)于局部配點(diǎn)法[2]，由于離散格式本身對(duì)于節(jié)點(diǎn)分布沒(méi)有限制，因而局部配點(diǎn)法的離散節(jié)點(diǎn)可以根據(jù)需要任意布置。用于局部配點(diǎn)法的節(jié)點(diǎn)細(xì)化方法有離散誤差分析[3]、小波分析[4-5]、多分辨率技術(shù)(Multiresolution techniques, MRT)[6-9]、密度函數(shù)法[10]等。對(duì)于全局配點(diǎn)法(又稱(chēng)偽光譜法)[11]，其離散節(jié)點(diǎn)是正交多項(xiàng)式的根，分布特點(diǎn)是中間稀兩端密，但是一旦節(jié)點(diǎn)數(shù)目給定，節(jié)點(diǎn)的分布情況也隨之確定。因此，全局配點(diǎn)法的離散節(jié)點(diǎn)不可隨意布置，其網(wǎng)格細(xì)化沒(méi)有局部配點(diǎn)法靈活。全局配點(diǎn)法的網(wǎng)格細(xì)化方法通常是通過(guò)添加結(jié)點(diǎn)(段與段之間的連接點(diǎn)稱(chēng)為結(jié)點(diǎn))將問(wèn)題分段優(yōu)化，利用結(jié)點(diǎn)附近的離散節(jié)點(diǎn)較密的特點(diǎn)達(dá)到加密網(wǎng)格的效果，比如結(jié)點(diǎn)偽譜法[12]，hp自適應(yīng)偽譜法[13]，ph偽譜法[14]，以及最新提出的既能添加又能減少離散節(jié)點(diǎn)的自適應(yīng)網(wǎng)格細(xì)化算法[15]。分段優(yōu)化面臨的主要問(wèn)題是在解出最優(yōu)控制問(wèn)題之前通常很難知道需要分多少段以及在哪分段，若通過(guò)優(yōu)化確定這些信息往往會(huì)存在分段過(guò)多、采用的離散節(jié)點(diǎn)較多等弊端[13-15]。

在各種網(wǎng)格細(xì)化算法中，多分辨率技術(shù)[6-7]比較有吸引力。該技術(shù)分析各個(gè)節(jié)點(diǎn)處的插值誤差，如果某個(gè)節(jié)點(diǎn)處的插值誤差較大，則在該節(jié)點(diǎn)附近進(jìn)行細(xì)化網(wǎng)格，否則不進(jìn)行細(xì)化。多分辨率技術(shù)具有原理簡(jiǎn)單、使用靈活、魯棒性好、采用的離散節(jié)點(diǎn)較少等多方面優(yōu)勢(shì)。但是，文獻(xiàn)[6-7]給出的多分辨率技術(shù)沒(méi)有對(duì)初始網(wǎng)格中的偶數(shù)節(jié)點(diǎn)進(jìn)行細(xì)化，由此可能會(huì)導(dǎo)致細(xì)化遺漏，并且在細(xì)化過(guò)程中由于非光滑區(qū)域的移動(dòng)可能會(huì)發(fā)生無(wú)法繼續(xù)細(xì)化的情況，即細(xì)化失敗。此外，文獻(xiàn)[6-7]給出的多分辨率技術(shù)基于二進(jìn)網(wǎng)格，其初始網(wǎng)格節(jié)點(diǎn)數(shù)目必須是2j+1(j為非負(fù)整數(shù))，使用不方便。文獻(xiàn)[9,16]將其推廣至廣義二分網(wǎng)格，克服了這一限制，但是優(yōu)化采用的初始網(wǎng)格的節(jié)點(diǎn)數(shù)目必須是基礎(chǔ)網(wǎng)格的二倍，這種不一致給使用帶來(lái)不便。文獻(xiàn)[16]通過(guò)引入節(jié)點(diǎn)檢測(cè)算法避免了細(xì)化遺漏或失敗問(wèn)題，但是該網(wǎng)格細(xì)化方法需要的迭代次數(shù)較多，細(xì)化效率不高。

針對(duì)多分辨率技術(shù)的不足之處，本文提出一種改進(jìn)多分辨率方法，既能夠避免多分辨率技術(shù)的細(xì)化遺漏問(wèn)題，又具有較高的細(xì)化效率，并且初始網(wǎng)格節(jié)點(diǎn)數(shù)可以取任意奇數(shù)而不限于2j+1。

1 最優(yōu)控制問(wèn)題數(shù)學(xué)描述

以Bolza型最優(yōu)控制問(wèn)題為例，可描述為：求解控制變量u(t)∈Rm，使得如下目標(biāo)函數(shù)最小化

(1)

式中：端點(diǎn)項(xiàng)M:Rm×R×Rn×Rm→R, 積分項(xiàng)L:Rn×Rm×R→R,x∈Rn，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,C:Rn×Rm×R→Rc。方程(1)～(4)描述的問(wèn)題稱(chēng)為Bolza型最優(yōu)控制問(wèn)題。

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

2.1 廣義二分網(wǎng)格

傳統(tǒng)二進(jìn)網(wǎng)格的初始節(jié)點(diǎn)數(shù)目必須是2j+1(j為非負(fù)整數(shù))，使用不方便。文獻(xiàn)[9,16]給出了一種廣義二分網(wǎng)格，但是用于優(yōu)化的初始網(wǎng)格節(jié)點(diǎn)數(shù)目必須是基礎(chǔ)網(wǎng)格的二倍，使用仍然不夠靈活。為此，本文重新定義一種廣義二分網(wǎng)格。

在區(qū)間[0, 1]內(nèi)，由N個(gè)均勻分布的初始子區(qū)間對(duì)分得到的二分網(wǎng)格為

Vj,N={τj,k∈[0,1]:τj,k=k/(2jN),
0≤k≤2jN}, -1≤j≤Jmax

(5)

式中：整數(shù)j表示網(wǎng)格分辨率，整數(shù)k表示節(jié)點(diǎn)位置，整數(shù)Jmax表示最高分辨率。由于j的最小值為-1，因而N必須為偶數(shù)。采用Wj, N表示屬于Vj+1, N但不屬于Vj, N的網(wǎng)格節(jié)點(diǎn)，即

(6)

因此，τj+1, k∈Vj+1, N滿(mǎn)足如下關(guān)系

(7)

二分網(wǎng)格的子空間Vj, N滿(mǎn)足嵌套關(guān)系，即V0, N?V1, N?…?VJmax, N，并且當(dāng)Jmax→∞時(shí)，VJmax, N=[0, 1]。子空間Wj, N滿(mǎn)足Wj, N∩Wl, N=?(j≠l)。由定義(5)和(7)可知，二分網(wǎng)格是將區(qū)間[0, 1]不斷細(xì)分獲得，并且當(dāng)k≥j時(shí)滿(mǎn)足Wk, N∩Vj, N=?。

圖 1給出一組廣義二分網(wǎng)格的示例(N=6,Jmax=5)。其中V0, N=V-1, N∪W-1, N通常作為初始網(wǎng)格。對(duì)于實(shí)際問(wèn)題，N可根據(jù)需要取任意偶數(shù)。

圖1 廣義二分網(wǎng)格的示例(N = 6, Jmax= 5)Fig.1 Generalized dyadic mesh with N = 6 and Jmax= 5

2.2 Runge-Kutta離散格式

(8)

并且滿(mǎn)足如下約束

(9)

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

(10)

(11)

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

(12)

式中：Δt=tf-t0，hi=τi+1-τi，fij=f(xij,uij,τij;t0,tf)，Lij=L(xij,uij,τij;t0,tf)，變量xij,uij和τij為區(qū)間[τi,τi+1]內(nèi)的中間變量。xij由下式給出

(13)

式中：τij=τi+hiρj,uij=u(τij)。ρj,βj,αjl均為已知常數(shù)并且滿(mǎn)足0≤ρ1≤ρ2≤…≤ρq≤1。當(dāng)αjl=0(l≤j)時(shí)，離散格式為顯式格式，否則為隱式格式。

變量q為離散格式的階數(shù)。常用的格式包括梯形格式(q= 2)，Hermite-Simpson格式(q= 3)，以及經(jīng)典四階Runge-Kutta格式(q= 4)。

3 改進(jìn)多分辨率技術(shù)

多分辨率技術(shù)的核心思想是根據(jù)插值誤差對(duì)網(wǎng)格進(jìn)行細(xì)化。對(duì)于一組離散最優(yōu)解，該算法在每個(gè)節(jié)點(diǎn)處通過(guò)其周?chē)?jié)點(diǎn)插值計(jì)算該節(jié)點(diǎn)的函數(shù)值。如果插值結(jié)果與離散最優(yōu)解差別較大，則在該節(jié)點(diǎn)附近細(xì)化網(wǎng)格，否則不進(jìn)行細(xì)化。原始多分辨率技術(shù)[6-7]沒(méi)有對(duì)初始網(wǎng)格的偶數(shù)節(jié)點(diǎn)進(jìn)行細(xì)化，可能存在細(xì)化遺漏，并且在細(xì)化過(guò)程中可能發(fā)生不能繼續(xù)細(xì)化的情況。文獻(xiàn)[16]對(duì)多分辨率技術(shù)進(jìn)行改進(jìn)，避免了細(xì)化遺漏問(wèn)題，但是細(xì)化效率較低。本節(jié)給出一種新的改進(jìn)多分辨率技術(shù)，既能避免細(xì)化遺漏或者失敗問(wèn)題，又具有較高的細(xì)化效率。

在進(jìn)行多分辨率優(yōu)化之前，結(jié)合軌跡優(yōu)化問(wèn)題的精度要求，選取初始節(jié)點(diǎn)參數(shù)N，網(wǎng)格細(xì)化誤差ε，最高分辨率Jmax，最大迭代次數(shù)Imax(Imax≥Jmax/2)。采用RK方法將軌跡優(yōu)化問(wèn)題轉(zhuǎn)化NLP。設(shè)置I= 1，Gold=V0, N，估計(jì)NLP優(yōu)化變量初值，將其記為Xold。那么，改進(jìn)多分辨率方法的流程如下：

1)以Gold為初始網(wǎng)格,Xold為初值求解NLP。

2)網(wǎng)格細(xì)化方法如下(步驟(1)～(4))：

令Φold={uj,k:τj,k∈Gold}，并將其改寫(xiě)為

Φold={φl(shuí)(τj,k):l=1,…,m;τj,k∈Gold}。

(1)初始化中間網(wǎng)格Gint=V0, N, 以及相應(yīng)的函數(shù)值Φint={φl(shuí)(τ0,k)∈Φold, 0≤k≤N;l=1, 2, …,m}, 設(shè)j=-1。那么，Gold的細(xì)化流程可以分為兩步：細(xì)化算法和檢測(cè)算法。

(2)細(xì)化算法。

(3)檢測(cè)算法。

(4)本次細(xì)化得到的非均勻網(wǎng)格為Gnew=Gint, 相應(yīng)的函數(shù)值Φnew=Φint。

3)置I′ =I+1。如果細(xì)化后網(wǎng)格保持不變，那么結(jié)束細(xì)化；否則，將步驟1)解出的NLP最優(yōu)解插值到新網(wǎng)格Gnew，并以此作為新的優(yōu)化初值Xold, 令Gold=Gnew, 返回步驟1)，繼續(xù)細(xì)化。

在上述改進(jìn)多分辨率技術(shù)的各個(gè)流程中，細(xì)化算法和檢測(cè)算法是關(guān)鍵，具體算法如下。

細(xì)化算法：

1)求Wj, N和Gold的交集

2)設(shè)置i=1，執(zhí)行如下步驟((1)～(6))。

(5)將新增加節(jié)點(diǎn)對(duì)應(yīng)的函數(shù)值添加到Φint。若新增節(jié)點(diǎn)對(duì)應(yīng)的函數(shù)值未知，則利用Gold和Φold通過(guò)p階ENO插值計(jì)算。

3)將j增加1。若j≤Jmax，返回步驟1)，否則轉(zhuǎn)入下一步。

4)結(jié)束網(wǎng)格細(xì)化算法。

檢測(cè)算法：

1)設(shè)置j為中間網(wǎng)格Gint所包含節(jié)點(diǎn)的次高分辨率，即j=min(2I,Jmax)-1。

2)求解Wj, N和Gint的交集

5)根據(jù)分辨率層次，將Gcheck分成以下兩組：

可見(jiàn)，G1是由Gcheck中的最低分辨率節(jié)點(diǎn)組成,G2由其它分辨率節(jié)點(diǎn)組成。為了細(xì)化Gcheck，首先檢查G2中的節(jié)點(diǎn)并且進(jìn)行必要的細(xì)化(步驟(1)～(7))，然后采用類(lèi)似的方法處理G1中的節(jié)點(diǎn)。

6)設(shè)置i= 1，執(zhí)行如下步驟((1)～(7))。

(6)將新增加節(jié)點(diǎn)處的函數(shù)值添加至Φint，若函數(shù)值未知，根據(jù)Gold和Φold插值得到。

7)采用類(lèi)似的方法(步驟(1)～(7))，檢查G1中的每個(gè)節(jié)點(diǎn)并進(jìn)行必要的細(xì)化。

4 優(yōu)化算例

本節(jié)應(yīng)用第3節(jié)描述的改進(jìn)多分辨率技術(shù)(簡(jiǎn)記為IMRT)求解公開(kāi)文獻(xiàn)中的優(yōu)化算例以驗(yàn)證其有效性，包括Bryson-Denham問(wèn)題[6-7]，小推力軌道轉(zhuǎn)移問(wèn)題[18]，航天器姿態(tài)調(diào)整問(wèn)題[19]。對(duì)于所有算例，采用Hermite-Simpson離散格式將其轉(zhuǎn)化為NLP，采用SNOPT 7[20]求解NLP，采用INTLAB 5.5[21]提供一階偏導(dǎo)數(shù)。采用的計(jì)算機(jī)為MacBook Air (處理器Intel Core i5-5250U 1.6 GHz，操作系統(tǒng)Windows 10企業(yè)版，編程語(yǔ)言MATLAB R2009a)。算例中給出的計(jì)算耗時(shí)為10次運(yùn)行的平均值。為了對(duì)比，每個(gè)算例中都給出了原始多分辨率技術(shù)MRT-I[6]和MRT-II[7]的優(yōu)化結(jié)果。二者的主要區(qū)別在于MRT-I在網(wǎng)格細(xì)化的每次迭代中添加一層高分辨率節(jié)點(diǎn)，而MRT-II則添加兩層高分辨率節(jié)點(diǎn)。在所有算例中，參數(shù)p= 2，N1= 1，N2= 1。

4.1 算例1：Bryson-Denham問(wèn)題

Bryson-Denham問(wèn)題[6-7]又稱(chēng)最小能量控制問(wèn)題，是典型的非光滑軌跡優(yōu)化算例。該問(wèn)題描述為求解最優(yōu)控制u(t)，使得如下目標(biāo)函數(shù)最小

(14)

并且滿(mǎn)足動(dòng)力學(xué)方程

(15)

初始條件x(0) = 0,v(0) = 1，終端條件x(1) = 0，v(1) = 1；以及如下路徑約束(狀態(tài)變量邊界約束)

x(t)≤l

(16)

式中：l為實(shí)數(shù)，本算例取l=0.04。

采用本文方法求解該問(wèn)題，取N=8，Jmax=7，ε=10-3。IMRT迭代5次停止，圖2給出優(yōu)化結(jié)果。圖2(a)為最優(yōu)控制變量以及離散點(diǎn)分布，圖中細(xì)實(shí)線為ENO插值結(jié)果，粗實(shí)線為解析解;圖2(b)狀態(tài)變量x隨時(shí)間變化曲線，圖中細(xì)實(shí)線為根據(jù)最優(yōu)控制采用數(shù)值積分得到，粗實(shí)線為解析解?？梢?jiàn)，本文方法從1025個(gè)均勻節(jié)點(diǎn)選取49個(gè)節(jié)點(diǎn)即可高精度逼近解析解。本文方法優(yōu)化的目標(biāo)函數(shù)11.11111101，與解析解(100/9)相差1.1×10-7。

對(duì)于該算例，若采用原始多分辨率技術(shù)求解(只細(xì)化控制變量的情況下)，則會(huì)發(fā)生細(xì)化不充分。圖3給出了采用MRT-I[6]優(yōu)化的結(jié)果,可以看出，MRT-I只細(xì)化至W1,N就停止細(xì)化。由圖3(a)可知，控制變量在切換點(diǎn)處的插值誤差顯然超過(guò)了誤差限(10-3)，因此MRT-I停止細(xì)化并不是因?yàn)榫冗_(dá)到要求，而是因?yàn)镸RT-I本身缺陷造成的。MRT-I優(yōu)化的目標(biāo)函數(shù)11.09858462，與解析解相差1.3×10-2。

對(duì)比圖2和圖3可以發(fā)現(xiàn)，盡管離散節(jié)點(diǎn)細(xì)化不夠充分引起的控制變量差異比較小，但是由此導(dǎo)致的狀態(tài)變量差異卻是顯著的。更為重要的是， 如圖3(b)中的局部放大圖所示，對(duì)于受約束的狀態(tài)變量x，盡管在離散點(diǎn)處均滿(mǎn)足路徑約束(式(16))，但是數(shù)值積分結(jié)果表明，在離散點(diǎn)之間x卻違反了路徑約束，x的峰值超出了路徑約束上限0.5%。

圖2 IMRT優(yōu)化的算例1最優(yōu)解(49個(gè)節(jié)點(diǎn))Fig.2 Solution for example 1 using IMRT (49 points)

圖3 MRT-I優(yōu)化的算例1最優(yōu)解(17個(gè)節(jié)點(diǎn))Fig.3 Solution for example 1 using MRT-I (17 points)

了MRT細(xì)化不充分問(wèn)題，但是由此帶來(lái)的潛在問(wèn)題包括：1)增加了算法復(fù)雜度，并且實(shí)際問(wèn)題的控制量和狀態(tài)量可能相差較大，需要分別選取誤差限；2)需要更多的離散節(jié)點(diǎn)，對(duì)比表中MRT-II和本文方法可證實(shí)這一點(diǎn)；3)這種處理方法并非一直有效，參見(jiàn)文獻(xiàn)[16]的論述和反例。此外，如前所述，文獻(xiàn)[6-7]的方法要求初始節(jié)點(diǎn)數(shù)目必須為2j+1。文獻(xiàn)[16]的方法突破了這一限制并且能夠避免細(xì)化失敗，但需要的網(wǎng)格細(xì)化迭代次數(shù)較多(見(jiàn)表1)。

表1 不同方法求解算例1的效果對(duì)比Table 1 Solutions from different methods for example 1

4.2 算例2：最短時(shí)間軌道轉(zhuǎn)移問(wèn)題

本算例研究的最短時(shí)間小推力軌道轉(zhuǎn)移問(wèn)題是由小推力能量最優(yōu)軌道轉(zhuǎn)移問(wèn)題簡(jiǎn)化而來(lái)，對(duì)于求解能量最優(yōu)問(wèn)題具有重要意義[18]。該問(wèn)題的特色之處是其最優(yōu)控制含有多個(gè)bang-bang切換結(jié)構(gòu)，能夠較好地檢驗(yàn)改進(jìn)多分辨率技術(shù)的有效性。

最短時(shí)間小推力軌道轉(zhuǎn)移問(wèn)題是通過(guò)優(yōu)化飛行器的切向加速度和法向加速度，使飛行器從初始軌道轉(zhuǎn)移至目標(biāo)軌道消耗的時(shí)間最短。對(duì)于初始軌道和目標(biāo)軌道均為圓軌道，問(wèn)題描述如下[18]：

MinJ[x(t),u(t),tf]=tfs.tr = vrθ·=vtvr=v2tr-μr2+uvt=-vrvtr+ut(r,θ,vr,vt)(t0) = (1,0,0,1)(r,vr,vt)(tf)=(4,0,0.5)

(17)

采用IMRT求解該問(wèn)題，取參數(shù)N= 30，Jmax= 6，ε=2×10-4。注意本算例中特意選取N使其不等于2j(j為非負(fù)整數(shù))。IMRT迭代4次中止，從最高分辨率節(jié)點(diǎn)V6, N的1921個(gè)均勻節(jié)點(diǎn)中選取153個(gè)離散節(jié)點(diǎn)，優(yōu)化耗時(shí)1.5 s。圖 4給出最優(yōu)控制變量以及離散節(jié)點(diǎn)分布。由圖4可知，本文方法在所有切換點(diǎn)處均達(dá)到預(yù)期的最高分辨率。本文方法優(yōu)化的最優(yōu)目標(biāo)函數(shù)為47.699776 TU(文獻(xiàn)[18]采用序列偽譜方法優(yōu)化的目標(biāo)函數(shù)為47.809 TU)。圖5給出最優(yōu)轉(zhuǎn)移軌道，其中圓圈為離散最優(yōu)解，細(xì)實(shí)線為根據(jù)圖 4(a)所示的最優(yōu)控制采用數(shù)值積分得到的結(jié)果。可見(jiàn)，飛行器雖然需要飛越多個(gè)周期才能達(dá)到目標(biāo)軌道，但是仍然能夠準(zhǔn)確進(jìn)入目標(biāo)軌道。

圖4 IMRT優(yōu)化的算例2最優(yōu)解(153個(gè)節(jié)點(diǎn))Fig.4 Solution for example 2 using IMRT (153 points)

圖5 IMRT優(yōu)化的算例2最優(yōu)轉(zhuǎn)移軌道Fig.5 Optimal transfer orbit for example 2 using IMRT

作為對(duì)比，圖6給出采用MRT-I優(yōu)化的控制變量和離散節(jié)點(diǎn)分布。從圖6(b)所示的節(jié)點(diǎn)分布可以發(fā)現(xiàn)，在控制變量的第2個(gè)和第5個(gè)切換點(diǎn)處的離散節(jié)點(diǎn)并沒(méi)有達(dá)到期望的最高分辨率。對(duì)比圖4(a)和圖 6(a)給出的控制變量局部放大插圖可以發(fā)現(xiàn)，圖6(a)中這兩個(gè)位置的節(jié)點(diǎn)分布比圖4(a)稀疏，并且二者的切換點(diǎn)位置也不完全一致，由此說(shuō)明了圖6(a)中這兩個(gè)位置的控制變量細(xì)化不夠充分。

圖6 MRT-I優(yōu)化的算例2最優(yōu)解(140個(gè)節(jié)點(diǎn))Fig.6 Solution for example 2 using MRT-I (140 points)

MRT-I優(yōu)化的最優(yōu)目標(biāo)函數(shù)為47.703654 TU，與IMRT的優(yōu)化結(jié)果相差+0.008%。因此，對(duì)于本算例，這種節(jié)點(diǎn)細(xì)化不夠充分對(duì)于目標(biāo)函數(shù)幾乎沒(méi)有影響。但是，進(jìn)一步研究表明，這種細(xì)化不充分卻對(duì)終端狀態(tài)變量的入軌誤差有著顯著影響。根據(jù)IMRT優(yōu)化的控制變量，采用四階龍格庫(kù)塔方法對(duì)動(dòng)力學(xué)系統(tǒng)進(jìn)行積分，得到終端矢徑、徑向速度、切向速度的入軌誤差分別為4.5×10-5，1.4×10-5，-3.0× 10-7；而根據(jù)MRT-I優(yōu)化的控制變量積分動(dòng)力學(xué)系統(tǒng)得到的終端入軌誤差分別為-1.1×10-3，3.7×10-5，1.2×10-4。可見(jiàn)，IMRT優(yōu)化的控制變量對(duì)應(yīng)的入軌誤差更小，即入軌精度更高。這是因?yàn)閿?shù)值積分動(dòng)力學(xué)系統(tǒng)時(shí)需要根據(jù)優(yōu)化的離散控制變量插值計(jì)算各個(gè)積分點(diǎn)處的控制變量，IMRT優(yōu)化的控制變量在所有切換點(diǎn)處都具有較高的分辨率，能夠更加準(zhǔn)確地插值出積分點(diǎn)處的控制變量，而MRT-I優(yōu)化的控制變量在部分切換點(diǎn)處的分辨率不夠高，在這些切換點(diǎn)處的插值誤差更大，插值誤差通過(guò)數(shù)值積分的累積效應(yīng)影響終端入軌精度。

為了進(jìn)一步分析本文的改進(jìn)多分辨率技術(shù)的性能，表2給出了采用文獻(xiàn)方法[6-7,16]求解該問(wèn)題的結(jié)果對(duì)比。其中，MRT-I[6]和MRT-II[7]均發(fā)生了細(xì)化遺漏，Zhao&Li[16]方法和本文方法均避免了細(xì)化遺漏。由表2可知，與前述結(jié)論類(lèi)似，是否發(fā)生細(xì)化遺漏對(duì)于該算例的目標(biāo)函數(shù)影響微小，但是對(duì)于終端約束誤差影響較大(未發(fā)生細(xì)化遺漏時(shí)終端狀態(tài)的約束誤差更小)。與MRT-I和MRT-II相比，本文方法避免了細(xì)化遺漏，終端約束誤差更小，并且需要的網(wǎng)格迭代次數(shù)更少，優(yōu)化耗時(shí)更短。與Zhao&Li方法相比，二者的終端約束誤差相當(dāng)，但是本文方法需要的網(wǎng)格迭代次數(shù)更少，因而耗時(shí)更短。

表2 不同方法求解算例2的效果對(duì)比Table 2 Solutions from different methods for example 2

4.3 算例3：航天器姿態(tài)調(diào)整問(wèn)題

本算例采用改進(jìn)多分辨率技術(shù)求解航天器的最優(yōu)姿態(tài)調(diào)整問(wèn)題。航天器為NASA的X射線計(jì)時(shí)探測(cè)器(XTE)，優(yōu)化目標(biāo)使得航天器從一個(gè)姿態(tài)調(diào)整到另一個(gè)姿態(tài)消耗的時(shí)間最短。由于控制力矩大小恒定，因此時(shí)間最短等價(jià)于能量最省。

將航天器假設(shè)為剛體，采用四元數(shù)法描述的剛體航天器姿態(tài)動(dòng)力學(xué)方程組如下[19]：

(18)

式中：q=[q1,q2,q3,q4]T為四元數(shù)，ω=[ω1,ω2,ω3]T為旋轉(zhuǎn)角速度，u=[u1,u2,u3]T為控制力矩，轉(zhuǎn)動(dòng)慣量Ix=5621 kg·m2,Iy=4547 kg·m2和Iz=2364 kg·m2。

路徑約束

(19)

控制變量約束

(20)

邊界約束

(21)

目標(biāo)函數(shù)

J[x(t),u(t),tf]=tf

(22)

采用IMRT求解該問(wèn)題，取N=20，Jmax=7，ε=0.1。IMRT迭代5次中止，從最高分辨率節(jié)點(diǎn)V7, N的2561個(gè)均勻節(jié)點(diǎn)中選取121個(gè)離散節(jié)點(diǎn)，優(yōu)化耗時(shí)4.3 s。圖 7給出最優(yōu)控制變量以及離散節(jié)點(diǎn)分布?？梢?jiàn)，本文方法在所有切換點(diǎn)處均達(dá)到預(yù)期的最高分辨率。本文方法優(yōu)化的目標(biāo)函數(shù)為28.630403 s。

圖7 IMRT優(yōu)化的算例3最優(yōu)解(121個(gè)節(jié)點(diǎn))Fig.7 Solution for example 3 using IMRT (121 points)

圖8 IMRT優(yōu)化的算例3最優(yōu)狀態(tài)變量Fig.8 State solution for example 3 using IMRT

圖 8給出了優(yōu)化的狀態(tài)變量，其中圓圈為離散最優(yōu)解，細(xì)實(shí)線為根據(jù)優(yōu)化的最優(yōu)控制變量采用四階龍格庫(kù)塔方法積分得到的結(jié)果。本算例的終端狀態(tài)變量全部受到約束，數(shù)值積分結(jié)果表明，終端狀態(tài)約束的誤差分別為-3.6×10-6，-2.5×10-5，-2.8×10-5，1.4× 10-5，-2.7×10-6，-8.6×10-6，1.4×10-6。

作為對(duì)比，圖 9給出了采用MRT-I優(yōu)化的最優(yōu)控制變量和節(jié)點(diǎn)分布，MRT-I優(yōu)化的目標(biāo)函數(shù)為28.630405 s。由圖 9(b)可知，在控制變量的第一個(gè)和最后一個(gè)切換點(diǎn)處離散節(jié)點(diǎn)沒(méi)有達(dá)到最高分辨率。對(duì)比圖 7(a)和圖 9(a)中的控制變量局部放大插圖可以發(fā)現(xiàn)，MRT-I細(xì)化的節(jié)點(diǎn)分布中存在遺漏。對(duì)于該算例，這種細(xì)化遺漏對(duì)目標(biāo)函數(shù)仍然幾乎沒(méi)有影響，但是根據(jù)MRT-I優(yōu)化的控制變量對(duì)動(dòng)力學(xué)系統(tǒng)積分，得到的終端狀態(tài)約束的誤差分別為2.1×10-5，5.8×10-4，-1.7×10-4，7.9×10-5，2.8×10-6，1.6×10-5，-1.3×10-4。可見(jiàn)，MRT-I得到的終端狀態(tài)約束最大誤差比IMRT的結(jié)果大一個(gè)量級(jí)以上。

表3給出采用本文方法和文獻(xiàn)方法[4,6-7,16]求解該問(wèn)題的結(jié)果對(duì)比。其中二代小波的結(jié)果來(lái)源于文獻(xiàn)[4]，最高分辨率為1921，其它方法的結(jié)果為本文采用文獻(xiàn)中的方法計(jì)算得出，最高分辨率均為2561。MRT-I[6]和MRT-II[7]方法發(fā)生了細(xì)化遺漏，Zhao& Li[16]方法和本文方法避免了細(xì)化遺漏，二代小波算法不確定是否發(fā)生細(xì)化遺漏。可見(jiàn)，與算例2類(lèi)似，是否發(fā)生細(xì)化遺漏對(duì)于目標(biāo)函數(shù)影響微小，但是對(duì)于終端約束誤差影響顯著。該算例再次表明，本文方法與其它多分辨方法相比需要的網(wǎng)格迭代次數(shù)和離散節(jié)點(diǎn)數(shù)目較少(優(yōu)化耗時(shí)更短)，并且能夠避免細(xì)化遺漏問(wèn)題(終端約束誤差更小)。

圖9 MRT-I優(yōu)化的算例3最優(yōu)控解(110個(gè)節(jié)點(diǎn))Fig.9 Solution for example 3 using MRT-I (110 points)

方法迭代次數(shù)節(jié)點(diǎn)數(shù)目耗時(shí)/s目標(biāo)函數(shù)終端約束最大誤差MRT-I81186.428.6304055.8×10-4MRT-II51334.728.6304011.4×10-4Zhao & Li81246.928.6304054.6×10-5二代小波8187—28.6304—本文方法51214.328.6304032.8×10-5

5 結(jié) 論

本文提出一種改進(jìn)多分辨率方法用于求解非光滑軌跡優(yōu)化問(wèn)題。針對(duì)原始多分辨率技術(shù)在網(wǎng)格細(xì)化過(guò)程中可能會(huì)發(fā)生細(xì)化遺漏或者失敗的缺陷，本文通過(guò)在網(wǎng)格細(xì)化過(guò)程中增加節(jié)點(diǎn)檢測(cè)算法確保在每個(gè)插值誤差較大的離散節(jié)點(diǎn)附近網(wǎng)格達(dá)到當(dāng)前迭代的最高分辨率。為了進(jìn)一步提高網(wǎng)格細(xì)化效率，該方法在每次細(xì)化時(shí)添加兩層高分辨率節(jié)點(diǎn)，從而減少網(wǎng)格細(xì)化迭代的次數(shù)。此外，針對(duì)原始多分辨率技術(shù)的初始網(wǎng)格節(jié)點(diǎn)數(shù)目必須是2j+1(j是非負(fù)整數(shù))的限制，定義了一種新的廣義二分網(wǎng)格，使得初始節(jié)點(diǎn)數(shù)目可以為任意奇數(shù)而不限于2j+1。采用多個(gè)非光滑軌跡優(yōu)化算例驗(yàn)證了方法的有效性和特色。結(jié)果表明，本文方法能夠在軌跡的非光滑區(qū)域持續(xù)細(xì)化網(wǎng)格，避免了細(xì)化遺漏或者失敗，并且具有較高的細(xì)化效率。此外，研究發(fā)現(xiàn)，本文方法避免細(xì)化遺漏對(duì)于目標(biāo)函數(shù)影響不大，但是卻能夠更加準(zhǔn)確地滿(mǎn)足路徑約束和終端狀態(tài)約束。

參 考 文 獻(xiàn)

[1] 雍恩米, 陳磊, 唐國(guó)金.飛行器軌跡優(yōu)化數(shù)值方法綜述[J]. 宇航學(xué)報(bào), 2008, 29(2): 397-406. [Yong En-mi, Chen Lei, Tang Guo-jin. A survey of numerical methods for trajectory optimization of spacecraft[J]. Journal of Astronautics, 2008, 29(2): 397-406. ]

[2] 趙吉松. 求解軌跡優(yōu)化問(wèn)題的局部配點(diǎn)法的稀疏性研究 [J]. 宇航學(xué)報(bào), 2017, 38(12): 1263-1272. [Zhao Ji-song. Exploiting sparsity in local collocation methods for solving trajectory optimization problems [J]. Journal of Astronautics, 2017, 38(12): 1263-1272.]

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

[4] 豐志偉, 張青斌, 唐乾剛, 等. 基于二代小波的軌跡優(yōu)化節(jié)點(diǎn)自適應(yīng)加密 [J]. 航空動(dòng)力學(xué)報(bào), 2013, 28(7): 1659-1665. [Feng Zhi-wei, Zhang Qing-bin, Tang Qian-gang, et al. Node adaptive refinement for trajectory optimization based on second-generation wavelets [J]. Journal of Aerospace Power, 2013, 28(7): 1659-1665.]

[5] Assassa F, Marquardt W. Dynamic optimization using adaptive direct multiple shooting [J]. Computers and Chemical Engine-ering, 2014, 60(2): 242-259.

[6] Jain S, Tsiotras P. Multiresolution-based direct trajectory optimization [C]. The 46th IEEE Conference on Decision and Control, Piscataway, USA, Dec 12-14, 2007.

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

[8] 陳小慶, 侯中喜, 劉建霞. 基于多分辨率技術(shù)的滑翔飛行器軌跡優(yōu)化算法 [J]. 宇航學(xué)報(bào), 2010, 31(8): 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-1950.]

[9] 趙吉松, 谷良賢, 佘文學(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. ]

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

[11] 張志國(guó), 余夢(mèng)倫, 耿光有, 等. 應(yīng)用偽譜法的運(yùn)載火箭在線制導(dǎo)方法研究 [J]. 宇航學(xué)報(bào), 2017, 38(3): 262-269. [Zhang Zhi-guo, Yu Meng-lun, Geng Guang-you, et al. Research on application of pseudo-spectral method in online guidance method for a launch vehicle [J]. Journal of Astronautics, 2017, 38(3): 262-269.]

[12] Ross I M, Fahroo F. Pseudospectral knotting methods for solving optimal control problems [J]. Journal of Guidance, Control, and Dynamics, 2004, 27(3): 397-405.

[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(3): 433-445.

[14] Patterson M A, Hager W W, Rao A V. A ph mesh refinement method for optimal control [J]. Optimal Control Applications and Methods, 2015, 36(4): 398-421.

[15] Liu F, Hager W W, Rao A V. Adaptive mesh refinement method for optimal control using nonsmoothness detection and mesh size reduction [J]. Journal of the Franklin Institute-Engineering and Applied Mathematics, 2015, 352(10): 4081-4106.

[16] Zhao J, Li S. Modified multiresolution technique for mesh refinement in numerical optimal control [J]. Journal of Guidance, Control, and Dynamics, 2017, 40(12): 3328-3338.

[17] Harten A, Engquist B, Osher S, et al. Uniformly high order accurate essentially non-oscillatory schemes III [J]. Journal of Computational Physics, 1997, 131(1): 3-47.

[18] Ross I M, Gong Q, Sekhavat P. Low-thrust, high-accuracy trajectory optimization [J]. Journal of Guidance, Control, and Dynamics, 2007, 30(4): 921-933.

[19] Fleming A, Sekhavat P, Ross I M. Minimum-time reorientation of a rigid body [J]. Journal of Guidance, Control, and Dyn-amics, 2010, 33(1): 160-170.

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

[21] Csendes T. Developments in reliable computing [M]. Dordrecht: Kluwer Academic Publishers, 1999.