朱 熠,李繼廣,郝向宇

(1.陸軍工程大學(xué) 野戰(zhàn)工程學(xué)院,南京 210007;2.西安航空學(xué)院 飛行器學(xué)院,西安 710077)

動(dòng)態(tài)滑翔也稱為梯度風(fēng)場(chǎng)翱翔。一些鳥類,如信天翁,可以長(zhǎng)時(shí)間利用梯度風(fēng)場(chǎng)沿著海岸滑翔,它們不會(huì)扇動(dòng)翅膀。自19世紀(jì)80年代,學(xué)者們開始研究動(dòng)態(tài)飆升的過(guò)程,取得了巨大成就,人們關(guān)注的焦點(diǎn)是最優(yōu)軌道的設(shè)計(jì)和風(fēng)能的獲取。

梯度風(fēng)場(chǎng)中的動(dòng)態(tài)滑翔對(duì)風(fēng)場(chǎng)的特性有一定的要求。通常,風(fēng)場(chǎng)的方向要求是水平的,高度方向上是梯度變化。動(dòng)態(tài)滑翔的典型飛行軌跡如圖1所示。一個(gè)完整的動(dòng)力翱翔過(guò)程一般包括四個(gè)階段,即逆風(fēng)爬升階段、高空轉(zhuǎn)彎階段、順風(fēng)潛水階段和低空轉(zhuǎn)彎階段。勢(shì)能在逆風(fēng)爬升階段增加,而動(dòng)能在順風(fēng)潛水階段增加[1]。部分文獻(xiàn)表明,風(fēng)能主要在逆風(fēng)爬升階段和高空轉(zhuǎn)彎早期獲得,也有一些在順風(fēng)潛水階段獲得能量;而部分文獻(xiàn)表明,只有逆風(fēng)爬升和順風(fēng)潛水才是風(fēng)能獲取的階段。

圖1 動(dòng)態(tài)滑翔的典型飛行軌跡

在梯度風(fēng)場(chǎng)中,可以設(shè)計(jì)無(wú)人機(jī)的飛行軌跡來(lái)優(yōu)化無(wú)人機(jī)系統(tǒng)的總能量,考慮到無(wú)人機(jī)的橫向運(yùn)動(dòng),最優(yōu)飛行軌跡較為復(fù)雜[2]。參考文獻(xiàn)[3]分析了無(wú)人機(jī)在無(wú)動(dòng)力飛行中產(chǎn)生最佳軌跡時(shí)受不同參數(shù)的影響。參考文獻(xiàn)[4]對(duì)有關(guān)動(dòng)態(tài)飆升的文獻(xiàn)進(jìn)行了綜述。參考文獻(xiàn)[5]研究了高空動(dòng)態(tài)滑翔的技術(shù)。

在眾多優(yōu)化算法中,偽譜方法是最近發(fā)展起來(lái)的一種求解最優(yōu)控制問(wèn)題的方法。該方法結(jié)合了間接方法和直接方法的優(yōu)點(diǎn):它不僅可以獲得最優(yōu)控制問(wèn)題的數(shù)值解,還可以提供協(xié)變量的精確信息[6-7]。通過(guò)這種方式,連續(xù)系統(tǒng)的最優(yōu)控制問(wèn)題被轉(zhuǎn)化為受一系列代數(shù)約束的非線性規(guī)劃(Nonlinear programming,NLP)問(wèn)題。最常用的偽譜方法包括:勒讓德(Legendre)偽譜方法、高斯(Gauss)偽譜方法和拉道(Radau)偽譜方法。

本文考慮非常規(guī)條件下的動(dòng)態(tài)滑翔軌跡優(yōu)化過(guò)程,采用一種大型無(wú)人機(jī)進(jìn)行動(dòng)態(tài)飛行軌跡優(yōu)化,并用NED坐標(biāo)系代替常用的ENU坐標(biāo)系推導(dǎo)了動(dòng)態(tài)騰空的運(yùn)動(dòng)方程,用Radau偽譜方法進(jìn)行軌跡優(yōu)化。無(wú)動(dòng)力飛行和有動(dòng)力飛行的最優(yōu)軌跡都是通過(guò)能量?jī)?yōu)化生成的。

1 動(dòng)態(tài)滑翔的數(shù)學(xué)描述

1.1 風(fēng)場(chǎng)模型和慣性坐標(biāo)系

(1)

因此,為了獲得與其他文獻(xiàn)相一致的等效表達(dá)式,假設(shè)風(fēng)場(chǎng)的方向?yàn)閥軸方向,即Vwx=Vwz=0.梯度風(fēng)場(chǎng)通常有三種形式的模型,即對(duì)數(shù)形式、指數(shù)形式和線性形式,Vwy滿足以下形式:

(1)對(duì)數(shù)模型

(2)

其中:Vwy(h)是在高度h時(shí)的風(fēng)速;Vref是參考高度href時(shí)的風(fēng)速;h0是校正系數(shù),代表地面的起伏程度。

(2)指數(shù)模型

(3)

(3)線性模型

在公式(3)中, 令p=1可得線性化模型:

(4)

1.2 動(dòng)態(tài)滑翔的運(yùn)動(dòng)方程

在本節(jié)中考慮了無(wú)動(dòng)力飛行的情況,使用了無(wú)人機(jī)在飛行軌跡坐標(biāo)系中的質(zhì)點(diǎn)模型,縱向和橫向運(yùn)動(dòng)都被考慮在內(nèi)。設(shè)γ是爬升俯沖角,φ是偏航角,μ是航跡滾轉(zhuǎn)角,它們?cè)陲w行軌跡坐標(biāo)系中。L是無(wú)人機(jī)的升力,D是無(wú)人機(jī)的阻力,m是無(wú)人機(jī)的質(zhì)量。根據(jù)牛頓第二定律和機(jī)身坐標(biāo)系與慣性坐標(biāo)系之間的轉(zhuǎn)換規(guī)則,可得:

(5)

(6)

(7)

因此,無(wú)人機(jī)在風(fēng)場(chǎng)中的運(yùn)動(dòng)方程可以表示為:

(8)

(9)

(10)

通過(guò)將公式(8)至公式(10)分別對(duì)時(shí)間t求導(dǎo),可得:

(11)

(12)

(13)

(14)

(15)

(16)

再對(duì)公式(14)和公式(15)兩邊分別乘以cosφ和sinφ, 得到:

(17)

(18)

然后,將公式(17)和公式(18)兩邊分別相加,可得:

(19)

公式(19)兩邊同乘以cosγ,可得:

(20)

公式(16)兩邊同乘以sinγ,得到:

(21)

公式(20)減去公式(21):

(22)

即:

(23)

公式(19)兩邊同乘以sinγ,公式(16)兩邊同乘以cosγ, 可得:

(24)

(25)

公式(24)和公式(25)相加:

(26)

即:

(27)

最后,把公式(23)和公式(27)替換到公式(14)和公式(15),可得:

(28)

因?yàn)槭窃贜ED坐標(biāo)系中,有h=-zi,因此公式(10)改寫為

(29)

公式(8)、(9)、(29)、(23)、(27)、(28)這六個(gè)微分方程就構(gòu)成了狀態(tài)空間方程,狀態(tài)變量為[xi,yi,h,Vt,γ′,φ]:

(30)

2 偽譜優(yōu)化算法的設(shè)計(jì)

Radau偽譜方法(RPM)求解最優(yōu)控制問(wèn)題的基本原理是:將未知狀態(tài)變量和控制變量離散在一系列Legendre-Gauss-Radau(LGR)點(diǎn)上;構(gòu)造全局插值多項(xiàng)式來(lái)逼近狀態(tài)變量和控制變量;然后,用狀態(tài)變量的導(dǎo)數(shù)代替動(dòng)態(tài)微分方程[8],具體設(shè)計(jì)過(guò)程如下。

2.1 時(shí)域區(qū)間轉(zhuǎn)換

偽譜方法最常見(jiàn)的三種點(diǎn)集配置方法是:Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), Legendre-Gauss-Lobatto (LGL)[9],其中LG、LGR和LGL分別對(duì)應(yīng)于開放區(qū)間τ∈(-1,1)、半閉區(qū)間τ∈[-1,1)和(-1,1]、閉區(qū)間τ∈[-1,1]。無(wú)論哪一組點(diǎn),都需要實(shí)現(xiàn)時(shí)域區(qū)間從t∈[t0,tf]到τ∈[-1,1]。將時(shí)域區(qū)間t∈[t0,tf]劃分成K個(gè)網(wǎng)格區(qū)間[tk-1,tk],k=1,…,K,其中,t0,…,tK是網(wǎng)格點(diǎn),并且滿足:

t0=t1<…

(31)

在每個(gè)網(wǎng)格區(qū)間t∈[tk-1,tk]中,令

(32)

由公式(32),區(qū)間t∈[tk-1,tk]被轉(zhuǎn)換到τ∈[-1,1],且有:

(33)

2.2 最優(yōu)控制策略

令x(k)(τ)和u(k)(τ)分別表示狀態(tài)變量和控制變量的第k個(gè)網(wǎng)格點(diǎn),且τ∈[-1,1]。代價(jià)函數(shù)的表達(dá)式如下:

(34)

具有動(dòng)態(tài)約束、不等式約束和邊界條件的方程為:

(35)

C[x(k)(τ),u(k)(τ),τ,tk-1,tk]≤0

(36)

Φ[x(1)(-1),t0,x(K)(1),tf]=0

(37)

2.3 多項(xiàng)式逼近

在處理RPM的連續(xù)時(shí)間最優(yōu)控制問(wèn)題時(shí),有必要根據(jù)以下規(guī)則使用拉格朗日多項(xiàng)式來(lái)近似每個(gè)網(wǎng)格區(qū)間k∈[1,…,K]上的連續(xù)最優(yōu)問(wèn)題:

(38)

(39)

(40)

2.4 多項(xiàng)式求導(dǎo)

在用全局插值多項(xiàng)式逼近狀態(tài)變量之后,狀態(tài)變量的導(dǎo)數(shù)可以用拉格朗日多項(xiàng)式關(guān)于參數(shù)τ的導(dǎo)數(shù)來(lái)近似。因此,動(dòng)態(tài)微分方程被轉(zhuǎn)換為代數(shù)約束。對(duì)公式(38)中的X(k)(τ)關(guān)于參數(shù)τ求導(dǎo):

(41)

動(dòng)態(tài)約束方程公式(35)通過(guò)公式(41)的Nk個(gè)LGR點(diǎn)進(jìn)行配置:

(42)

(43)

(44)

為了確保網(wǎng)格點(diǎn)的連續(xù)性,還需要滿足以下等式:

(45)

由代數(shù)約束方程公式(42)至公式(45),利用Radau偽譜近似方法,將連續(xù)最優(yōu)控制問(wèn)題轉(zhuǎn)化為一個(gè)最小化代價(jià)函數(shù)公式(41)的非線性規(guī)劃問(wèn)題。

3 仿真分析

3.1 最小風(fēng)場(chǎng)無(wú)動(dòng)力滑翔

在無(wú)動(dòng)力飛行中,無(wú)人機(jī)的運(yùn)動(dòng)方程如公式(30)。狀態(tài)變量和控制變量分別為[xi,yi,h,Vt,γ,φ]和[CL,μ]。阻力系數(shù)CD沒(méi)有額外考慮,因?yàn)镃D和升力系數(shù)CL之間的存在如下關(guān)系:

(46)

其中,CD0是寄生阻力系數(shù)。Rmax是最大升阻比,K是誘導(dǎo)阻力系數(shù),它是CD0和Rmax的函數(shù):

(47)

(48)

表1 樣本無(wú)人機(jī)的參數(shù)值

表1中,xi(tf)、yi(tf)和φ(tf)的邊界條件并未給出,這是因?yàn)橥ㄟ^(guò)這三個(gè)約束條件的不同設(shè)置,可以獲得以下三種特殊的飛行狀態(tài)。

(1)飛行狀態(tài)1

邊界約束:φ(tf)=φ(t0)+π/2,xi(tf)和yi(tf)無(wú)約束。飛行軌跡被稱為動(dòng)態(tài)翱翔的基本軌跡,航向的變化是90°,φ的范圍是[-π/2,π/2]。限制φ范圍的目的是設(shè)置搜索區(qū)間,這個(gè)范圍的合理選擇將大大增加找到最佳軌跡的機(jī)會(huì)。

(2)飛行狀態(tài)2

邊界約束:φ(tf)=φ(t0)+π/2,xi(tf)=xi(t0),yi(tf)=yi(t0),飛行軌跡近似為圓形,也稱為游蕩軌跡。φ的范圍是[-π,π].

(3)飛行狀態(tài)3

邊界約束:φ(tf)=φ(t0)+π,xi(tf)和yi(tf)無(wú)約束。 這個(gè)飛行軌跡是基本軌跡的一個(gè)變體,初始和最終航向方向正好相反,即航向方向的變化180°,φ的范圍是[-π/2,3π/2]。

三種飛行狀態(tài)的最佳飛行軌跡如圖2所示。圖2中的軌跡與三種飛行狀態(tài)的邊界條件完全一致。從圖2中可以看出,高度的范圍很大,這是因?yàn)橛幸环N假設(shè),即風(fēng)場(chǎng)的強(qiáng)度隨著海拔高度的增加而增加。由于無(wú)人機(jī)體積較大,它需要飛得更高才能獲得更多的風(fēng)能。無(wú)人機(jī)的初始速度也很大,因此前方飛行距離也比通常情況下大。

圖2 三種飛行狀態(tài)的最佳飛行軌跡

在上述三種飛行狀態(tài)中,系統(tǒng)的狀態(tài)變量為[xi,yi,h,Vt,γ,φ],控制變量為[C1,μ]。[xi,yi,h,Vt]、[γ,φ]、[C1,μ]的變化曲線和系統(tǒng)的能量變化分別如圖3(a)～5(a)、3(b)～5(b)、3(c)～5(c)和3(d)～5(d)所示。通過(guò)比較圖3至圖5中的曲線,可以看到真空速Vt在三種狀態(tài)下變化趨勢(shì)相同, 變化范圍為40～130 m·s-1。爬升角/下滑角γ在-50°到50°之間,γ的變化趨勢(shì)相同。φ的范圍為邊界約束,分別是π/2,2π和π。升力系數(shù)CL在飛行狀態(tài)1中達(dá)到最大值,大約為1.3,而CL在飛行狀態(tài)3中小于0.9,飛行狀態(tài)2中的CL大小在其他兩種狀態(tài)之間。μ的變化趨勢(shì)和CL相似,μ大部分時(shí)間保持正值,這意味著無(wú)人機(jī)以正方向飛行,這與圖2一致,圖2中的三個(gè)軌跡都沿著x軸的正方向增加。此外μ被限幅在±75°之間且μ在三種飛行狀態(tài)下都達(dá)到極限值并保持一段時(shí)間。 三種飛行狀態(tài)的能量變化趨勢(shì)大致相同,重力勢(shì)能先增大后減小,動(dòng)能先減小后增大,總能量趨于恒定。

圖3 飛行狀態(tài)1的變量曲線

圖4 飛行狀態(tài)2的變量曲線

圖5 飛行狀態(tài)3的變量曲線

3.2 固定風(fēng)場(chǎng)下的最小功率滑翔

考慮具有恒定推力的動(dòng)力飛行,公式(23)可改寫為

(49)

圖6 和最小推力T的關(guān)系

圖7 三種取值情況下的最優(yōu)軌跡

圖8 恒定推力下的各變量曲線

4 結(jié)論

本文提出了一種無(wú)人機(jī)動(dòng)態(tài)翱翔飛行軌跡優(yōu)化方法,給出了獲得無(wú)人機(jī)動(dòng)態(tài)翱翔運(yùn)動(dòng)方程的具體步驟;選擇了梯度風(fēng)場(chǎng)模型和慣性坐標(biāo)系,使這些方程與大多數(shù)文獻(xiàn)的表達(dá)式一致;采用Radau偽譜方法對(duì)飛行軌跡進(jìn)行優(yōu)化。由于所選無(wú)人機(jī)比普通無(wú)人機(jī)大,因此在仿真中需要正確選擇參數(shù)。仿真結(jié)果表明,該方法可以成功地生成無(wú)人機(jī)動(dòng)態(tài)飛行的最優(yōu)軌跡。此外,通過(guò)對(duì)無(wú)動(dòng)力飛行和有動(dòng)力飛行的比較,發(fā)現(xiàn)最佳發(fā)動(dòng)機(jī)推力和風(fēng)場(chǎng)強(qiáng)度近似成反比。