張笑妍，程昊宇，韓博，閆杰

（1.西北工業(yè)大學(xué)無人系統(tǒng)技術(shù)研究院，西安 710072；2.中國人民解放軍93525部隊，日喀則 857000）

0 引言

高超聲速飛行器（Hypersonic flight vehicle，HFV）由于具備快速的全球打擊能力和靈活的機動性等特點，受到研究人員的廣泛關(guān)注［1-2］。HFV一般在距地表約100～120 km 的高度進入大氣層，整個彈道通常在約20～30 km 的高度擴散到終端區(qū)域［3］。再入軌跡的優(yōu)化在HFV 性能中起著重要作用，其產(chǎn)生的控制指令能夠引導(dǎo)飛行器從初始狀態(tài)飛向特定的終端狀態(tài)，確保HFV 能夠精確和快速地到達目標點［4］。

HFV 再入段有大包線、多約束飛行等特點，使得再入?yún)⒖架壽E的設(shè)計優(yōu)化尤其具有挑戰(zhàn)性。HFV除了會在復(fù)雜環(huán)境中受到大量約束，如終端條件約束、狀態(tài)邊界約束、熱流速率約束、動壓約束以及法向過載約束等［5］，其再入軌跡優(yōu)化還需要根據(jù)任務(wù)要求滿足終端時間最短、能量最小或終端速度誤差最小等性能指標。因此，再入軌跡優(yōu)化實際上是一個多約束優(yōu)化問題［6］。在現(xiàn)有文獻中，求解高超聲速飛行器再入軌跡優(yōu)化問題的方法一般可分為2 大類：數(shù)值方法和啟發(fā)式優(yōu)化算法。近年，學(xué)者通過數(shù)值法中的直接方法對作為最優(yōu)控制問題的HFV 的再入軌跡優(yōu)化進行了大量研究［7］。直接方法，特別是偽譜法，可以將高超聲速飛行器軌跡優(yōu)化問題轉(zhuǎn)化為非線性規(guī)劃（Nonlinear programming problem，NLP）問題，直接通過傳統(tǒng)方法解決，而不需要推導(dǎo)出一階必要條件。相較于其他直接方法，偽譜法求解速度快且求解精度高，被廣泛應(yīng)用于復(fù)雜多約束優(yōu)化問題的求解中。Patterson 提出了高斯偽譜法（Gauss pseudospectral method，GPM），用于解決多階段最優(yōu)控制問題［8］。利用GPM 解決HFV 再入段軌跡優(yōu)化問題，可以快速獲得滿足約束條件的軌跡。但實際解決HFV 再入段軌跡優(yōu)化問題時，約束條件多、優(yōu)化問題的復(fù)雜性強，由于GPM 全局搜索能力差，其性能對初始估計值的依賴性會增強［9］。初始值選擇不當會造成算法陷入局部最優(yōu)，甚至發(fā)散，這一缺點使得調(diào)參難度大，軌跡優(yōu)化結(jié)果受設(shè)計者經(jīng)驗影響。

針對以上問題，廣大學(xué)者利用啟發(fā)式智能優(yōu)化算法簡單、靈活、全局搜索能力強等優(yōu)點，開展軌跡優(yōu)化算法研究，以克服GPM 算法過于依賴初始值、全局搜索能力不足的缺點。對于復(fù)雜優(yōu)化問題，啟發(fā)式算法在全局搜索方面具有顯著優(yōu)勢，但在全局最優(yōu)附近的收斂速度和收斂精度相對于數(shù)值方法較差［10-11］。為了充分發(fā)揮數(shù)值方法和啟發(fā)式優(yōu)化算法各自的優(yōu)點，大量學(xué)者對混合算法進行了研究。文獻［12］將鯨魚優(yōu)化算法和高斯偽譜法相結(jié)合，為偽譜法提供更好的初值。文獻［13］改進了麻雀優(yōu)化算法并結(jié)合控制變量參數(shù)化方法設(shè)計混合算法，應(yīng)用于高超聲速飛行器的軌跡優(yōu)化。上述方法在處理約束時均將約束轉(zhuǎn)化為懲罰函數(shù)的形式，參數(shù)復(fù)雜、對研究人員的經(jīng)驗依賴性強，且上層中算法為GPM 提供的初始解效能較差，使混合算法收斂精度不足，魯棒性不強。

基于擴散限制聚集理論，可將利希滕貝格圖（Lichtenberg figure，LF）［14］應(yīng)用于利希滕貝格算法（Lichtenberg algorithm，LA）中。與許多文獻中的算法不同，LA 是一種成功地利用了分形的混合算法，它由基于群體和軌跡的搜索方法組成，這種組成模式能夠提高迭代過程的穩(wěn)定性，在探索和利用之間產(chǎn)生良好的平衡，具有更高的收斂精度和魯棒性。該算法在復(fù)雜多約束優(yōu)化問題中具有優(yōu)越的求解性能［15］，適合應(yīng)用在高超聲速飛行器再入軌跡優(yōu)化問題中。

為提升高超聲速飛行器再入軌跡優(yōu)化問題的求解效率和精度，本文首先提出了一種自適應(yīng)分段利希滕貝格算法（Adaptive piecewise Lichtenberg algorithm，APLA），其在解決帶約束的優(yōu)化問題方面具有明顯優(yōu)勢。算法針對所解決問題的特點，基于拉丁超立方抽樣（Latin hypercube sampling，LHS）［16］改進了初始化方法，最大化提高算法前期收斂速度和全局搜索能力。在迭代階段采用了全局至局部的自適應(yīng)分層策略分階段維護多樣性，提升搜索精度。提出的策略使APLA 在保證收斂速度的前提下提升算法的全局搜索能力。接著，將APLA 和GPM結(jié)合提出再入軌跡優(yōu)化混合算法即APLA_GPM 算法，有效減少了參數(shù)的數(shù)目，降低了參數(shù)的調(diào)節(jié)難度，提高了算法的收斂速度、求解精度以及魯棒性。首先，利用GPM 算法將軌跡優(yōu)化問題轉(zhuǎn)化為NLP 問題，繼而利用具備較好全局搜索能力的APLA 算法作為初始化搜索器，求出再入過程的近似最優(yōu)解序列。最后，將前一搜索階段中，APLA 獲得的近似最優(yōu)解作為GPM 中NLP 求解器的初始解。利用GPM優(yōu)越的收斂速度和更高的搜索精度，在近似最優(yōu)解附近尋找精確的最優(yōu)解。

1 高超聲速飛行器再入軌跡優(yōu)化問題建模

1.1 高超聲速飛行器動力學(xué)與約束

考慮高超聲速飛行器再入段動力學(xué)模型：

式 中：狀態(tài)量x=[r，?，θ，v，γ，ψ]T；控制量u=[α，σ]T；t表示時間；m是飛行器的質(zhì)量；g為引力常量；r為地心距離；γ為航跡傾角；v為速率；?為經(jīng)度；θ為緯度；ψ為方位角；α為攻角；σ為傾側(cè)角；D和L分別是氣動阻力和升力大小，可以表示為：

式中：ρ為大氣密度；ρ0為海平面空氣密度；Re代表地球半徑；hd是大氣密度標高；CD和CL分別為阻力系數(shù)和升力系數(shù)；M是馬赫數(shù)；S是HFV的參考面積。

式中：kQ是一個常數(shù)，與HFV的結(jié)構(gòu)特性相關(guān)。

在建立HFV 模型的基礎(chǔ)上，選取攻角α及傾側(cè)角σ作為控制量即u=[α，σ]T，需要滿足控制約束（式（5））。狀態(tài)量x=[r，?，θ，v，γ，ψ]T需滿足邊界約束、初值約束和終端約束（式（6））：

式中：“≤”表示每個對應(yīng)元素都滿足“≤”關(guān)系；Umin=[αmin，σmin]T和Umax=[αmax，σmax]T分別代表控制量的下邊界和上邊界；x0=[r0，?0，θ0，v0，γ0，ψ0]T為狀態(tài)量的初值；xf=[rf，?f，θf，vf，γf，ψf]T為狀態(tài)的終值；Xmin和Xmax分別為狀態(tài)的下邊界和上邊界；t0和tf分別表示初始時間和終端時間。

1.2 高超聲速飛行器再入任務(wù)代價函數(shù)

考慮如下高超聲速飛行器再入段軌跡優(yōu)化問題，將最大化縱程軌跡及軌跡平滑作為任務(wù)，則可將問題轉(zhuǎn)化為如下所示的NLP問題：

式中：τ為權(quán)重因子；Gmax=[，Pmax，Nmax]T為約束量的上邊界矩陣；Gmin=[，Pmin，Nmin]T為約束量的下邊界矩陣；g(x)=[，Pd(x)，nL(x)]T。

2 自適應(yīng)分段利希滕貝格算法

2.1 利希滕貝格算法

LF 以矩陣F的形式表示。在LA 中構(gòu)建LF 需要3個重要的參數(shù)，粒子數(shù)量Np、粘性系數(shù)Sp和創(chuàng)建半徑Rc。其中Rc決定矩陣的行列。粒子在整個矩陣中隨機釋放，將其行走隨機地、徑向地繪制出來。集群中的每個粒子都可以轉(zhuǎn)化為笛卡爾平面上的位置，并將位置四舍五入為一個具有行號和列號的矩陣元素。Sp對圖的密度有很大影響，Sp越小，粒子粘在某個集群上的概率越小，集群的密度就會增加。

LA用有限數(shù)量的點來映射搜索空間，以便在目標函數(shù)中進行評估。其中，前一次迭代的最佳點Pb總是當前迭代的觸發(fā)點Xs，但它不一定是要再次評估的點。在第一次迭代中，它是隨機選擇的LF點中的一個。每次迭代時在整個LF 中選擇代表群體的LF 特征點，稱為利希滕貝格點（Lichtenberg point，LP），并保證這些LP 一定在規(guī)定的搜索空間內(nèi)。這種形式使LA 兼?zhèn)淙后w和軌跡2種搜索方式，成為一種混合算法。LP的生成公式如下所示：

2.2 改進利希滕貝格算法

LA算法存在兩方面缺點，一方面是算法初始化階段粒子在空間內(nèi)分布不均勻，導(dǎo)致全局搜索能力差且收斂速度慢。另一方面是迭代搜索過程中，算法前期全局搜索能力強但收斂速率慢，后期粒子多樣性差，全局搜索性能降低。針對LA 算法的不足，及高超聲速飛行器再入軌跡優(yōu)化問題多約束、強耦合的特點對LA 算法進行了改進，提出了APLA 算法。通過優(yōu)化初始解的效能加速算法前期的收斂速度并提高前期全局搜索能力。高超聲速飛行器再入軌跡優(yōu)化問題約束多且復(fù)雜性較高，算法易陷入局部最優(yōu)，使求解精度降低。本文提出了全局-局部交叉自適應(yīng)更新方法，在保證算法全局搜索性能的前提下，提高算法收斂速度。

圖1 是APLA 算法的流程。其中，行向量LjP為矩陣LP的第j行；Pb(i)為第i次迭代中的最優(yōu)位置。

圖1 APLA流程圖Fig.1 APLA flowchart

2.2.1 拉丁超立方抽樣初始化

初始觸發(fā)點需要在算法迭代前選出，它的質(zhì)量影響著APLA 后續(xù)迭代的效率。良好的初始觸發(fā)點能夠加速算法收斂，提升算法前期解的質(zhì)量。本文通過引入LHS 改進粒子位置初始化，如式（14）～（16）所示，提升觸發(fā)點備選解在搜索空間中分布的均勻性，從而提升初始觸發(fā)點的全局優(yōu)越性。LHS是一種分層抽樣技術(shù)，它使用多變量參數(shù)分布的近似隨機抽樣方法來創(chuàng)建初始位置，適用于本文研究的高超聲速飛行器再入軌跡優(yōu)化問題。

2.2.2 全局——局部交叉自適應(yīng)更新

在保證算法全局搜索效率的同時，為了加快算法的收斂速度，設(shè)計自適應(yīng)參數(shù)，收縮率ξ∈[0，1]，用以加速局部搜索。通過ξ創(chuàng)建副圖，其與主圖具有相同的觸發(fā)點，使算法能夠在一個較小的局部區(qū)間內(nèi)快速搜索尋優(yōu)，從而加速算法的收斂。為了在提高收斂速率的基礎(chǔ)上防止算法早熟，在迭代在前期ξ是一個較大的值，后期粒子分布較集中，收縮率將隨著迭代的進行不斷減小，增加主圖所占比例。ξ的更新公式如下所示：

式中：M表示算法最大迭代次數(shù)；i表示當前迭代數(shù)；Δ ∈(0，0.2)是正常數(shù)，根據(jù)實際問題設(shè)置。

在算法迭代的過程中，用于計算目標函數(shù)的粒子由取自主圖的粒子（LPG）和取自副圖的粒子（LPL）組成。LPG 主要用于提升粒子多樣性，防止算法陷入局部最優(yōu)。LPL用于在當前最優(yōu)點附近搜索，加速算法收斂。前期多樣性較為豐富且ξ較大，LPL 應(yīng)占主導(dǎo)地位。后期粒子多樣性減少，應(yīng)增加LPG 所占比例。綜上，需要對2 類粒子的數(shù)量進行設(shè)計：

式中：Pg為從主圖選取的粒子數(shù)；i表示當前迭代數(shù)；τ∈[1，1.7]是正常數(shù)，根據(jù)實際問題設(shè)置。

綜上，經(jīng)過改進的粒子更新公式如下所示：

圖2為APLA迭代中粒子分布示意圖，x1，x2為粒子在分布空間中的坐標分量。對比圖2（a）與圖2（b）可知：算法前期局部LF 占比大，LP 在空間中分布均勻，算法全局搜索能力強；后期局部LF 占比逐漸縮小，增強算法局部搜索能力并加快收斂，LP 分布在全局LF 中的比例增加，保證種群多樣性，防止算法陷入局部最優(yōu)。

圖2 利希滕貝格圖中的LP分布Fig.2 LP distribution in the Lichtenberg figures

3 APLA_GPM 高超聲速飛行器再入軌跡優(yōu)化方法

3.1 高斯偽譜法

本文通過GPM［17］將時間區(qū)間為[t0，tf]的高超聲速飛行器再入軌跡優(yōu)化問題（式（1）～（7））轉(zhuǎn)化為NLP 問題進行求解，首先需要將其時間區(qū)間轉(zhuǎn)化到GPM的時間區(qū)間[-1，1]上：

進而，微分狀態(tài)為：

式中：D∈R(N-1)×N為微分矩陣，k=1，2，…，N-1。

將原問題的動力學(xué)微分約束轉(zhuǎn)換為在配點上完全等價的代數(shù)約束：

利用高斯積分近似離散化的終端狀態(tài)及代價函數(shù)得到：

另外，終端約束條件及路徑約束條件分別表示為式（28）～（29）：

綜上，GPM 將高超聲速飛行器再入軌跡優(yōu)化問題轉(zhuǎn)化為如式（26）～（29）所示的NLP問題。

3.2 APLA約束處理

經(jīng)典LA 算法最初提出的目的即求解多約束優(yōu)化問題，因此它是求解具有線性和非線性約束實際問題的強大工具。APLA 在經(jīng)典LA 對約束的處理方法的基礎(chǔ)上，根據(jù)高超聲速飛行器再入軌跡優(yōu)化問題的特點優(yōu)化處理參數(shù)，提升算法針對本問題的求解效能。轉(zhuǎn)化后高超聲速飛行器再入軌跡優(yōu)化的目標函數(shù)為：

式中：k為不等式約束的總數(shù)；Hi(x)表示第i個不等式的被約束量的值；Bi為當前約束的邊界值；λ=10為懲罰因子；εi為懲罰標記，當不滿足當前約束時εi=1，滿足約束則εi=0。

3.3 APLA_GPM優(yōu)化方法

本文將高超聲速飛行器再入軌跡優(yōu)化問題轉(zhuǎn)化為NLP 問題。GPM 由于具備求解精度高和搜索速度快等卓越性能，被廣泛用于求解NLP 問題。但其對初始值的依賴性強，調(diào)試難度大，容易陷入局部最優(yōu)甚至發(fā)散無法求解。針對這一缺點，提出APLA_GPM，利用APLA 算法強大的全局搜索能力及優(yōu)秀的收斂速度為GPM 提供良好的初始解，而GPM 則彌補了APLA 在最優(yōu)解附近求解精度較低的缺點。將兩種算法的優(yōu)點結(jié)合，極大地增強了問題的求解效率。圖3 為APLA_GPM 軌跡優(yōu)化算法框架，詳細過程如下：

圖3 軌跡優(yōu)化算法框架Fig.3 Framework of trajectory optimization algorithm

1) 高超聲速飛行器再入軌跡優(yōu)化問題建模（式（1）～（7））；

2) 設(shè)置APLA 所需的參數(shù)，包括Np、Sp、Rc、p、M以及優(yōu)化誤差eA；

3) 初始化LPs求出初始觸發(fā)點，并生成LF；

4) 執(zhí)行APLA 更新LPs 得到當前時刻下的最佳控制量；

5) 保存當前狀態(tài)量及控制量并判斷當前狀態(tài)是否滿足‖x(t) -xf‖2

6) 以APLA 獲取的最優(yōu)解序列作為GPM 的初始解；

7) 利用GPM 的NLP 求解器在初始解附近搜索最佳解，若NLP求解精度達到eG，終止GPM；

8) 輸出最優(yōu)狀態(tài)量及控制量，完成高超聲速飛行器再入軌跡優(yōu)化。

4 仿真校驗

4.1 APLA數(shù)值測試分析

利用5 個帶約束基準優(yōu)化問題（G3 函數(shù)（F1）、Martin 和Gaddy 函數(shù)（F2）、用立方體和直線約束的Rosenbrock 函 數(shù)（F3）、Bird 函 數(shù)（F4）、G2 函 數(shù)（F5））［18］以及焊接梁問題（F6）［19］進行數(shù)值測試，將APLA 的性能與WOA［20］、SSA［21］、PSO［22］以及LA 的性能進行對比分析，以驗證本研究提出算法優(yōu)越的全局搜索能力與收斂精度。APLA中τ越大，Pg的初始值和最大值越大，Pg變化越平緩，Pg的變化范圍越小，本文中τ=1.3，M=100，p=10d。圖4 為各算法應(yīng)用于相應(yīng)測試函數(shù)的收斂曲線。表1為將上述算法應(yīng)用于各個測試函數(shù)并進行蒙特卡洛仿真的結(jié)果，其中最優(yōu)解、平均最優(yōu)解和標準差為各算法獨立運行1 500次得出的解所求。

表1 算法應(yīng)用于各測試函數(shù)中的結(jié)果Table 1 Results of the algorithm applied to each test function

圖4 算法在各測試函數(shù)上的收斂曲線對比Fig.4 Comparison of convergence curves of algorithms on various test functions

如圖4 所示，對比5 種算法在各帶約束問題中的收斂曲線，可以得出：APLA 具有良好的初始解，加速了算法的收斂過程。采用全局-局部交叉自適應(yīng)更新策略，使得APLA 不僅以更小的迭代次數(shù)收斂至最優(yōu)解，且全局搜索性能明顯強于WOA、SSA、PSO以及LA。

進一步，通過分析表1中的數(shù)據(jù)可以看出，本文所提算法在測試實例中的最優(yōu)解和平均最優(yōu)解均優(yōu)于其他算法，且標準差很小。這一結(jié)果證明了APLA 魯棒性強且收斂精度高，適合解決帶約束的優(yōu)化問題。

4.2 APLA_GPM仿真分析

本節(jié)通過將APLA_GPM與傳統(tǒng)GPM、WOAGPM、SSAGPM 以及LAGPM 進行對比，驗證本文所提方法在解決高超聲速飛行器再入軌跡優(yōu)化問題上的優(yōu)越性。

在當前高超聲速飛行器再入軌跡優(yōu)化問題中，算法參數(shù)設(shè)置為τ=1.3，Δ=0.06，M=100，p=100?？刂屏康纳舷藓拖孪薹謩e設(shè)置為Umin=[10，-70]T，Umax=[20，70]T。路徑約束中Qmax=3 000 kW/m2，Pmax=75 kPa，Nmax=2.5 g。高超聲速飛行器再入軌跡參數(shù)設(shè)置見表2，飛行器初始條件與終端條件如表3所示。

表2 高超聲速飛行器再入軌跡參數(shù)Table 2 Re-entry trajectory parameters of hypersonic vehicle

表3 飛行器初始和終端條件設(shè)置Table 3 Initial and terminal condition setting of vehicle

圖5給出了高超聲速飛行器初始化過程中的收斂曲線?？梢钥闯?，APLA 可以用更少的迭代次數(shù)獲得更小的適應(yīng)度值。這一結(jié)果表明相對于其他算法，APLA 可以用更少的時間完成上層中的初始化搜索階段，并切換到下層，為APLA_GPM 提供了尋找全局最優(yōu)解的最佳初始解集。

圖5 再入軌跡優(yōu)化初始化收斂曲線Fig.5 Convergence curves of initialization for the re-entry trajectory optimization

圖6為各算法求解高超聲速飛行器再入軌跡優(yōu)化問題的對比仿真結(jié)果?？梢钥闯?，GPM、WOAGPM、SSAGPM、LAGPM 以及APLA_GPM 5 種算法所求得的解均滿足高度、速度、航跡傾角和方位角終端約束。

圖6 再入軌跡優(yōu)化結(jié)果對比Fig.6 Re-entry trajectory optimization results comparison

進一步由圖6（b）可知，GPM 算法優(yōu)化得到的縱程軌跡最短，而由APLA_GPM 算法求解得到的縱程軌跡最長，LAGPM 次之。證明APLA_GPM 在最大化縱程軌跡任務(wù)中表現(xiàn)優(yōu)越。觀察圖6（d）可得，由APLA_GPM 算法求解的航跡傾角隨時間的變化曲線波動最小，最為平滑，而由WOAGPM 和SSAGPM等算法求解的航跡傾角隨時間的變化曲線波動劇烈，觀察地心距離、速度以及方位角隨時間的變化曲線得到的結(jié)論相同。

通過以上分析可以得出如下結(jié)論：相同條件下，由APLA_GPM 求得的解，在平滑度、精度上都優(yōu)于GPM、WOAGPM、SSAGPM 以 及LAGPM 所 得的解。

為進一步驗證本文提出算法的優(yōu)越性，對GPM、WOAGPM、SSAGPM、LAGPM 以及APLA_GPM進行了蒙特卡羅仿真。每個算法獨立運行1 000次。為了表征各算法所求解的平滑度、精度以及算法的求解速度，統(tǒng)計各算法所得高度、速度、航跡傾角以及方位角隨時間變化函數(shù)的導(dǎo)數(shù)（r′，v′，γ′和ψ′）的標準差平均值、終端經(jīng)度平均值、目標函數(shù)的平均值和最優(yōu)值以及算法平均運行時間，統(tǒng)計數(shù)據(jù)如表4所示。

表4 蒙特卡洛仿真數(shù)據(jù)統(tǒng)計Table 4 Monte Carlo simulation data statistics

分析數(shù)據(jù)，可知APLA_GPM 所解得的高度、速度、航跡傾角以及方位角函數(shù)導(dǎo)數(shù)的標準差均小于其他方法，證明其所得軌跡更加平滑。從各個算法所求得的終端經(jīng)度平均值可得出，本研究提出的算法縱程軌跡最大化能力優(yōu)于其他算法。最后，APLA_GPM 目標函數(shù)的平均值最小且平均運行時間最短，證明其不僅在優(yōu)化精確性和魯棒性方面的表現(xiàn)明顯優(yōu)于GPM、WOAGPM、SSAGPM、LAGPM，且在收斂速度方面也優(yōu)于其他算法。

綜上，通過仿真分析得出，APLA_GPM 算法用于解決本文中提出的高超聲速飛行器再入軌跡優(yōu)化問題具有優(yōu)越性。

5 結(jié)論

針對具有復(fù)雜約束的高超聲速飛行器再入軌跡優(yōu)化問題，本文提出了APLA_GPM 混合方法。該方法利用APLA 快速生成最佳初始解，結(jié)合GPM 強大的局部搜索能力和優(yōu)越的收斂速度在初始解周圍搜索最優(yōu)解。數(shù)值仿真分析結(jié)果表明，APLA 算法不僅收斂速度快且具有優(yōu)越的全局搜索性能，彌補了傳統(tǒng)偽譜法對初值的依賴性強以及全局搜索性能差的缺點。高超聲速飛行器再入軌跡優(yōu)化仿真結(jié)果表明，APLA_GPM 以更高的搜索精度、更快的搜索速度完成優(yōu)化指標且魯棒性優(yōu)于其他算法，解得的軌跡更平滑且縱程軌跡更長。