admin 针对采用直接法求解轨迹优化问题中精度和效率之间的矛盾，我们提出了一种基于二代小波轨迹优化节点自适应加密。该方法采用RK( Runge-Kutta)离散方法将原轨迹优化问题转化为非线性规划问题，并采用成熟的非线性规划算法求解，对控制或状态函数进行小波变换得到小波系数。
一、轨迹优化问题与离散方法
目前求解轨迹优化问题的数值方法主要有2种：间接法和直接法。间接法基于变分原理和Pontryagin极大值原理，将原轨迹优化问题转化为Hamilton边值问题(HBVP)，然后通过数值方法求解，间接法的主要优点是得到的解较精确并满足1阶最优条件。但是由于其收敛半径非常小，协（伴随）状态的初始值不易给出，限制了其工程应用．而直接法则将状态变量和控制变量对时间离散化，将原始优化问题转化为非线性规划问题。
RK离散方法由于其精度较高，在常微分方程的数值求解中应用广泛，因而也用于早期的间接和直接打靶法中并逐渐用于控制／状态参数化方法中。由于RK法不要求节点固定，因而易于进行节点自适应调整。
1、轨迹优化问题
考虑如下轨迹优化问题：确定控制函数u(t)∈Rm（和时间t0，tf），使得如下性能函数最小：
满足状态方程：
边界条件：
路径约束：
其中初始状态x0=x(to)、结束状态Xf= X(tf)。
这里假设初始时间to固定，结束时间tf固定或自由。
对于性能函数含有积分指标的问题和自由时间问题，均可以通过引入附加状态转化为固定时间问题。
为了便于程序实现，通过变换t=(t -to)／(tf-t0)将时间区间归一化，即t∈[0，1]。变换后式(1)～式(4)分别转化为：
其中△t=tf-t0。
2、Runge-Kutta离散方法与非线性规划问题
不失一般性，假设状态量和控制量在一组非均匀时间点T={ti：ti∈[0，1]）（下标i表示时间点序号）上进行离散，则状态方程离散的q级Runge-Kutta离散方法如下：
其中：
这里1≤m≤q，i=0，1，…，N-1，N为时间节点数．αml，βm，ρm为常数，满足O≤ρ1≤…≤1。
时间节点上的变量(xi，ui)必须同时在整个时间区间[O，1]上确定，称为全局变量；而节点之间的变量(xim，uim)位于时间区间[t，Ti+l]内，称为局部变量，为了使最优解满足相容性条件，本文采用经典4阶Runge-Kutta离散方法。对于4阶Runge-Kutta离散方法，设计变量中还需加入（ti+ti+1）／2时刻对应的局部变量Ui2=Ui3=Ui。定义如下集合：
通过以上离散处理，可将轨迹优化问题转化为如下NLP问题：
确定变量：
使得如下性能函数最优：
且满足约束条件：
对于式(1 9)～式(22)描述的非线性规划问题，可采用序列二次规划(SQP)算法求解，如求解大规模优化问题的软件包SNOPT。为了提高计算速度，计算过程中利用自动微分软件Intlab计算雅克比矩阵，并采用稀疏矩阵表示法。
二、二代小波节点自适应算法
利用上述Runge-Kutta离散方法可将轨迹优化问题转换为非线性规划问题求解，其优点是节点位置、时间步长可以变化，若可以在优化过程中动态调整节点则可满足一定精度情况下减小节点的数量，提高计算效率。由于小波系数幅值反应了函数（高阶导数）的局部光滑性，可以基于小波系数幅值的大小进行节点自动调整。
1、二代小波
二代小波是提升插值小波的统称，为某些函数空间上的Reisz基，同时具有时间和空间局部性和多项式消失矩。与经典小波相比，可以直接在区间上构造。由于构造第二代小波的数学推导较为复杂，为了节约篇幅，本文仅给出部分必要的结论。
函数空间L上的二代小波多分辨分析M由一组闭子空间M=（vjC|li∈J）（j表示分辨水平）组成，并满足如下条件：
1)vj∈Vj+l；
2) Ui∈Jvj在L中稠密；
3)对每个j∈J，尺度函数{φkj∈Kj）是驴的Reisz基（其中Kj是指标集）。
Sweldens给出了二代小波变换(second gen-eration wavelets transform，SGWT)的构建方法，并构造了二代小波函数。以一维小波为例，二代小波的小波系数碰和尺度系数cjk具有如下正变化关系式：
其逆变换关系式为：
式中的vklj，vkl-j是构造提升小波的插值系数，可以通过插值细分原理求取。
2、节点自适应算法
Vasilyev等提出了求解偏微分方程的二代小波配点法，包括基于小波的节点自适应方法和自适应节点上的导数计算方法，下面结合轨迹优化问题的求解，简要给出节点自适应方法。
考虑定义在闭区间Ω=[O，1]上的函数f(t)（可以是控制量或状态量），并在如下二分节点上构造插值小波：
其中节点tkj满足嵌套关系：
为了不产生歧义，gj中的节点称为背景节点，只有在小波变换中用到。根据二代小波构造方法，控制函数U(T)可以在分辨率J的每个水平上用尺度函数φkj(t)（k∈Kj）和小波函数φlj(t)(l∈lj)表示为：
小波系数绝对值的大小表示了函数的局部光滑性质，函数越光滑，则小波系数越小。对于局部突变的函数，其大部分小波系数是足够小的，仅在函数突变处较大，因此即使舍去大部分较小的小波系数，也可以保持对原函数保持较好的近似。由于尺度函数与背景节点tkj唯一对应，而小波φlj(t)与背景节t2l+1j唯一对应。所以进行小波分解后，在最精细分辨层，上的每个节点要么与小波函数对应要么与尺度函数对应，因此，如果相应的小波系数从近似关系式(28)中省略，则相应的背景节点可以删除，可得到计算节点ti。
将式(28)改写为两项之和：
其中：
表示小波系数幅值大于某个阀值ε的部分。
表示小波系数幅值小于某个阀值ε的部分。
Donoho给出了方程的截断误差：
其中系数Cl与函数uj(t)相关。
根据小波系数与节点的分层对应关系可以确定最少的自适应节点。但是根据这些节点上的函数值无法还原全部节点的函数值，因此还需要稍微扩大自适应节点的范围，以使得从自适应节点能够恢复函数uj≥(t)。一般来讲，上述过程增加的节点是非常少的。
利用基于小波的节点自适应算法获得新的自适应节点的过程如下：
1)获得控制量ukj(t)在节点gj≥的函数值；
2)根据式(23)和式(24)执行正向小波变换，获得尺度系数Ck0（k∈K）和每个分辨率水平上的小波系数dlj（l∈L，O≤j≤J-1）；
3)分析小波系数dlj，建立所有背景节点tkj的指标集L，L中包含小波系数|dlj|≥ε的所有节点及其相邻的小波系数(包括同一分辨水平j和下一分辨水平j+1)；
4)在L中增加所有尺度系数c2（k∈K）对应的节点；
5)构建新的计算节点GL>+l（即L），用于下一步的优化计算。
在自适应节点Gj≥上执行小波变换可以保证所有小波系数与uj≥(t)在背景节点Gr进行小波变换然后从自适应节点删除等于0的小波系数得到的结果一致。
基于二代小波的节点自适应算法输入参数包括：Jo，初始分辨率水平（控制初始节点个数）；Jmax，最大分辨率水平（控制背景节点个数，即最小节点间隔）；ε为小波系数幅值的阀值，终止条件：前后两次节点位置相同或达到预定的迭代次数。
计算流程如下：
1)根据初始分辨率Jo生成均匀节点gjo(或k))，以及该均匀节点上的状态函数x和控制函数u的初始猜测；
2)利用非线性规划算法SNOPT优化式(18)～式(22)给出的非线性规划问题，得到状态函数x和最优控制函数H；
3)根据节点自适应算法确定下一步优化的自适应节点，若存在多个控制函数，则将这些控制函数分别得到的自适应节点进行合并；
4)根据前一次优化计算的解x和M，利用小波插值计算新的自适应节点上的初始猜测，并进行循环迭代，直到满足终止条件。
需要说明的是：小波系数阀值的取值通常凭经验选取，一般可取为ε=α（Umax- Umin），对于间断函数，01可取0.005～0.O1，对于连续函数，可取0.001～0. 005。
三、数值算例
1、Breakwell问题
考虑如下Breakwell问题：
目标函数为：
状态方程为：
边界条件为：
对该问题采用小波方法进行节点自适应调整，取小波系数阀值ε=0．OOl（umax- umin）。优化的状态初始猜测为初末状态之间的线性函数，控制量的初始值为0。
表1给出了不同的Z取值时，优化的节点数目、目标函数误差|J-J*|以及控制函数误差|| u-ux||ω max ui一u*（ti）||（为数值解与解析解误差函数的最大值范数）。这里与高斯伪谱法Matlab软件包(GPOPS)c14]进行了对比，计算采用CPU为Intel双核1.86GHz，内存为1.5GB的台式机。从表中可以看出：目标函数误差的量级约为l0-7～l0-9，控制函数误差最大值范数的量级约为l0-3～l0-4。结果显示：与GPOPS相比，本文算法的节点数大约节省1O%，最优指标精度大约提高1个数量级，而计算时间相当。
l=0.1时的节点自适应过程如图1所示。
状态函数和控制函数变化曲线分别如图2和图3所示。图1中垂直线为解析解中控制函数导数不连续时间点，迭代所采用的节点个数依次为33，45，57，63和72。从图中可以看出：每次节点调整后更加集中于控制函数导数不连续处，且在两侧都有加密。图2中实线为状态量解析解，“*”和“+"分别为本文算法计算得到的状态x和状态v。从结果可知，最优解与解析解非常吻合。
给出了基于密度函数的节点自适应算法(density function-based mesh refinement al-gorithm，DENMRA)，并给出了该算例的结果。与DENMRA相比，基于小波的算法精度稍差，但是比SOCS( sparse optimal control software)的结果要好。虽然DENMRA精度较高，但是需要针对具体问题设计密度函数，对于复杂问题，具有一定的难度，而本文基于小波的节点自适应算法的优势在于只需简单设定小波系数阀值，即可自动在控制函数（或高阶导数）不连续处加密节点。
2、航天器姿态机动问题
算例1通过与高斯伪谱法对比验证了算法的有效性和精度，同时计算效率与高斯伪谱法相当。为了验证算法处理不连续控制函数的性能以及在工程问题是的适用性，下面利用本文算法求解给出的刚体航天器姿态机动问题。采用四元数描述的航天器姿态动力学方程如下：
其中四元数定义如下：
其中||q ||=1（这里为2范数）。ωT[ω1，ω2，ω3]为角速度。刚体惯量参数使用NASA XTE (X-ray tlming experiment)航天器参数，如下：
控制力矩约束为：
机动目标是使用最短的时间使航天器绕x轴转动150°，因此取如下边界条件：
采用基于二代小波的节点自适应算法进行优化，计算中小波系数阀值取ui=0.=005(uimax -Uimin)，其中i=l，2，3.经过8次节点调整，最终保留187个节点。计算硬件与算例1相同。最优机动时间为28. 6304 s，与文献结果一致，图4是控制量曲线，图5是四元数曲线，图6是角速度曲线，图7是节点自适应过程中节点位置图。需要说明的是：对于该问题，控制量维数为3，且相互独立，在进行节点自适应过程中，对3个控制量分别计算需要的节点，最后将节点合并，因此会出现在第1个控制量光滑处存在节点加密情况。
从优化结果可以看出：基于二代小波的节点自适应算法在控制函数突变处进行了合理加密，提高了对控制函数逼近的精度，由于大部分节点分布在控制函数间断处，可提高控制开关处的精度，而在其他区域则可用较少节点。结果说明本文算法可用于具有多个控制函数、多个开关的复杂Bang-Bang控制问题的节点自适应。对于该类问题，若节点个数较多，可适当增大小波系数阀值。
小知识之小波小波（Wavelet)这一术语，顾名思义，“小波”就是小区域、长度有限、均值为0的波形。所谓“小”是指它具有衰减性；而称之为“波”则是指它的波动性，其振幅正负相间的震荡形式。与Fourier变换相比，小波变换是时间（空间）频率的局部化分析，它通过伸缩平移运算对信号（函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，解决了Fourier变换的困难问题，成为继Fourier变换以来在科学方法上的重大突破。有人把小波变换称为“数学显微镜”。
