UKF濾波算法在兩點邊值問題求解中的應用

藏潔，劉升剛

(1.北京空間飛行器總體設計部，北京 100094；2.北京航空航天大學 宇航學院，北京 100191)

0 引言

在最優化問題的間接法求解算法中，根據最優點的必要性條件，原最優化問題轉化為一個兩點邊值問題(two point boundary value problem,TPBVP)，兩點邊值問題的解就是原問題的最優解[1]。通常兩點邊值問題所描述的是一個有限維動力系統，在其有限的單維自變量區間的2個端點上存在與系統維數相同個數的若干約束，要求求解滿足上述約束的系統狀態的時間軌跡。

約束是分別施加于2個端點，故這2個端點的狀態都存在一定的自由度。兩點邊值問題求解的關鍵就是尋找2個端點狀態之間精確的定量關系，或者具有足夠精度的近似關系。不過除了十分簡單的情況，一般來說從原最優化問題轉化得到的兩點邊值問題所對應的系統都是比較復雜的非線性系統，自變量區間2個端點處的系統狀態之間的關系都是具有高敏度、高非線性特性的，且很難得到解析表達式。現在主流的兩點邊值問題求解算法都是數值求解算法，基本可分為2類:打靶法(shooting method)和轉錄/配置法(transcription/collocation)[2]。

打靶法的設計變量是考慮始端約束條件下，始端自由狀態變量，通過數值方法計算得到終端系統狀態變量以及終端狀態變量相對始端狀態變量的雅可比矩陣，然后根據終端狀態約束和雅可比矩陣計算始端狀態的修正量，由此迭代重復進行上述計算，直至終端約束滿足，得到兩點邊值問題的精確解[3-5]。打靶法中終端狀態和始端狀態之間定量關系是對原系統的一階線性近似，對于高敏度高非線性特性的系統模型，要求打靶迭代的狀態初始猜測值具有較高的精度，否則難以收斂到精確解。為了降低系統敏度對打靶法的影響，提高打靶法的魯棒性，在此基礎上發展了多重打靶法[6]，將系統狀態軌跡分為若干段，相鄰段之間需要滿足連續性約束條件[7]，每一段內采用打靶法求解。多重打靶法降低了模型敏度，提高了魯棒性，但本質上并沒有解決打靶法對狀態初值敏感的問題，并且引入額外的設計變量和約束。

轉錄/配置法[8-9]將自變量在取值區間上離散化，以所有離散節點所對應的系統狀態變量為設計變量，節點之間的狀態變量通過多項式插值得到。在自變量區間上選擇配置點，要求在配置點上滿足系統的動力學方程，這實際上是有限維的等式約束，由此將動力系統的動力學方程約束轉化為代數方程約束，同時考慮原問題的邊值約束，將兩點邊值問題轉化為一個非線性規劃問題[7,10]。為了保證在配置點上的導數計算精度，一般采取Hermite-Simpson插值方法[11]計算系統狀態變量的導數。與打靶法相比，打靶法在求解過程中始終保證滿足系統動力學方程，而轉錄/配置法在求解過程中不能保證滿足系統動力學方程，一般不需要進行大計算量的數值積分，但是其設計變量維數大大增加，并且同樣存在對狀態初值敏感的問題。

針對目前兩點邊值問題求解算法對初值十分敏感、難以收斂的問題，本文提出基于UKF(unscented Kalman filter)濾波算法的新的兩點邊值問題求解方法。其基本思路是將原動力系統隨機化，假設系統狀態變量都滿足高斯分布，將原兩點邊值問題轉換為一個最優參數估計問題。采用Unscented變換，可以根據系統始端狀態變量的概率分布，得到經過非線性變換后系統終端狀態變量概率分布的二階以上精度近似，由此得到系統始端狀態變量和終端狀態變量之間二階以上精度的定量關系。采用UKF濾波算法，可以遞推修正系統始端狀態，直至濾波收斂，得到原兩點邊值問題的解。由于此方法對系統的近似描述具有二階以上精度，直觀上可以確定其收斂域要遠大于上述僅僅具有一階精度的傳統求解方法，對狀態初值敏感的程度大大降低。

1 Kalman濾波和Unscented Kalman濾波算法簡介

目前，估計理論主要有三大類方法:最小二乘估計、極大似然估計和貝葉斯估計。從最小二乘法和貝葉斯估計出發，在一定假設條件下，都可以推導出Kalman濾波。從最小二乘法理論理解，Kalman濾波本質上是線性系統的線性最小方差估計；從貝葉斯估計理論理解，Kalman濾波本質上是高斯線性系統的序列Bayes估計[12-13]。通常貝葉斯估計需要相關變量完整的概率分布描述，在高斯分布的假設下，只需要變量的均值、方差和協方差即可。最小二乘估計(最小方差估計)只需要變量的均值、方差和協方差，但精確的最小方差估計求解比較困難，以線性最小方差為目標，可以采用Hilbert空間投影定理直接求解。

考慮在狀態空間中描述的一個離散系統，

(1)

(2)

若式(1)是線性系統的情況，上述遞推線性最小方差估計的實現(很容易計算隨機變量經線性變換后的均值和協方差)就是Kalman filter(KF)。對于式(1)是非線性系統的情況，產生的困難是準確計算隨機變量經非線性變換后的均值和協方差。Extended Kalman filter(EKF)是將非線性系統近似為一階線性系統，但精度較低。Unscented Kalman filter(UKF)由Julier首先提出[14]，是在Kalman filter基礎上增加了UT變換(unscented transformation)，專門處理非線性系統的最優估計問題。

UT變換可以得到隨機變量經非線性變換后二階精度以上的均值和協方差，基本思想是近似非線性函數的概率分布比近似非線性函數本身更容易，它采取一定規則的Sigma采樣點集合S來表征隨機變量的基本統計特征，然后對所有Sigma點進行非線性變換，通過加權計算非線性變換后的隨機變量的基本統計特征。這種方法不依賴于非線性函數的具體形式，關鍵在于采樣策略，即Sigma點的個數，分布和相應的權值。

(3)

(4)

將所有的Sigma點進行非線性變換G處理，

(5)

(6)

上面公式中一些常數的含義如下：

(2)常量κ≥0，以保證方差矩陣的半正定性，一般取為0或者3-n。

(3)β是與x的先驗分布相關的常量，對于高斯分布，β=2是最優的。對于非高斯分布，這個參數還可以用來控制誤差。

(7)

不考慮系統噪聲，經式(1)非線性變換，

(8)

xk經過(1)作用后的均值和方差為

(9)

(10)

(11)

不考慮觀測噪聲，經式(1)非線性變換，

(12)

根據UT變換，

(13)

(14)

而xk+1和yk+1的協方差矩陣為

(15)

對于線性最小方差估計，

(16)

(17)

KF和UKF本質上都是高斯系統的線性最小方差估計(只需要隨機變量概率密度分布的一階矩和二階矩即可計算最優估計)，而KF是線性系統的精確實現，UKF是非線性系統二階精度以上的實現。濾波過程的核心是計算相關變量的方差矩陣，只有正確計算方差矩陣才能保證濾波的收斂性和無偏性。系統狀態變量的方差分為2部分，之前時刻變量方差的傳遞和當前系統噪聲的影響。之前時刻的狀態變量方差體現了對狀態變量先驗知識的了解，系統的噪聲體現了對模型先驗知識的了解。在實際應用中，系統模型是不可能完全精確建模，為了彌補模型和真實情況的偏差，必須給定合適的系統噪聲，以體現這部分的影響，否則濾波器沒有足夠的增益對狀態變量進行更新修正，導致濾波發散，具體的體現可能就是狀態變量的方差很小，但濾波結果遠偏離實際值。適當增加系統噪聲的方差，可以增加濾波器的穩定性[15]，降低對狀態初值的敏感性，但是這會降低模型的預測精度，增加濾波結果的誤差。從另一角度來理解此結論，Kalman濾波發散是由于模型的偏差使得只能在短時間內保證預測正確性，而Kalman濾波算法是根據全部的測量預測結果來估計當前的狀態。增加系統噪聲方差，使得在濾波過程中，更早的觀測數據會更快地被舍棄掉，而給最近的觀測數據賦予更大的權值，由此可以避免濾波發散。

2 UPE算法簡介及其在兩點邊值問題求解中的應用

UPE，即UKF parameter estimation。假設存在非線性映射，

yk=G(xk,w),k=1,2,…,∞,

(18)

式中：(xk,yk)為已知的映射對序列；w為一組未知的參數，參數估計即根據上述非線性映射和映射對序列，估計未知的參數w。將非線性映射進行如下變換：

(19)

這2種方法思路都是在濾波初段設置較大的系統噪聲方差，保證濾波收斂，在濾波后段使得系統噪聲方差減小，保證濾波精度。算法流程如下：

(1)初始化

(20)

(2)狀態更新

對于k∈{1,2,…,∞}狀態更新和Sigma點計算

(21)

Yk|k-1=G(xk,Wk|k-1),

(22)

(3)觀測更新

(23)

兩點邊值問題也可以看作一個參數估計問題。存在如下動力系統，

(24)

式中：x為n維狀態變量；u為給定驅動，在時間區間t∈[t0,tf]上分析系統響應。時間區間2個端點的狀態變量記為x0=x(t0)和xf=x(tf)，當給定x0即可根據式(24)確定xf，xf可記為x0的函數xf=g(x0)。假設對x0存在k個約束，對xf存在n-k個約束，形式為

(25)

上述即為常用的兩點邊值問題完整描述。由式(25)可以消除x0的k個自由度，假設余下的n-k個自由度記為w，x0可記為w的函數x0=h(x)。式(25)的終端約束可以變換為

0=G(w)=Gf(g(h(w))).

(26)

式(26)可看作一個參數估計問題，和式(18)相比，xk隱含在w中，yk恒等于0，可以采用UPE算法得到兩點邊值問題的解。

3 算例

考慮高超聲速飛行器上升段飛行過程，動力為固體火箭發動機，采用縱向平面內動力學方程，不考慮地球自轉和弧度，忽略控制系統延遲，根據瞬時平衡假設，歸一化的質心動力學方程為

(27)

(28)

氣動力為

(29)

(30)

Cx=(a2M2+a1M+a0)+(b2M2+b1M+b0)α2，

(31)

Cy=(c2M2+c1M+c0)+(d2M2+d1M+d0)α,

(32)

式中：M為馬赫數；α為攻角；ai，bi，ci，di為常系數。

初始邊界為

(33)

終端邊界為

(34)

優化指標為

(35)

(36)

協態方程為

(37)

式(37)右端下標代表對應偏導數，具體形式略。最優控制α*使得沿最優軌跡哈密爾頓函數最小，滿足必要條件，

(38)

將式(29)～(30)代入式(38)，得

(39)

式(39)可能存在多解，需要確定使得式(36)最小的值。將式(39)再次對α求偏導

Hαα=k1sinα+k2cosα+k3=

(40)

式中：k1，k2，k3對應相應的系數。

如果式(40)無解，那么式(39)只有一個解，如果式(40)有解，將其記為

(41)

式中:β定義如下：

(42)

橫截條件為

(43)

組成了兩點邊值問題的描述。算例初末狀態見表1。

表1 初末軌道要素和計算參數Table 1 Elementary and final orbital elements and calculation parameters

圖1 考慮不同大氣影響下狀態變量時間歷程Fig.1 Considering the time history of state variables under different atmospheric influences

圖2 考慮不同大氣影響下協態變量時間歷程Fig.2 Considering the time history of co-state variables under different atmospheric influences

圖3 考慮不同大氣影響下控制變量時間歷程Fig.3 Considering the time history of control variables under different atmospheric influences

根據仿真結果，所有的邊界條件和最優必要條件都滿足。對于本文給定末端條件的顯著特點是，有大氣影響時末速最大的彈道并不是直接達到給定的15 km高度，而是盡快達到較高的高度，然后下滑拉平。從物理意義來看，在具有足夠的爬高能力時，快速達到較高的高度，可以顯著減少氣動阻力損失，增大最終末狀態的能量。

4 結束語

UPE是基于UKF的參數估計方法，適用于非線性模型的參數估計問題，具有較大的收斂域。兩點邊值問題可以轉化為一個非線性參數估計問題，因此可以采用UPE進行求解。與傳統的兩點邊值問題求解算法相比，其收斂域更大，對初值精度要求不高，因此大大降低了兩點邊值問題的求解難度。

兩點邊值問題在諸多科學問題中都有涉及，特別是間接法最優化問題中，最終將最優化問題轉化為一個兩點邊值問題來求解，這在航空航天器軌跡優化與最優制導中應用廣泛。本文給出的算例，驗證了UPE算法求解兩點邊值問題應用在軌跡優化中的可行性和優越性，但對于大氣層內的軌跡優化問題，仍存在一些不足之處需要改進。

后續研究方向主要分為3個方向：①UKF濾波的方差初值和根據濾波情況自適應調整系統方差的算法；②UPE算法在直接法最優化問題中的應用；③粒子濾波算法在兩點邊值問題中的應用。