轨迹优化的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，使得如下目标函数最小：

并满足约束条件

式中，阶数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。可见，本文的非均匀网格能以较少的节点获得较高精度的解。

图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