帶有控制變量變化率約束的偽譜軌跡優(yōu)化方法

趙吉松，尚 騰，張金明，程曉明，劉蛟龍

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

0 引 言

軌跡優(yōu)化對(duì)于導(dǎo)彈等飛行器的制導(dǎo)控制系統(tǒng)設(shè)計(jì)和總體設(shè)計(jì)具有重要意義[1-4]。軌跡優(yōu)化問題本質(zhì)上屬于最優(yōu)控制問題，其求解方法主要分為間接法和直接法[5]。間接法借助變分法或者最大值原理，把泛函優(yōu)化問題轉(zhuǎn)化為兩點(diǎn)邊值問題求解；直接法通過對(duì)控制變量和/或狀態(tài)變量進(jìn)行離散把泛函優(yōu)化轉(zhuǎn)化為非線性規(guī)劃(Non-linear programming，NLP)問題，然后采用各種非線性規(guī)劃算法求解，比如基于序列二次規(guī)劃(Sequential quadratic progra-mming，SQP)算法的SNOPT[6]和基于內(nèi)點(diǎn)法的IPOPT[7]。直接法中的配點(diǎn)法[5]由于不需要推導(dǎo)最優(yōu)性必要條件，并且對(duì)初值的敏感性較低，容易收斂，近年來得到廣泛研究和應(yīng)用。特別是其中的偽譜法[8-11]由于具有快速收斂性和較高優(yōu)化精度，近年來成為軌跡優(yōu)化和最優(yōu)控制領(lǐng)域一種非常常用的數(shù)值優(yōu)化方法。其中，開源軟件GPOPS[11]更是推動(dòng)了偽譜法的發(fā)展和應(yīng)用，逐漸發(fā)展為行業(yè)標(biāo)準(zhǔn)。

對(duì)于軌跡優(yōu)化而言，控制變量的變化率約束和端點(diǎn)約束具有重要工程價(jià)值，甚至影響軌跡方案的可行性。以導(dǎo)彈為例，初始攻角通常由發(fā)射條件給定，有時(shí)為了達(dá)到更好的攻擊效果對(duì)末端攻角也有限制。在實(shí)際飛行過程中，還需要對(duì)攻角變化率進(jìn)行限制。對(duì)于流量可調(diào)的超燃沖壓發(fā)動(dòng)機(jī)，不僅攻角變化率有約束要求，發(fā)動(dòng)機(jī)的流量變化率也有嚴(yán)格的約束要求[12]。在軌跡優(yōu)化時(shí)如果不考慮這些約束可能會(huì)導(dǎo)致優(yōu)化的控制變量變化率過大甚至發(fā)生突變，制導(dǎo)控制系統(tǒng)難以對(duì)其進(jìn)行跟蹤。為了考慮控制變量的變化率約束和端點(diǎn)約束，目前常用的處理方式是引入控制變量的導(dǎo)數(shù)作為虛擬控制變量，將實(shí)際控制變量轉(zhuǎn)換為狀態(tài)變量[12-14]。這種處理方式的優(yōu)點(diǎn)是形式簡單，只需要增加1個(gè)狀態(tài)變量和1個(gè)狀態(tài)方程，對(duì)軌跡優(yōu)化方法的繼承性比較好，幾乎不需要對(duì)軌跡優(yōu)化方法進(jìn)行改動(dòng)。但是，這種引入控制變量導(dǎo)數(shù)作為虛擬控制變量的處理方法會(huì)帶來一些問題[5]。一方面，引入新的狀態(tài)變量顯然增加了狀態(tài)變量和狀態(tài)方程的個(gè)數(shù)，進(jìn)而增加了離散后的NLP的約束和優(yōu)化變量的個(gè)數(shù)。另一方面，增加的控制變量微分方程為線性微分方程，很可能會(huì)產(chǎn)生奇異弧進(jìn)而引起虛擬控制變量振蕩甚至發(fā)散。這兩方面原因都會(huì)增加NLP的求解難度(需要更多的迭代次數(shù)才能收斂，甚至不收斂)。

最優(yōu)控制領(lǐng)域的著名學(xué)者Betts[5]在研究局部配點(diǎn)法時(shí)提出一種直接施加控制變量變化率約束的方法，其思路是采用有限差分法逼近控制變量的變化率，從而將控制變量變化率約束轉(zhuǎn)化為相鄰離散節(jié)點(diǎn)的控制變量的線性約束。Betts將該方法成功用于局部配點(diǎn)法的梯形離散格式(二階格式)和Hermite-Simpson格式(四階格式)，避免了傳統(tǒng)方法在處理控制變量變化率約束時(shí)需要引入虛擬控制變量以及由此導(dǎo)致的虛擬控制變量振蕩問題，容易收斂。該方法能夠方便地施加控制變量變化率約束，并且添加的控制變量變化率約束全部為線性約束，計(jì)算量比較小。與局部配點(diǎn)法相比，偽譜軌跡優(yōu)化方法[8-11]在計(jì)算效率和精度方面具有更好的綜合性能，近年來得到了快速發(fā)展和廣泛應(yīng)用，但是目前公開文獻(xiàn)中還沒有基于偽譜法開展直接施加控制變量變化率約束的研究。如果將Betts針對(duì)局部配點(diǎn)法提出的直接施加控制變量變化率約束的方法與偽譜法相結(jié)合，有可能使偽譜法能夠方便、快速地求解帶有控制變量變化率約束的軌跡優(yōu)化問題，這對(duì)于拓展偽譜法的功能和性能顯然具有重要意義。

本文在傳統(tǒng)多區(qū)間偽譜法的基礎(chǔ)上，引入一種控制變量變化率約束和端點(diǎn)約束的處理方法，以空空導(dǎo)彈軌跡優(yōu)化問題和臨近空間高超聲速滑翔軌跡優(yōu)化問題為例，通過數(shù)值仿真和對(duì)比驗(yàn)證了所建立的控制變量變化率約束處理方法的有效性。

1 帶有控制變量變化率約束的軌跡優(yōu)化方法

1.1 多區(qū)間Radau偽譜法

與全局偽譜法相比，多區(qū)間偽譜法[11]具有更高的計(jì)算效率。軌跡優(yōu)化本質(zhì)上是最優(yōu)控制問題。為了便于描述，本節(jié)以一般形式的Bolza型最優(yōu)控制問題為例，介紹多區(qū)間Radau偽譜法。

一般形式的Bolza型問題可描述為：求解控制變量u(t)∈Rm，使得如下目標(biāo)函數(shù)最小化

(1)

式中：M：Rn×R×Rn×R→R，L：Rn×Rm×R→R，x∈Rn，u∈Rm，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×Rn×R→Re；C:Rn×Rm×R→Rc。方程(1)～(4)所描述的連續(xù)最優(yōu)控制問題稱為Bolza型最優(yōu)控制問題。

為了便于應(yīng)用偽譜法對(duì)上述一般形式的最優(yōu)控制問題進(jìn)行離散，需要將最優(yōu)控制問題的時(shí)間區(qū)間t∈[t0,tf]變換至τ∈[-1,+1]，變換方式如下：

(5)

多區(qū)間偽譜法將時(shí)間區(qū)間τ∈[-1,+1]劃分為K個(gè)子區(qū)間Ik=[Tk-1,Tk]，k=1,…,K，其中-1=T0

(6)

(7)

將方程(6)中X(k)(τ)關(guān)于τ求微分得到

(8)

將方程(2)和(8)相結(jié)合，可以得到狀態(tài)方程在區(qū)間Ik(k=1,…,K)上的離散形式為

i=1,…,Nk

(9)

目標(biāo)函數(shù)(1)的離散形式為

(10)

端點(diǎn)條件(3)的離散形式為

(11)

路徑約束(4)的離散形式為

(12)

(13)

作為多區(qū)間偽譜法的特例，全局配點(diǎn)法可以看作只包含1個(gè)子區(qū)間的多區(qū)間偽譜法，并且在該子區(qū)間應(yīng)用高階多項(xiàng)式離散最優(yōu)控制問題。

為了提高計(jì)算效率和優(yōu)化精度，本文應(yīng)用文獻(xiàn)[10]給出的hp自適應(yīng)節(jié)點(diǎn)細(xì)化技術(shù)在優(yōu)化過程中調(diào)整子區(qū)間分布以及每個(gè)子區(qū)間內(nèi)的節(jié)點(diǎn)數(shù)量。

1.2 控制變量斜率約束

假設(shè)一般形式的控制變量變化率約束為

(14)

式中：u(t)表示原始軌跡優(yōu)化問題的控制變量；rL和rU分別為控制變量變化率的下限和上限。

目前常用的處理方法是引入實(shí)際控制變量的導(dǎo)數(shù)作為虛擬控制變量，從而將原始控制變量轉(zhuǎn)換為狀態(tài)變量，其變化規(guī)律滿足以下狀態(tài)方程

(15)

式中：σ(t)為引入的虛擬控制變量。方程(15)與原始動(dòng)力學(xué)方程(2)組成擴(kuò)展形式的微分方程組。

虛擬控制變量的邊界約束如下

rL≤σ(t)≤rU

(16)

如果實(shí)際控制變量存在初始約束和終端約束，可以直接通過如下方式施加約束

u(t0)=u0,u(tf)=uf

(17)

式中：u0和uf分別為控制變量的初值和終值。

這種引入控制變量導(dǎo)數(shù)作為虛擬控制變量的處理方法雖然簡單，但是會(huì)帶來一些問題。一方面，引入了新狀態(tài)變量，顯然增加了狀態(tài)變量和狀態(tài)方程的個(gè)數(shù)，進(jìn)而增加了NLP的約束和優(yōu)化變量的個(gè)數(shù)。盡管方程是線性微分方程，但是其離散形式是非線性的(與初始時(shí)間t0和終端時(shí)間tf存在耦合)。另一方面，方程為線性微分方程，容易產(chǎn)生奇異弧進(jìn)而引起虛擬控制變量振蕩甚至發(fā)散。這兩方面原因都會(huì)增加NLP的求解難度和計(jì)算量。

為了克服上述控制變量斜率約束處理方法的弊端，本文引入有限差分法對(duì)約束進(jìn)行離散

(18)

式中：uk=u(τk); Δτk=τk+1-τk。對(duì)其進(jìn)行變換和整理，可以得到如下兩個(gè)不等式約束

0≤uk+1-uk-[Δτk(tf-t0)]rL

(19)

0≥uk+1-uk-[Δτk(tf-t0)]rU

(20)

將其進(jìn)一步整理成標(biāo)準(zhǔn)形式得到

uk-uk+1-ΔτkrLt0+ΔτkrLt0≤0

(21)

-uk+uk+1+ΔτkrUt0-ΔτkrUtf≤0

(22)

對(duì)于控制變量的初始約束和終端約束，可以直接通過以下方式進(jìn)行限制

u(τ0)=u0,u(τf)=uf

(23)

式中：u0和uf分別為控制變量的初值和終值。

方程(21)、(22)和(23)全部為線性約束(離散控制變量、初始時(shí)間和終端時(shí)間的線性函數(shù))，其特點(diǎn)是一階偏導(dǎo)數(shù)為常數(shù)、二階偏導(dǎo)數(shù)為零，在求解NLP時(shí)不需要在每步迭代中反復(fù)計(jì)算偏導(dǎo)數(shù)，因而不會(huì)明顯增加計(jì)算量。在方程(9)-(13)描述的傳統(tǒng)偽譜軌跡法的基礎(chǔ)上，增加方程(21)、(22)和(23)描述的控制變量變化率約束，即可得到考慮控制變量變化率約束的軌跡優(yōu)化方法。圖1給出該方法的流程圖。相對(duì)于傳統(tǒng)偽譜法，在非線性規(guī)劃、歸一化處理和求解器模塊需要擴(kuò)展以考慮控制變量變化率約束，在稀疏型分析和偏導(dǎo)數(shù)計(jì)算模塊需要考慮控制變量變化率約束的常值偏導(dǎo)數(shù)，其它模塊相同。

圖1 帶有控制量變化率約束的軌跡優(yōu)化方法流程圖Fig.1 Procedure of trajectory optimization method with the constraint on the change rate of control variables

2 仿真算例與結(jié)果

本節(jié)采用兩個(gè)軌跡優(yōu)化算例展示本文方法的有效性，并與GPOPS軟件進(jìn)行了對(duì)比。其中，計(jì)算平臺(tái)為桌面計(jì)算機(jī)，處理器為Intel Core i3-12100 3.3 GHz，內(nèi)存為DDR4 8 GB，操作系統(tǒng)為Windows 11，編程語言為MATLAB 2017a。采用SNOPT求解NLP，采用稀疏有限差分法[15]提供一階偏導(dǎo)數(shù)。文中給出的優(yōu)化耗時(shí)為10次運(yùn)行的平均耗時(shí)。

2.1 雙脈沖空空彈軌跡優(yōu)化

本算例優(yōu)化空空導(dǎo)彈的攻角和推力脈沖時(shí)間間隔，使得導(dǎo)彈飛行時(shí)間最短或者射程最大。在縱向剖面內(nèi)，描述導(dǎo)彈質(zhì)心運(yùn)行的微分方程組為

(24)

式中：x為水平位置；h為飛行高度；v為飛行速度；γ為航跡角；m為質(zhì)量；T為發(fā)動(dòng)機(jī)推力，Isp為燃料比沖，g0為海平面重力加速度；L和D分別為升力和阻力；α為攻角，為需要優(yōu)化的控制變量。

升力和阻力表達(dá)式如下：

(25)

式中：ρ為大氣密度；Sref為氣動(dòng)參考面積；CL和CD分別為升力系數(shù)和阻力系數(shù)；Ma為馬赫數(shù)。表1給出升力系數(shù)和阻力系數(shù)隨馬赫數(shù)和攻角的變化特性，其中參考面積Sref=0.02488 m2。該氣動(dòng)力數(shù)據(jù)是通過對(duì)美國空空導(dǎo)彈AIM-120C進(jìn)行氣動(dòng)外形建模，然后采用CFD方法計(jì)算得出。與解析形式的氣動(dòng)數(shù)據(jù)相比，離散表格形式的氣動(dòng)數(shù)據(jù)在實(shí)際工程應(yīng)用中更為常見，更具有代表性。

導(dǎo)彈采用雙脈沖發(fā)動(dòng)機(jī)提供推力，推力剖面如圖2所示。其中，兩個(gè)脈沖推力的大小和持續(xù)時(shí)間均為已知量，兩個(gè)脈沖之間的時(shí)間間隔為設(shè)計(jì)變量，本文通過軌跡優(yōu)化確定最優(yōu)的時(shí)間間隔。

狀態(tài)變量的初始條件為

(26)

表1 導(dǎo)彈的氣動(dòng)力系數(shù)Table 1 Aerodynamic coefficients of the missile

式中：r0,h0,v0,γ0和m0分別為狀態(tài)變量的初值。

圖2 雙脈沖發(fā)動(dòng)機(jī)推力剖面Fig.2 Thrust profile of the dual pulse engine

控制變量的初始條件為

α(t0)=α0

(27)

式中：α0為控制變量的初值。

為了充分發(fā)揮導(dǎo)彈的攻擊效果，導(dǎo)彈命中目標(biāo)時(shí)刻的速度需要滿足如下約束

v(tf)≥vf

(28)

式中：tf為終端時(shí)刻；vf為最低攻擊速度。

導(dǎo)彈的終端高度為目標(biāo)高度，即

h(tf)=hf

(29)

式中：hf為命中點(diǎn)處導(dǎo)彈的高度要求。

控制變量攻角的變化范圍約束如下：

αmin≤α(t)≤αmax

(30)

式中：αmin和αmax分別為攻角的下限和上限。

攻角的變化率需要滿足如下約束

(31)

式中：umin和umax分別為攻角變化率的下限和上限。

為了防止導(dǎo)彈速度過大帶來嚴(yán)重的氣動(dòng)加熱問題，需要對(duì)導(dǎo)彈的飛行速度進(jìn)行約束，即

Ma(t)≤Mamax

(32)

式中：Mamax為允許的最大飛行馬赫數(shù)。

本節(jié)研究兩類不同目標(biāo)函數(shù)的軌跡優(yōu)化問題。第一類為最短時(shí)間軌跡優(yōu)化問題(簡記為P1)，即導(dǎo)彈飛行的時(shí)間最短，目標(biāo)函數(shù)為

minJ2=tf

(33)

對(duì)于P1問題，通常還需要限定射程，即

x(tf)=xf

(34)

式中：xf為導(dǎo)彈的射程要求。

本節(jié)研究的第二類軌跡優(yōu)化問題為最大射程問題(簡記為P2)，目標(biāo)函數(shù)為

minJ2=-x(tf)

(35)

該軌跡優(yōu)化問題描述為：確定最優(yōu)攻角α(t)和脈沖時(shí)間間隔tg，使得目標(biāo)函數(shù)(33)或者(35)最小化，并且滿足動(dòng)力學(xué)方程組(24)，初始條件(26)和(27)，終端約束(28)和(29)，路徑約束(30)～(32)。

在射程給定的情況下，飛行時(shí)間更短意味著能夠更快地?fù)糁心繕?biāo)，因而P1問題具有工程實(shí)際價(jià)值。本文首先對(duì)P1問題進(jìn)行仿真。根據(jù)雙脈沖發(fā)動(dòng)機(jī)的工作特性，空空導(dǎo)彈軌跡優(yōu)化問題分為4段，分別為：第1脈沖工作段，脈沖間隔段，第2脈沖工作段，以及無動(dòng)力滑翔段。采用的仿真參數(shù)為：r0=0 m，h0=6000 m，v0=400 m/s，γ0=0，m0=157 kg，α0=0，αmin=-15°，αmax=15°，umin=-2 (°)/s，umax=2 (°)/s，Mαmax=4.6，hf=6000 m，vf=600 m/s，xf=80000 m，T1=22775 N，T2=18980 N，t1=3.5 s，t2=4.0 s，Isp=247.23 s，g0=9.80665 m/s2。其中，攻角變化率參照文獻(xiàn)[13]設(shè)置(文獻(xiàn)[13]給出的攻角變化率為1 (°)/s，本文將其適當(dāng)放大至2 (°)/s。雖然目前工程上允許設(shè)置更大的攻角變化率，但是變化緩和的攻角顯然更容易跟蹤。攻角變化率約束值越大，意味著攻角變化率約束越寬松(在極限情況下，攻角變化率為無窮大意味著攻角變化率實(shí)際上不受約束)；攻角變化率約束值越小，更能體現(xiàn)出攻角變化率約束的影響。對(duì)于本文方法，不同的攻角變化率約束值對(duì)于展示方法的功能沒有本質(zhì)影響。

作為對(duì)比，本文首先采用開源軟件GPOPS求解該問題。由于GPOPS不能直接施加控制變量初值和變化率約束，需要引入攻角導(dǎo)數(shù)作為虛擬控制變量，將攻角轉(zhuǎn)換為狀態(tài)變量。圖3給出GPOPS優(yōu)化的虛擬控制變量(攻角導(dǎo)數(shù))和作為狀態(tài)變量解算出的攻角隨時(shí)間變化曲線。其中h方法為多區(qū)間偽譜法(采用均勻分布的子區(qū)間，每個(gè)子區(qū)間包含4個(gè)離散節(jié)點(diǎn))，p方法為全局偽譜法，即單區(qū)間偽譜法，hp方法為采用hp節(jié)點(diǎn)自適應(yīng)方法調(diào)整節(jié)點(diǎn)分布的多區(qū)間偽譜法。對(duì)于h方法和p方法，各段軌跡采用的配點(diǎn)數(shù)目依次為12, 12, 12和32。對(duì)于hp方法，節(jié)點(diǎn)數(shù)目由節(jié)點(diǎn)細(xì)化算法迭代確定?？梢?，無論采用何種離散節(jié)點(diǎn)分布方式，GPOPS優(yōu)化的虛擬控制變量(攻角導(dǎo)數(shù))都發(fā)生了明顯的振蕩，并且在軌跡最右端的值超出了允許的變化范圍，不是可行解。盡管如此，在GPOPS中作為狀態(tài)變量的攻角沒有表現(xiàn)出明顯振蕩，只是在軌跡最右端存在一定程度的振蕩。但是對(duì)于GPOPS而言，虛擬控制變量在求解過程中被當(dāng)作實(shí)際控制變量，是需要優(yōu)化算法通過迭代求解的控制變量，該變量振蕩會(huì)導(dǎo)致NLP收斂緩慢，甚至收斂困難乃至不收斂，由此帶來的影響是增加優(yōu)化耗時(shí)，降低了方法的可靠性。

圖3 GPOPS優(yōu)化的P1問題的控制變量Fig.3 Control solution to P1 using GPOPS

圖4給出采用本文方法處理斜率約束時(shí)優(yōu)化的攻角隨時(shí)間變化曲線。其中，h方法和p方法的離散點(diǎn)分布方式和數(shù)量與GPOPS相同，hp方法的離散點(diǎn)數(shù)量和分布方式由節(jié)點(diǎn)細(xì)化算法[10]確定?？梢姡瑹o論采用何種離散點(diǎn)分布方式，本文方法優(yōu)化的攻角都不發(fā)生振蕩。圖5給出本文優(yōu)化的離散軌跡和數(shù)值積分軌跡的對(duì)比，采用的離散點(diǎn)為hp自適應(yīng)節(jié)點(diǎn)?？梢?，離散軌跡和數(shù)值積分軌跡非常一致，說明了優(yōu)化精度高。從圖5還可以看出，本文方法優(yōu)化的終端高度、射程和飛行馬赫數(shù)都嚴(yán)格滿足約束要求。導(dǎo)彈飛行80 km需要的最短飛行時(shí)間為82.67 s，雙脈沖推力的最佳時(shí)間間隔為6.35 s。

對(duì)比本文方法和GPOPS軟件可知：GPOPS軟件為了施加控制變量變化率約束，需要引入虛擬控制變量作為控制變量，而且在迭代求解過程中虛擬控制變量會(huì)振蕩；本文方法不需要引入虛擬控制變量，因而避免了虛擬控制變量振蕩問題，更容易收斂，能夠減少優(yōu)化耗時(shí)，可靠性更高(因?yàn)檎袷幱锌赡軙?huì)導(dǎo)致不收斂)。表2給出本文方法與GPOPS的優(yōu)化結(jié)果對(duì)比?？梢?，本文方法優(yōu)化耗時(shí)明顯少于GPOPS，對(duì)于不同的離散點(diǎn)分布方式，本文方法的優(yōu)化耗時(shí)約為GPOPS的1/4～1/5。在目標(biāo)函數(shù)和終端誤差方面，本文方法和GPOPS沒有明顯差異。其中，終端誤差定義為離散軌跡與積分軌跡在終端時(shí)刻的差異，Δr表示位置誤差。在評(píng)估GPOPS結(jié)果的終端誤差時(shí)，由于圖3(a)所示的攻角導(dǎo)數(shù)振蕩嚴(yán)重，據(jù)此積分帶有方程(15)的擴(kuò)展動(dòng)力學(xué)系統(tǒng)會(huì)導(dǎo)致終端誤差過大?？紤]到方程(15)只與攻角直接相關(guān)，與其它方程的關(guān)系是通過攻角間接聯(lián)系的，本文直接采用圖3(b)所示的GPOPS優(yōu)化的攻角作為控制輸入，積分方程(24)所示的動(dòng)力學(xué)方程組評(píng)估GPOPS結(jié)果的終端誤差。對(duì)于該算例，GPOPS結(jié)果的攻角導(dǎo)數(shù)振蕩并沒有引起明顯的攻角振蕩，因而不影響其它狀態(tài)方程的離散誤差。對(duì)于求解器而言，控制變量振蕩意味著收斂性差，增加了優(yōu)化耗時(shí)。

圖4 本文方法優(yōu)化的P1問題的控制變量Fig.4 Control solution to P1using the proposed method

圖5 本文方法優(yōu)化的P1問題的狀態(tài)變量Fig.5 State solution to P1 using the proposed method

表2 本文方法與GPOPS的優(yōu)化結(jié)果對(duì)比Table 2 Comparison of results obtained using GPOPS and the proposed method

在導(dǎo)彈概念設(shè)計(jì)和方案論證階段，最大射程是一個(gè)重要性能指標(biāo)。對(duì)于P2問題，仿真條件與P1問題相同(需要移除射程約束)。同樣為了對(duì)比，本文首先采用GPOPS軟件求解該問題。為了施加攻角變化率約束，GPOPS需要引入虛擬控制變量(攻角導(dǎo)數(shù))作為控制變量。圖6給出采用不同分布方式的離散點(diǎn)時(shí)GPOPS優(yōu)化的攻角導(dǎo)數(shù)和攻角隨時(shí)間變化曲線?？梢?，無論采用何種離散點(diǎn)分布方式，虛擬控制變量(攻角導(dǎo)數(shù))都發(fā)生了明顯振蕩，并且在軌跡最右端超出了限定的變化率約束范圍。這些特性與求解P1問題時(shí)相同，進(jìn)一步驗(yàn)證了GPOPS通過引入虛擬控制變量施加斜率約束容易發(fā)生振蕩是一種普遍現(xiàn)象而非特例。因?yàn)槊枋龉ソ亲兓实奈⒎址匠淌蔷€性微分方程，從最優(yōu)控制的角度而言容易產(chǎn)生奇異弧進(jìn)而引起攻角振蕩甚至發(fā)散，Betts[5]在其專著中對(duì)于這一問題進(jìn)行了論述。

圖7為在不同離散節(jié)點(diǎn)分布情況下采用本文方法求解該問題得到的最優(yōu)攻角隨時(shí)間變化曲線。以其中的hp節(jié)點(diǎn)為例，導(dǎo)彈的最大射程為100.24 km，雙脈沖推力的最佳時(shí)間間隔為6.63 s。

對(duì)比圖6和圖7可知，本文方法優(yōu)化的攻角隨時(shí)間變化更加光滑，其變化率嚴(yán)格滿足約束要求，在軌跡末端也不發(fā)生突變。更重要的是，本文方法直接以攻角作為控制變量，能夠直接施加攻角變化率約束，而GPOPS為了施加攻角變化率約束，需要引入虛擬控制變量(攻角導(dǎo)數(shù))作為控制變量，在迭代求解過程中面臨著虛擬控制變量振蕩問題。這些因素使得本文方法得到的NLP規(guī)模更小，而且避免了GPOPS面臨的虛擬控制變量振蕩問題，計(jì)算量更小，收斂速度更快，可靠性更高。以采用hp節(jié)點(diǎn)為例，GPOPS求解P2問題耗時(shí)11.65 s，本文方法求解該問題耗時(shí)2.24 s，約為GPOPS的1/5。在軌跡終端精度方面，本文方法與GPOPS沒有明顯差異。這些特性與求解P1問題時(shí)類似。

圖6 GPOPS優(yōu)化的P2問題的控制變量Fig.6 Control solution to P2 using GPOPS

圖7 本文方法求解的P2問題的控制變量Fig.7 Control solution to P2 using the proposed method

2.2 臨近空間滑翔軌跡優(yōu)化

臨近空間高超聲速滑翔飛行器[16]在快速打擊領(lǐng)域具有重要應(yīng)用價(jià)值。參數(shù)化剖面控制由于簡單、工程上容易實(shí)現(xiàn)，近年來逐漸受到重視[17]。本節(jié)研究基于參數(shù)化剖面控制的滑翔軌跡優(yōu)化問題。

如圖1所示，臨近空間高超聲速滑翔軌跡的飛行剖面可分成3段：第1段從起始點(diǎn)G0下滑到G1點(diǎn)，第2段從G1點(diǎn)持續(xù)滑翔到G2點(diǎn)(滑翔高度在hmin和hmax之間)，第3段從G2點(diǎn)持續(xù)到G3點(diǎn)。假設(shè)每段攻角為常值參數(shù)，那么總共有3個(gè)參數(shù)。本文的參數(shù)化剖面優(yōu)化就是對(duì)每段軌跡的攻角進(jìn)行優(yōu)化。由于攻角為常數(shù)意味著攻角變化率為0，因而本文建立的帶有控制變量變化率約束的軌跡優(yōu)化方法通過施加攻角變化率約束可以求解該問題。

圖8 分段常值攻角飛行剖面示意圖Fig.8 Sketch of sectional constant angle of attack profile

在二維縱向平面內(nèi)，描述臨近空間高超聲速滑翔飛行器質(zhì)心運(yùn)動(dòng)的微分方程組為

(36)

式中：r是飛行器距地心距離；θ是射程角；v是飛行器速度；γ是航跡角；m是飛行器的質(zhì)量；g是重力加速度，g=μ/r2，μ是地球引力常數(shù)；L和D分別為升力和阻力，定義與前述P1問題中相同。

再入滑翔軌跡的初始條件為

(37)

式中：t0為初始時(shí)間；h0,θ0,v0,γ0和m0分別為狀態(tài)變量的初值；高度h=r-Re，Re為地球半徑。

為了命中定點(diǎn)目標(biāo)，終端條件為

h(tf)=hf,θ(tf)=θf

(38)

式中：tf為終端時(shí)間；hf和θf分為目標(biāo)位置坐標(biāo)。

為了充分發(fā)揮武器的打擊效果，命中點(diǎn)的速度大小和方向需要滿足如下約束條件

v(tf)≥vf,γ(tf)≤γf

(39)

式中：vf為和γf為終端速度大小和方向約束邊界。

在軌跡的第2段，軌跡跳躍的高度范圍需要滿足如下約束(不包括初始段和末端下壓段)

hmin≤h(t)≤hmax

(40)

由于每段軌跡的攻角參數(shù)均為常數(shù)，因而在每段內(nèi)，攻角變化率需要滿足如下約束

(41)

在最短時(shí)間內(nèi)擊中給定的目標(biāo)顯然具有重要實(shí)戰(zhàn)意義。本文以時(shí)間最短為目標(biāo)函數(shù)，即

J=tf

(42)

臨近空間高超聲速滑翔軌跡優(yōu)化問題(簡記為P3)可描述為：求解3個(gè)常值攻角，使得目標(biāo)函數(shù)(42)最小化，并且滿足狀態(tài)方程(36)，初始條件(37)，終端條件(38)和(39)，路徑約束(40)和(41)。

對(duì)于P3問題，仿真參數(shù)為：h0=70000 m，θ0=0，v0=2835.9 m/s，γ0=0，m0=907.2 kg，hf=1000 m，θf=0.1570 rad，vf=1000 m/s，γf=-1.40 rad，hmin=35 km，hmax=45 km。飛行器的氣動(dòng)數(shù)據(jù)采用文獻(xiàn)[16]給出的CAV氣動(dòng)數(shù)據(jù)插值計(jì)算，大氣密度模型采用美國1976版本標(biāo)準(zhǔn)大氣模型插值計(jì)算。

圖9給出本文方法優(yōu)化的3段常值攻角，圖10給出以常值攻角作為輸入控制，采用數(shù)值積分方法得到的軌跡狀態(tài)變量曲線。其中，終端射程、終端速度、終端航跡角以及第2段軌跡的飛行高度范圍全部滿足約束要求?？梢?，本文方法通過限制控制變量的變化率為0，實(shí)現(xiàn)了對(duì)攻角參數(shù)的優(yōu)化。本算例優(yōu)化耗時(shí)為0.31 s，最優(yōu)目標(biāo)函數(shù)即飛行器滑翔給定射程需要的最短時(shí)間為444.26 s。與連續(xù)變化的攻角相比，分段常值攻角方案的優(yōu)勢是更容易跟蹤(實(shí)際工程應(yīng)用時(shí)還需要考慮相鄰常值攻角之間的連續(xù)過渡，本文暫未考慮這一細(xì)節(jié))。

圖11給出采用連續(xù)變化攻角方案優(yōu)化的攻角和滑翔軌跡曲線。其中，最短飛行時(shí)間為432.09 s，比分段常值攻角方案減少12.17 s(2.74%)。由圖11(a)可知，連續(xù)變化攻角方案的最優(yōu)攻角曲線比較復(fù)雜，顯然會(huì)增加制導(dǎo)跟蹤的難度。由圖11(b)可知，連續(xù)變化攻角方案的飛行軌跡與常值攻角方案的軌跡比較接近，說明了常值攻角方案的合理性。

圖9 本文方法優(yōu)化的P3問題的控制變量Fig.9 Control solution to P3 using the proposed method

圖10 本文方法優(yōu)化的P3問題的狀態(tài)變量Fig.10 State solution to P3 using the proposed method

3 結(jié) 論

本文提出一種考慮控制變量變化率約束的軌跡優(yōu)化方法。該方法在傳統(tǒng)偽譜法的基礎(chǔ)上，引入有限差分法將控制變量變化率約束轉(zhuǎn)化為計(jì)算量較小的線性約束，然后采用非線性規(guī)劃算法求解。該方法能夠方便地施加控制變量變化率約束，避免了傳統(tǒng)偽譜法為了施加控制變量變化率約束需要引入虛擬控制變量以及由此導(dǎo)致的虛擬控制變量振蕩問題，因而該方法具有更好的性能。采用帶有控制變量變化率約束的空空導(dǎo)彈軌跡優(yōu)化問題和臨近空間高超聲速飛行器滑翔軌跡優(yōu)化問題驗(yàn)證了方法的有效性。仿真結(jié)果表明，與引入控制量導(dǎo)數(shù)作為虛擬控制量的控制變量變化率約束處理方法相比，本文方法能夠?qū)④壽E優(yōu)化效率提高大約5倍。在可靠性方面，該方法由于避免了傳統(tǒng)方法引入虛擬控制變量導(dǎo)致的虛擬控制量振蕩問題，更容易收斂，更加可靠。本文的研究工作拓展了偽譜法的功能。

圖11 連續(xù)變化攻角剖面情況下P3問題的最優(yōu)解Fig.11 Optimal solution to P3 using continuously varying control profile