分析插值誤差和斜率的軌跡優化網格細化方法

趙吉松，尚 騰

(1. 南京航空航天大學航天學院，南京 210016；2. 北京航天自動控制研究所，北京 100854)

0 引 言

軌跡優化對于飛行器設計有著重要意義和工程實際價值。軌跡優化屬于最優控制問題，其求解方法主要分為間接法和直接法[1]。間接法借助變分法或者龐特里亞金最大值原理，把泛函優化轉化為邊值問題求解；直接法通過對控制變量/狀態變量進行離散把泛函優化轉化為非線性規劃(Non-linear programming，NLP)求解。直接法中的配點法由于不需要推導最優性必要條件，并且對優化初值的敏感性較低，容易收斂，近年來得到廣泛研究。

如果需要高精度求解軌跡優化問題，直接增加節點數量不是一種較好的方案。因為這樣處理會導致計算量顯著增加，需要更多的計算資源，并且隨著節點數目的增加，算法的收斂性通常變差。這種情況下有必要引入網格細化技術對節點進行局部加密或者調整。目前，國內外研究者在網格細化方面已經做了不少工作[2-17]。對于局部配點法，由于離散格式對于節點分布沒有限制，因而局部配點法的節點可以根據需要任意布置。用于局部配點法的網格細化方法主要有離散誤差分析法[2]、小波分析法[3-4]、多分辨率技術(MRT)[5-10]、密度函數法(DENMRA)[11]等。對于全局配點法(又稱偽譜法)，其離散節點是正交多項式的根，在離散區間的分布特點是中間稀兩端密，但是一旦節點數目給定，節點的分布情況也隨之確定。因此，全局配點法的離散節點不可以隨意布置，在網格細化方面沒有局部配點法靈活。全局配點法的網格細化方法通常是通過添加結點(段與段之間的連接點稱為結點)將問題分段優化，利用結點附近的離散節點較密的特點達到加密網格的效果，比如結點偽譜法[12]，hp自適應偽譜法[13-14]，ph偽譜法[15]，以及最新提出的既能添加節點又能減少節點的自適應網格細化算法[16]。分段優化面臨的主要問題是通常很難知道需要分多少段以及在哪分段，若通過優化迭代確定這些信息往往會存在分段過多、采用的離散節點較多等弊端[12-16]。

在眾多網格細化算法中，多分辨率技術[5-10]非常有吸引力。該技術分析各個節點處的插值誤差，如果某個節點的插值誤差較大，則在該節點附近細化網格，否則不進行細化。多分辨率技術原理簡單、魯棒性好，采用的節點較少。但是，文獻[9]發現，原始多分辨率技術[5-6]可能會出現細化遺漏或者細化失敗的問題。另外，多分辨率技術基于二進網格，其初始網格的節點數必須是2j+1(j為非負整數)。文獻[9]通過引入廣義二分網格并增加檢測算法解決了這些問題，但是給出的細化算法較為復雜并且需要較多的網格迭代次數。文獻[10]給出一種改進多分辨率技術(IMRT)，細化效率較高，但是仍然需要檢測算法才能避免細化遺漏或者失敗。文獻[5-10]本質上都屬于多分辨體系。多分辨率技術發生細化不充分的根本原因是由于這種方法在細化過程中將離散節點分成不同的分辨率層次造成的。此外，這種多層次結構的分析方式沒有充分利用現有網格中的信息，通常會添加一些不必要的節點。

本文通過直接分析離散節點的插值誤差和斜率提出一種新的網格細化方法，能夠以較少的網格迭代次數和節點數目高精度求解軌跡優化問題。

1 軌跡優化問題數學描述

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

(1)

式中:端點項M:Rm×R×Rn×Rm→R，積分項L:Rn×Rm×R→R,x∈Rn，t∈[t0,tf]?R。

狀態方程為

(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，C:Rn×Rm×R→Rc。方程(1)～(4)描述的問題稱為Bolza型最優控制問題。

2 基于Runge-Kutta格式的配點法離散

首先利用積分變換τ= (t-t0)/(tf-t0)將最優控制問題(方程(1)～(4))變換至時間區間τ∈[0, 1]。假設單位區間[0, 1]上的N個離散區間的節點為

G={τi:τi∈[0,1],i=0,1,…,N;τ0=0,

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

(5)

式中：τi稱為節點或網格點，τi在[0, 1]上可以均勻分布，也可以非均勻分布。

記xi=x(τi),ui=u(τi), 對于狀態方程(2)，采用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均為已知常數并且滿足0≤ρ1≤…≤ρq≤1。當αjl=0 (l≥j)時，方程(6)為顯式格式；否則，方程(6)為隱式格式。采用類似的方法，可將目標函數離散化。常用的離散格式包括梯形格式(q=2)，Hermite-Simpson格式(q=3)，以及經典四階Runge-Kutta格式(q=4)。

(8)

并且滿足如下約束

(9)

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

(10)

(11)

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

(12)

式中:

3 網格細化算法

由于狀態方程決定了狀態變量隨控制變量的變化規律，狀態變量的非光滑特性通常在控制變量中有所反應(對應于控制變量的劇烈變化甚至不連續)，因而本文只對控制變量進行細化。當然，如果有必要，也可以同時對狀態變量和控制變量進行細化。對于一組離散節點G={τi:i=0,1,…,N}和相應的離散最優控制變量φ={ui:i=0,1,…,N}，網格細化是根據軌跡的非光滑特性在某些區域插入新的離散節點或者刪除多余的離散節點，從而能夠采用較少的離散節點高精度求解軌跡優化問題。

3.1 節點插入方法

(13)

這種細化方法使得節點τM比τi-1和τi高一個分辨率(步長減半)，節點τQ比τM和τi高一個分辨率。這種分辨率層次關系同樣適用于節點τM ′和τQ ′。

對于每一個節點進行上述處理，便實現對網格的細化。如果優化問題有多個控制變量，那么在某個節點處只要有任意一個控制變量分量的插值誤差超過了預定的誤差限，就需要在該節點周圍細化網格。與多分辨率技術[5-6]相比，上述節點插入方法具有如下優勢。首先，這種方法非常簡單、易操作，并且對初始網格的節點數目沒有限制。其次，這種方法能夠在軌跡的任意插值誤差較大的區域持續細化網格，直到滿足精度要求或者達到預定的最高分辨率，避免了細化遺漏或者細化失敗，也避免了文獻[9-10]采用的較為復雜的檢測算法。最后，這種方法在計算某個節點的插值誤差時從其周圍的所有節點中選取最優插值節點，而文獻[5-10]只從其周圍的同等分辨率節點和更低分辨率節點中選取最優插值節點。可見，本文方法利用了當前網格中的更多信息，因而需要插入的節點數量可能會更少。

3.2 節點刪除方法

前述節點插入方法可能會添加一些多余的離散節點，若能夠將其刪除，顯然有利于減小優化計算量。本文只刪除左右斜率都為零的節點。以圖 1(a)所示的節點τi為例，其左、右斜率分別定義為

(14)

由于上述節點刪除方法只能刪除左右斜率均為零的節點，而實際上多余節點的斜率不一定為零，因而上述方法是一種比較保守的方法。未來研究中可以設計更加高效的節點刪除方法。

3.3 網格細化算法

基于前述節點插入算法和節點刪除算法，可設計出軌跡優化網格細化算法。在進行網格細化之前，需要選取初始離散區間數目N，網格細化誤差限ε，最高分辨率Jmax，最大迭代次數Imax(Imax≥Jmax/2+1)。采用RK離散方法將軌跡優化問題轉化NLP。設置I= 1，Gold=G0，提供NLP優化變量的初值并記為Xold。其中G0為區間[0, 1]上均勻分布的N+1個離散節點。那么，網格細化算法流程如下：

1)在網格Gold上以Xold為初值求解NLP。若I≥Imax，那么結束細化；否則轉入下一步。

2)節點插入算法。

3)節點刪除算法。

4)本次細化得到的非均勻網格Gnew=Gint，相應的函數值Φnew=Φint。

5)如果細化后的新網格Gnew與舊網格Gold相同，那么結束細化；否則，置I=I+1，將步驟1)解出的NLP最優解插值到Gnew，此作為新的初值Xold，令Gold=Gnew，返回步驟1)，繼續細化。

在上述網格細化技術的流程中，節點插入算法和節點刪除算法是關鍵，下文分別給出。

節點插入算法：

1)令Φold={uk:τk∈Gold,k=1,…,Nold}，將其改寫為Φold={φl(τk):l=1,…,m,τk∈Gold}。

2)初始化網格Gint=Gold以及函數值Φint=Φold。

3)設置i=2，執行如下步驟3(a)～3(c)。

(c)將i增加1。若i≤Nold-1，返回步驟3(a)；否則轉入下一步。

4)設置i= 3，執行如下步驟(4(a)和4(b))。

(a)若節點τi-1∈Gold和τi+1∈Gold的周圍都已經進行了細化，那么在節點τi周圍進行補充細化，即將以下節點添加到Gint中

其中,τl和τr分別為節點τi∈Gold在Gint中的左側相鄰節點和右側相鄰節點。

(b)將i增加1。若i≤Nold-2，返回步驟4(a)；否則轉入下一步。

5)網格Gint為節點插入算法得到的新網格，相應的函數值為Φint。若新增節點的函數值未知，則利用Gold和Φold通過p階ENO插值計算。

節點刪除算法：

1)找出Gint中達到當前最高分辨率的節點，即距離小于1/(22IN)的相鄰節點，記為Gmin。

2)用Gcheck表示屬于Gint但是不屬于G0也不屬于Gmin的節點的集合，即

Gcheck={τki:τki∈Gint(Gbase∪Gmin),i=1,…,Ncheck}

集合Gcheck表示需要分析是否可以刪除的節點。

3)設置i= 1，執行如下步驟3(a)～3(c)。

(c)將i增加1。若i≤Ncheck，返回步驟3(a)；否則，結束節點刪除算法。

4 優化算例

本節應用第3節描述的網格細化方法求解優化算例以驗證其有效性和特色，包括Bryson-Denham問題[6]、月球著落問題[13]和航天器姿態調整問題[18]。對于每個算例，采用Hermite-Simpson格式將其轉化為NLP，采用SNOPT 7[19]求解，采用INTLAB 5.5[20]提供一階偏導數。計算平臺為MacBook Air(處理器Intel Core i5-5250U 1.6 GHz，操作系統Windows 10企業版，編程語言MATLAB R2009a)。算例中的優化耗時為10次運行的平均值。另外，在每個算例中都給出了其它方法的優化結果作為對比。

4.1 算例1：Bryson-Denham問題

Bryson-Denham問題[5]又稱最小能量控制問題，是非光滑軌跡優化的常用算例。該問題具有解析最優解，能夠檢驗本文方法的精度。該問題描述為求解最優控制u(t)，使得如下目標函數最小

(15)

并且滿足動力學方程

(16)

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

x(t)≤l

(17)

式中：l為實數,本算例取l= 0.04。

采用本文的網格細化方法求解該問題，取N= 8，Jmax= 8，ε=10-3。網格細化算法迭代5次中止，從2049個均勻節點中選用43個節點高精度逼近解析解，總耗時1.13 s。圖 2給出優化的控制變量，為了對比，圖中給出了解析解。本文優化的目標函數為11.11111101，與解析最優解(100/9)的前6位小數相同。圖 3給出網格細化的迭代過程，其中豎直虛線為理論切換點位置。可見，在網格細化迭代過程中，離散點在切換點附近逐漸加密。

需要說明的是，表1中的hp方法2的耗時與本文方法也比較接近，盡管其采用的節點數目將近是本文方法的2倍。這是因為優化耗時除了與網格迭代次數和節點數目相關，還與NLP的稀疏特性以及偏導數的計算方法等因素有關，hp方法[13-14]集成了這方面的研究成果[21]。本文方法若應用這方面的研究成果，其優化效率也可以進一步提高。

4.2 算例2：月球著落問題

本算例研究月球軟著陸最優控制問題[13]。該問題描述如下：最小化目標函數

(18)

并且滿足狀態方程

(19)

邊界條件

[h(0),v(0),h(tf),v(tf)]=(10,-2,0,0)

(20)

控制變量約束

umin≤u≤umax

(21)

式中:umin= 0，umax= 3，g= 1.5，tf自由。該問題同樣具有解析最優解，具體參見文獻[13]。

采用本文方法求解該問題，取參數N= 8，Jmax= 10，ε=10-3。網格細化算法迭代6次中止，從8193個均勻節點中選用26個節點逼近解析解，總耗時0.15 s。圖 4給出優化的控制變量。為了對比，圖中給出了解析解。由圖4中的局部放大圖可知，本文方法準確捕捉住控制變量切換點位置。圖 5給出了網格細化過程，其中豎直虛線為解析切換點的位置。可見，隨著網格迭代細化，位于平坦區域的多余節點被及時刪除，新增節點逐漸向切換點聚集。

為了對比，表2統計了文獻中不同網格細化方法求解該問題的結果。其中，hp方法和h方法的結果與參數設置有關，表中只列出了最佳設置的結果。可見，與hp方法或h方法相比，本文方法的優化精度更高，并且網格迭代次數和節點數目顯著減少。因而本文方法的優化耗時更短，盡管本文沒有利用NLP的稀疏性并且采用的計算機主頻更低(文獻[13-14]采用的計算機主頻分別為2.0 GHz和2.5 GHz)。另外，表中還給出了本文采用IMRT[10]求解該算例的結果。可見，與IMRT相比，本文的細化方法將節點數目減小31.2%，從而進一步減小優化耗時。該算例再次驗證了本文方法的優勢。

方法狀態變量最大誤差迭代次數節點數目耗時/shp方法1[13]1.82×10-617453.71hp方法2[14]6.97×10-411860.29h方法[14]7.20×10-48720.20IMRT[10]9.10×10-86380.19本文方法3.22×10-86260.15

4.3 航天器姿態調整問題

本算例研究航天器最短時間姿態調整問題。假設航天器為剛體，采用四元數法描述的航天器姿態動力學方程組(即狀態方程)如下[18]

(22)

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

路徑約束

(23)

控制變量約束

(24)

邊界約束(初始條件和終端約束)

(25)

式中：φ為歐拉特征軸旋轉角，φ=150°。

目標函數

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

(26)

采用本文的網格細化方法求解該問題，取參數N= 20，Jmax= 8，ε= 10-3。本算例中特意選取N使其不等于2j(j為非負整數)。網格細化算法迭代5次中止，從5121個均勻節點中選用102個節點求解最優控制，最優目標函數為28.630387 s，優化耗時2.79 s。圖6給出優化的控制變量和相應的節點分布。由圖6(b)可知(注意圖中的局部放大圖)，本文方法在每個切換點附近都達到了預期的最高分辨率(節點間距越小，分辨率越高)，避免了細化遺漏或者細化失敗，從而準確捕捉住控制變量的切換點位置(圖中豎直虛線所示)。圖 7給出了網格迭代過程，圖中豎直虛線為最優控制變量的切換點位置。

為了對比，表3給出本文方法和文獻方法[3,10]求解該問題的結果。其中，二代小波算法[3]從2049個均勻節點中選取187個節點，IMRT[10]從2561個均勻節點中選取121個節點。顯然，本文方法采用的節點數目最少但是分辨率最高(用于選取離散節點的均勻節點基數越大，分辨率越高)。研究表明[10]，提高分辨率能夠減小終端約束的誤差。與二代小波算法相比，本文方法的網格迭代次數和節點數都顯著減少，并且對初始網格的節點數沒有限制(二代小波算法要求初始網格的節點數為2j+1)。與IMRT相比，本文的網格細化方法更簡單，并且將網格節點數目減少約16%，從而減少優化耗時。

從上述三個算例可以看出，本文的網格細化方法能夠采用較少的離散節點高精度求解軌跡優化問題。本文方法原理簡單、易操作，特別是提出的節點插入算法不需要采用多分辨率技術的多層次結構就能夠持續細化網格直至滿足精度或者達到最高分辨率，避免了細化遺漏或者失敗，并且由于在計算插值誤差時充分利用了網格信息，使得插入的節點數量較少，減小了優化計算量。此外，該方法對初始網格的節點數目沒有限制。本文方法的不足之處是目前采用的節點刪除算法只能刪除斜率為零的節點。對于最優控制變量含有常值分量的軌跡優化問題(比如Bang-Bang控制問題，本文中的算例都屬于這一類)，本文方法的優化效率較高。與之對比，多分辨率技術及其改進形式[5-10]更適合求解最優解不含常值分量的軌跡優化問題。對于最優解不含常值分量的問題，雖然本文方法也能夠求解，但是由于不能刪除多余節點，可能會導致節點數量較多。當然，這種缺陷是由于節點刪除算法的局限性造成的，未來若能設計出高效的節點刪除算法，那么將有效擴展本文的節點插入算法的適用性。

方法迭代次數節點數目耗時/s目標函數終端約束最大誤差二代小波[3]8187—28.6304—IMRT[10]51214.328.630 4032.8×10-5本文方法51022.7928.630 3871.2×10-5

5 結 論

本文提出了一種基于插值誤差和斜率分析的軌跡優化自適應網格細化方法。該方法采用局部配點法將軌跡優化問題轉化為NLP，通過分析每個離散節點處的控制變量插值誤差判斷是否需要在該節點附近加密網格，通過計算每個離散節點處的控制變量斜率判斷是否可以刪除相應節點。采用三個軌跡優化算例驗證了方法的有效性和特色。仿真結果表明，與文獻中的網格細化方法相比，本文方法生成的網格規模更小，甚至需要更少的網格迭代次數，能夠快速、高精度求解最優控制變量含有常值分量的軌跡優化問題。本文方法的不足之處是采用的節點刪除算法只能刪除斜率為零的節點。未來研究中需要設計更加通用高效的節點刪除算法。