李 臻，許冰青，李慶波，鄢雄偉，李博雅

（1.上海機(jī)電工程研究所，上海 201109；2.陸軍駐上海地區(qū)軍事代表室，上海201109）

0 引 言

助推滑翔式飛行器在飛行時，長時間處于稠密大氣層內(nèi)，熱力學(xué)環(huán)境復(fù)雜惡劣，為了保證彈載設(shè)備及機(jī)身結(jié)構(gòu)正常，飛行器必須滿足熱流、動壓、過載等多方面的約束［1-2］。在工程優(yōu)化領(lǐng)域，針對這類非線性規(guī)劃問題，一般基于直接法求解，如序列二次規(guī)劃等方法。例如，黃魯豫等［3］基于序列二次規(guī)劃方法研究了飛行器復(fù)合控制技術(shù)；李征等［4］基于序列二次規(guī)劃解決了偽譜法轉(zhuǎn)換后的無人機(jī)集群規(guī)劃問題；Cui等［5］利用高斯偽譜法求解了基于線性協(xié)方差的彈道優(yōu)化方法，解決火星再入問題。然而，這類方法在使用中常存在初值敏感、收斂較慢的問題。

近年來，凸優(yōu)化方法以其獨(dú)特的理論優(yōu)勢逐漸受到青睞。對于一個凸問題而言，基本性質(zhì)明確了其得到的局部最優(yōu)解即為全局最優(yōu)解，算法解集的最優(yōu)性得到了保證［6］。然而，凸優(yōu)化方法在一般非線性問題中的應(yīng)用存在一定限制，要求原問題保持凸性，而一般的非線性問題具有高度的非凸性［7］。為此，需要將原非線性規(guī)劃問題轉(zhuǎn)化為一系列近似的線性規(guī)劃子問題，通過迭代逐步逼近原問題的解，這一過程稱為序列凸優(yōu)化［8］。本文以助推滑翔式飛行器為應(yīng)用背景，設(shè)計(jì)了在序列凸優(yōu)化算法下的非線性規(guī)劃方法，具有初值不敏感性、快速收斂性。

1 規(guī)劃問題與模型

在半速度系下，建立無量綱的飛行器動力學(xué)模型，可表達(dá)為［9］

式中：r為單位地心距；θ、?分別為經(jīng)度和緯度；v為單位速度；γ、ψ分別為航跡角和航向角；α、σ分別為攻角與傾側(cè)角；ω為地球自轉(zhuǎn)速度；D、L為歸一化后阻力和升力。

式（1）是關(guān)于6 個狀態(tài)變量與2 個控制變量的強(qiáng)非線性方程，本文通過文獻(xiàn)［10］所述方法，設(shè)計(jì)新的控制變量來對運(yùn)動方程進(jìn)行解耦處理，將原控制量攻角、傾側(cè)角引入狀態(tài)量中，擴(kuò)維狀態(tài)為

取攻角導(dǎo)數(shù)、傾側(cè)角導(dǎo)數(shù)為新的控制變量u，新的擴(kuò)維動力學(xué)方程可表達(dá)為

式中：動力學(xué)列向量fp(x)∈R8；控制量系數(shù)陣B∈R8×2；角速度相關(guān)列向量fω(x)∈R8。解耦后動力學(xué)模型中控制系數(shù)陣B為一常量，獨(dú)立于狀態(tài)參數(shù)。

飛行器規(guī)劃問題中還存在一系列約束，如：①禁飛區(qū)約束，用Nf表示；②熱流Q?、動壓q、過載n過程約束；③控制系統(tǒng)性能約束；④初末值約束等常見過程約束形式?？梢杂洖?/p>

式中：rn、θn、?n為禁飛區(qū)的半徑、中心點(diǎn)經(jīng)緯度；ρ為大氣密度；Q?max、qmax、nmax分別為熱流、動壓、過載的約束極值；umax為控制的約束極值；xf為末狀態(tài)。

基于上述分析，軌跡優(yōu)化問題可以表述為：確定容許控制u（t）∈R2以及容許狀態(tài)x（t）∈R8，使得由一個微分方程組（式（3））確定的系統(tǒng)，從給定的初始狀態(tài)過渡到給定的終端狀態(tài)，在滿足規(guī)定約束（式（4））的同時，使性能指標(biāo)函數(shù)J達(dá)到最小。軌跡優(yōu)化問題P1的數(shù)學(xué)形式可表達(dá)為

Problem 1：

2 參數(shù)化方法

軌跡優(yōu)化問題P1具有連續(xù)非線性非凸的形式，本文采用分段高斯偽譜法，將上述非凸問題離散。將原問題時域均分為N個區(qū)間，每段區(qū)間上選擇N1個高斯點(diǎn)。原動力學(xué)方程在每段區(qū)間上可以離散為式（5）所示形式，其推導(dǎo)參見文獻(xiàn)［11］。

式中：τ為［-1，1］區(qū)間上的時間變量；ti為第i段區(qū)間的時間起點(diǎn)；由Djk構(gòu)成的系數(shù)矩陣稱為微分近似陣，j＝1，…，N1，k＝0，…，N1，可用D來表示 ，D∈RN1×(N1＋1)。Djk的形式為

式（5）對應(yīng)的是一個區(qū)間段上的微分等式約束，而該問題在全區(qū)間段[t0，tf]上共分為了N段，因此在整個規(guī)劃問題中需要同時考慮這N個微分等式約束。為簡化表達(dá)，以式（7）代表這一系列的等式約束。

高斯偽譜法下系統(tǒng)微分方程只在高斯點(diǎn)處進(jìn)行了離散近似，對于區(qū)間端點(diǎn)還需利用高斯求積公式約束，用x(τN1＋1)、x(τ0)表示[ti，ti＋1]區(qū)間端點(diǎn)狀態(tài)

式中：wj為每一節(jié)點(diǎn)處的高斯權(quán)重，可表達(dá)為

同樣，在全區(qū)間上需考慮N次端點(diǎn)約束，采用相同方式，式（8）簡記為

新的非線性規(guī)劃問題P2形式為

Problem 2：

3 約束轉(zhuǎn)化

對于擴(kuò)維動力學(xué)方程（式（3）），以x*表示上一步子問題解，一階線性展開后新的動力學(xué)方程可表達(dá)為

式中：A(x*)為動力學(xué)項(xiàng)fp(x)的偏導(dǎo)系數(shù)陣。

動力學(xué)方程中未展開地球自轉(zhuǎn)項(xiàng)，原因在于該項(xiàng)在動力學(xué)方程中量級較小，近似處理不會產(chǎn)生較大誤差。

采用相似的方式處理式（4）中的禁飛區(qū)約束、熱流約束、動壓約束、過載約束，形式為

線性化只有在基準(zhǔn)值附近才能保證近似精度，因此需要在原問題中附加一個置信域約束

式中：δx表示各歸一化后狀態(tài)量的極限偏差取值范圍。

該問題中還存在以下約束：①對無動力滑翔飛行器而言，飛行過程伴隨著能量耗散，因而速度曲線呈現(xiàn)出遞減的趨勢，滿足Δv≤0；②飛行器作平衡滑翔飛行，滿足δγ≤γ≤0；③飛行器高度會一直下降，滿足Δr≤0。

至此，已完成了對問題P2中所有非線性約束的轉(zhuǎn)化，非線性規(guī)劃問題已經(jīng)轉(zhuǎn)化為有限維的線性規(guī)劃問題P3。

4 帶罰函數(shù)的序列凸優(yōu)化算法

在一般性初值下，直接迭代求解問題P3是不可行的［12］。為了使算法具有初值不敏感性，并能以較快速度收斂到可行解。本文給出一種帶罰函數(shù)的序列凸優(yōu)化算法。設(shè)計(jì)如下：

1）第一步，在原問題P3 中對所有轉(zhuǎn)化后的非線性約束項(xiàng)引入松弛變量ξ與ζ，并放棄置信域約束，該部分目標(biāo)函數(shù)設(shè)置為式（14）所述形式，式中c1、c2、p為各項(xiàng)的常系數(shù)。

以第（k-1）次迭代值為基準(zhǔn)，第k次迭代求解子問題SP1。子問題中用(·)k-1的形式代替前文(·)*，作為展開基值。(·)k元素為第k次迭代時待規(guī)劃量。

式中：ξk1、ξk2為微分等式約束對應(yīng)的松弛變量；ζ ki為不等式約束下的松弛變量。所有松弛變量趨于0 時，第一步規(guī)劃子問題SP1 與問題P3 在約束上是相同的。該步目的在于，在一個較差的初值下能尋找到合適可行域以初步滿足約束條件。

由于子問題SP1 中不等式約束比微分約束更容易滿足，令e0＝不等式約束判斷條件為：若e0＜10-5，認(rèn)為當(dāng)前解在不等式約束上已經(jīng)足夠完備，轉(zhuǎn)入下一步規(guī)劃。

2）第二步，舍棄不等式約束相應(yīng)的松弛因子，求解新的過渡子問題SP2，該部分目標(biāo)函數(shù)設(shè)置為

該步規(guī)劃采用嚴(yán)格不等式約束，保留松弛微分約束，加速尋找初值過程。

從第二步規(guī)劃過渡到第三步規(guī)劃時，需考慮微分誤差的求解情況。令，微分約束判別條件為：當(dāng)eˉ≤10-8時，認(rèn)為當(dāng)前解對應(yīng)的子問題SP2與問題P3在約束意義上足夠接近，過渡到下一步。

3）第三步，回歸問題P3，子問題SP3 目標(biāo)函數(shù)重設(shè)為最小化置信域誤差

式中：q為置信域系數(shù)；Δxkp由k次迭代時原控制變量Δαk與Δσk構(gòu)成。

求解子問題SP3 時，若某一次計(jì)算結(jié)果使|Δxk|≤ε成立，則認(rèn)為當(dāng)前規(guī)劃結(jié)果是滿足原非線性規(guī)劃問題P1的一組可行解。

帶罰函數(shù)的序列凸優(yōu)化算法流程如圖1所示。

圖1 帶罰函數(shù)的序列凸優(yōu)化算法流程Fig.1 The flow of sequential convex optimization with penalty function

5 仿真實(shí)驗(yàn)

對于求解單個凸優(yōu)化問題的過程，本文采用的是SDTP3工具包。在應(yīng)用時，要求輸入的所有等式約束是仿射的，所有的不等式約束是凸的，即只解決數(shù)學(xué)意義上滿足凸性的問題。因此，在配合本文所設(shè)計(jì)算法的情況下，可以應(yīng)用于非線性規(guī)劃問題求解。

該部分設(shè)計(jì)了一個仿真算例來驗(yàn)證上述序列凸優(yōu)化算法的可行性。某型飛行器啟滑段狀態(tài)參數(shù)固定，根據(jù)飛行任務(wù)，要求在1 500s 后抵達(dá)西南方向1 500 km外的固定目標(biāo)點(diǎn)，飛行中途存在一個禁飛區(qū)?？刂萍s束與路徑約束見表1。

表1 控制約束與路徑約束Tab.1 Control constraint and path constraint

算法中置信域約束與收斂條件取值為

松弛因子與置信域誤差的常系數(shù)分別取c1＝100、c2＝200、p＝10、q＝0.2。

為了更好展現(xiàn)本節(jié)算法的性能，將其與非線性規(guī)劃領(lǐng)域已廣泛應(yīng)用的求解工具包GPOPS-2 進(jìn)行對比。在數(shù)值仿真階段，GPOPS 軟件與序列凸優(yōu)化方法輸入的原問題相同，為P1的形式。并且求解動力學(xué)模型均采用式（3）所述的擴(kuò)維模型。偽譜法參數(shù)如表2所示。

表2 偽譜法參數(shù)Tab.2 The parameters in the pseudo-spectral method

仿真結(jié)果如圖2～11 所示。為體現(xiàn)算法的初值不敏感性，迭代的初值彈道直接按常值攻角、常值傾側(cè)角的控制規(guī)律提供，不滿足絕大部分約束條件。

圖2 迭代中水平面軌跡變化Fig.2 The trajectories during iteration

在該初值下，序列凸優(yōu)化算法經(jīng)過7 次迭代運(yùn)算解收斂，耗時142.2 s。迭代中每次運(yùn)算的軌跡變化如圖2所示，其中紅色加粗實(shí)線段代表最后收斂的軌跡，初值軌跡則是圖中最下方穿過禁飛區(qū)的那一條曲線。可以明顯看出，求解過程中，解對應(yīng)的軌跡不斷向禁飛區(qū)左側(cè)挪動，最后收斂。運(yùn)算中目標(biāo)函數(shù)變化如圖3 所示，隨迭代進(jìn)行趨于0，因此最后得到的解是可信的。

圖3 迭代中目標(biāo)函數(shù)值變化Fig.3 Objective function values during iteration

圖4 兩種規(guī)劃水平面軌跡Fig.4 The horizontal trajectory using two methods

圖5 兩種規(guī)劃空間軌跡Fig.5 The space trajectory using two methods

圖6 兩種規(guī)劃攻角規(guī)律Fig.6 The profile of attack angle using two methods

圖7 兩種規(guī)劃傾側(cè)角規(guī)律Fig.7 The profile of bank angle using two methods

兩種優(yōu)化方法下的可行解結(jié)果對比見圖4～11，其中紅色虛線代表凸優(yōu)化仿真結(jié)果，藍(lán)色點(diǎn)劃線代表凸優(yōu)化初值，黃色實(shí)線代表GPOPS 方法仿真結(jié)果。

圖8 兩種規(guī)劃速度變化Fig.8 The profile of velocity using two methods

圖9 兩種規(guī)劃航跡角變化Fig.9 The profile of flight path angel using two methods

圖10 兩種規(guī)劃熱流密度變化Fig.10 The density of heat flow using two methods

由于未額外約束可行解過程參數(shù)，因此從結(jié)果上看二者有差異，但均滿足了所有約束要求。凸優(yōu)化方法找到的可行解是從禁飛區(qū)左側(cè)繞飛后抵達(dá)目標(biāo)，而GPOPS 方法的解是從禁飛區(qū)右側(cè)繞飛后抵達(dá)目標(biāo)，因此控制量規(guī)律也會產(chǎn)生相應(yīng)變化。

總體而言，序列凸優(yōu)化方法與GPOPS 軟件均找到了合適可行解，結(jié)果對比如表3 所示。在相同仿真環(huán)境下，GPOPS方法經(jīng)過了6次迭代，耗時465.5 s，慢于本文設(shè)計(jì)的序列凸優(yōu)化算法的速度；求解精度為10-3，也不及本文設(shè)計(jì)的算法。

圖11 兩種規(guī)劃動壓變化Fig.11 The dynamic pressure using two methods

表3 方法結(jié)果對比Tab.3 Comparison of Results

6 結(jié)束語

本文以助推滑翔式飛行器在多約束條件下軌跡規(guī)劃為應(yīng)用背景，設(shè)計(jì)了針對一般非線性規(guī)劃問題尋找可行解的序列凸優(yōu)化算法。通過算例仿真，驗(yàn)證了設(shè)計(jì)算法可行性，能克服初值敏感并以較快速度收斂，滿足多約束要求。