再入飛行器深度確定性策略梯度制導(dǎo)方法研究

郭冬子, 黃 榮, 許河川, 孫立偉, 崔乃剛,*

(1. 哈爾濱工業(yè)大學(xué)航天學(xué)院, 黑龍江 哈爾濱 150001; 2. 北京控制與電子工程研究所, 北京 100038;3. 中國兵器工業(yè)集團(tuán)航空彈藥研究院, 黑龍江 哈爾濱 150030)

0 引 言

從1911年Thorndike提出效用法則,闡述了生物學(xué)上的情景增強(qiáng)與動作再現(xiàn)之間的關(guān)系,到1954年Minsky首次實(shí)現(xiàn)了“試錯(cuò)學(xué)習(xí)”,并在1961年提出“reinforcement learning”概念,再到1957～1960年Bellman和Howard提出強(qiáng)化學(xué)習(xí)框架基礎(chǔ)——馬爾可夫決策過程(Markov dec process, MDP)及其策略迭代方法,為強(qiáng)化學(xué)習(xí)奠定了扎實(shí)的現(xiàn)代理論研究基礎(chǔ);1956～1994年,動態(tài)規(guī)劃、時(shí)間差分、Q-learning、Saras等算法相繼被提出,進(jìn)一步豐富了強(qiáng)化學(xué)習(xí)的內(nèi)涵和應(yīng)用手段,強(qiáng)化學(xué)習(xí)逐漸進(jìn)入了理論研究與實(shí)際應(yīng)用的“高并發(fā)”熱潮;1992～2018年,不僅有經(jīng)典的神經(jīng)動態(tài)規(guī)劃、置信上限樹、Deep-Q-Network、異步優(yōu)勢執(zhí)行器評價(jià)器(asynchronous advantage actor-crific, A3C)、分布式策略優(yōu)化(distributed proximal policy optimization, DPPO)和深度確定性策略(deep deterministic policy gradient, DDPG)等算法,還有世人盡知的Backgammon游戲、AlphaGo、AlphaGo Zero和人工智能(artificial intelligence, AI)參與的DOTA比賽等,強(qiáng)化學(xué)習(xí)已經(jīng)逐步深入到了社會發(fā)展的各個(gè)行業(yè)[1]。

再入飛行器在稠密大氣層內(nèi)飛行時(shí),由高馬赫引起的氣動加熱、氣動過載和動壓影響嚴(yán)重[2],傳統(tǒng)跟蹤制導(dǎo)方法難以充分利用飛行器的剩余飛行能力以滿足終端約束;可采用強(qiáng)化學(xué)習(xí)方法,將再入飛行器視作智能體agent,通過與飛行環(huán)境的不斷交互,尋找最佳的制導(dǎo)系統(tǒng)模型,以在線確定最佳制導(dǎo)指令,從而有效克服再入飛行過程中內(nèi)外擾動與模型不確定性引起的偏差,滿足終端約束要求。

國內(nèi)外相關(guān)學(xué)者已經(jīng)就強(qiáng)化學(xué)習(xí)方法在飛行器軌跡優(yōu)化與制導(dǎo)控制領(lǐng)域中的應(yīng)用問題,進(jìn)行了廣泛地研究。文獻(xiàn)[3]提出了一種基于強(qiáng)化學(xué)習(xí)與神經(jīng)網(wǎng)絡(luò)的軌跡優(yōu)化智能決策算法,解決了再入軌跡優(yōu)化算法在建模偏差或者擾動下任務(wù)適應(yīng)性較差的問題。文獻(xiàn)[4-10]針對飛行器姿態(tài)穩(wěn)定控制問題,采用了actor-critic、PILCO (probabilistic infe-rence for learning control)等方法進(jìn)行了復(fù)雜干擾條件下的控制器設(shè)計(jì)與控制參數(shù)優(yōu)化研究。文獻(xiàn)[11]針對傳統(tǒng)預(yù)測校正算法迭代預(yù)測再入軌跡時(shí)占用大量計(jì)算資源的問題,提出了一種基于actor-critic強(qiáng)化學(xué)習(xí)的跨周期迭代再入飛行器預(yù)測修正制導(dǎo)方法,對初始條件和飛行參數(shù)的不確定性具有良好的適應(yīng)性,但其制導(dǎo)參數(shù)存在周期性相互影響,并且攻角剖面給定,難以最大程度地利用再入飛行器的剩余飛行能力。文獻(xiàn)[12]提出了基于進(jìn)化算法與強(qiáng)化學(xué)習(xí)算法的路徑規(guī)劃方法,解決了有效探索和超參數(shù)敏感等動態(tài)規(guī)劃問題,但針對的是簡單的“小車爬坡”和“月面末段著陸”等簡單動力學(xué)問題,并不涉及三維空間中受復(fù)雜約束條件的再入飛行軌跡規(guī)劃問題。文獻(xiàn)[13]提出了一種基于深度強(qiáng)化學(xué)習(xí)的航跡規(guī)劃策略自學(xué)習(xí)方法,能夠引導(dǎo)無人飛行器在未知飛行環(huán)境中安全無碰地到達(dá)目的地,但并未針對飛行器動力學(xué)模型和復(fù)雜飛行約束條件進(jìn)行充分考慮。文獻(xiàn)[14]研究了采用深度強(qiáng)化學(xué)習(xí)方法生成過載指令以最短時(shí)間到達(dá)目標(biāo),但忽略了引力的影響,飛行器動力學(xué)模型較為簡單。

本文針對再入飛行器復(fù)雜動力學(xué)模型和飛行約束問題,利用DDPG強(qiáng)化學(xué)習(xí)方法,在三維連續(xù)空間中將再入飛行器制導(dǎo)規(guī)劃問題轉(zhuǎn)化為已知當(dāng)前狀態(tài)和動作的后續(xù)最佳動作與策略的尋找問題,通過大量偏差樣本下的離線網(wǎng)絡(luò)學(xué)習(xí),構(gòu)建滿足在線性能指標(biāo)要求的剖面參數(shù)調(diào)整網(wǎng)絡(luò),最后通過極限拉偏仿真,證明了制導(dǎo)方法的有效性。

1 數(shù)學(xué)模型

1.1 動力學(xué)模型

本文研究對象為大氣層內(nèi)無動力再入飛行器,為更方便地表征飛行器運(yùn)動狀態(tài),在位置坐標(biāo)系(以地心為原點(diǎn),Xp軸在飛行器軌道平面內(nèi)并指向飛行器位置方向,Yp軸與飛行器所處位置處的緯線相切)下建立動力學(xué)方程[15],如下所示:

(1)

(2)

式中:CL和CD為升力系數(shù)和阻力系數(shù),可由攻角α和馬赫數(shù)Ma插值氣動文件進(jìn)行求取;q為動壓;Sref為飛行器氣動參考面積。

(3)

1.2 約束條件

再入飛行器飛行過程受到熱流密度、過載、動壓等復(fù)雜過程約束限制[16]。

(1) 熱流密度約束

采用熱流密度計(jì)算公式:

(4)

式中:Rh為飛行器頭部曲率半徑;Cf為與飛行器特征相關(guān)的常數(shù);海平面大氣密度常量ρ0=1.225 kg/m3。

(2) 過載約束

為保證飛行器結(jié)構(gòu)設(shè)備穩(wěn)定工作,氣動力總過載需要限定在一定范圍內(nèi):

(5)

式中：g為重力加速度；m為質(zhì)量。

(3) 動壓約束

動壓對控制系統(tǒng)橫側(cè)向的穩(wěn)定具有重要影響,需小于限定閾值:

(6)

(4) 終端約束

根據(jù)能量管理或末制導(dǎo)要求,終端目標(biāo)狀態(tài)需滿足如下約束:

(7)

1.3 控制量剖面

攻角α和傾側(cè)角σ為再入飛行器的控制量,攻角設(shè)計(jì)[17]為與最大熱流密度約束和過載約束相關(guān)的隨速度變化的剖面形式:

(8)

式中:αmax為與氣動特性和熱流密度相關(guān)的最大攻角；α2為滑翔飛行段的攻角,通常取為最大升阻比對應(yīng)攻角；V1、V2為分段處對應(yīng)速度大小。

設(shè)計(jì)傾側(cè)角大小|σ|為速度V的函數(shù):

(9)

式中:V0為當(dāng)前速度大小;σ0為當(dāng)前狀態(tài)下滿足準(zhǔn)平衡畫像條件的傾側(cè)角大小;Vf為飛行終端速度大小;σf為終端速度狀態(tài)下滿足準(zhǔn)平衡滑翔條件的傾側(cè)角大小;σmid為速度中點(diǎn)狀態(tài)對應(yīng)的待規(guī)劃傾側(cè)角大小。

傾側(cè)角符號則按照經(jīng)典的航向誤差走廊翻轉(zhuǎn)形式給定,后續(xù)將基于本節(jié)描述的攻角剖面和傾側(cè)角剖面形式,采用強(qiáng)化學(xué)習(xí)方法調(diào)整剖面參數(shù)。

2 制導(dǎo)系統(tǒng)設(shè)計(jì)與關(guān)鍵算法

2.1 制導(dǎo)系統(tǒng)設(shè)計(jì)

再入飛行過程中,攻角/傾側(cè)角剖面對飛行器的制導(dǎo)性能影響關(guān)系較為復(fù)雜,傳統(tǒng)制導(dǎo)方法通常根據(jù)專家經(jīng)驗(yàn)離線調(diào)整設(shè)計(jì)并固定飛行剖面,難以適應(yīng)強(qiáng)擾動下的飛行環(huán)境[18]。本文采用強(qiáng)化學(xué)習(xí)方法,根據(jù)飛行狀態(tài)自適應(yīng)調(diào)整制導(dǎo)剖面參數(shù)。

考慮到各剖面參數(shù)均為在一定區(qū)間內(nèi)的連續(xù)量,若采用深度Q網(wǎng)絡(luò)方法,則需要將參數(shù)離散化,產(chǎn)生維度災(zāi)難[19]。為此,本文采用深度確定性策略梯度方法[20],將剖面參數(shù)調(diào)整視為連續(xù)動作空間強(qiáng)化學(xué)習(xí)問題。建立actor-critic架構(gòu)的智能體,actor部分根據(jù)狀態(tài)st進(jìn)行決策輸出動作μ(st),critic部分則根據(jù)狀態(tài)st及動作at估計(jì)Q值。對于每一部分,分別設(shè)置兩個(gè)結(jié)構(gòu)相同,參數(shù)不同的神經(jīng)網(wǎng)絡(luò),即online網(wǎng)絡(luò)與target網(wǎng)絡(luò)。

在訓(xùn)練過程中,智能體的online-actor網(wǎng)絡(luò)生成剖面參數(shù),并添加動作噪聲,完成對攻角及傾側(cè)角剖面參數(shù)的調(diào)整;飛行器基于調(diào)整后的剖面生成攻角、傾側(cè)角指令并輸入動力學(xué)模型中;動力學(xué)模型返回飛行器下一步狀態(tài)及反映制導(dǎo)結(jié)果的剩余航程、速度誤差。將狀態(tài)轉(zhuǎn)換過程信息,即原狀態(tài)-動作-獎勵-轉(zhuǎn)換狀態(tài)[st,at,rt,st+1]存入經(jīng)驗(yàn)池。按照設(shè)定的訓(xùn)練頻次,隨機(jī)抽取經(jīng)驗(yàn)池中的樣本,訓(xùn)練online網(wǎng)絡(luò)并對target網(wǎng)絡(luò)參數(shù)進(jìn)行軟更新。訓(xùn)練過程如圖1所示。

訓(xùn)練完畢后,提取智能體的online-actor網(wǎng)絡(luò)參數(shù)并固定下來。在飛行過程中根據(jù)當(dāng)前狀態(tài),預(yù)測輸出調(diào)整后攻角/傾側(cè)角剖面參數(shù),飛行器即可根據(jù)剖面輸出攻角/傾側(cè)角指令,進(jìn)行后續(xù)飛行。此外,為進(jìn)一步提高飛行控制指令的合理性與飛行軌跡的平滑性,在由智能體網(wǎng)絡(luò)得到攻角指令αnet的基礎(chǔ)上,進(jìn)一步加入速度與飛行路徑角反饋,并加入角速率限幅:

αh=αnet-kVγ

(10)

(11)

式中:k為調(diào)整系數(shù),通常為一較小的非負(fù)常數(shù);αh為速度與飛行路徑角反饋得到的攻角指令;αup為上一步實(shí)際攻角指令;εα為角速率限幅,可取5°/s。

2.2 再入制導(dǎo)問題的強(qiáng)化學(xué)習(xí)要素

(2) 強(qiáng)化學(xué)習(xí)的動作量為傾側(cè)角剖面參數(shù)σmid及攻角剖面參數(shù)α2,V2,根據(jù)前文確定的飛行走廊,動作空間的取值范圍為σmid∈[6°,70°],α2∈[8°,16°],V2∈[3 500 m/s,4 500 m/s]。

(3) 對于再入制導(dǎo)問題而言,調(diào)整剖面參數(shù)的主要目標(biāo)是減少終端剩余航程誤差及速度誤差。因此,當(dāng)飛行器滑翔段飛行結(jié)束時(shí),可將終端剩余航程誤差及速度誤差的平方和的相反數(shù)作為懲罰項(xiàng):

(12)

為了彌補(bǔ)獎勵函數(shù)過于稀疏的問題,在其他時(shí)刻,過程獎勵函數(shù)如下:

(13)

2.3 強(qiáng)化學(xué)習(xí)關(guān)鍵算法

在強(qiáng)化學(xué)習(xí)訓(xùn)練過程中,智能體采用online網(wǎng)絡(luò)產(chǎn)生動作與環(huán)境交互,考慮到DDPG方法中actor網(wǎng)絡(luò)為確定性策略,若智能體直接采用actor網(wǎng)絡(luò)輸出的動作,將減弱智能體探索動作空間的能力,易陷入局部最優(yōu)。為此,一般在訓(xùn)練過程中,需對actor網(wǎng)絡(luò)產(chǎn)生的動作添加噪聲作為智能體的動作輸出[21]。對于飛行器再入制導(dǎo)這樣的連續(xù)控制問題,若動作噪聲完全隨機(jī),易導(dǎo)致各步動作差別較大,不符合真實(shí)物理過程。為此,在進(jìn)行強(qiáng)化學(xué)習(xí)訓(xùn)練時(shí),采用具有自相關(guān)特性的Ornstein-Uhlenbeck過程[22]作為動作噪聲:

dxt=κo(ηo-xt)dt+σodWt

(14)

式中:κo,ηo,σo均為常數(shù);Wt為標(biāo)準(zhǔn)隨機(jī)維納過程。自相關(guān)特性的動作噪聲可有效提高網(wǎng)絡(luò)的探索效果[23]。則訓(xùn)練過程中實(shí)際輸出的動作為

at=xt+μθoa(st)

(15)

通過帶噪聲的動作,智能體不斷與環(huán)境交互,并將狀態(tài)轉(zhuǎn)換過程信息[st,at,rt,st+1]存入經(jīng)驗(yàn)池。每當(dāng)智能體與環(huán)境交互sonline步時(shí),從經(jīng)驗(yàn)池中抽取ntrain個(gè)樣本進(jìn)行訓(xùn)練,采用梯度下降法[24]更新online-critic網(wǎng)絡(luò)的參數(shù)θoc及online-actor網(wǎng)絡(luò)參數(shù)θoa:

(16)

式中:rt為當(dāng)前步的獎勵值,若當(dāng)前步為終端步,則rt=r終端,否則rt=r過程;αc為online-critic網(wǎng)絡(luò)的學(xué)習(xí)率,αa為online-actor網(wǎng)絡(luò)的學(xué)習(xí)率,ζ為獎勵折扣率。每當(dāng)online網(wǎng)絡(luò)更新starget次時(shí),軟更新各target網(wǎng)絡(luò):

(17)

通過以上步驟,在與環(huán)境交互過程中,完成對各網(wǎng)絡(luò)的參數(shù)更新。在訓(xùn)練結(jié)束后,取出online-actor網(wǎng)絡(luò)μθoa(s),即可用于再入飛行制導(dǎo)剖面的在線自適應(yīng)調(diào)整。

3 實(shí)驗(yàn)與結(jié)果

3.1 仿真條件設(shè)定

仿真采用通用航空飛行器(Common Aero Vehicle, CAV)再入飛行器,總體參數(shù)如表1所示。

表1 飛行器總體參數(shù)

為驗(yàn)證方法性能,設(shè)計(jì)飛行器從赤道上格林尼治零點(diǎn)開始,向正東飛行,相關(guān)的任務(wù)條件如表2所示。

表2 飛行器任務(wù)參數(shù)列表

仿真步長為0.01 s。仿真過程中,輸出傾側(cè)角、攻角指令的制導(dǎo)周期為0.1 s。

再入制導(dǎo)強(qiáng)化學(xué)習(xí)中,更改剖面參數(shù)與環(huán)境交互的周期為1 s。actor網(wǎng)絡(luò)及critic網(wǎng)絡(luò)隱含層均設(shè)置為[200,200,20]的結(jié)構(gòu),強(qiáng)化學(xué)習(xí)訓(xùn)練超參數(shù)如表3所示。

表3 強(qiáng)化學(xué)習(xí)訓(xùn)練超參數(shù)設(shè)置

本文進(jìn)行訓(xùn)練及仿真的計(jì)算機(jī)配置為intel core i7-10700+Nvidia RTX 2060。

3.2 再入制導(dǎo)智能體的強(qiáng)化學(xué)習(xí)訓(xùn)練

基于python3.8+tensorflow2.4編寫強(qiáng)化學(xué)習(xí)訓(xùn)練程序。為了提高強(qiáng)化學(xué)習(xí)智能體的魯棒性,對環(huán)境引入正態(tài)分布隨機(jī)偏差如表4所示。

表4 訓(xùn)練偏差范圍

通過1 500次訓(xùn)練,如圖3所示,獎勵函數(shù)逐漸上升,最終完成收斂。

3.3 制導(dǎo)仿真結(jié)果

在強(qiáng)化學(xué)習(xí)訓(xùn)練完成后,獲得的智能體online-actor網(wǎng)絡(luò)可用于飛行器制導(dǎo)的攻角、傾側(cè)角剖面調(diào)整。為了驗(yàn)證本文提出的強(qiáng)化學(xué)習(xí)再入制導(dǎo)方法的精度及計(jì)算效率,給定極限拉偏條件:

(18)

分別采用線性二次型調(diào)節(jié)(linear quadratic regulator,LQR)制導(dǎo)律及強(qiáng)化學(xué)習(xí)方法進(jìn)行再入制導(dǎo)仿真,結(jié)果如圖4～圖11所示。

在極限偏差情況下,按照終端高度跳出時(shí),采用LQR制導(dǎo)方法最大終端速度誤差為99.326 4 m/s,最大剩余航程誤差為1 997.43 m;采用強(qiáng)化學(xué)習(xí)制導(dǎo)方法最大終端速度誤差為34.223 3 m/s,最大剩余航程誤差為212.408 m。

在仿真過程中,LQR方法單次輸出指令的計(jì)算耗時(shí)為0.001 1 s,而本文提出的確定性策略梯度方法單次輸出指令的計(jì)算耗時(shí)為0.000 9 s,單次調(diào)整剖面參數(shù)的計(jì)算耗時(shí)為0.01 s,可滿足計(jì)算實(shí)時(shí)性要求。

4 結(jié) 論

本文針對再入飛行器復(fù)雜動力學(xué)模型和飛行約束問題,利用深度確定性策略梯度方法,在三維連續(xù)空間中將再入飛行器制導(dǎo)問題轉(zhuǎn)化為已知當(dāng)前狀態(tài)和動作的后續(xù)最佳動作與策略的尋找問題,通過大量偏差樣本下的離線網(wǎng)絡(luò)學(xué)習(xí),構(gòu)建滿足在線性能指標(biāo)要求的剖面調(diào)整動作網(wǎng)絡(luò)。實(shí)驗(yàn)結(jié)果表明:強(qiáng)化學(xué)習(xí)再入制導(dǎo)方法較傳統(tǒng)LQR跟蹤制導(dǎo)方法精度更高,能夠在初始偏差和過程干擾條件下達(dá)到終端速度誤差小于35 m/s,終端剩余航程誤差小于500 m,滿足終端約束要求,并且具有較好的實(shí)時(shí)性。