基于偽譜法的滑翔軌跡多級迭代優化策略

(西南科技大學 計算機科學與技術學院，四川 綿陽 621010)

0 引言

自中國進入“超級航空2018”，就預示著在未來幾十年中我國將在航空領域迎來更大的突破，而高超聲速技術則是突破航空領域的重要研究點。飛行器最優軌跡規劃是高超聲速技術研究的關鍵點之一，其求解方法主要有兩大類：間接法和直接法[1]。偽譜法以較少的參數、較高的計算精度、指數級的收斂速度，同時兼具了直接法和間接法的優勢，被認為是目前最具潛力的軌跡最優控制求解方法[2]。

直接求解飛行器軌跡最優控制問題難度較大，利用偽譜法將軌跡最優控制問題進行離散，對離散化得到的非線性規劃(Nonlinear Programming，NLP)問題進行插值獲得原問題的解，可將求解最優軌跡控制問題上的難度轉移到NLP問題的求解上，轉化得到的NLP問題表現出優化空間大、約束控制量多等特點[3-4]。近幾年，在選擇求解該類NLP問題的策略上，很多研究者偏向使用序列二次規劃(Sequential Quadratic Programming,SQP)算法[5]來求解。由于SQP算法可以將一些約束添加到變量中去，對于強約束問題的求解有很大優勢。近兩年有研究者使用SQP算法對高超聲速飛行器周期巡航軌跡[6]、爬升軌跡[7]進行求解，取得了較好的結果。

作為求解強約束NLP問題的主流方法，SQP算法具有計算速度快，求解精度高等優勢，但當NLP問題規模過大時，使用SQP算法求解存在初值敏感，收斂速度減慢的情況。因此求解較大規模的NLP問題時，選擇一個可靠的初值是十分必要的。K.Okuda等人[8]在求解帶翼火箭彈道優化時，提出使用具有全局最優性的遺傳算法為SQP算法提供初始值，有效的解決初值選取問題。但遺傳算法等全局優化算法在求解約束性較強的NLP問題時，存在精度低，收斂速度慢，結果隨機抖動等問題。

本文結合偽譜法與SQP算法各自的優勢，提出基于偽譜法的滑翔軌跡多級迭代優化策略。以X-51A[9]相似飛行器模型為研究對象，利用增量法與查表插值建立飛行器縱向氣動力模型。使用偽譜法與SQP算法對滑翔段軌跡進行求解，提出多級迭代優化策略彌補軌跡規劃過程中，SQP算法求解大規模NLP問題的不足之處，并與傳統方法求解出的狀態量與控制量仿真飛行狀態進行對比。仿真結果表明本文提出的多級迭代優化策略穩定性好，計算精度高，收斂速度快，在實際工程應用中起到了很好的支撐作用。

1 滑翔段建模

1.1 氣動模型

X-51A相似飛行器是氣動推進一體化飛行器，如圖1所示，模型總長4.26 m，最大寬度0.58 m，最大高度0.566 m，翼展最大0.927 m，兩舵面與水平面成45°，內流道寬0.225 m，進氣道總收縮比為4.9，內收縮比為2.4。

圖1 一體化飛行器數模

本文借鑒意大利航空航天中心(CIRA)[10]利用增量法實現PRORA-USV1高超聲速飛行器風洞試驗數據與數值計算數據的融合建模思路，使用數值方法模擬X-51A相似飛行器的數據集進行建模，為仿真實驗提供數據基礎。結合增量建模理論，在基準數據集與舵效數據集的基礎上構建增量模型，流程圖如圖2所示。

圖2 增量模型建模流程圖

(1)

(2)

式中，ρ∞為大氣密度，L為升力，D為阻力，m表示飛行器質量，Sref=1m2,bref=1m2，V∞表示飛行速度。

1.2 滑翔段軌跡描述

滑翔段軌跡規劃是X-51A相似飛行器軌跡優化的重要階段，滑翔飛行能延長飛行器整體航程，實現遠距離投送或打擊。由于高超聲速飛行器的動力學系統較為復雜，其建模所涉及的因素繁多，為研究需要，將地球視作一個均勻球體，認為滑翔飛行過程中飛行器是一個質點，質量保持不變，并忽略外界環境的干擾。滑翔段的運動學方程形式含多個狀態變量與非線性等式約束，且方程中隱含有控制變量。本文以最大滑翔距離作為優化目標，其余指標作為約束條件，構建最優軌跡控制問題，以分析X-51A相似飛行器的滑翔彈道軌跡，運動模型表達式如下：

(3)

式中，x=[H,V,γ]為軌跡狀態變量，分別表示飛行高度、飛行速度、飛行航跡角；u=α為軌跡的控制變量，α為攻角；tf為終端時間；D為阻力方向的合力，L為升力，m為飛行器質量，g為重力加速度，δ為最優控制變量α對應飛行姿態下保持縱向平衡的俯仰舵偏。J為目標函數，表示滑翔段的最大飛行距離。

飛行器滑翔段為無動力飛行，為保證平衡飛行，在滑翔段軌跡設計過程中要考慮相應的約束條件，本文的參數約束如表1所示。

表1 滑翔段約束條件

2 偽譜法參數化

偽譜法源于求解微分方程的譜方法，其思想是在離散的狀態變量與控制變量基礎上，構造逼近式，求得微分矩陣與積分算子，將最優控制問題轉換為NLP問題,其核心技術包括：離散網格劃分策略的選擇；離散點插值函數對連續函數的逼近；插值函數的微積分對連續函數的微積分的近似。偽譜法參數化過程細化為以下幾步。

1)時域變換：

本文使用Legendre多項式作為插值基，其系數{Ln(τ)}是τ∈[-1,1]上的正交多項式。當n≥1時滿足以下遞推關系：

(4)

由于插值多項式的作用域為τ∈[-1,1]，所以需要將t∈[0,tf]時域變換到τ∈[-1,1]：

(5)

2)變量離散化：

通過離散網格劃分策略選取合適的配點，能有效避免逼近可能出現的Runge現象，通過Gauss-Lobato求積，取n個節點將問題離散，獲得離散時間：Tn=[τ1,τ2…τn+1]。離散時間對應的狀態變量xn=[Hn,Vn,γn]，Hn=[h1,h2…hn+1]，Vn=[v1,v2…vn+1]，γn=[γ1,γ2…γn+1]，mn=[m1,m2…mn+1]；控制變量un=[αn]，αn=[α1,α2…αn+1]。

3)LGL配點近似逼近狀態與控制變量：

(6)

其中，Xi∈Rnx,Ui∈Rnu，φi(τ)為n階Lagrange插值基數：

(7)

4)最優控制轉化為NLP問題：

通過前幾個步驟，可得到參數化后的NLP問題，具體表達式如下所示：

(8)

其中: [Hn,Vn,γn]為軌跡狀態變量的離散向量，n為偽譜配點數；D,L分別為各離散狀態量下的阻力與升力；D為微分矩陣；wi為積分算子在對應第i個狀態變量下的積分權重：

(9)

3 非線性規劃問題求解

3.1 多級迭代優化策略

初始值的選取對于NLP問題的求解來說具有很大影響。不可靠的初值，不僅會影響收斂速度，甚至有可能導致無法收斂，找不到可行解或最優解的情況。并且隨著NLP問題規模的增大，約束性的強增，初值對求解結果的影響也會不斷增大。偽譜法求解飛行器最優軌跡的精度是由配點數量的多少所決定的，一般情況下配點數越多，求解精度越高。但是轉化得到的NLP問題的規模也會變得更大，在使用SQP算法求解時，會出現初值敏感、收斂速度減慢、陷入局部最優的情況。而配點數過少，則會導致偽譜法求解精度無法滿足要求。

本文在滑翔軌跡問題離散過程中取n個配點時，會形成n+1個離散變量，對應決策向量的維度為4n+5，其中狀態變量與控制變量構成4(n+1)維決策變量，最后一維決策變量為終端時間tf，采用一個多級迭代優化策略，對含高維度決策變量的NLP問題求解初值進行選取。

首先設置較少偽譜配點，保證問題規模不至于影響SQP算法的求解，對于得到的結果進行插值，然后將插值后的結果作為增加配點后的大規模NLP問題的初始值。詳細過程如下：

1)在變量上下邊界條件之間進行初始猜測插值，獲得一個初始值x0；

2)設置少量偽譜配點，通過偽譜法參數化構建中小規模NLP問題，以x0作為初始值，使用SQP算法求解NLP問題得到粗略解x1；

3)增加配點數，對獲得的粗略解x1進行插值，將插值結果作為增大規模后的NLP問題的初始值，求解NLP問題，重復，直到收斂。

多級迭代優化策略流程圖如圖3所示。

圖3 多級迭代優化策略流程圖

在初始猜測插值過程中，由于控制變量不光滑，如果與狀態變量一樣選擇使用Lagrange多項式插值，其變化過大時會出現Gipps現象，導致控制變量超出邊界。分段三次Hermite插值多項式不僅要求在節點上的值要等于目標函數值，而且要求在節點上的導數值也要與目標函數的導數值相等，光滑性更好。選擇使用分段三次Hermite插值對離散的控制變量進行插值，可有效避免出現Gipps現象。

3.2 NLP問題求解

SQP算法是求解強約束NLP問題最有效的方法之一，其思想是將原NLP問題拆分為多個二次規劃子問題，通過逐步求解二次規劃子問題的解，為主迭代過程提供搜索方向，保證主迭代逐漸收斂于NLP問題的最優解。SQP算法的簡要計算步驟如下：

1)初值設置：選擇初始點x0，初始矩陣H1是n×n的正定矩陣，令k=1，計算精度0<ε<Δ；

2)求解二次規劃子問題，確定下次迭代搜索方向d：

(10)

得到滿足原問題Kuhn-Tucker條件的近似解(d(k),λ(k));

3)若‖d(k)‖≤ε，則x(k+1)=x(k)+d(k)為求得的最優解，停止計算；

4)如果P(x(k)+d(k),r)

5)一維搜索：對minP(x(k)+λd(k),r)進行求解，得到最優步長λ(k)，令x(k+1)=x(k)+λd(k)；

6)更新矩陣H(k)：根據擬牛頓公式更新Hessen矩陣H(k)，確定新的正定矩陣H(k+1)，令k=k+1，轉到第2)步。

4 實驗結果與分析

通過MATLAB進行仿真驗證，仿真實驗過程中設置X-51A相似飛行器質量m=1 600 kg，考慮到滑翔飛行的終點與投射、打擊目標之間還存在一定的距離，因此需要對高度與速度設置相應的終端條件。設置初始速度V0=6 mach，初始高度H0=27 000 m，終端速度Vf=3 mach，終端高度Hf=1 000 m，為得到最遠滑翔航程，不對終端時間做限制。

經過反復計算，整理結果如表2～3所示。表2中n為配點數，z為對應的決策變量數，t(s)為使用傳統方法求解對應配點數規模的NLP問題計算時間，titer(s)為使用多級迭代優化策略下的計算時間。通過對比可以看出隨著偽譜法配點個數的增加，計算所需要的時間也在快速增加，多級迭代優化策略可有效的減少計算時間。配點數越多，節省的計算時間越多，體現出的效果越好。

仿真實驗過程中設置的最大迭代次數為1 000，表3中Kopt為最優迭代次數，L為計算求得的目標值。通過分析比較，本文提出的多級迭代優化策略能大幅度減少迭代次數，加快收斂速度，獲得更精確的解。圖4、圖5為配點數n=40時，傳統方法與使用多級迭代優化策略計算下，目標值的迭代情況。可以看到，使用多級迭代優化策略時，SQP算法能更快的收斂得到最優值，減少計算時間。

表2 計算時間對比情況

表3 目標值與迭代次數

圖4 傳統方法目標值迭代情況

圖5 多級迭代優化策略目標值迭代情況

SQP算法對偽譜法的離散配點數量敏感，配點數越多決策變量維度越高，算法收斂時間越長。但配點數過少，則會影響結果精度。對不同配點數求得的結果，進行反解常微分方程，結果表明配點數過少時，求得的結果會出現很大偏差，情況如圖6所示。

圖6 配點數對解的影響

由于控制變量的不光滑性，使用Lagrange多項式插值，在變量變化較大的情況下會出現Gipps現象，造成控制變量超出邊界，情況如圖7所示。采用分段三次Hermite插值方法對控制變量進行插值，能保證光滑性，有效避免Gipps現象與Runge現象。

圖7 PCHIP插值克服Gipps現象

通過仿真實驗，對比本文方法與傳統方法求解的狀態變量與俯仰舵偏的飛行狀態。結果表明多級迭代優化策略相比于傳統方法，在求解精度和收斂速度上有很大的優勢，圖8～11對比了兩種方法下配點數n=40時，狀態變量和俯仰舵偏角的變化情況，虛線為使用傳統方法的變化情況。傳統方法求解下，飛行高度與終端條件設置的高度相差大約1818 m，速度相差大約0.3 mach。多級迭代優化策略方法求解下，飛行高度與終端條件設置的高度相差大約86 m，速度相差大約0.01 mach。本文方法與傳統方法相比，求解偏差降低了很多，所得結果更加接近最優軌跡。并且在多級迭代優化策略求解下，飛行器滑翔時間更長，得到的滑翔水平距離比傳統方法遠大約220 m，能獲得更遠航程，實現更遠距離的投送或打擊。通過以上各項指標，證明了本文提出的多級迭代優化策略的有效性與高效性。

圖8 滑翔軌跡對比圖

圖9 速度變化對比圖

圖10 航跡角變化對比圖

圖11 俯仰舵偏變化對比圖

5 結束語

本文提出的多級迭代優化策略有效的結合了偽譜法參數化最優軌跡控制問題與SQP算法求解強約束非線性規劃問題的優勢。偽譜法參數化后的NLP問題具有大規模、強約束等特征，直接使用SQP算法求解時存在對初值敏感、收斂速度減慢等問題。多級迭代優化策略，為求解大規模NLP問題提供一個較為有效的初始值，可有效解決SQP算法在求解過程中存在的不足。仿真驗證表明，在多級迭代優化策略下，NLP求解算法求解精度更高，計算時間更少。且該方法在實際工程應用中取得了良好效果。

隨著高超聲速技術的發展，全軌跡優化的研究也在火熱進行中，所面對的問題規模將更大，約束性將更強，多級迭代策略體現的價值也將更大。