梯度風場中無人機動態滑翔飛行軌跡優化

朱 熠,李繼廣,郝向宇

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

動態滑翔也稱為梯度風場翱翔。一些鳥類,如信天翁,可以長時間利用梯度風場沿著海岸滑翔,它們不會扇動翅膀。自19世紀80年代,學者們開始研究動態飆升的過程,取得了巨大成就,人們關注的焦點是最優軌道的設計和風能的獲取。

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

圖1 動態滑翔的典型飛行軌跡

在梯度風場中,可以設計無人機的飛行軌跡來優化無人機系統的總能量,考慮到無人機的橫向運動,最優飛行軌跡較為復雜[2]。參考文獻[3]分析了無人機在無動力飛行中產生最佳軌跡時受不同參數的影響。參考文獻[4]對有關動態飆升的文獻進行了綜述。參考文獻[5]研究了高空動態滑翔的技術。

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

本文考慮非常規條件下的動態滑翔軌跡優化過程,采用一種大型無人機進行動態飛行軌跡優化,并用NED坐標系代替常用的ENU坐標系推導了動態騰空的運動方程,用Radau偽譜方法進行軌跡優化。無動力飛行和有動力飛行的最優軌跡都是通過能量優化生成的。

1 動態滑翔的數學描述

1.1 風場模型和慣性坐標系

(1)

因此,為了獲得與其他文獻相一致的等效表達式,假設風場的方向為y軸方向,即Vwx=Vwz=0.梯度風場通常有三種形式的模型,即對數形式、指數形式和線性形式,Vwy滿足以下形式:

(1)對數模型

(2)

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

(2)指數模型

(3)

(3)線性模型

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

(4)

1.2 動態滑翔的運動方程

在本節中考慮了無動力飛行的情況,使用了無人機在飛行軌跡坐標系中的質點模型,縱向和橫向運動都被考慮在內。設γ是爬升俯沖角,φ是偏航角,μ是航跡滾轉角,它們在飛行軌跡坐標系中。L是無人機的升力,D是無人機的阻力,m是無人機的質量。根據牛頓第二定律和機身坐標系與慣性坐標系之間的轉換規則,可得:

(5)

(6)

(7)

因此,無人機在風場中的運動方程可以表示為:

(8)

(9)

(10)

通過將公式(8)至公式(10)分別對時間t求導,可得:

(11)

(12)

(13)

(14)

(15)

(16)

再對公式(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)

因為是在NED坐標系中,有h=-zi,因此公式(10)改寫為

(29)

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

(30)

2 偽譜優化算法的設計

Radau偽譜方法(RPM)求解最優控制問題的基本原理是:將未知狀態變量和控制變量離散在一系列Legendre-Gauss-Radau(LGR)點上;構造全局插值多項式來逼近狀態變量和控制變量;然后,用狀態變量的導數代替動態微分方程[8],具體設計過程如下。

2.1 時域區間轉換

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

t0=t1<…

(31)

在每個網格區間t∈[tk-1,tk]中,令

(32)

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

(33)

2.2 最優控制策略

令x(k)(τ)和u(k)(τ)分別表示狀態變量和控制變量的第k個網格點,且τ∈[-1,1]。代價函數的表達式如下:

(34)

具有動態約束、不等式約束和邊界條件的方程為:

(35)

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

(36)

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

(37)

2.3 多項式逼近

在處理RPM的連續時間最優控制問題時,有必要根據以下規則使用拉格朗日多項式來近似每個網格區間k∈[1,…,K]上的連續最優問題:

(38)

(39)

(40)

2.4 多項式求導

在用全局插值多項式逼近狀態變量之后,狀態變量的導數可以用拉格朗日多項式關于參數τ的導數來近似。因此,動態微分方程被轉換為代數約束。對公式(38)中的X(k)(τ)關于參數τ求導:

(41)

動態約束方程公式(35)通過公式(41)的Nk個LGR點進行配置:

(42)

(43)

(44)

為了確保網格點的連續性,還需要滿足以下等式:

(45)

由代數約束方程公式(42)至公式(45),利用Radau偽譜近似方法,將連續最優控制問題轉化為一個最小化代價函數公式(41)的非線性規劃問題。

3 仿真分析

3.1 最小風場無動力滑翔

在無動力飛行中,無人機的運動方程如公式(30)。狀態變量和控制變量分別為[xi,yi,h,Vt,γ,φ]和[CL,μ]。阻力系數CD沒有額外考慮,因為CD和升力系數CL之間的存在如下關系:

(46)

其中,CD0是寄生阻力系數。Rmax是最大升阻比,K是誘導阻力系數,它是CD0和Rmax的函數:

(47)

(48)

表1 樣本無人機的參數值

表1中,xi(tf)、yi(tf)和φ(tf)的邊界條件并未給出,這是因為通過這三個約束條件的不同設置,可以獲得以下三種特殊的飛行狀態。

(1)飛行狀態1

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

(2)飛行狀態2

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

(3)飛行狀態3

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

三種飛行狀態的最佳飛行軌跡如圖2所示。圖2中的軌跡與三種飛行狀態的邊界條件完全一致。從圖2中可以看出,高度的范圍很大,這是因為有一種假設,即風場的強度隨著海拔高度的增加而增加。由于無人機體積較大,它需要飛得更高才能獲得更多的風能。無人機的初始速度也很大,因此前方飛行距離也比通常情況下大。

圖2 三種飛行狀態的最佳飛行軌跡

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

圖3 飛行狀態1的變量曲線

圖4 飛行狀態2的變量曲線

圖5 飛行狀態3的變量曲線

3.2 固定風場下的最小功率滑翔

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

(49)

圖6 和最小推力T的關系

圖7 三種取值情況下的最優軌跡

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

4 結論

本文提出了一種無人機動態翱翔飛行軌跡優化方法,給出了獲得無人機動態翱翔運動方程的具體步驟;選擇了梯度風場模型和慣性坐標系,使這些方程與大多數文獻的表達式一致;采用Radau偽譜方法對飛行軌跡進行優化。由于所選無人機比普通無人機大,因此在仿真中需要正確選擇參數。仿真結果表明,該方法可以成功地生成無人機動態飛行的最優軌跡。此外,通過對無動力飛行和有動力飛行的比較,發現最佳發動機推力和風場強度近似成反比。