張　　松，侯　明　善

（1.西北工業(yè)大學(xué)自動化學(xué)院，陜西　西安 710072；2.中航工業(yè)第一飛機(jī)設(shè)計研究院，陜西　西安 710089）

軌跡優(yōu)化的LASSO網(wǎng)格自適應(yīng)加密方法

張松1，2，侯明善1

（1.西北工業(yè)大學(xué)自動化學(xué)院，陜西西安 710072；2.中航工業(yè)第一飛機(jī)設(shè)計研究院，陜西西安 710089）

針對軌跡優(yōu)化直接方法，提出了以控制變量曲率為基礎(chǔ)的最小絕對收縮與選擇算子（least absolute shrinkage and selection operator，LASSO）網(wǎng)格自適應(yīng)加密策略，用于提高優(yōu)化精度。以高分辨率二分網(wǎng)格節(jié)點為中心，構(gòu)造徑向基函數(shù)逼近控制曲線，利用LASSO方法估計徑向基函數(shù)系數(shù)，并自動篩選出位于控制曲線曲率極大區(qū)間的高分辨率節(jié)點加密當(dāng)前網(wǎng)格。本文方法不需要進(jìn)行狀態(tài)和控制誤差估計，適應(yīng)性和通用性強(qiáng)。兩組典型算例驗證了方法的有效性。

軌跡優(yōu)化；網(wǎng)格加密；最小絕對收縮與選擇；徑向基函數(shù)

網(wǎng)址：www.sys-ele.com

0　引　言

軌跡優(yōu)化問題廣泛存在于航空航天領(lǐng)域中，如導(dǎo)彈中制導(dǎo)、飛機(jī)航路設(shè)計和航天器軌道轉(zhuǎn)移等，其數(shù)值解法主要分為間接法和直接法［1 2］。間接法根據(jù)最優(yōu)控制原理將優(yōu)化問題轉(zhuǎn)換為Hamilton邊值問題進(jìn)行求解，直接法則在一定網(wǎng)格上離散化狀態(tài)、控制，將連續(xù)最優(yōu)控制問題轉(zhuǎn)化為有限維非線性規(guī)劃問題迭代求解。直接法的離散網(wǎng)格可根據(jù)節(jié)點的分布情況分為均勻網(wǎng)格和非均勻網(wǎng)格兩類。使用均勻網(wǎng)格進(jìn)行離散狀態(tài)或控制會引入大量冗余節(jié)點，優(yōu)化計算效率較低且穩(wěn)定性較差。因此，通常需要根據(jù)問題特性生成非均勻網(wǎng)格。文獻(xiàn)［3］通過解一組關(guān)于離散誤差的整數(shù)規(guī)劃對舊網(wǎng)格進(jìn)行局部加密。文獻(xiàn)［4-7］針對偽譜法給出了不同的分段網(wǎng)格加密方法，基本思想是根據(jù)狀態(tài)偏差確定可能的不連續(xù)點并以此為依據(jù)調(diào)整或增加分段網(wǎng)格和離散節(jié)點數(shù)目。文獻(xiàn)［8-9］基于二分網(wǎng)格提出了多分辨率網(wǎng)格加密技術(shù)，能夠以較少的離散節(jié)點獲得較高的精度，文獻(xiàn)［10］從細(xì)化效率、易用性和適應(yīng)性等角度出發(fā)對該方法進(jìn)行了改進(jìn)。以上方法都是基于局部積分或插值誤差來配置節(jié)點，當(dāng)最優(yōu)解具有周期特性時往往存在求解困難。此外，文獻(xiàn)［11］引入二代小波技術(shù)對控制或狀態(tài)進(jìn)行小波變換，并根據(jù)小波系數(shù)確定迭代中所使用的節(jié)點，有效降低了所需節(jié)點數(shù)，但小波轉(zhuǎn)換過程極為復(fù)雜。文獻(xiàn)［12］則提出密度函數(shù)法，根據(jù)控制變量的連續(xù)性來重新配置節(jié)點，但由于采用差分方法估算曲率，算法穩(wěn)定性欠佳。

簡言之，軌跡優(yōu)化問題中的網(wǎng)格加密就是在軌跡變化劇烈部分增加節(jié)點，以準(zhǔn)確捕獲軌跡特性，且變化越劇烈增加節(jié)點應(yīng)越多。對于一般軌跡優(yōu)化問題，軌跡的特性主要由控制變量決定，因而可以直接根據(jù)控制曲線的曲率特性來加密網(wǎng)格。筆者研究發(fā)現(xiàn)，若能構(gòu)造一組合適的基函數(shù)對控制變量進(jìn)行逼近，那么可以利用最小絕對收縮與選擇（least absolute shrinkage and selection operator，LASSO）方法自動選取加密節(jié)點加密網(wǎng)格。LASSO方法是一種變量選擇方法，基本思想是用逼近函數(shù)系數(shù)的絕對值作為罰函數(shù)來壓縮逼近函數(shù)的系數(shù)，使絕對值較小的系數(shù)自動壓縮為零，從而實現(xiàn)顯著性基函數(shù)的選擇和對應(yīng)參數(shù)的估計［13］。這里，選出的顯著性基函數(shù)通常位于曲線曲率較大或不連續(xù)處。本文將LASSO方法與徑向基函數(shù)逼近方法應(yīng)用于軌跡優(yōu)化問題中的網(wǎng)格加密，在每次優(yōu)化迭代完成后生成高分辨率節(jié)點并構(gòu)造徑向基函數(shù)逼近控制曲線，利用LASSO方法估計徑向基函數(shù)系數(shù)，并將系數(shù)非零基函數(shù)對應(yīng)的高分辨率節(jié)點加入當(dāng)前網(wǎng)格，形成新的多分辨率離散網(wǎng)格。與現(xiàn)有方法相比，本文方法不需要對狀態(tài)和控制變量誤差進(jìn)行估計，且所需設(shè)定的參數(shù)少，適應(yīng)性和通用性強(qiáng)。兩組典型算例驗證了本文方法的有效性。

1　優(yōu)化問題描述

1.1連續(xù)優(yōu)化問題

一般軌跡優(yōu)化問題可表述為：確定控制量u（t）∈Rm，使得Bolza型代價函數(shù)

式（1）～式（4）為連續(xù)Bolza問題的數(shù)學(xué)描述。

1.2Runge-Kutta離散

不失一般性，假設(shè)狀態(tài)和控制變量在一組時間點

上進(jìn)行離散。記xi=x（ti），ui=u（ti），則對狀態(tài)方程（4）的q階Runge-Kutta離散格式為

ρj，βj，αjl均為已知常數(shù)并且滿足0≤ρ1≤ρ2≤…≤ρq≤1，當(dāng)αjl=0（l≥j）時，離散格式為顯式，否則為隱式格式。

目標(biāo)函數(shù)（）可離散化為

則由上述連續(xù)Bolza問題離散化得到的NLP可描述為：確定離散狀態(tài)變量X，離散控制變量U和ˉU，初始時間t0，終端時間tf，使得如下目標(biāo)函數(shù)最?。?/p>

并滿足約束條件

式中，階數(shù)q取不同值可以得到不同的離散格式；中間控制變量uij可通過插值求得而不作為優(yōu)化變量。本文采用經(jīng)典4階Runge-Kutta離散格式。

2　網(wǎng)格加密算法

在給出網(wǎng)格加密算法之前，本節(jié)將先引入二分網(wǎng)格和徑向基函數(shù)的概念。

2.1二分網(wǎng)格與徑向基函數(shù)

二分網(wǎng)格是指［0，1］區(qū)間上的一組均勻劃分，具體形式如下：

式中，j表示網(wǎng)格的分辨率層級；k為空間位置。當(dāng)時間區(qū)間端點t0和tf已知時，與二分網(wǎng)格Gj對應(yīng)的時間區(qū)間網(wǎng)格可表示為

可見，分辨率層級j越高，網(wǎng)格Tj就越細(xì)密。

徑向基函數(shù)是指函數(shù)值只依賴自變量與中心間距離的實值函數(shù)，其一般表達(dá)式如下：

式中，z為自變量；c為中心。常見徑向基函數(shù)包括高斯函數(shù)和多面函數(shù)等。徑向基函數(shù)可用于函數(shù)逼近：給定徑向基函數(shù)φ：R+→R，可以將待逼近函數(shù)f（z）表示成徑向基函數(shù)和的形式，即

式中，qj為逼近系數(shù)。逼近系數(shù)可以根據(jù)已知數(shù)據(jù)進(jìn)行估算。

如果f（z）表示控制變量，中心取為比當(dāng)前節(jié)點分辨率更高的節(jié)點，那么可以通過LASSO方法選出用于加密網(wǎng)格的高分辨率節(jié)點。

2.2加密節(jié)點選取方案

首先，設(shè)當(dāng)前離散網(wǎng)格為Tc，構(gòu)造高分辨率節(jié)點集合：

式中，Tj如式（15）所示；jmin和jmax分別表示集合^T 的最小和最大分辨率，jmax＞jmin，且jmin大于當(dāng)前離散網(wǎng)格Tc的分辨率層級。

其次，構(gòu)造控制曲線逼近函數(shù)。根據(jù)式（18）構(gòu)造徑向基函數(shù)集

式中，qk為與式（18）中節(jié)點對應(yīng)的待定系數(shù)。

最后，利用LASSO方法估計逼近函數(shù)式（21）的系數(shù)，并將qk（k=1，2，…，m）中非零項對應(yīng)的高分辨率節(jié)點加入當(dāng)前網(wǎng)格形成下一輪優(yōu)化所使用的網(wǎng)格。LASSO系數(shù)估計的具體步驟如下：

步驟1產(chǎn)生估計樣本。給定采樣時刻集合

根據(jù)當(dāng)前時間網(wǎng)格Tc及其對應(yīng)的最優(yōu)控制解，利用插值方法采集n個最優(yōu)控制樣本

式中，n為一較大整數(shù)。

得到系數(shù)向量Q，式中‖·‖1和‖·‖2分別為l1和l2范數(shù)。式（24）前半部分是常規(guī)的最小二次型估計，而后半部分是關(guān)于系數(shù)向量Q的罰函數(shù)項。通過解算式（24）估計系數(shù)向量Q的方法即為LASSO方法。LASSO方法能將非顯著性基函數(shù)的系數(shù)準(zhǔn)確地縮減到零，即高分辨率網(wǎng)格^T中位于曲線平滑處的節(jié)點對應(yīng)基函數(shù)的系數(shù)將被縮減到零。因此，與系數(shù)向量Q中非零元對應(yīng)的節(jié)點即為位于曲線變化劇烈處的細(xì)化節(jié)點，利用此類節(jié)點細(xì)化網(wǎng)格將極大提高解算精度。問題（24）可以采用小角回歸方法［14］進(jìn)行計算。

需要指出的是，式（24）中參數(shù)λ的值對節(jié)點選取結(jié)果具有直接影響，下面給出具體算例進(jìn)行說明?？紤]函數(shù)

令jmin=4和jmax=7，根據(jù)式（18）生成包含129個相異節(jié)點的高分辨率節(jié)點集合^T，并構(gòu)造形如式（19）的徑向基函數(shù)集。采樣節(jié)點取為區(qū)間［0.1，10］內(nèi)100個等距分布點。采用第2.2節(jié)方法，解回歸分析問題（25）得系數(shù)向量Q。參數(shù)λ取9.905 0×10-3、9.905 0×10-5和9.905 0× 10-7時的節(jié)點選取結(jié)果如圖1所示，圖中實線為函數(shù)g（t），圓圈為篩選出的節(jié)點，節(jié)點數(shù)分別為16、66和124。由圖1可知，當(dāng)λ較大時，選出的加密節(jié)點數(shù)量較少，且均位于曲線曲率變化劇烈處，隨著λ值的減小，加密節(jié)點數(shù)量不斷增加，且在曲線較為平坦處也有節(jié)點分布。因此，需要對參數(shù)λ的值進(jìn)行折中選取，以保證加密節(jié)點數(shù)量在合理水平。本文采用十折交叉校驗方法確定參數(shù)λ的值［15］。

圖1　參數(shù)λ對加密節(jié)點的影響

2.3多分辨軌跡優(yōu)化算法

本文多分辨率軌跡優(yōu)化算法如下：

初始化首先，設(shè)定滿足最低精度要求的基準(zhǔn)分辨率層級Jmin，并令J=Jmin，生成初始離散網(wǎng)格

式中，GJ如式（14）所述。

其次，設(shè)置最大分辨率層級Jmax，每次迭代分辨率升級參數(shù)ΔJ，以及樣本數(shù)p。

步驟1令Told=Titer，在網(wǎng)格Titer上利用1.2節(jié)方法解軌跡優(yōu)化問題。記優(yōu)化控制變量為Uopt，優(yōu)化末端時間為tfopt。產(chǎn)生采樣時間序列Tsopt

步驟2加密節(jié)點選擇：

（1）計算高分辨率節(jié)點集合

（4）選擇高分辨率節(jié)點集合^T中與Qopt非零元對應(yīng)且不屬于網(wǎng)格Told的節(jié)點為加密節(jié)點，記為Tadd。

步驟3網(wǎng)格加密：將選出的節(jié)點加入當(dāng)前網(wǎng)格得Tnew= Told∪Tadd，令J=J+ΔJ。

步驟4如果J＜Jmax，則令Titer=Tnew，插值產(chǎn)生與Tnew對應(yīng)的優(yōu)化變量初值，轉(zhuǎn)入步驟1；否則轉(zhuǎn)入步驟5。

步驟5終止算法，輸出優(yōu)化結(jié)果。

圖2給出了多分辨率軌跡優(yōu)化流程圖。

（2）構(gòu)造與節(jié)點集合式（27）或式（28）對應(yīng)的徑向基函數(shù)集

（3）解回歸分析問題

圖2　多分辨率軌跡優(yōu)化流程圖

3　算法驗證

3.1雙積分最小能量控制問題

雙積分最小能量問題如下［16］：

確定控制量u（t）使得目標(biāo)函數(shù)

最小，并滿足微分方程

邊界條件為

路徑約束為p=50。圖3～圖6給出了最后一次迭代的優(yōu)化結(jié)果，其中，圖3和圖4顯示狀態(tài)變量變化規(guī)律，圖5顯示控制變量變化規(guī)律，圖中實線為解析解，圓圈為本文數(shù)值解。圖6則顯示了多分辨率網(wǎng)格加密過程，其中，第1次迭代為均勻初始網(wǎng)格，第2次迭代起為不斷加密的網(wǎng)格，各次迭代節(jié)點數(shù)分別為17、27、37和51。如圖5和圖6所示，隨著迭代的進(jìn)行，本文方法不斷選出位于控制變量尖角處的高分辨率節(jié)點入到當(dāng)前網(wǎng)格，而在控制變量平滑處的節(jié)點保持不變。4次迭代中目標(biāo)函數(shù)數(shù)值解與解析間的絕對偏差為0.003 8、5.2× 10-5、1.9×10-7和1.3×10-7。若采用均勻網(wǎng)格，節(jié)點數(shù)取27、37和51時目標(biāo)函數(shù)偏差分別為0.001 0、5.36×10-4和7.18×10-5?？梢?，本文的非均勻網(wǎng)格能以較少的節(jié)點獲得較高精度的解。

圖3　狀態(tài)變量x

圖4　狀態(tài)變量υ

圖5　控制變量u

圖6　網(wǎng)格加密過程

3.2飛機(jī)最短時間爬升問題

鉛垂面內(nèi)飛機(jī)的動力學(xué)方程［17］為

式中，距離x、高度h、速度V和飛行路徑角γ為狀態(tài)變量，載荷系數(shù)為控制變量：

式中，L為升力；質(zhì)量m和重力常數(shù)g分別為

阻力D和推力T滿足如下關(guān)系：

式中，a（h）為聲速；θ為開氏溫度；高度h和ˉh的單位分別為m和km。

優(yōu)化目標(biāo)為飛機(jī)爬升時間最短，即最小化

簡便起見，優(yōu)化中不考慮式（33）。邊界條件給定為

狀態(tài)及控制約束分別為

網(wǎng)格加密參數(shù)取為：Jmin=5，Jmax=8，ΔJ=1，樣本數(shù)p=30。圖7～圖11為飛機(jī)最短時間爬升問題最后一次迭代的優(yōu)化結(jié)果。圖7～圖9分別為高度、速度和飛行路徑角時間歷程，圖10為控制變量載荷系數(shù)的時間歷程。圖11顯示了網(wǎng)格加密過程，其中，第一次迭代為均勻初始網(wǎng)格，第二次迭代起為不斷加密的網(wǎng)格，各次迭代節(jié)點數(shù)分別為33，43，62和103。經(jīng)優(yōu)化，飛機(jī)最短爬升時間為172.9s。由圖10和圖11可知，本文方法能夠準(zhǔn)確地捕獲控制變量變化劇烈區(qū)域并在迭代中不斷加密相應(yīng)區(qū)域的網(wǎng)格，以提高解算精度。

圖7　高度時間歷

圖8　速度時間歷程

圖9　飛行路徑角時間歷程

圖10　載荷系數(shù)時間歷程

圖11　網(wǎng)格加密過程

4　結(jié)　論

本文將LASSO方法與徑向基函數(shù)逼近方法用于軌跡優(yōu)化中的網(wǎng)格自適應(yīng)加密，將加密節(jié)點的選擇問題轉(zhuǎn)化為對控制變量逼近函數(shù)系數(shù)的估計問題，并給出了完整的多分辨率軌跡優(yōu)化流程。由于不需要進(jìn)行誤差估計，且所需設(shè)定的參數(shù)較少，本文方法適應(yīng)性和通用性強(qiáng)。多組典型算例仿真表明，本文方法能夠準(zhǔn)確地捕獲軌跡優(yōu)化問題中的非光滑特性，有效提高優(yōu)化精度和計算效率。

［1］Betts J T.Survey of numerical methods for trajectory optimization［J］.Journal of Guidance，Control and Dynamic，1998，21（2）：193-206.

［2］Yong E M，Chen L，Tang G J.A survey of numerical methods for trajectory optimization of spacecraft［J］.Journal of Astronuatics，2008，29（2）：7-16.（雍恩米，陳磊，唐國金.飛行器軌跡優(yōu)化數(shù)值方法綜述［J］.宇航學(xué)報，2008，29（2）：7-16.）

［3］Betts J T，Huffman W.Mesh refinement in direct transcription methods for optimal control［J］.Optimal Control Application and Methods，1998，19（1）：1-21.

［4］Darby C L，Rao A V.A state approximation-based mesh refinement algorithm for solving optimal control problems［C］//Proc. of the AIAA Guidance，Naυigation and Control Conference，2009.

［5］Darby C L，Hager WW，Rao A V.An hp-adaptive pseudospectral method for solving optimal control problems［J］.Optimal Control Application and Methods，2011，32（4）：476-502.

［6］Darby C L，Hager WW，Rao A V.Direct trajectory optimization using a variable low-order adaptive pseudospectral method［J］.Journal of Spacecraft and Rocket，2011，48（3）：433-445.

［7］Liu Y B，Zhu H W，Huang X N，et al.Optimal mesh segmentation algorithm for pseudospectral methods for non-smooth optimal control problems［J］.Systems Engineering and Electronics，2013，35（11）：2396-2399.（劉淵博，朱恒偉，黃小念，等.偽譜法求解非光滑最有控制問題的網(wǎng)格優(yōu)化［J］.系統(tǒng)工程與電子技術(shù)，2013，35（11）：2396-2399.）

［8］Jain S，Tsiotras P.Trajectory optimization using multiresolution techniques［J］.Journal of Guidance，Control，and Dynamics，2008，31（5）：1425-1436.

［9］Jain S，Tsiotras P.Multiresolution-based direct trajectory optimization［C］//Proc.of the 46th IEEE Conference on Decision and Control，2008：5991-5996.

［10］Zhao J S，Gu L X，She W X.Application of lacal collocation method and mesh refinement to non-smooth trajectory optimization［J］.Journal of Astronuatics，2013，34（11）：1442 1450.（趙吉松，谷良賢，佘文學(xué).配點法和網(wǎng)格加密技術(shù)用于非光滑軌跡優(yōu)化［J］.宇航學(xué)報，2013，34（11）：1442-1450.）

［11］Feng Z W，Zhang Q B，Tang Q G，et al.Node adaptive refinement for trajectory optimization based on second-generation wavelets［J］.Journal of Aerospace Power，2013，28（7）：1659 1665.（豐志偉，張青斌，唐乾剛，等.基于二代小波的軌跡優(yōu)化節(jié)點自適應(yīng)加密［J］.航空動力學(xué)報，2013，28（7）：1659-1665.）

［12］Zhao Y M，Tsiotras P.Density function for mesh refinement in numerical optimization［J］.Journal of Guidance，Control，and Dynamics，2011，34（1）：271-277.

［13］Tibshirani R.Regression shrinkage and selection via the LASSO［J］. Journal of the Royal Statistical Society，Series B（Methodological），1996，58（1）：267-288.

［14］Efron B，Hastie T，Johnstone I，et al.Least angle regression［J］. The Annals of Statistical，2004，32（2）：407-499.

［15］Prabir B.A comparative study of ordinary cross-validation，r-fold cross-validation and the repeated learning-testing methods［J］.Biometrika，1989，76（3）：503-514.

［16］Bryson A E，Ho Y C.Applied optimal control［M］.New York：Hemisphere Publishing，1975.

［17］Seywald H，Cliff EM，Well K H.Range optimal trajectories for an aircraft flying in the vertical plane［J］.Journal of Guidance，Control，and Dynamics，1994，17（2）：389-398.

LASSO-based node adaptive refinement in trajectory optimization

ZHANG Song1，2，HOUMing-shan1
（1.School of Automation，Northwestern Polytechnical Uniυersity，Xi'an 710072，China；
2.AVIC the First Aircraft Institute，Xi'an 710089，China）

A least absolute shrinkage and selection operator（LASSO）based node adaptive refinement approach for the direct method in trajectory optimization is proposed.Firstly，a higher resolution grid and its associated radial basis function set are created.Sequentially，the control variables are approximated using the resulting radial basis function set，and its sampling sequence is generated by interpolation.Finally，the coefficients of the formulated approximation function are estimated based on the statistical variable selection method-LASSO. The higher multi-resolution nodes associated with radial basis functions with non-zero coefficient are selected as new nodes.The proposed method refines the mesh without estimation of states and/or errors controls，and few extra parameters are involved.Therefore，the formulated trajectory optimization algorithm behaves strong adaptability and generality.The validity of this method is demonstrated by several typical examples.

trajectory optimization；mesh refinement；least absolute shrinkage and selection operator（LASSO）；radial basis function

V 19

A

10.3969/j.issn.1001-506X.2016.05.34

1001-506X（2016）05-1195-06

2014-11-04；

2015-10-10；網(wǎng)絡(luò)優(yōu)先出版日期：2016-01-12。

網(wǎng)絡(luò)優(yōu)先出版地址：http://www.cnki.net/kcms/detail/11.2422.TN.20160112.1737.006.html

張松（1985-）男，博士研究生，主要研究方向為飛行器軌跡優(yōu)化與制導(dǎo)技術(shù)。

E-mail：zhangsong.gz@outlook.com

侯明善（1959-）男，教授，博士，主要研究方向為飛行器導(dǎo)航、制導(dǎo)與控制以及控制理論與應(yīng)用。

E-mail：mingshan@nwpu.edu.cn