劉曉慧，聶萬勝，楊新壘

（裝備學院，北京，101416）

高超聲速飛行器多約束再入滑翔機動彈道優(yōu)化設計

劉曉慧，聶萬勝，楊新壘

（裝備學院，北京，101416）

建立了升力體再入滑翔飛行器的氣動模型和多約束模型。多約束模型除了包括熱流密度、氣動過載、動壓和終端約束等典型約束外，還建立了更符合實際任務的路徑點和禁飛區(qū)約束模型，并利用路徑點、禁飛區(qū)和終端約束劃分彈道，在各段分別使用高斯偽譜法進行彈道求解，將多段多約束的最優(yōu)控制問題轉換為非線性規(guī)劃問題。改進的準平衡滑翔條件保證了彈道平緩。最后通過Matlab仿真計算驗證了所用分段高斯偽譜法規(guī)劃彈道比傳統(tǒng)的高斯偽譜法具有更精確的優(yōu)化結果和更高的優(yōu)化效率。

高超聲速飛行器；機動彈道；再入滑翔；分段高斯偽譜法

0 引 言

實際任務往往要求高超聲速飛行器能夠平緩再入大氣層和大范圍滑翔機動飛行。再入段彈道平緩可便于控制系統(tǒng)操控[1]，而滑翔機動飛行則可以規(guī)避攔截和進行防御，因此有必要根據(jù)實際需求進行高超聲速飛行器再入滑翔機動彈道的優(yōu)化設計。高超聲速飛行器的彈道優(yōu)化涉及眾多非線性運動方程和約束條件，是一種復雜的非線性多約束最優(yōu)控制問題[2]，需要通過數(shù)值方法進行計算。目前數(shù)值方法主要有兩類：

a）間接優(yōu)化目標函數(shù)，通過變分法或最小值原理將最優(yōu)控制問題轉換為漢密爾頓邊界值問題（Hamilton Boundary Value Problem，HBVP）進行求解。HBVP計算精確度高，但對未知邊界條件極其敏感而且求解復雜困難[3]。

b）直接優(yōu)化目標函數(shù)，如直接打靶法和偽譜法，通過在特定時間節(jié)點對微分方程進行離散，將最優(yōu)控制問題轉換為非線性規(guī)劃（Nonlinear Programming，NLP）問題進行數(shù)值求解。NLP操作簡易，對復雜問題有更高的優(yōu)化效率[4]，但因未計算共態(tài)變量，求解精確度相對會下降。

很多學者就此類問題進行了研究。Shen Zuojun使用準平衡滑翔條件，將包括熱流密度、氣動過載和動壓在內的路徑約束轉換為傾側角約束[5]。Reddien首次使用高斯偽譜法（Gauss Pseudospectral Method，GPM）求解最優(yōu)控制問題，使用全局正交多項式在特定節(jié)點同時離散狀態(tài)量和控制量[4]。Benson詳細推導了積分和微分形式的GPM[6]。Huntington改進了GPM，提出了計算邊界控制量的新方法[7]。Jorris和 Cobb使用GPM 首次研究了路徑點和禁飛區(qū)約束下高超聲速飛行器的三維彈道優(yōu)化[8]。

在實際任務中，飛行器往往需要經過特定任務路徑點并順利避開敵方攔截或勿入區(qū)域。因此，除了典型的熱流密度、氣動過載、動壓、終端條件等約束外，還要考慮路徑點和禁飛區(qū)約束，這給傳統(tǒng)的GPM求解帶來一定難度：

a）到達路徑點的時刻未知，不利于狀態(tài)量和控制量離散節(jié)點的選取；

b）離散節(jié)點的分布在整個積分的兩端密集，中間稀疏，擬合結果會降低禁飛區(qū)附近的精確度；

c）多階次的插值多項式也會導致狀態(tài)量和控制量出現(xiàn)較大的擬合誤差。

本文提出的分段GPM，利用指定路徑點和禁飛區(qū)擬接觸點來劃分彈道，能夠很好地解決上述難點，有效優(yōu)化高超聲速飛行器再入滑翔彈道并滿足包括典型約束和路徑點、禁飛區(qū)約束在內的多種復雜約束。

1 數(shù)學模型

1.1 飛行器氣動模型

飛行器采用升力體構型，氣動外形如圖1所示。

對于高超聲速飛行器，其入口條件取來流狀態(tài)參數(shù)（速度、溫度、壓強），出口條件采用外推方式獲得，飛行器表面滿足無穿透、無滑移的壁面邊界條件[9]。

查標準大氣屬性表得到 35 km處的大氣壓強為558.87 Pa，溫度為237 K，密度為8.2×10-3kg/m3，選用標準k-ε湍流模型分別計算得到Ma為6、8、10和12時飛行器升力系數(shù)CL、阻力系數(shù)CD和升阻比L/D隨攻角α的變化曲線，如圖2所示。

由圖 2可知，最大升阻比出現(xiàn)在α為 6°左右且Ma在10以上，CL、CD和L/D隨Ma變化極小，可近似認為只與α有關，與高馬赫數(shù)無關理論一致[10]。由于彈道再入滑翔段的飛行速度是Ma為10左右，因此可將CL和CD分別擬合為α的一次函數(shù)和二次函數(shù)：

同時采用美國1976標準大氣模型來計算隨高度變化的地球重力加速度g、大氣密度ρ、溫度T和聲速c。

1.2 飛行器滑翔動力學模型

由于高升阻比飛行器其氣動力和地球重力分別約為科氏慣性力和地球自轉引起的慣性離心力的100倍和1 000倍，因此可以假設地球為非旋轉橢球體，坐標系原點O位于飛行器質心，OX軸沿飛行速度方向，OY軸位于包含速度向量的縱向對稱平面內并垂直于OX軸指向上，OZ軸與OXY平面垂直構成右手準則，飛行器的三自由度質點動力學模型如下[11]：

式中V為飛行器速度；m為飛行器質量，m=907 kg；D為飛行器氣動阻力，D= 0.5ρV2SmCD；L為飛行器氣動升力，L= 0.5ρV2SmCL；Sm為飛行器參考面積，Sm= 0.48 m2；θ為飛行路徑角；?為飛行航向角（定義從北向南順時針方向為正）；β為飛行傾側角；λ和ψ分別為地球經度和緯度；H為飛行高度；Re為標準地球平均半徑，Re= 6 478 km；r為飛行器到地心的徑向距離，r=H+Re；LV為飛行航程。對于在臨近空間飛行的高超聲速飛行器，其飛行高度H遠小于地球平均半徑Re，因此航程的導數(shù)可進一步簡化為[12]

1.3 約束模型

1.3.1 控制約束

以飛行攻角α和傾側角β為控制變量，為使飛行器平穩(wěn)飛行，有如下約束方程：

1.3.2 終端約束

為保證在一定的落角（終端路徑角）和速度下將飛行器引導至指定終端目標位置，約束方程為

式中LVf表示指定航程；Hf為終端高度；λf和ψf分別為終端目標的經度和緯度；θf為指定終端路徑角；Vf為終端速度。

為滿足末制導階段初始條件，需要對滑翔飛行的終端高度和速度進行約束，本文假設高度約束為35 km以上，速度約束為3 km/s。

1.3.3 路徑約束

a）物理量約束。

熱流密度氣動過載n和動壓q的約束方程如下：

式中k為常數(shù)，k= 5.21×10-8；為最大熱流密度，根據(jù)文獻[13]，取= 300 W/cm2；nmax為最大氣動過載，nmax= 4g；qmax為最大動壓，qmax= 40 kPa。

b）路徑點約束。

設飛行器在時刻ti到達第i個路徑點時的坐標為（x(ti),y(ti)），指定路徑點的坐標為（xi,yi），則約束方程為[11]

c）禁飛區(qū)約束。

假設禁飛區(qū)為圓柱型區(qū)域，第j個禁飛區(qū)中心位置的坐標為（xcj,ycj），其半徑為Rj，飛行器到第j個禁飛區(qū)中心的距離為 ?xj=x(t) ?xcj， ?yj=y(t) ?ycj。則規(guī)避禁飛區(qū)的約束方程為

2 最優(yōu)彈道求解方法

2.1 準平衡滑翔條件的改進

由于彈道再入滑翔段飛行路徑角較小且變化相對較緩，因此可假設cosθ= 1，，將式（3）的第2個方程化簡為

式（10）稱為準平衡滑翔條件（Quasi-Equilibrium Glide Condition，QEGC）[14]。理論上，只要氣動升力足夠，式（10）能夠使彈道平緩沿直線前行。但實際上，飛行器再入后的速度會隨著時間變小，為了保持準平衡滑翔，就需要更大的攻角，這樣氣動升力才能夠抵消掉重力的影響。然而，攻角達到其上限后將不能再增大，致使可用的升力不足以保持準平衡滑翔，使得，θ迅速減小，這將導致嚴重不穩(wěn)定[15]。因此對式（10）QEGC進行了一定改進后得到：

式中ε為一極小負預選值。這樣，可以使θ在整個滑翔段穩(wěn)定平緩地減小，即使攻角α達到上限，也能夠保證θ減小得更少。

由于CL是α的函數(shù)，QEGC建立了α和β兩個控制變量間的關系。已知β，則α可用β表示，反之亦然。很明顯，QEGC減少了優(yōu)化變量，提高了優(yōu)化效率。

2.2 高斯偽譜法及其改進

2.2.1 傳統(tǒng)高斯偽譜法

彈道優(yōu)化屬于最優(yōu)控制問題，可通過數(shù)值模擬的方法確定近似解，將連續(xù)性最優(yōu)控制問題轉換為NLP問題。其中一種常用方法就是GPM，傳統(tǒng)的GPM求解優(yōu)化問題的基本思路如下[4]。

首先定義狀態(tài)量為x(t) ∈Rn，控制量為u(t) ∈Rm，初始時間為t0，最終時間為tf，則在廣義區(qū)間τ∈ [τ0,τf] =[? 1,1]上，有如下轉換關系：

最優(yōu)控制問題要求通過控制變量使目標函數(shù)最小。選取N個勒讓德-高斯（LG）點τk（k= 1,…,N），加上2個邊界點τ0和τf，總共N+2個離散點。因此通過N+1個拉格朗日插值多項式Ls來逼近狀態(tài)量，同時通過N個拉格朗日插值多項式Lc來逼近控制量：

將式（13）微分求導，得到：

每個拉格朗日多項式在 LG點的導數(shù)都可表示為微分形式的近似矩陣D∈RN(N+1)，該矩陣的各項為

通過該矩陣將動力學約束改寫為代數(shù)約束：

終端狀態(tài)Xf可通過高斯積分用Xk≡X(τk) ∈Rn，Uk≡U(τk)∈Rm和初始狀態(tài)X0表示：

目標函數(shù)J的積分部分也用高斯積分近似：

式中ωk為高斯權重。

邊界約束為

LG點的路徑約束為

至此，目標函數(shù)和代數(shù)約束轉換為NLP問題。

2.2.2 分段高斯偽譜法

分段高斯偽譜法在傳統(tǒng)GPM的基礎上加以改進，以適用于路徑點和禁飛區(qū)約束的彈道優(yōu)化，基本思路為：首先根據(jù)指定路徑點和禁飛區(qū)擬接觸點將整個彈道規(guī)劃分為多個階段；在每個階段使用傳統(tǒng)GPM將狀態(tài)量和控制量分別在 LG點離散，轉換為代數(shù)約束；終端狀態(tài)由初始狀態(tài)和高斯積分表示；相鄰兩個階段之間滿足連續(xù)性連接條件和路徑點及禁飛區(qū)約束條件，保證狀態(tài)量和控制量的連續(xù)性；將所有階段組合共同確定目標函數(shù)達到全局最優(yōu)。

3 仿真計算與分析

設飛行器初始及滑翔段終端條件、對應的路徑點和禁飛區(qū)約束條件分別如表1、表2和表3所示。

表1 飛行器初始端和滑翔段終端條件

表2 路徑點約束條件

表3 禁飛區(qū)約束條件

3.1 無路徑點和禁飛區(qū)約束的彈道

作為對比，先忽略路徑點和禁飛區(qū)約束，目標函數(shù)為J= max(tf)，采用30個LG點來離散最優(yōu)控制問題，使用傳統(tǒng)GPM求解得到的三維彈道如圖3所示。

計算結果顯示，使用傳統(tǒng)GPM能夠得到滿足終端約束條件的再入滑翔段彈道：滑翔終端位置為[56.331 2°W，29.680 3°N]，高度為35.012 km，速度為 3 007.9 m/s；落點位置為[25.095 2°W，25.001 3°N]，速度為800.4 m/s，飛行路徑角65.22°，航向角-31.26°，總飛行時間為2 369 s。從圖3的θ曲線可以看出，在改進的QEGC下，再入和滑翔段θ變化很小，彈道基本保持平緩，便于控制系統(tǒng)操控，利于減少熱燒蝕。

3.2 滿足路徑點約束的彈道

在第3.1節(jié)的算例中，飛行器的初始端和滑翔終端條件不變，在其間增加2個路徑點約束，進而將整個彈道分為初始端至路徑點1、路徑點1至路徑點2、路徑 點2至滑翔終端、滑翔終端至目標4段，相鄰兩段中，前一階段的最終狀態(tài)值即為后一階段的初始狀態(tài)值，并通過路徑點約束將相鄰2個階段的狀態(tài)連接起來。

先使用傳統(tǒng)GPM求解彈道，采用30個LG點進行離散；再使用分段GPM，將彈道分為4段，每段采用7個LG點進行離散，得到的彈道如圖4和圖5所示。

兩種方法均能得到經過路徑點且滿足終端約束條件的再入滑翔段彈道。使用傳統(tǒng) GPM的彈道的路徑點1實際位置為[99.987 0°W，33.485 9°N]，路徑點2實際位置為[70.008 3°W，33.989 0°N]；滑翔終端位置為[58.618 5°W，34.263 7°N]，高度為34.959 km，速度為 2 993.1 m/s；落點位置為[25.078 2°W，25.021 5°N]，速度為 797.0 m/s，飛行路徑角22.41°，航向角-19.11°，總飛行時間為2518 s。使用分段GPM的彈道的路徑點 1實際位置為[99.999 8°W，33.500 1°N]，路徑點 2實際位置為[69.999 8°W，34.100 2°N]；滑翔終端位置為[56.250 0°W，33.976 7°N]，高度為35.006 km，速度為3 092.7 m/s；落點位置為[25.340 2°W，25.011 8°N]，速度為809.9 m/s，飛行路徑角63.28°，航向角-50.36°，總飛行時間為2 476 s。

相比于傳統(tǒng)GPM求解的彈道，分段GPM選擇了更少的 LG點，且只在初始位置、終端位置和路徑點位置分布密集，使得彈道能夠更加嚴格的滿足約束條件，且滑翔終端速度更大，落點誤差更小，飛行時間更短，因此有效的提高了優(yōu)化效率。另外，從圖5中θ與?曲線的對比可以看出，使用分段GPM求解得到的飛行路徑角變化更小，彈道更為平緩。

3.3 滿足路徑點約束和禁飛區(qū)約束的彈道

在第 3.2節(jié)的算例中，分別在初始端至路徑點 1和路徑點2至滑翔終端2段中加入2個禁飛區(qū)。同樣，先使用傳統(tǒng)GPM求解彈道，采用30個LG點進行離散；再使用分段GPM，將彈道分為6段，每段采用5個LG點進行離散，得到結果如圖6所示。

計算結果顯示，傳統(tǒng)GPM未能避開禁飛區(qū)，無法 求解到滿足禁飛區(qū)約束條件的彈道；分段GPM則能夠順利求解到符合各種約束條件的彈道，飛行經過2個路徑點，順利避開了2個禁飛區(qū)并緊貼其邊緣以減少額外能耗。禁飛區(qū)1的接觸點位置為[107.100 3°W，33.010 9°N]，禁飛區(qū)2的接觸點位置為[77.155 5°W，33.823 7°N]；路徑點1實際位置為[100.001 5°W，33.499 6°N]，路徑點 2實際位置為[70.000 9°W，34.099 5°N]；滑翔終端位置為[55.328 1°W，33.910 9°N]，高度為35.015 km，速度為3176.1 m/s；落點位置為[西經25.035 6°W，北緯24.994 2°N]，速度為 814.4 m/s，飛行路徑角 53.31°，航向角為-58.03°，飛行總時間為2 597s。

圖7為使用分段GPM求解得到的路徑點和禁飛區(qū)約束彈道的飛行路徑角θ、航向角?、飛行速度V、熱流密度氣動過載n及動壓q的變化情況。

圖7中熱流密度在飛行初始段和飛行末段均控制在較低水平，這是因為，初始高度處大氣極其稀薄，ρ很小，而飛行末段V大大降低，在滑翔段出現(xiàn)最大值，接近，即達到最大值；相應地，由圖2可知，在飛行初始段先增大到最大值16°后迅速減小并保持在最大升阻比攻角6°附近，就是為了保證彈道平緩，減少熱燒蝕，控制熱流密度。氣動過載n一直控制在較低水平，直到飛行末段接近最大值，這是由末制導段飛行高度下降快且大氣密度急劇增大導致的，動壓q亦同理為逐漸增大并在飛行末段達到最大值；分析對應攻角變化，在飛行末段，α先迅速增大以滿足n和q的約束，之后迅速減小至0.76°，以壓低彈道，保證攻擊目標時具有一定的落角。通過分析可以發(fā)現(xiàn)，熱流密度對飛行器滑翔段的約束更強，氣動過載和動壓則是對飛行末段的約束更強。

4 結 論

本文在準平衡滑翔假設下，提出了一種分段高斯偽譜法來優(yōu)化高超聲速彈道，實現(xiàn)平緩再入和機動滑翔，并且滿足包括典型約束、路徑點約束和禁飛區(qū)約束在內的多種約束條件。該方法優(yōu)點如下：

a）在改進的準平衡滑翔條件假設下，再入段彈道平緩，避免了高度振蕩，便于控制系統(tǒng)操控。

b）使用分段 GPM，根據(jù)指定路徑點將彈道劃分為多段，無需考慮到達每個路徑點的時刻，而是將路徑點視作每段的終端約束條件，再使用GPM求解。相比傳統(tǒng)GPM，分段GPM可以更少的離散節(jié)點完成彈道求解，降低了插值多項式的階次，進而減小了擬合誤差，優(yōu)化結果更精確，優(yōu)化效率更高。

c）使用分段 GPM，根據(jù)禁飛區(qū)的擬接觸點來劃分彈道，在接觸點附近的離散節(jié)點更密集，彈道平滑，可順利避開禁飛區(qū)并通過路徑點，且各種典型約束都能控制在較低水平，而傳統(tǒng)GPM則無法得到可行解。

[1] 和爭春, 何開鋒, 朱國林. 過載限制條件下升力體再入大氣層彈道優(yōu)化設計[J]. 彈道學報, 2007,19(2): 8-12.

[2] 張燁琛, 包為民, 閔昌萬, 張寧寧. 高超聲速無動力滑翔飛行器最優(yōu)彈道研究[J]. 導彈與航天運載技術, 2014(4): 5-8.

[3] 李佳峰, 周浩, 陳萬春. 高超聲速飛行器滑行段最優(yōu)彈道的間接算法[J].飛行力學, 2009,27(3): 41-49.

[4] Reddien G W. Collocation at Gauss points as a discre-tization in optimal control[J]. Control Optim, 1979, 17(2): 298-306.

[5] Shen Zuojun. Onboard generation of three-dimensional constrained entry trajectories[J]. Guidance, Control and Dynamics, 2003,26(1): 111-121.

[6] Benson D A. Direct trajectory optimization and costate estimation via an orthogonal collocation method[J]. Guid Control Dyn, 2006,29(6): 1435-1440.

[7] Huntington G T. Advancement and analysis of a Gauss pseudospectral transcription for optimal control prob-lems dissertation[D]. Massachusetts Institute of Technology, 2007.

[8] Jorris T R, Cobb G C. Three-dimensional trajectory op-timization satisfying waypoint and no-fly zone con-straints[J]. Guidance，Control and Dynamics, 2009,32(2): 543-553.

[9] 丁洪波, 蔡洪, 張士峰, 李安梁. 高起聲速滑翔式再入飛行器最大航程飛行軌跡分析[J]. 國防科技大學學報, 2009,31(6): 67-72.

[10] Philips T H. A common aero vehicle (CAV) model, de-scription, and employment guide[R]. Schafer Corporation for AFRL and AFSPC, 2003.

[11] 黃育秋, 何麟書. 升力式再入飛行器多約束多階段彈道優(yōu)化設計[J].哈爾濱工業(yè)大學學報, 2011,43(7): 144-148.

[12] 方群, 李新三. 臨近空間高超聲速無動力滑翔飛行器最優(yōu)軌跡設計及制導研究[J]. 宇航學報, 2008,29(5): 1485-1491.

[13] 劉欣, 李建林, 葛健全, 楊濤. 滑翔式飛行器再人彈道設計[J]. 彈箭與制導學報, 2011,31(6): 161-164.

[14] 張科南, 周浩, 陳萬春. 高超聲速飛行器多約束多種機動突防模式彈道規(guī)劃[J]. 彈道學報, 2012,23(3): 85-90.

[15] 陳小慶, 侯中喜, 劉建震. 高超聲速滑翔飛行器彈道特性分析[J]. 導彈與航天運載技術, 2011(2): 5-9.

Optimization and Design of Reentry Gliding Maneuvering Trajectory with Multi-Constraints for Hypersonic Vehicle

Liu Xiao-hui, Nie Wan-sheng, Yang Xin-lei
(Equipment Academy, Beijing, 101416)

The aerodynamics model of the reentry gliding vehicle and the multi-constraints model including heating rate, overload, dynamic pressure, way point and no-fly zone constraints are built. The study is based on the assumption of quasi-equilibrium glide condition to keep the trajectory smooth. The trajectory is divided into multi-phases according to the waypoints, no-fly zones and terminal constraints and the optimal trajectory is solved by employing multi-phases gauss pseudospectral method, which transfers the multi-constraints and multi-phases trajectory optimization problem into a nonlinear programming problem. The simulations can well demonstrate that the improved method in this paper is more precise and efficient for the optimization than the traditional method.

Hypersonic vehicle; Maneuvering trajectory; Reentry gliding; Multi-phases Gauss pseudospectral method

V421

A

1004-7182(2017)02-0006-06

10.7654/j.issn.1004-7182.20170202

2016-07-08；

2017-01-10

劉曉慧（1991-），女，碩士研究生，主要研究方向為高超聲速飛行器總體技術