王麗英，張友安，趙國榮

（海軍航空工程學院1.系統(tǒng)科學與數(shù)學研究所；2.控制工程系，山東 煙臺264001）

近年來，偽譜法（p法）和網(wǎng)格細化法（h法）在求解非線性最優(yōu)控制問題的數(shù)值解上得到了廣泛的應用[1－3]。文獻[1]提出了一種p法網(wǎng)格細化策略，缺點是要提前獲知不連續(xù)點和奇點個數(shù)及位置。文獻[2－3]分別提出了一種基于多分辨率的自適應h法軌跡優(yōu)化算法和基于密度函數(shù)的網(wǎng)格點分布算法，優(yōu)點是可以較好地捕捉狀態(tài)變量、控制變量的不連續(xù)性和高階非平滑性；缺點是計算效率低。起源于有限元法的自適應hp法，結合了p法和h法的優(yōu)點，在流體力學和偏微分方程的求解中得到廣泛應用[4－7]。為同時提高偽譜法求解非線性最優(yōu)控制問題的精度和計算效率，文獻[8]將hp思想引入最優(yōu)控制問題的求解中，提出了hp自適應偽譜法，在實施過程中忽略了軌跡的平滑性；文獻[9－10]則給出了一種變階自適應偽譜法，以軌跡的曲率作為改進精度方式的標準，相比文獻[8]，該算法考慮了軌跡的平滑性，但在執(zhí)行h法細化時，是根據(jù)軌跡曲率密度函數(shù)定義新增網(wǎng)格點的個數(shù)和位置，對于無約束或約束較少的情況，該算法很實用，但對于多約束條件下的復雜優(yōu)化問題來說，卻要耗費大量的優(yōu)化時間。

本文在文獻[8－9]的基礎上，基于Radau偽譜法給出了一種改進的hp自適應網(wǎng)格細化算法。以高超再入飛行器軌跡優(yōu)化為例，將該算法和全局偽譜法及局部偽譜法在解的精度、優(yōu)化計算時間、迭代次數(shù)等方面進行了比較，結果表明hp自適應網(wǎng)格細化算法能以較少的計算代價得到較高精度的解。

1 多區(qū)間軌跡優(yōu)化問題闡述

多約束條件下的軌跡優(yōu)化問題可描述為下列非線性最優(yōu)控制問題：定義狀態(tài)－控制變量，使Bolza型代價函數(shù)最小化：

滿足約束條件：

多區(qū)間非線性最優(yōu)控制問題是指將問題（1）、問題（2）的時間區(qū)間[t0，tf]劃分為K個子區(qū)間，用t0＜t1＜…＜tK＝tf表示K＋1個網(wǎng)點，Nk（k＝1，…，K－1）表示第k個子區(qū)間[tk－1，tk]內的配點數(shù)（這里，定義網(wǎng)點為相鄰區(qū)間的連接點；配點為單個區(qū)間內的節(jié)點；整個區(qū)間的節(jié)點數(shù)等于網(wǎng)點數(shù)與配點數(shù)之和）。進行下述時域變換：

將時間t∈[tk－1，tk]轉變成τ∈[－1，＋1]。

對于轉換后的多區(qū)間最優(yōu)控制問題，采用hp－LGR偽譜法作為基本的離散化方法，將其轉化為非線性規(guī)劃問題，具體步驟參見文獻[10]。

2 改進的hp自適應網(wǎng)格細化算法

2.1 誤差評估準則

在數(shù)值求解過程中，一般不能得到原連續(xù)時間最優(yōu)控制問題的確切解，但若給出的數(shù)值求解方法能夠精確地近似原問題，則微分－代數(shù)約束方程應該在任意點上被滿足。因此，以微分－代數(shù)約束方程在配點間的滿足程度作為解的近似誤差評估準則。

取相鄰配點間等間距分布的3個點作為誤差評估準則的采樣點，即

為敘述方便，將[tk－1，tk]內的采樣點定義為，于是微分－代數(shù)約束方程在采樣點上的殘差表示為

式中：l＝1，…，n，表示狀態(tài)變量的個數(shù)；j＝1，…，s，表示路徑約束的個數(shù)；p＝1，…，3（Nk－1），表示采樣點表示第l個狀態(tài)在第p個采樣點處的殘差絕對值表示第j個路徑約束在第p個采樣點處的殘差。式（4）表示動態(tài)約束殘差矩陣；式（5）表示路徑約束殘差矩陣。對于路徑約束來說，當式（5）中的元素全為負值時，表示滿足設定的要求；當存在大于零的元素時，則表示不能滿足設定的要求。

設ε為給定的誤差門限值。于是，當式（4）中的所有元素都小于設定的ε且式（5）中的元素全為負值時，認為達到求解精度要求；反之則需要進一步提高求解精度。下面以第k個子區(qū)間[tk－1，tk]為例，具體介紹改進求解精度要求的措施。

2.2 提高求解精度的策略

2.2.1 路徑約束殘差起主要作用

2.2.2 動態(tài)約束殘差起主要作用

令（xl，p）n×3（Nk－1）表示n個狀態(tài)變量在p個采樣點上的離散值所構成的狀態(tài)矩陣，為矩陣（xl，p）n×3（Nk－1）中相應 于最大 殘差的狀態(tài)所在的行向量，即，則第l個狀態(tài)在第p個采樣點處的曲率為

設表示曲率的算術平均值，即

令

若式（6）中的每個元素都比ρ?。é褳橐辉O定值），則認為該區(qū)間內的曲率具有一致形式，此時，通過增加區(qū)間內插值多項式的維數(shù)來提高求解精度。設表示區(qū)間[tk－1，tk]在第i次迭代時的配點數(shù)，表示第i＋1次迭代時的配點總數(shù)，N（＋）表示新增加的配點數(shù)，則

若式（6）中存在比ρ大的元素，則認為該區(qū)間內存在一個（或多個）較大的曲率，即軌跡相對來說是非平滑的，此時，需要將該區(qū)間作進一步的細化以提高求解精度。

①確定新網(wǎng)點的位置。

取式（6）中元素的局部最大值對應的采樣點作為新增加的網(wǎng)點。

②確定新增加區(qū)間內配點的個數(shù)。

與2.2.1節(jié)相同，在每個新增加的區(qū)間內設定N（0）個配點。

2.3 算法具體步驟

具體算法步驟總結如下。

①初始化。給定誤差門限值ε和初始狀態(tài)x0；選取K個子區(qū)間，每個區(qū)間內設定N（0）個配點。

②在給定的狀態(tài)估計值和節(jié)點分布下求解NLP，若所有區(qū)間的＜ε，則停止迭代，若存在某個ε，則繼續(xù)③；

④在路徑約束的殘差最大處進行網(wǎng)格細化，令N（0）為新增區(qū)間內的配點數(shù)，轉⑥；

⑤若式（6）具有一致形式，則增加該區(qū)間內的配點數(shù)，如式（7）所示；若式（6）具有非一致形式，則進行網(wǎng)格細化，取式（6）中元素的局部最大值對應的采樣點作為新增的網(wǎng)點位置，在每個新增加的區(qū)間內設定N（0）個配點；

⑥在所有的網(wǎng)格區(qū)間都被更新后，返回②。

3 應用實例

下面以再入飛行器軌跡優(yōu)化問題為例，對本文算法進行驗證。

3.1 再入軌跡優(yōu)化問題

再入軌跡優(yōu)化的數(shù)學模型，即詳細的3－D運動方程參見文獻[11]。同時將升力參數(shù)CL和阻力參數(shù)CD表示為攻角α和馬赫數(shù)的函數(shù)[11]：

為限制執(zhí)行機構的變化速率，令＝u1，＝u2，于是將3－D運動方程改寫為向量形式的8狀態(tài)運動方程，其狀態(tài)空間X1和控制空間U1定義如下：

式中：r0為飛行器距地心的距離；θ，φ分別為經(jīng)度、緯度；v為飛行速度；γ，ψ，α，σ分別為飛行路徑角、航向角、攻角、傾側角；u1，u2分別為攻角變化率、傾側角變化率。

以再入航程最大為優(yōu)化目標，即J＝－φ（tf），其中，tf為末端時刻。同時考慮以下的約束條件：

式中：FL為升力，F(xiàn)D為阻力；c＝34 882.985 5，為常數(shù)＝2 453kW/m2；nZ，max＝4，qmax＝335.2kPa。

②終端條件約束。

3.2 仿真分析

仿真中狀態(tài)變量、控制變量的邊界限制如下：

仿真中涉及到的各種參數(shù)設置見表1。

取N（0）＝2，N（＋）＝5，ρ＝2；N0，Nc分別表示不同算法初始化時選取的區(qū)間個數(shù)和每個區(qū)間內的配點數(shù)。為簡化表達形式，以下表述中的p代表全局Radau偽譜法；h代表分段Radau偽譜法；hp代表改進的自適應 hp－Radau偽譜法。E，tc，Cpt，MI，IT，OF分別表示優(yōu)化結束時的最大微分－代數(shù)約束的近似誤差、優(yōu)化計算時間、優(yōu)化結束時的配點個數(shù)、網(wǎng)格區(qū)間個數(shù)、迭代次數(shù)和目標函數(shù)值。整個求解過程是在普通PC機上進行的，其中，CPU為1.73G/Pentimu Dual，內存為DDR 1.0G，操作系統(tǒng)為Windows XP，編程環(huán)境為 Matlab R2009a。

表2給出了3種偽譜方法在不同精度要求下的最大近似誤差、優(yōu)化時間、優(yōu)化結束時所需的配點總數(shù)、子區(qū)間個數(shù)、迭代次數(shù)及代價函數(shù)絕對值間的比較結果，從中可得到以下結論：①不同精度標準下，p法在優(yōu)化過程中得到的微分－代數(shù)約束誤差均最大；②10－3和10－4的精度要求下，h法和hp法在優(yōu)化速度、求解精度上相差不大，相對來說hp法的計算精度最高；③隨著精度要求的提高（10－5），優(yōu)化所需的配點總數(shù)、子區(qū)間的個數(shù)、迭代次數(shù)均逐漸增多，相應地也導致了優(yōu)化時間的增長，符合偽譜法求解最優(yōu)控制問題的特點。同時可以看出，hp法在求解精度和優(yōu)化時間上最具優(yōu)勢。

表1 仿真參數(shù)設置

表2 不同偽譜方法的優(yōu)化結果比較

圖1給出了10－4精度下不同偽譜法的節(jié)點分布特點（圖中縱坐標無實際物理意義，僅為圖示方便，將采用不同方法優(yōu)化得到的節(jié)點分別投影于y＝1，y＝2和y＝3所代表的3條直線上），從圖中可以看出：①p法采用的是一個區(qū)間，通過增加區(qū)間內節(jié)點個數(shù)提高求解精度要求，保持了偽譜法節(jié)點兩端密、中間疏的分布特點；②h法通過細化子區(qū)間的方式提高求解精度要求，所有子區(qū)間內的配點數(shù)是相同的，可看出每3個點構成一個區(qū)間；③hp法的節(jié)點分布明顯不同于其它2種偽譜方法，不具備統(tǒng)一節(jié)點分布標準。

圖1 不同偽譜法的節(jié)點分布

上述分析過程表明：hp法與其他2種優(yōu)化方法相比，能夠以較少的計算代價得到較高精度的解。

圖2、圖3分別給出了本文算法與文獻[8－9]中的算法在求解3.1節(jié)軌跡優(yōu)化問題中得到的部分狀態(tài)軌跡和路徑約束變化曲線；3種hp算法優(yōu)化求解過程所需的優(yōu)化時間、節(jié)點總數(shù)及迭代次數(shù)間的比較結果如表3所示。

圖2 高度變化曲線

取誤差門限值ε＝10－4，由圖可以看出，從最優(yōu)軌跡的平滑性上考慮，文獻[9]的算法最好，其次是本文算法，文獻[8]最差；而由表3則可看出，文獻[9]所需的優(yōu)化時間最長。由于研究目的是實現(xiàn)再入軌跡的快速優(yōu)化，因此，綜合衡量軌跡精度和優(yōu)化時間，本文所給算法結合了文獻[8]和文獻[9]的優(yōu)點，更適合快速求解高超飛行器的再入軌跡優(yōu)化問題。

圖3 動壓變化曲線

圖4和圖5分別為優(yōu)化得到的最優(yōu)控制信號及其隨時間的變化速率曲線。結合兩圖可以看出：①優(yōu)化過程中控制變量的大小及變化速率均在要求的設定范圍內；②攻角大部分時間保持在15°～20°之間，即保持大攻角飛行，最大程度地減少熱載的作用，符合高超熱防護的需求。優(yōu)化得到的三維最優(yōu)軌跡如圖6所示。

圖4 控制信號

圖5 控制信號變化率

圖6 三維最優(yōu)軌跡

為分析hp－LGR偽譜法計算最優(yōu)軌跡的精度，將計算得到的最優(yōu)控制信號帶入原微分方程，用Matlab中的ode45命令進行數(shù)值積分，狀態(tài)變量的積分軌跡與最優(yōu)軌跡間的誤差曲線如圖7，可以看出兩者間的誤差較小；表4為末端狀態(tài)變量實際值與期望值間的比較結果，從中可以看出終端約束均得到滿足，期望的落點位置與實際位置間的誤差很小，表明了hp－LGR偽譜法具有較高的求解精度。

圖7 誤差曲線

表4 末端狀態(tài)誤差

圖8為路徑約束隨時間變化曲線，從圖中可以看出：熱流密度、動壓和過載約束均在設定的范圍內；且熱流密度的峰值出現(xiàn)在再入初期，隨著高度的降低、大氣密度的增大，熱流密度逐漸減小，動壓和過載則逐漸變大，符合高超再入軌跡的特點。

圖8 路徑約束隨時間變化曲線

4 結論

基于hp－LGR偽譜法給出了一種求解非線性最優(yōu)控制問題的hp自適應網(wǎng)格細化算法，通過增加插值多項式維數(shù)或細化網(wǎng)格的方式提高解的精度，優(yōu)化結束時的微分－代數(shù)約束最大誤差為0.0802×10－4，滿足設定的精度要求，優(yōu)化時間為4.197 4s。仿真結果表明：

①與全局偽譜法和局部偽譜法相比，hp偽譜法充分結合了p法的快速收斂性和h法的計算稀疏性，能夠以較少的計算代價得到較高精度的解；

②與已有hp算法相比，本文方法更適合求解多約束條件下軌跡優(yōu)化問題，具有實時應用的潛力，可作為實時制導的基礎。

[1]ROSS I M，F(xiàn)AHROO F.Pseudospectral knotting methods for solving optimal control problems[J].Journal of Guidance，Control，and Dynamics，2004，27（3）：397－405.

[2]JAIN D，TSIOTRAS P.Trajectory optimization using multiresolution techniques[J].Journal of Guidance，Control，and Dynamics，2008，31（5）：1 424－1 436.

[3]ZHAO Y P，TSIOTRAS A.Density－function based mesh refinement algorithm for solving optimal control problems[C]//Infotech and Aerospace Conference.Seattle，Washington：AIAA，2009.

[4]GUI W，BABUSKA I.The h，p，and hp versions of the finite element method in 1dimension.Part I：the error analysis of the p version[J].Numerische Mathematik，1986，49：577－612.

[5]GUI W，BABUSKA I.The h，p，and hp versions of the finite element method in 1dimension.Part II：the error analysis of the h and h－p versions[J].Numerische Mathematik，1986，49：613－657.

[6]GUI W，BABUSKA I.The h，p，and hp versions of the finite element method in 1dimension.Part III：the adaptive h－p version[J].Numerische Mathematik，1986，49：659－683.

[7]GALVAO A，GERRITSMA M，MAERSCHALK B D.Hp－adaptive least－squares spectral element method for hyperbolic partial differential equations[J].Journal of Computational and Applied Mathematics，2008，215（2）：409－418.

[8]DARBY C L，HAGER W W，ANIL V R.An hp－adaptive pseudospectral method for solving optimal control problems[J].Optimal control applications and methods，2011，32（4）：476－502.

[9]DARBY C L，HAGER W W，ANIL V R.A preliminary analysis of a variable－order approach to solving optiaml control problems using pseudospectral methods[C]//AIAA/AAS Astrodynamics Specialist Conference.Toronto， Canada：AIAA，2010.

[10]DARBY C L，HAGER W W，ANIL V R.Direct trajectory optimization using a variable low－order adaptive pseudospectral method[J].Journal of Spacecraft and Rockets，2011，48（3）：433－445.

[11]BETTS J T.Practical methods for optimal control and estimation using nonlinear programming[M].Washington：Society for Industrial and Applied Mathematics，2010.