(西南科技大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽(yáng) 621010)

0 引言

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

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

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

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

1 滑翔段建模

1.1 氣動(dòng)模型

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

圖1 一體化飛行器數(shù)模

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

圖2 增量模型建模流程圖

(1)

(2)

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

1.2 滑翔段軌跡描述

滑翔段軌跡規(guī)劃是X-51A相似飛行器軌跡優(yōu)化的重要階段，滑翔飛行能延長(zhǎng)飛行器整體航程，實(shí)現(xiàn)遠(yuǎn)距離投送或打擊。由于高超聲速飛行器的動(dòng)力學(xué)系統(tǒng)較為復(fù)雜，其建模所涉及的因素繁多，為研究需要，將地球視作一個(gè)均勻球體，認(rèn)為滑翔飛行過(guò)程中飛行器是一個(gè)質(zhì)點(diǎn)，質(zhì)量保持不變，并忽略外界環(huán)境的干擾?；瓒蔚倪\(yùn)動(dòng)學(xué)方程形式含多個(gè)狀態(tài)變量與非線性等式約束，且方程中隱含有控制變量。本文以最大滑翔距離作為優(yōu)化目標(biāo)，其余指標(biāo)作為約束條件，構(gòu)建最優(yōu)軌跡控制問(wèn)題，以分析X-51A相似飛行器的滑翔彈道軌跡，運(yùn)動(dòng)模型表達(dá)式如下：

(3)

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

飛行器滑翔段為無(wú)動(dòng)力飛行，為保證平衡飛行，在滑翔段軌跡設(shè)計(jì)過(guò)程中要考慮相應(yīng)的約束條件，本文的參數(shù)約束如表1所示。

表1 滑翔段約束條件

2 偽譜法參數(shù)化

偽譜法源于求解微分方程的譜方法，其思想是在離散的狀態(tài)變量與控制變量基礎(chǔ)上，構(gòu)造逼近式，求得微分矩陣與積分算子，將最優(yōu)控制問(wèn)題轉(zhuǎn)換為NLP問(wèn)題,其核心技術(shù)包括：離散網(wǎng)格劃分策略的選擇；離散點(diǎn)插值函數(shù)對(duì)連續(xù)函數(shù)的逼近；插值函數(shù)的微積分對(duì)連續(xù)函數(shù)的微積分的近似。偽譜法參數(shù)化過(guò)程細(xì)化為以下幾步。

1)時(shí)域變換：

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

(4)

由于插值多項(xiàng)式的作用域?yàn)棣印蔥-1,1]，所以需要將t∈[0,tf]時(shí)域變換到τ∈[-1,1]：

(5)

2)變量離散化：

通過(guò)離散網(wǎng)格劃分策略選取合適的配點(diǎn)，能有效避免逼近可能出現(xiàn)的Runge現(xiàn)象，通過(guò)Gauss-Lobato求積，取n個(gè)節(jié)點(diǎn)將問(wèn)題離散，獲得離散時(shí)間：Tn=[τ1,τ2…τn+1]。離散時(shí)間對(duì)應(yīng)的狀態(tài)變量xn=[Hn,Vn,γn]，Hn=[h1,h2…h(huán)n+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配點(diǎn)近似逼近狀態(tài)與控制變量：

(6)

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

(7)

4)最優(yōu)控制轉(zhuǎn)化為NLP問(wèn)題：

通過(guò)前幾個(gè)步驟，可得到參數(shù)化后的NLP問(wèn)題，具體表達(dá)式如下所示：

(8)

其中: [Hn,Vn,γn]為軌跡狀態(tài)變量的離散向量，n為偽譜配點(diǎn)數(shù)；D,L分別為各離散狀態(tài)量下的阻力與升力；D為微分矩陣；wi為積分算子在對(duì)應(yīng)第i個(gè)狀態(tài)變量下的積分權(quán)重：

(9)

3 非線性規(guī)劃問(wèn)題求解

3.1 多級(jí)迭代優(yōu)化策略

初始值的選取對(duì)于NLP問(wèn)題的求解來(lái)說(shuō)具有很大影響。不可靠的初值，不僅會(huì)影響收斂速度，甚至有可能導(dǎo)致無(wú)法收斂，找不到可行解或最優(yōu)解的情況。并且隨著NLP問(wèn)題規(guī)模的增大，約束性的強(qiáng)增，初值對(duì)求解結(jié)果的影響也會(huì)不斷增大。偽譜法求解飛行器最優(yōu)軌跡的精度是由配點(diǎn)數(shù)量的多少所決定的，一般情況下配點(diǎn)數(shù)越多，求解精度越高。但是轉(zhuǎn)化得到的NLP問(wèn)題的規(guī)模也會(huì)變得更大，在使用SQP算法求解時(shí)，會(huì)出現(xiàn)初值敏感、收斂速度減慢、陷入局部最優(yōu)的情況。而配點(diǎn)數(shù)過(guò)少，則會(huì)導(dǎo)致偽譜法求解精度無(wú)法滿足要求。

本文在滑翔軌跡問(wèn)題離散過(guò)程中取n個(gè)配點(diǎn)時(shí)，會(huì)形成n+1個(gè)離散變量，對(duì)應(yīng)決策向量的維度為4n+5，其中狀態(tài)變量與控制變量構(gòu)成4(n+1)維決策變量，最后一維決策變量為終端時(shí)間tf，采用一個(gè)多級(jí)迭代優(yōu)化策略，對(duì)含高維度決策變量的NLP問(wèn)題求解初值進(jìn)行選取。

首先設(shè)置較少偽譜配點(diǎn)，保證問(wèn)題規(guī)模不至于影響SQP算法的求解，對(duì)于得到的結(jié)果進(jìn)行插值，然后將插值后的結(jié)果作為增加配點(diǎn)后的大規(guī)模NLP問(wèn)題的初始值。詳細(xì)過(guò)程如下：

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

2)設(shè)置少量偽譜配點(diǎn)，通過(guò)偽譜法參數(shù)化構(gòu)建中小規(guī)模NLP問(wèn)題，以x0作為初始值，使用SQP算法求解NLP問(wèn)題得到粗略解x1；

3)增加配點(diǎn)數(shù)，對(duì)獲得的粗略解x1進(jìn)行插值，將插值結(jié)果作為增大規(guī)模后的NLP問(wèn)題的初始值，求解NLP問(wèn)題，重復(fù)，直到收斂。

多級(jí)迭代優(yōu)化策略流程圖如圖3所示。

圖3 多級(jí)迭代優(yōu)化策略流程圖

在初始猜測(cè)插值過(guò)程中，由于控制變量不光滑，如果與狀態(tài)變量一樣選擇使用Lagrange多項(xiàng)式插值，其變化過(guò)大時(shí)會(huì)出現(xiàn)Gipps現(xiàn)象，導(dǎo)致控制變量超出邊界。分段三次Hermite插值多項(xiàng)式不僅要求在節(jié)點(diǎn)上的值要等于目標(biāo)函數(shù)值，而且要求在節(jié)點(diǎn)上的導(dǎo)數(shù)值也要與目標(biāo)函數(shù)的導(dǎo)數(shù)值相等，光滑性更好。選擇使用分段三次Hermite插值對(duì)離散的控制變量進(jìn)行插值，可有效避免出現(xiàn)Gipps現(xiàn)象。

3.2 NLP問(wèn)題求解

SQP算法是求解強(qiáng)約束NLP問(wèn)題最有效的方法之一，其思想是將原NLP問(wèn)題拆分為多個(gè)二次規(guī)劃子問(wèn)題，通過(guò)逐步求解二次規(guī)劃子問(wèn)題的解，為主迭代過(guò)程提供搜索方向，保證主迭代逐漸收斂于NLP問(wèn)題的最優(yōu)解。SQP算法的簡(jiǎn)要計(jì)算步驟如下：

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

2)求解二次規(guī)劃子問(wèn)題，確定下次迭代搜索方向d：

(10)

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

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

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

5)一維搜索：對(duì)minP(x(k)+λd(k),r)進(jìn)行求解，得到最優(yōu)步長(zhǎng)λ(k)，令x(k+1)=x(k)+λd(k)；

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

4 實(shí)驗(yàn)結(jié)果與分析

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

經(jīng)過(guò)反復(fù)計(jì)算，整理結(jié)果如表2～3所示。表2中n為配點(diǎn)數(shù)，z為對(duì)應(yīng)的決策變量數(shù)，t(s)為使用傳統(tǒng)方法求解對(duì)應(yīng)配點(diǎn)數(shù)規(guī)模的NLP問(wèn)題計(jì)算時(shí)間，titer(s)為使用多級(jí)迭代優(yōu)化策略下的計(jì)算時(shí)間。通過(guò)對(duì)比可以看出隨著偽譜法配點(diǎn)個(gè)數(shù)的增加，計(jì)算所需要的時(shí)間也在快速增加，多級(jí)迭代優(yōu)化策略可有效的減少計(jì)算時(shí)間。配點(diǎn)數(shù)越多，節(jié)省的計(jì)算時(shí)間越多，體現(xiàn)出的效果越好。

仿真實(shí)驗(yàn)過(guò)程中設(shè)置的最大迭代次數(shù)為1 000，表3中Kopt為最優(yōu)迭代次數(shù)，L為計(jì)算求得的目標(biāo)值。通過(guò)分析比較，本文提出的多級(jí)迭代優(yōu)化策略能大幅度減少迭代次數(shù)，加快收斂速度，獲得更精確的解。圖4、圖5為配點(diǎn)數(shù)n=40時(shí)，傳統(tǒng)方法與使用多級(jí)迭代優(yōu)化策略計(jì)算下，目標(biāo)值的迭代情況?？梢钥吹剑褂枚嗉?jí)迭代優(yōu)化策略時(shí)，SQP算法能更快的收斂得到最優(yōu)值，減少計(jì)算時(shí)間。

表2 計(jì)算時(shí)間對(duì)比情況

表3 目標(biāo)值與迭代次數(shù)

圖4 傳統(tǒng)方法目標(biāo)值迭代情況

圖5 多級(jí)迭代優(yōu)化策略目標(biāo)值迭代情況

SQP算法對(duì)偽譜法的離散配點(diǎn)數(shù)量敏感，配點(diǎn)數(shù)越多決策變量維度越高，算法收斂時(shí)間越長(zhǎng)。但配點(diǎn)數(shù)過(guò)少，則會(huì)影響結(jié)果精度。對(duì)不同配點(diǎn)數(shù)求得的結(jié)果，進(jìn)行反解常微分方程，結(jié)果表明配點(diǎn)數(shù)過(guò)少時(shí)，求得的結(jié)果會(huì)出現(xiàn)很大偏差，情況如圖6所示。

圖6 配點(diǎn)數(shù)對(duì)解的影響

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

圖7 PCHIP插值克服Gipps現(xiàn)象

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

圖8 滑翔軌跡對(duì)比圖

圖9 速度變化對(duì)比圖

圖10 航跡角變化對(duì)比圖

圖11 俯仰舵偏變化對(duì)比圖

5 結(jié)束語(yǔ)

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

隨著高超聲速技術(shù)的發(fā)展，全軌跡優(yōu)化的研究也在火熱進(jìn)行中，所面對(duì)的問(wèn)題規(guī)模將更大，約束性將更強(qiáng)，多級(jí)迭代策略體現(xiàn)的價(jià)值也將更大。