軌跡優化的LASSO網格自適應加密方法

張　　松，侯　明　善

（1.西北工業大學自動化學院，陜西　西安 710072；2.中航工業第一飛機設計研究院，陜西　西安 710089）

軌跡優化的LASSO網格自適應加密方法

張松1，2，侯明善1

（1.西北工業大學自動化學院，陜西西安 710072；2.中航工業第一飛機設計研究院，陜西西安 710089）

針對軌跡優化直接方法，提出了以控制變量曲率為基礎的最小絕對收縮與選擇算子（least absolute shrinkage and selection operator，LASSO）網格自適應加密策略，用于提高優化精度。以高分辨率二分網格節點為中心，構造徑向基函數逼近控制曲線，利用LASSO方法估計徑向基函數系數，并自動篩選出位于控制曲線曲率極大區間的高分辨率節點加密當前網格。本文方法不需要進行狀態和控制誤差估計，適應性和通用性強。兩組典型算例驗證了方法的有效性。

軌跡優化；網格加密；最小絕對收縮與選擇；徑向基函數

網址：www.sys-ele.com

0　引　言

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

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

1　優化問題描述

1.1連續優化問題

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

式（1）～式（4）為連續Bolza問題的數學描述。

1.2Runge-Kutta離散

不失一般性，假設狀態和控制變量在一組時間點

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

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

目標函數（）可離散化為

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

并滿足約束條件

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

2　網格加密算法

在給出網格加密算法之前，本節將先引入二分網格和徑向基函數的概念。

2.1二分網格與徑向基函數

二分網格是指［0，1］區間上的一組均勻劃分，具體形式如下：

式中，j表示網格的分辨率層級；k為空間位置。當時間區間端點t0和tf已知時，與二分網格Gj對應的時間區間網格可表示為

可見，分辨率層級j越高，網格Tj就越細密。

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

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

式中，qj為逼近系數。逼近系數可以根據已知數據進行估算。

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

2.2加密節點選取方案

首先，設當前離散網格為Tc，構造高分辨率節點集合：

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

其次，構造控制曲線逼近函數。根據式（18）構造徑向基函數集

式中，qk為與式（18）中節點對應的待定系數。

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

步驟1產生估計樣本。給定采樣時刻集合

根據當前時間網格Tc及其對應的最優控制解，利用插值方法采集n個最優控制樣本

式中，n為一較大整數。

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

需要指出的是，式（24）中參數λ的值對節點選取結果具有直接影響，下面給出具體算例進行說明。考慮函數

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

圖1　參數λ對加密節點的影響

2.3多分辨軌跡優化算法

本文多分辨率軌跡優化算法如下：

初始化首先，設定滿足最低精度要求的基準分辨率層級Jmin，并令J=Jmin，生成初始離散網格

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

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

步驟1令Told=Titer，在網格Titer上利用1.2節方法解軌跡優化問題。記優化控制變量為Uopt，優化末端時間為tfopt。產生采樣時間序列Tsopt

步驟2加密節點選擇：

（1）計算高分辨率節點集合

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

步驟3網格加密：將選出的節點加入當前網格得Tnew= Told∪Tadd，令J=J+ΔJ。

步驟4如果J＜Jmax，則令Titer=Tnew，插值產生與Tnew對應的優化變量初值，轉入步驟1；否則轉入步驟5。

步驟5終止算法，輸出優化結果。

圖2給出了多分辨率軌跡優化流程圖。

（2）構造與節點集合式（27）或式（28）對應的徑向基函數集

（3）解回歸分析問題

圖2　多分辨率軌跡優化流程圖

3　算法驗證

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

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

確定控制量u（t）使得目標函數

最小，并滿足微分方程

邊界條件為

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

圖3　狀態變量x

圖4　狀態變量υ

圖5　控制變量u

圖6　網格加密過程

3.2飛機最短時間爬升問題

鉛垂面內飛機的動力學方程［17］為

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

式中，L為升力；質量m和重力常數g分別為

阻力D和推力T滿足如下關系：

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

優化目標為飛機爬升時間最短，即最小化

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

狀態及控制約束分別為

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

圖7　高度時間歷

圖8　速度時間歷程

圖9　飛行路徑角時間歷程

圖10　載荷系數時間歷程

圖11　網格加密過程

4　結　論

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

［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.（雍恩米，陳磊，唐國金.飛行器軌跡優化數值方法綜述［J］.宇航學報，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.（劉淵博，朱恒偉，黃小念，等.偽譜法求解非光滑最有控制問題的網格優化［J］.系統工程與電子技術，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.（趙吉松，谷良賢，佘文學.配點法和網格加密技術用于非光滑軌跡優化［J］.宇航學報，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.（豐志偉，張青斌，唐乾剛，等.基于二代小波的軌跡優化節點自適應加密［J］.航空動力學報，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；網絡優先出版日期：2016-01-12。

網絡優先出版地址：http://www.cnki.net/kcms/detail/11.2422.TN.20160112.1737.006.html

張松（1985-）男，博士研究生，主要研究方向為飛行器軌跡優化與制導技術。

E-mail：zhangsong.gz@outlook.com

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

E-mail：mingshan@nwpu.edu.cn