1、求解轨迹优化问题的局部配点法的稀疏性研究 赵吉松 南京航空航天大学航天学院 摘 要： 直接配点法通过对控制变量和状态变量都进行离散将轨迹优化问题转化为非线性规划 (NLP) 问题。为了提高 NLP 的求解效率, 需要利用其偏导数的稀疏特性并建立偏导数的高效计算方法。本文研究了局部配点法离散得到的 NLP 的一阶偏导数的稀疏特性, 建立了一阶偏导数的高效计算方法。推导了 NLP 的目标函数梯度和约束雅克比矩阵的数学表达式, 得到了 NLP 偏导数的稀疏型, 并且将NLP 的偏导数分解为原始轨迹优化问题的偏导数。由于原始轨迹优化问题的约束和变量的数量远少于 NLP 的约束和变量的数量, 从而显著减
2、小了 NLP 的一阶偏导数的计算量。含有离散气动力和推力数据的仿真算例验证了本文方法的有效性。仿真结果表明, 与有限差分法直接计算 NLP 的偏导数相比, 本文方法能够将优化耗时减小至 4%以内, 随着离散节点数目的增加, 计算效率的提升更为显著。关键词： 轨迹优化; 局部配点法; 非线性规划; 一阶偏导数; 稀疏特性; 作者简介：赵吉松 (1984-) , 男, 博士, 南京航空航天大学航天学院讲师, 主要从事飞行器总体设计与轨迹优化等方面的研究。通信地址:南京市秦淮区御道街 29 号 (210016) 电话:18260412336E-mail:收稿日期：2017-07-06基金：国家自然科
3、学基金 (11602107) Exploiting Sparsity in Local Collocation Methods for Solving Trajectory Optimization ProblemsZHAO Ji-song College of Astronautics, Nanjing University of Aeronautics and Astronautics; Abstract： In a direct local collocation method, a trajectory optimization problem is transcribed into
4、a nonlinear programming (NLP) problem. Solving this NLP as efficiently as possible requires that the sparsity of the NLP derivatives should be exploited and the derivatives should be efficiently calculated. In this paper, a computational efficient method is developed for computing the first derivati
5、ves of the NLP functions arising from a local discretization of a trajectory optimization problem. Specifically, the expressions are derived for the NLP objective function gradient and constraint Jacobian. It is shown that the NLP derivatives can be reduced to the first derivatives of the functions
6、in the trajectory optimization problem. As a result, the method derived in this paper reduces significantly the amount of computation required to obtain the first-derivatives required by a NLP solver. The approach derived in this paper is demonstrated by an example with discrete aerodynamic data and
7、 thrust data where it is forund that the time required to solve the NLP is reduced to less than 4% compared with the direct differentiation of the NLP functions using a finite difference method, and the efficiency improvement is more significant as the number of the grid points increases.Keyword： Tr
8、ajectory optimization; Local collocation method; NLP; First derivatives; Sparsity; Received： 2017-07-060 引言轨迹优化对于飞行器设计有着十分重要的意义和工程实际价值1。轨迹优化本质上属于最优控制问题, 其求解方法主要分为间接法和直接法。其中, 直接法中的配点法通过对控制变量和状态变量进行离散将轨迹优化问题转化为非线性规划 (NLP) 问题, 降低了对初值的敏感性并且具有很好的收敛性, 近年来得到广泛的研究和应用2-7。虽然配点法具有多方面优势, 但是其具体实施方法对于提高优化效率具有重要影响
9、。目前, 以序列二次规划 (SQP) 为代表的基于梯度算法的 NLP 求解器需要提供 NLP 的目标函数和约束的一阶偏导数, 甚至二阶偏导数。其中, 一阶偏导数 NLP 求解器需要提供一阶偏导数, 然后采用拟牛顿法 (DFP 法或者 BFGS 法等) 构造近似的二阶偏导数, 比如 SNOPT8;二阶偏导数 NLP 求解器除了需要一阶偏导数, 还需要准确的二阶偏导数, 比如 IPOPT9。但是, 偏导数的计算量通常比较大, 甚至超过优化算法本身。因此, 提高 NLP 的一阶/二阶偏导数的计算效率对于提高轨迹优化效率有重要意义。研究发现, 由轨迹优化离散得到的 NLP 是非常稀疏的, 即 NLP
10、的一阶偏导数和二阶偏导数含有大量零元素2,10-12。Betts 等10较早地研究了局部配点法的稀疏特性, 得到了梯形格式、Hermite-Simpson 格式和 Runge-Kutta 格式的状态方程离散残差约束的偏导数的稀疏型, 并将其中的非零元素转化为最优控制问题的偏导数, 有效地减小了计算量。但是, Betts 等10的研究存在一些不足之处:一是只给出了状态方程离散残差约束的偏导数稀疏特性和非零元素的计算方法, 没有给出目标函数、路径约束和端点约束的偏导数的稀疏型和计算方法;二是对于最常用的 Hermite-Simpson 格式, 没有完全探索出紧凑形式的状态方程离散残差的偏导数的稀疏
11、型, 而是通过在离散节点中间位置添加离散格式约束和状态变量将其转换为分离形式才得到完全的稀疏型。在 Betts 等的研究工作的基础上, Patterson 等11研究了 Radau 伪谱法的一阶和二阶偏导数的稀疏性, 得到了完整的稀疏型, 并且推导出非零元素的高效计算方法 (将NLP 偏导数分解为最优控制问题的偏导数) 。因此, 如果能够基于上述研究, 完全探索出局部配点法的 NLP 偏导数的稀疏型并建立偏导数的高效计算方法, 对于提高局部配点法的优化效率具有重要意义。此外, 与全局配点法 (又称伪谱法, 离散节点是正交多项式的根) 相比, 局部配点法的离散节点可以根据需要任意布置, 在网格细
12、化方面具有更好的灵活性13-19, 适合求解非光滑轨迹优化问题。因而, 提高局部配点法的优化效率还有能够促进局部配点法在非光滑轨迹优化领域的应用。在实际应用中, 一阶偏导数 NLP 求解器比二阶偏导数 NLP 求解器更为常用, 因为二阶偏导数的计算通常比较繁琐并且计算量较大。以 NLP 的一阶偏导数为例, 其计算方法可分为两大类。一类方法是直接计算 NLP 的偏导数, 包括自动微分法20、复变量微分法21、有限差分法等。其中自动微分法计算量小, 精度高, 得到了广泛应用, 其局限性在于要求优化模型解析可导, 不适用于带有离散数据的轨迹优化问题。对于带有离散数据的优化模型, 只能采用有限差分法。
13、但是, 采用有限差分直接计算 NLP 偏导数的效率较低, 需要耗费大量的机时。另一类方法是将 NLP 的偏导数分解为最优控制问题的偏导数, 然后采用各种微分算法计算最优控制问题的偏导数, 并组装得到 NLP 的偏导数。因为与 NLP 相比, 最优控制的约束和变量的数量大幅减少, 因而这样处理可以显著提高偏导数的计算效率。本文在 Betts 等10和 Patterson 等11的研究工作的基础上, 以 Hermite-Simpson 格式为例, 研究了局部配点法离散得到的 NLP 的一阶偏导数 (目标函数梯度和约束雅克比矩阵) 的稀疏性, 建立非零元素的高效计算方法。采用带有离散参数模型的优化算
14、例验证了所述方法的有效性。仿真结果表明, 与采用有限差分法直接计算 NLP 的偏导数相比, 本文方法能够将优化耗时减小至 4%以下, 并且随着离散节点数量的增加, 计算效率的提升更为显著。1 轨迹优化问题数学描述轨迹优化问题本质上属于最优控制问题, 以 Bolza 型最优控制问题为例, 可描述为:求解控制变量 u (t) R, 使得如下目标函数最小化状态方程为端点条件为路径约束为式中 f:RRRR, E:RRRRR, 1:RRRR。方程 (1) - (4) 所描述的问题称为连续 Bolza 型最优控制问题。2 基于 Runge-Kutta 格式的配点法离散首先利用积分变换 = (t-t 0)
15、/ (tf-t0) 将轨迹优化问题 (方程 (1) - (4) ) 变换至时间区间 (0, 1。假设单位区间0, 1上的 N 个离散区间的节点为式中 i称为节点或网格点, i在0, 1上可以均匀分布, 也可以非均匀分布。记 xi=x ( i) , ui=u ( i) , 对于状态方程, 基于 q 阶 Runge-Kutta (RK) 方法的离散格式为式中:t=t f-t0, hi= i+1- i, fij=f (xij, uij, ij;t0, tf) , xij, uij和 ij为中间变量, x ij由下式给出式中: ij= i+hi j, uij=u ( ij) . j, j, jl均为已
16、知常数并且满足0 1 2 q1。当 jl=0 (lj) 时, 离散格式为显式格式, 否则为隐式格式。采用类似的方法, 可将目标函数可离散化。常用的离散格式包括梯形格式 (q=2) , HermiteSimpson 格式 (q=3) , 以及经典四阶 Runge-Kutta格式 (q=4) 。由此离散得到的非线性规划问题是求解变量 X, U, U, t0和 tf使得如下目标函数最小并且满足如下约束式中常用的离散格式有梯形格式 (q=2) , HermiteSimpson 格式 (q=3, 简记 HS 格式) , 以及经典四阶 Runge-Kutta 格式 (q=4, 简记 RK 格式) 。以 H
17、S 格式为例, 该格式需要用到区间中点的变量和函数值, 为此将区间中点的控制变量作为优化变量, 并且在区间中点添加路径约束, 即HS 格式得到的 NLP 的优化变量为 , 目标函数为约束条件为其中在数值优化时, 为了使问题具有实际物理意义, 还需要添加时间差约束3 NLP 偏导数计算方法3.1 依赖关系矩阵在推导 NLP 一阶偏导数的稀疏特性时, 需要用到原始轨迹优化问题对自变量的依赖关系。由于状态方程、路径约束和目标函数 Lagrange 积分项都定义在整个时域区间, 因而本文将这三项对自变量的偏导数定义在一起, 其中 G1的每一项仍为矩阵, 以 为例, 易知, G 1为 (n+c+1) (
18、n+m+1) 维矩阵。通常情况下, G 1是稀疏矩阵。为了描述 G1的稀疏型, 定义如下 struct 函数记式中 struct (G1) 表示对 G1的每个元素进行 struct 运算。S 1表示 G1的稀疏型。为了得到 S1, 不需要计算 G1的每个元素的具体值, 只需要判断是否为 0。类似可以定义端点约束和目标函数的 Mayer 项对自变量的依赖关系矩阵和稀疏型易知, G 2为 (e+1) 2 (n+1) 维矩阵。3.2 变量记法前述离散格式将同一个节点处的变量或约束记为一个列向量, 这种记法与数值积分格式的形式一致, 但是不利于推导偏导数矩阵的稀疏特性。为此, 本文定义一种新的变量记法
19、, 将变量或约束的同一个分量在不同节点的值记为一个新的向量。以状态方程离散残差为例, 定义式中 i, j的下标表示第 i 节点, 第 2 个下标表示第 j 分量。易知, :, j为N1 向量。类似地可定义 :。利用这种新的变量记法, 可以可将 HS 格式离散得到的 NLP 问题重新描述为:求解优化变量 z, 使得如下目标函数最小化并且满足约束条件式中:目标函数 J (z) 的表达式参见方程, 优化变量 z 和约束函数 F (z) 的定义如下3.3 目标函数梯度目标函数梯度是指目标函数对优化变量的偏导数, 具体定义如下将目标函数写成矩阵乘积形式可得到式中:D=1, 1, , 1 1N, :, i
20、的表达式如下对角阵 h=diag (h0, h1, , hN-1) , 其中 hi (i=0, 1, , N-1) 为积分步长, 矩阵 D1和 D2定义如下D1和 D2均为 N (N+1) 矩阵, 其中空白元素为 0。应用向量链式求导规则, 可推导出目标函数对 , t0和 tf的偏导数分别为式中可见, NLP 的目标函数梯度可以分解为轨迹优化问题的目标函数和状态方程的偏导数。3.4 雅克比矩阵NLP 的雅克比矩阵定义为 NLP 的约束对优化变量的偏导数矩阵, 对于 HS 格式, 形式如下式中向量 F 和 z 的定义参见方程 (29) 。雅克比矩阵 GF的展开形式遵循向量求偏导数运算规则 (参见
21、方程 (21) ) 。下文推导 GF的数学表达式。3.4.1 状态方程离散残差约束的偏导数结合前述定义的变量记法, 将状态方程离散残差约束即方程写成如下形式矩阵 D1和 D2的定义与方程中的相同。将方程分别对 , t0和 tf求偏导数导得到式中3.4.2 路径约束的偏导数应用向量链式求导, 可以得到节点路径约束的偏导数。其中 已经是轨迹优化问题的偏导数, , 其余项为类似地, 可以推导出得到节点中点处的路径约束对优化变量 z 的偏导数, 其中直接从轨迹优化问题的偏导数中提取, 其余项如下3.4.3 端点约束的偏导数将端点约束对变量 z 求偏导可得到相应的表达式, 其中 已经是轨迹优化问题的偏导
22、数, , 其余项如下3.4.4 时间差约束的偏导数时间差约束为线性约束, 易知 , 其余项为可见, NLP 的雅克比矩阵可分解为轨迹优化问题的状态方程、路径约束、端点约束和时间约束的偏导数。计算出这些约束在离散节点和区间中点的偏导数 (对于 f 和 C) 以及端点处的偏导数 (对于 E 和 t) 之后, 采用本节的方法组装得到雅克比矩阵。对于 HS 格式, 以 N=4 为例, 其雅克比矩阵的稀疏型如图 1所示。其中, 空白元素表示恒为零, “”表示非零元素, “”表示可能不为零的元素。对于可能不为零的元素, 其具体取值与依赖关系矩阵有关, 以 为例4 仿真算例飞行器最短时间爬升问题最初由 Bryson 等22提出, 此后得到广泛研究2。该算例的特色之处是推力和气动力数据以离散表格形式给出, 与实际工程问题
