【附源码和详细的公式推导】Minimum Snap轨迹生成，闭式求解Minimum Snap问题，机器人轨迹优化，多项式轨迹路径生成与优化

该篇博客内容较多，工作量也很大，难免在理解上表达上有错误，如有发现烦请指教。如有问题在博客中留言，或者github的仓库中提Issues都行，看到后我会尽快回复。

Minimum Snap源代码：Minimum Snap(如果代码有用，别忘记给我点Star呀)

主要内容：

介绍Minimum Snap轨迹生成方法的由来以及用处；
详细推导Minimum Snap轨迹生成的数学原理；
对Minimum Snap问题进行闭式求解；
使用C++代码进行算法复现，达到甚至超过论文中的实验结果。

1.什么是Minimum Snap？

考虑一下这样的运动规划问题，上图的两个五角星分别是起点和终点，我们通过某一种路径搜索方式(A star、RRT star等等)，找到了一段轨迹路径，搜索出来的路径是不考虑运动学约束的，因此轨迹上会出现明显的不光滑点，可以想象，机器人不可能在某一点出现运动突变，如果是万向轮的机器人要想按照这个轨迹走，它必须每走完一段路，就要停下来，然后旋转一个角度对着路径然后再加速往前走。这样很明显浪费很多时间和效率。如果是阿克曼的机器人那么它就无法按照这个路径进行运动。

但是如果我们以轨迹中的一些航迹点为基础，画出一条光滑的连接每一个航迹点曲线，那么这条曲线明显就非常的适合机器人去运动。

因此，基于上述这个思路，学者们提出了一些方法去解决这个问题，其中文献【1】提出了一种通过多项式来拟合这个轨迹的方法，可以使得轨迹是光滑并且满足机器人动力学约束。该方法首先在无人机中提出，但是已经被扩展到了很多机器人的运动规划中，后续基于他们工作的研究也非常的多，因此这篇文章的方法十分值得去深入学习。

本篇博客就是总结一些以这篇论文为基础的一些工作，主要关注的是基于Minimum Snap的轨迹生成。

2.一段多项式轨迹

对于空间中的两个点，我们可以使用直线进行拟合，对于三个点我们可以使用二次函数进行拟合。更广义一点来说，一次函数可以拟合两个参数点，二次函数可以拟合三个参数点，NNN次函数也可以拟合N+1N+1N+1参数点。

想象一下如果我们知道了一段一维轨迹，并且我们以一个五次多项式去拟合它：
P(t)=p0t0+p1t1+⋯+p5t5(1)P(t) = p_0t^0+p_1t^1+\cdots+p_5t^5 \tag{1}P(t)=p0t0+p1t1+⋯+p5t5(1)

上式（1）中有6个未知的系数，如果我们要求解出它，我们至少需要6个方程。一个很明显的思路是我们可以取6个位置点进行拟合。多数情况下我们不光关心机器人的位置，对于它的速度和加速度我们同样很关心，因此我们如果选取轨迹头和尾的位置以及速度和加速度，那么我们也可以唯一确定这个多项式，而且机器人还能以我们想要的速度和加速度开始和结束，这样不是更好吗？而后一种思路正是论文中所采用的。

3.多段多项式轨迹

在第2小节，介绍了一段轨迹使用多项式去拟合的问题，那么如果是多段轨迹，我们又该如何操作呢？

也许你可能会想到多种方式去拟合这整条轨迹,在此我就不做扩展.

直接引用论文中思路进行介绍，论文中作者提出对于一个经过多个中间航迹点的轨迹，分别使用多个轨迹段进行拟合，也就是对于上述T0−T1−T2−T3T_0-T_1-T_2-T_3T0−T1−T2−T3这条轨迹，我们将其分为三段分别进行轨迹拟合，分别拟合轨迹段T0−T1,T1−T2,T2−T3T_0-T_1,T_1-T_2,T_2-T_3T0−T1,T1−T2,T2−T3。每一段轨迹的拟合沿用第2小节中对于一段轨迹拟合的思路。

但是与一段轨迹拟合有一些区别，多数情况下，我们很难去给定机器人在经过中间航迹点的速度、加速度、Jerk等信息，于是我们就认为中间点的速度、加速度等信息是一个可以自由变化的状态，那么对于这种中间状态不确定的问题，我们可以构建一个优化问题，使得中间状态取得某一个值时，我们制定的代价函数达到最小值，那么我们一方面就不去设计机器人以什么样的状态通过中间点，另一方面还可以通过优化代价函数得到一个最优的轨迹。

上面的文字读起来可能会有一些难理解，以上图所示的轨迹为例，假设我们已知T0T_0T0和T3T_3T3的位置、速度、加速度等运动信息，并且我们还知道T1T_1T1和T2T_2T2的位置信息，而对于T2T_2T2和T3T_3T3的速度、加速度等信息，我们认为它是一个变量，可以任意变化的量，但是也不是毫无目的的变化，一般我们会设置一个目标代价函数，找到这几个变量在取一定值时使得我们设计的代价函数取得最小值。更通俗的说：就是这几个变量取一定值的时候可以使机器人能量耗费最少，或者机器人运动的转动更少等等。

听起来这个想法非常好，那它是否可行呢？答案是肯定的！

以上就是Minimum Snap算法的最核心思路！

由于四旋翼无人机以及双轮差速轮机器人(不准确的简化)，它们轨迹在各个坐标轴上是独立的，因此我们可以对其分别进行轨迹拟合。也就是说我们可以分别对它们在x,y,zx,y,zx,y,z轴上进行路径生成，然后直接将三个轴合成就可以得到一个完整的空间轨迹。
因此后续算法细节介绍中，我将不再强调轨迹是一维、二维还是三维，均以一维来进行介绍，你可以很容易将其推广到任意维度。

4. Minimum Snap

后续的数学表达中，有时出现向量和标量字母混淆的问题，但是从文章前后逻辑的角度而言并不会引起歧义，我就不再细致的修改了，哈哈，还请留意！

4.1 Minimum Snap的代价函数

我们已经对基于多项式的轨迹生成有了一些感性的认识，下面我们就从数学的层面对其进行具体建模。

类似于上图的轨迹，它是一个一维的轨迹，横轴是时间ttt，纵轴是位置Pm(t)P_m(t)Pm(t)，我们以NNN次多项式来拟合每一个分段轨迹，那么对于第mmm段轨迹:
Pm(t)=pm,0t0+pm,1t1+⋯+pm,NtN=[pm,0pm,1⋯pm,N][t0t1⋮tN]P_m(t) = p_{m,0}t^0+p_{m,1}t^1+\cdots+p_{m,N}t^N \\ = \left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,N} \end{matrix} \right] \left[ \begin{matrix} t^0 \\ t^1 \\ \vdots\\ t^N \end{matrix} \right] Pm(t)=pm,0t0+pm,1t1+⋯+pm,NtN=[pm,0pm,1⋯pm,N]⎣⎢⎢⎢⎡t0t1⋮tN⎦⎥⎥⎥⎤

那么，这个多项式就有N+1N+1N+1个未知系数。

对其进行一阶求导，相当于得到了它的速度函数：
dP(t)dt=[pm,0pm,1⋯pm,N][01⋮NtN−1]\frac{dP(t)}{dt} = \left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,N} \end{matrix} \right] \left[ \begin{matrix} 0 \\ 1 \\ \vdots \\ Nt^{N-1} \end{matrix} \right]dtdP(t)=[pm,0pm,1⋯pm,N]⎣⎢⎢⎢⎡01⋮NtN−1⎦⎥⎥⎥⎤

对其进行二阶求导，相当于得到的了它的加速度函数：
d2P(t)dt2=[pm,0pm,1⋯pm,N][01∗02∗1⋮N(N−1)tN−2]\frac{d^2P(t)}{dt^2} = \left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,N} \end{matrix} \right] \left[ \begin{matrix} 0 \\ 1* 0 \\ 2 * 1 \\ \vdots \\ N(N-1)t^{N-2} \end{matrix} \right] dt2d2P(t)=[pm,0pm,1⋯pm,N]⎣⎢⎢⎢⎢⎢⎡01∗02∗1⋮N(N−1)tN−2⎦⎥⎥⎥⎥⎥⎤

三阶求导，得到Jerk函数：
d3P(t)dt3=[pm,0pm,1⋯pm,N][0∗0∗01∗0∗02∗1∗0⋮N(N−1)(N−2)tN−3]\frac{d^3P(t)}{dt^3} = \left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,N} \end{matrix} \right] \left[ \begin{matrix} 0 * 0 * 0\\ 1* 0 *0 \\ 2 * 1 * 0 \\ \vdots \\ N(N-1)(N-2)t^{N-3} \end{matrix} \right] dt3d3P(t)=[pm,0pm,1⋯pm,N]⎣⎢⎢⎢⎢⎢⎡0∗0∗01∗0∗02∗1∗0⋮N(N−1)(N−2)tN−3⎦⎥⎥⎥⎥⎥⎤

四阶求导，则得到了Snap函数：
d4P(t)dt4=[pm,0pm,1⋯pm,N][0∗0∗0∗01∗0∗0∗02∗1∗0∗03∗2∗1∗0⋮N(N−1)(N−2)(N−3)tN−4]\frac{d^4P(t)}{dt^4} = \left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,N} \end{matrix} \right] \left[ \begin{matrix} 0 * 0 * 0 *0\\ 1* 0 *0 *0 \\ 2 * 1 * 0 *0 \\ 3*2*1*0 \\ \vdots \\ N(N-1)(N-2)(N-3)t^{N-4} \end{matrix} \right] dt4d4P(t)=[pm,0pm,1⋯pm,N]⎣⎢⎢⎢⎢⎢⎢⎢⎡0∗0∗0∗01∗0∗0∗02∗1∗0∗03∗2∗1∗0⋮N(N−1)(N−2)(N−3)tN−4⎦⎥⎥⎥⎥⎥⎥⎥⎤

以上我们就写出了一段轨迹的多项式函数，以及它们对应的导数。

还记得前面铺垫吗？为了得到中间点的运动状态，我们构建一个代价函数，然后求解最优中间运动状态。

在论文【1】中，对于最优的轨迹生成问题，它构建了一个如下的代价含数JmJ_mJm，对应的是一段多项式轨迹的四阶导数(Snap)平方的积分。

后续的过程中我会交叉使用最小化Jerk平方的积分以及最小化Snap平方的积分，看到后不要感到困惑！你甚至还可以最小化七阶导数平方的积分！

Jm=∫Tm−1Tm(d4P(t)dt4)2=∫Tm−1Tm[pNpN−1⋮p0]T[N(N−1)(N−2)(N−3)tN−4⋮3∗2∗1∗02∗1∗0∗01∗0∗0∗00∗0∗0∗0][N(N−1)(N−2)(N−3)tN−4⋮3∗2∗1∗02∗1∗0∗01∗0∗0∗00∗0∗0∗0]T[pNpN−1⋮p0]=[pNpN−1⋮pi⋮p0]T[N(N−1)(N−2)(N−3)N(N−1)(N−2)(N−3)(Tm−Tm−1)N+N−7N+N−7N(N−1)(N−2)(N−3)(N−1)(N−2)(N−3)(N−4)(Tm−Tm−1)N+N−1−7N+N−1−7⋯⋮⋮⋮⋯i(i−1)(i−2)(i−3)l(l−1)(l−2)(l−3)(Tm−Tm−1)i+l−7i+l−7⋯⋮⋮⋮00⋯][pNpN−1⋮pl⋮p0]=PmTQmPmJ_m = \int_{T_{m-1}}^{T_{m}}(\frac{d^4P(t)}{dt^4})^2 \\ = \int_{T_{m-1}}^{T_{m}} \left[ \begin{matrix} p_N \\ p_{N-1} \\ \vdots \\ p_0 \end{matrix} \right]^T \left[ \begin{matrix} N(N-1)(N-2)(N-3)t^{N-4} \\ \vdots \\3*2*1*0\\ 2 * 1 * 0 *0\\ 1* 0 *0 *0 \\0 * 0 * 0 *0 \end{matrix} \right] \left[ \begin{matrix} N(N-1)(N-2)(N-3)t^{N-4} \\ \vdots \\3*2*1*0\\ 2 * 1 * 0 *0\\ 1* 0 *0 *0 \\0 * 0 * 0 *0 \end{matrix} \right]^T \left[ \begin{matrix} p_N \\ p_{N-1} \\ \vdots \\ p_0 \end{matrix} \right] \\= \left[ \begin{matrix} p_N \\ p_{N-1} \\ \vdots \\ p_i \\ \vdots\\ p_0 \end{matrix} \right]^T \left[ \begin{matrix} \frac{N(N-1)(N-2)(N-3)N(N-1)(N-2)(N-3)(T_m-T_{m-1})^{N+N-7}}{N+N-7} & \frac{N(N-1)(N-2)(N-3)(N-1)(N-2)(N-3)(N-4)(T_m-T_{m-1})^{N+N-1-7}}{N+N-1-7} & \cdots \\ \vdots & \vdots & \vdots \\ \cdots & \frac{i(i-1)(i-2)(i-3)l(l-1)(l-2)(l-3)(T_m-T_{m-1})^{i+l-7}}{i+l-7} & \cdots \\ \vdots& \vdots&\vdots\\ 0 & 0 &\cdots \end{matrix} \right] \left[ \begin{matrix} p_N \\ p_{N-1} \\ \vdots \\ p_l \\ \vdots\\ p_0 \end{matrix} \right] \\= {P_m}^TQ_m{P_m} Jm=∫Tm−1Tm(dt4d4P(t))2=∫Tm−1Tm⎣⎢⎢⎢⎡pNpN−1⋮p0⎦⎥⎥⎥⎤T⎣⎢⎢⎢⎢⎢⎢⎢⎡N(N−1)(N−2)(N−3)tN−4⋮3∗2∗1∗02∗1∗0∗01∗0∗0∗00∗0∗0∗0⎦⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎡N(N−1)(N−2)(N−3)tN−4⋮3∗2∗1∗02∗1∗0∗01∗0∗0∗00∗0∗0∗0⎦⎥⎥⎥⎥⎥⎥⎥⎤T⎣⎢⎢⎢⎡pNpN−1⋮p0⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡pNpN−1⋮pi⋮p0⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤T⎣⎢⎢⎢⎢⎢⎢⎢⎡N+N−7N(N−1)(N−2)(N−3)N(N−1)(N−2)(N−3)(Tm−Tm−1)N+N−7⋮⋯⋮0N+N−1−7N(N−1)(N−2)(N−3)(N−1)(N−2)(N−3)(N−4)(Tm−Tm−1)N+N−1−7⋮i+l−7i(i−1)(i−2)(i−3)l(l−1)(l−2)(l−3)(Tm−Tm−1)i+l−7⋮0⋯⋮⋯⋮⋯⎦⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡pNpN−1⋮pl⋮p0⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=PmTQmPm

上式推导思路：上式思路：如果两个向量相乘为常数，那么它们乘积等于它们乘积的转置。如a⃗∗b⃗=c\vec{a} * \vec{b} = ca∗b=c，则(a⃗∗b⃗)T=b⃗T∗a⃗T=c(\vec{a} * \vec{b})^T = \vec{b}^T * \vec{a}^T = c(a∗b)T=bT∗aT=c

注意：上式中(Tm−Tm−1)N+N−7(T_m-T_{m-1})^{N+N-7}(Tm−Tm−1)N+N−7实际就是这段轨迹运动是所用的时长，为了方便，后续不再采用每段轨迹的绝对时间，而是采用相对时间，也就是对于任意一段轨迹认为它们的时间都是(0，Tm)(0，T_m)(0，Tm)，那么就可以将(Tm−Tm−1)N+N−7(T_m-T_{m-1})^{N+N-7}(Tm−Tm−1)N+N−7简记为TN+N−7T^{N+N-7}TN+N−7，这一点请理解并留意。并且假设我们事先已经知道了各段轨迹的时间。

小例子：

为了大家更容易理解，我这里举一个实际的例子向大家展示一下QQQ矩阵的实际结果：论文中给出的是最小Snap，为了形成差异性我这里最小化Jerk，我们使用5次多项式，则N=5N=5N=5，它一共有6个未知数，那么它的QQQ矩阵如下：
Q=[720360120000360192720001207236000000000000000000000]Q = \left[ \begin{matrix} 720 & 360 & 120 & 0 & 0 & 0\\ 360 &192 & 72 & 0 & 0 & 0\\ 120 & 72 & 36 & 0 & 0 & 0\\ 0 &0 &0 & 0 & 0 & 0\\ 0 &0& 0& 0 & 0 & 0\\ 0 &0 &0 &0 &0 &0 \end{matrix} \right]Q=⎣⎢⎢⎢⎢⎢⎢⎡720360120000360192720001207236000000000000000000000⎦⎥⎥⎥⎥⎥⎥⎤

跳转到生成Q矩阵的代码

总的代价函数
假设我们有MMM段轨迹，那么每一段轨迹都是一个高阶的多项式函数，写出代价函数和如下：
J=J0+J1+⋯+JM=∑m=0MPm⃗TQmPm⃗=[P0P1⋮PM]T[Q00⋯⋯00Q10⋯0⋮⋮⋮⋯⋮000⋯QM][P0P1⋮PM]J = J_0 + J_1 + \cdots + J_M \\ = \sum_{m=0}^M \vec{P_m}^TQ_m\vec{P_m} \\ = \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right]^T \left[ \begin{matrix} Q_0 & 0 & \cdots & \cdots & 0 \\ 0 & Q_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots &\cdots & \vdots \\ 0 & 0 & 0 & \cdots & Q_M \end{matrix} \right] \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right] J=J0+J1+⋯+JM=m=0∑MPmTQmPm=⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤T⎣⎢⎢⎢⎡Q00⋮00Q1⋮0⋯0⋮0⋯⋯⋯⋯00⋮QM⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤

到此我们以多项式为基础推导出了代价函数的具体形式，如果我们直接对上式进行优化得到其最小值，看似确实完成了我们的目标，但是这其中忽略了一个问题，那就是可能优化出来的结果是，路标点之间的速度或者加速度不连续，显然这是可能的，因为我们没有限制航迹点在前后两端轨迹中的状态是否一致，也就没有限制要连续。但是我们的机器人只能平滑的运动，前后状态一定是要连续的，因此我们需要加一些额外的约束去限制连续性。

4.2 连续性约束

首先，我们需要限制中间的航迹点在两段轨迹的接合处是左右两边状态要连续。如果我们要最小化的目标是Jerk，那么我们必须保证接合处三阶导数是连续的。因为我们要对Jerk进行积分，那么多项式在求导3次之后必须要连续，所以要确保0阶、1阶、2阶导数都是光滑的，也就是保证求导到加速度函数依然是光滑的，这一点就体现了多项式轨迹的优势，它天然就是一个高阶光滑的函数。

写成数学的形式，那么就是：
Pm(k)(Tm)=Pm+1(k)(Tm)P^{(k)}_m(T_m) = P^{(k)}_{m+1}(T_m) Pm(k)(Tm)=Pm+1(k)(Tm)
上式的意义就是轨迹段PmP_mPm和Pm+1P_{m+1}Pm+1在TmT_{m}Tm时刻具有相同的kkk阶导数。

对上式进一步推导：
Pm(k)(Tm)=Pm+1(k)(Tm)⇒∑i≥ki!(i−k)!Tmi−kpm,i−∑l≥kl!(l−k)!Tml−kpm+1,l=0⇒[Am−Am+1][PmPm+1]=0P^{(k)}_m(T_m) = P^{(k)}_{m+1}(T_m) \\ \Rightarrow \sum_{i\ge k}\frac{i!}{(i-k)!}T_m^{i-k}p_{m,i}-\sum_{l\ge k}\frac{l!}{(l-k)!}T_m^{l-k}p_{m+1,l} = 0 \\ \Rightarrow \left[ \begin{matrix} A_m & -A_{m+1} \end{matrix} \right] \left[ \begin{matrix} P_m \\ P_{m+1} \end{matrix} \right] = 0 Pm(k)(Tm)=Pm+1(k)(Tm)⇒i≥k∑(i−k)!i!Tmi−kpm,i−l≥k∑(l−k)!l!Tml−kpm+1,l=0⇒[Am−Am+1][PmPm+1]=0

4.3 微分约束

连续性约束保证了轨迹的光滑性。我们可以想一下虽然限制了轨迹的光滑性，但是我们并没有限制轨迹必须通过某一些点，或者以一个确定的状态通过某一些中间的航迹点。也就是上述的这些多项式系数中，有一些是已知的。但是大多数情况下，我们都是对轨迹的位置有一些限制的，并不能随意生成，所以还需要给中间的航迹点加一些限制。包括位置、速度、加速度、Jerk、Snap等等的限制。

首先我们考虑很实际的问题，机器人的起点和终点我们肯定事先是知道的，并且机器人以什么样的速度或者加速度开始，我们也是可以设置的，机器人的终点速度和加速度也是可以事先给定的。机器人中间航迹点的位置也是需要事先给定的，但是机器人中间点的速度和加速度等信息，我们没法去约束它，因为很难去确定机器人到达某一个点的速度或者加速度，所以对于中间点的速度和加速度，或者更高阶的Jerk甚至Snap等，我们不对其进行微分限制，让它们只具有连续性约束，然后让优化器给它们分配最优的状态。

那么用数学方法去表述上面的内容，如下：

起始点和终点的位置、速度、加速度等使用确定值限制：
P0(0)=确定值P0(1)(0)=确定值P0(2)(0)=确定值PM(TM)=确定值PM(1)(TM)=确定值PM(2)(TM)=确定值⋮P_0(0)=确定值 \\ P^{(1)}_0(0)=确定值 \\ P^{(2)}_0(0)=确定值 \\ P_M(T_M)=确定值 \\ P_M^{(1)}(T_M)=确定值 \\ P_M^{(2)}(T_M)=确定值 \\ \vdots P0(0)=确定值P0(1)(0)=确定值P0(2)(0)=确定值PM(TM)=确定值PM(1)(TM)=确定值PM(2)(TM)=确定值⋮

按照之前的思路，如果我们最小化Snap，那我们最少需要给起点终点位置、速度、加速度、Jerk一个确定值，但是实际工程中，对于高于加速度之后的导数，我们都将其设置为0，也就是你最小化Snap，但是我只给位置、速度、加速度设置值，而Jerk就默认设置为0。更高阶的情况也如此。

中间航迹点使用位置限制：
P0(T0)=确定值P1(T0)=确定值P1(T1)=确定值⋮P_0(T_0) = 确定值 \\ P_1(T_0) = 确定值\\ P_1(T_1) = 确定值\\ \vdots P0(T0)=确定值P1(T0)=确定值P1(T1)=确定值⋮
画图表示

把上面的过程写成数学等式(以任意一段多项式轨迹为例)：
[Tm−1NTm−1N−1⋯Tm−10TmNTmN−1⋯Tm0NTm−1N−1(N−1)Tm−1N−2⋯0⋮⋮⋮⋮⋯i!(i−k)!Tm−1i−k⋯⋯⋯i!(i−k)!Tmi−k⋯⋯][pNpN−1⋮p0]=[ρTm−1ρTmνTm−1νTmαTm−1αTm⋮]⇒AmP⃗m=d⃗m\left[ \begin{matrix} T_{m-1}^N & T_{m-1}^{N-1} &\cdots & T_{m-1}^0 \\ T_m^N & T_m^{N-1} & \cdots & T_m^0 \\ NT_{m-1}^{N-1} & (N-1)T_{m-1}^{N-2} & \cdots & 0 \\ \vdots & \vdots & \vdots &\vdots \\ \cdots & \frac{i!}{(i-k)!}T_{m-1}^{i-k} & \cdots & \cdots \\ \cdots & \frac{i!}{(i-k)!}T_{m}^{i-k} & \cdots & \cdots \end{matrix} \right] \left[ \begin{matrix} p_N \\ p_{N-1} \\ \vdots \\ p_0 \end{matrix} \right] =\left[ \begin{matrix} \rho_{T_{m-1}} \\ \rho_{T_m} \\ \nu_{T_{m-1}} \\ \nu_{T_m} \\ \alpha_{T_{m-1}} \\ \alpha_{T_{m}} \\ \vdots \end{matrix} \right] \\ \Rightarrow A_m \vec{P}_m = \vec{d}_m ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡Tm−1NTmNNTm−1N−1⋮⋯⋯Tm−1N−1TmN−1(N−1)Tm−1N−2⋮(i−k)!i!Tm−1i−k(i−k)!i!Tmi−k⋯⋯⋯⋮⋯⋯Tm−10Tm00⋮⋯⋯⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡pNpN−1⋮p0⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ρTm−1ρTmνTm−1νTmαTm−1αTm⋮⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⇒AmPm=dm
以上就是对一段轨迹的头和尾进行微分限制，对于中间的航迹点。通常做法，对于它们高于0阶之后的导数不进行微分限制。

如果将连续性约束和微分约束进行合并，则可以得到如下式子：
[A0000⋯0A0−A100⋯00A100⋯00A1−A20⋯0⋮⋮⋮⋮⋮⋮00000AM][P0P1P2P3⋮PM]=[d00d10⋮dM]\left[ \begin{matrix} A_0 & 0 & 0 & 0 & \cdots & 0 \\ A_0 & -A_1 & 0 & 0 &\cdots &0 \\ 0 & A_1 & 0 & 0 & \cdots & 0 \\ 0 & A_1 & -A_2 & 0 &\cdots &0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 &0 &0&0&0&A_M \end{matrix} \right] \left[ \begin{matrix} P_0 \\ P_1 \\P_2\\ P_3 \\ \vdots \\ P_M \\ \end{matrix} \right] = \left[ \begin{matrix} d_0 \\ 0 \\d_1\\ 0 \\ \vdots \\ d_M \\ \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎢⎡A0A000⋮00−A1A1A1⋮0000−A2⋮00000⋮0⋯⋯⋯⋯⋮00000⋮AM⎦⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎡P0P1P2P3⋮PM⎦⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎡d00d10⋮dM⎦⎥⎥⎥⎥⎥⎥⎥⎤
最终，我们构建了一个带约束的QP问题
min⁡[P0P1⋮PM]T[Q00⋯⋯00Q10⋯0⋮⋮⋮⋱⋮000⋯QM][P0P1⋮PM]=PTQPs.t.[A0000⋯0A0−A100⋯00A100⋯00A1−A20⋯0⋮⋮⋮⋮⋮⋮00000AM][P0P1P2P3⋮PM]=[d00d10⋮dM]\min \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right]^T \left[ \begin{matrix} Q_0 & 0 & \cdots & \cdots & 0 \\ 0 & Q_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 0 & 0 & 0 & \cdots & Q_M \end{matrix} \right] \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right] = P^TQP \\ s.t. \left[ \begin{matrix} A_0 & 0 & 0 & 0 & \cdots & 0 \\ A_0 & -A_1 & 0 & 0 &\cdots &0 \\ 0 & A_1 & 0 & 0 & \cdots & 0 \\ 0 & A_1 & -A_2 & 0 &\cdots &0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 &0 &0&0&0&A_M \end{matrix} \right] \left[ \begin{matrix} P_0 \\ P_1 \\P_2\\ P_3 \\ \vdots \\ P_M \\ \end{matrix} \right] = \left[ \begin{matrix} d_0 \\ 0 \\d_1\\ 0 \\ \vdots \\ d_M \\ \end{matrix} \right] min⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤T⎣⎢⎢⎢⎡Q00⋮00Q1⋮0⋯0⋮0⋯⋯⋱⋯00⋮QM⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤=PTQPs.t.⎣⎢⎢⎢⎢⎢⎢⎢⎡A0A000⋮00−A1A1A1⋮0000−A2⋮00000⋮0⋯⋯⋯⋯⋮00000⋮AM⎦⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎡P0P1P2P3⋮PM⎦⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎡d00d10⋮dM⎦⎥⎥⎥⎥⎥⎥⎥⎤
构建好上述QP问题，然后将其带入第三方求解库就可以求到一个全局最优解。

推荐使用的QP问题求解器：OOQP、Mosek

5. 闭式求解Minimum Snap

对于上面构造的QP问题，它的约束是一个等式约束，实际上可以通过一些数学的带入方法将约束条件转换到代价函数中，那么就变成了无约束的优化问题，那么就可以进行闭式求解。闭式求解的好处的求解速度更快，数值稳定性更高。

上面的代价函数中，我们优化的变量是多项式的系数，我们可以观察一下其中一段多项式的某一项pitip_it^ipiti，由于ppp没有实际的物理意义，所以在优化的时候可以会出现一些特别小的值，此时可能就会计算精度的影响，导致数值不稳定。

为了解决这个系数不稳定的问题，在论文【2】中提出了将系数通过一个映射矩阵转换成轨迹端点的导数，也就是说我们不再优化轨迹的系数，而是对轨迹在端点处的位置、速度、加速度、Jerk等进行优化，由于这些量都是有实际物理含义，不会出现太离谱的数值，所以在一定程度上是比较稳定的。

对于某一段轨迹，它的多项式表示如下：
Pm(t)=pm,0t0+pm,1t1+⋯+pm,NtN=[pm,0pm,1⋯pm,N][t0t1⋮tN]P_m(t) = p_{m,0}t^0+p_{m,1}t^1+\cdots+p_{m,N}t^N =\left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,N} \end{matrix} \right] \left[ \begin{matrix} t^0 \\ t^1 \\ \vdots\\ t^N \end{matrix} \right] Pm(t)=pm,0t0+pm,1t1+⋯+pm,NtN=[pm,0pm,1⋯pm,N]⎣⎢⎢⎢⎡t0t1⋮tN⎦⎥⎥⎥⎤
根据我们之前说的多项式个数与阶数的关系，我们知道对于上式未知数一共有N+1N+1N+1个，对应的两个端点出可以提供N+1N+1N+1个微分约束和连续性约束。

以最小化Snap为例，由于我们是对多项式进行了4阶求导，那么可以写出多项式的向量形式如下：
Pm(t)=pm,0t0+pm,1t1+⋯+pm,Nt7=[pm,0pm,1⋯pm,7][t0t1⋮t7]P_m(t) = p_{m,0}t^0+p_{m,1}t^1+\cdots+p_{m,N}t^7 =\left[ \begin{matrix} p_{m,0} & p_{m,1} & \cdots & p_{m,7} \end{matrix} \right] \left[ \begin{matrix} t^0 \\ t^1 \\ \vdots\\ t^7 \end{matrix} \right] Pm(t)=pm,0t0+pm,1t1+⋯+pm,Nt7=[pm,0pm,1⋯pm,7]⎣⎢⎢⎢⎡t0t1⋮t7⎦⎥⎥⎥⎤

每段多项式有8个未知的系数，因此我们需要对轨迹的端点提供位置、速度、加速度、jerk约束(包括微分约束和连续性约束)，那么我们8个未知系数就可以使用两端的位置、速度、加速度、Jerk运动量来表示，每段轨迹可以提供8个运动量，因此我们需要找到某一种方式将系数映射到成运动量(位置、速度、加速度等)。

我们不妨来探索一下，这次为了表述方便和公式写法的简化，我以最小化Jerk为例，那么根据前面的结论，我们需要构建次数为5次的多项式，那么它将有6个未知的系数，我们需要提供轨迹两端的位置、速度、加速度约束。下面我们对轨迹进行多次求导，看看是否可以从中找到这种映射关系：
[Pm(t)Pm(1)(t)Pm(2)(t)]=[pm,5t5pm,4t4pm,3t3pm,2t2pm,1t1pm,0t05pm,5t44pm,4t33pm,3t22pm,2t1pm,1t005∗4pm,5t34∗3pm,4t23∗2pm,3t12∗1pm,2t000]=[t5t4t3t2t1t05t44t33t22t1t005∗4t34∗3t23∗2t12∗1t000][pm,5pm,4pm,3pm,2pm,1pm,0]\left[ \begin{matrix} P_m(t) \\ P_m^{(1)}(t) \\ P_m^{(2)}(t) \end{matrix} \right] = \left[ \begin{matrix} p_{m,5}t^5 & p_{m,4}t^4 & p_{m,3}t^3 & p_{m,2}t^2 & p_{m,1}t^1 & p_{m,0}t^0 \\ 5p_{m,5}t^4 & 4p_{m,4}t^3 & 3p_{m,3}t^2 & 2p_{m,2}t^1 & p_{m,1}t^0 & 0 \\ 5*4p_{m,5}t^3 & 4*3p_{m,4}t^2 & 3*2p_{m,3}t^1 & 2*1p_{m,2}t^0 & 0 & 0 \end{matrix} \right] \\= \left[ \begin{matrix} t^5 & t^4 & t^3 & t^2 & t^1 & t^0 \\ 5t^4 & 4t^3 & 3t^2 & 2t^1 & t^0 & 0 \\ 5*4t^3 & 4*3t^2 & 3*2t^1 & 2*1t^0 & 0 & 0 \end{matrix} \right] \left[ \begin{matrix} p_{m,5}\\ p_{m,4}\\ p_{m,3}\\ p_{m,2} \\ p_{m,1} \\ p_{m,0} \end{matrix} \right] ⎣⎢⎡Pm(t)Pm(1)(t)Pm(2)(t)⎦⎥⎤=⎣⎡pm,5t55pm,5t45∗4pm,5t3pm,4t44pm,4t34∗3pm,4t2pm,3t33pm,3t23∗2pm,3t1pm,2t22pm,2t12∗1pm,2t0pm,1t1pm,1t00pm,0t000⎦⎤=⎣⎡t55t45∗4t3t44t34∗3t2t33t23∗2t1t22t12∗1t0t1t00t000⎦⎤⎣⎢⎢⎢⎢⎢⎢⎡pm,5pm,4pm,3pm,2pm,1pm,0⎦⎥⎥⎥⎥⎥⎥⎤

上式中的Pm(t),Pm(1)(t),Pm(2)(t)P_m(t),P_m^{(1)}(t),P_m^{(2)}(t)Pm(t),Pm(1)(t),Pm(2)(t)其实分别就是对应轨迹在ttt时刻的位置、速度、加速度

似乎已经很接近答案了，我们只要将分段轨迹两端的时间都带入上式，不就找到了多项式系数到运动微分的映射关系了吗！(与之前一样我们使用相对时间)
[P(0)P(1)(0)P(2)(0)P(T)P(1)(T)P(2)(T)]=[0504030201005∗044∗033∗022∗010005∗4∗034∗3∗023∗2∗012∗1∗0000T5T4T3T2T1T05∗T44∗T33∗T22∗T1T005∗4∗T34∗3∗T23∗2∗T12∗1∗T000][p5p4p3p2p1p0]⇒dm=AmPm⇒Pm=Am−1dm\left[ \begin{matrix} P(0) \\ P^{(1)}(0) \\ P^{(2)}(0) \\ P(T) \\ P^{(1)}(T) \\ P^{(2)}(T) \end{matrix} \right] =\left[ \begin{matrix} 0^5 & 0^4 & 0^3 & 0^2 & 0^1 & 0^0 \\ 5*0^4 & 4*0^3 & 3*0^2 & 2*0^1 & 0^0 & 0 \\ 5*4*0^3 & 4*3*0^2 & 3*2*0^1 & 2*1*0^0 & 0 & 0 \\ T^5 & T^4 & T^3 & T^2 & T^1 & T^0 \\ 5*T^4 & 4*T^3 & 3*T^2 & 2*T^1 & T^0 & 0 \\ 5*4*T^3 & 4*3*T^2 & 3*2*T^1 & 2*1*T^0 & 0 & 0 \end{matrix} \right] \left[ \begin{matrix} p_5\\ p_4\\ p_3\\ p_2 \\ p_1 \\ p_0 \end{matrix} \right] \\ \Rightarrow d_m = A_m P_m \\ \Rightarrow P_m = A_m^{-1}d_m ⎣⎢⎢⎢⎢⎢⎢⎡P(0)P(1)(0)P(2)(0)P(T)P(1)(T)P(2)(T)⎦⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎡055∗045∗4∗03T55∗T45∗4∗T3044∗034∗3∗02T44∗T34∗3∗T2033∗023∗2∗01T33∗T23∗2∗T1022∗012∗1∗00T22∗T12∗1∗T001000T1T000000T000⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎡p5p4p3p2p1p0⎦⎥⎥⎥⎥⎥⎥⎤⇒dm=AmPm⇒Pm=Am−1dm
因此找到了映射矩阵，类似的可以推导任意阶多项式的映射矩阵：
Am=[0N0N−1⋯020100N∗0N−1(N−1)∗0N−2⋯2∗01000N(N−1)∗0N−2(N−1)(N−2)∗0N−3⋯2∗1∗0000⋮⋮⋮⋮⋮⋮TNTN−1⋯T2T1T0N∗TN−1(N−1)∗TN−2⋯2∗T1T00N(N−1)∗TN−2(N−1)(N−2)∗TN−3⋯2∗1∗T000⋮⋮⋮⋮⋮⋮]A_m = \left[ \begin{matrix} 0^N & 0^{N-1} & \cdots & 0^2 & 0^1 & 0^0 \\ N*0^{N-1} & (N-1)*0^{N-2} & \cdots & 2*0^1 & 0^0 & 0 \\ N(N-1)*0^{N-2} & (N-1)(N-2)*0^{N-3} & \cdots & 2*1*0^0 & 0 & 0 \\ \vdots&\vdots& \vdots & \vdots & \vdots & \vdots \\ T^N & T^{N-1} & \cdots & T^2 & T^1 & T^0 \\ N*T^{N-1} & (N-1)*T^{N-2} & \cdots & 2*T^1 & T^0 & 0 \\ N(N-1)*T^{N-2} & (N-1)(N-2)*T^{N-3} & \cdots & 2*1*T^0 & 0 & 0 \\ \vdots&\vdots& \vdots & \vdots & \vdots & \vdots \\ \end{matrix} \right] Am=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡0NN∗0N−1N(N−1)∗0N−2⋮TNN∗TN−1N(N−1)∗TN−2⋮0N−1(N−1)∗0N−2(N−1)(N−2)∗0N−3⋮TN−1(N−1)∗TN−2(N−1)(N−2)∗TN−3⋮⋯⋯⋯⋮⋯⋯⋯⋮022∗012∗1∗00⋮T22∗T12∗1∗T0⋮01000⋮T1T00⋮0000⋮T000⋮⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

跳转到生成A矩阵的代码

将代价函数中的系数向量通过映射矩阵AmA_mAm进行替换，则得到下式：
min⁡[P0P1⋮PM]T[Q00⋯⋯00Q10⋯0⋮⋮⋮⋱⋮000⋯QM][P0P1⋮PM]⇒min⁡[d0d1⋮dM]T[A00⋯⋯00A10⋯0⋮⋮⋮⋱⋮000⋯AM]−T[Q00⋯⋯00Q10⋯0⋮⋮⋮⋱⋮000⋯QM][A00⋯⋯00A10⋯0⋮⋮⋮⋱⋮000⋯AM]−1[P0P1⋮PM]⇒min⁡dTA−TQA−1d\min \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right]^T \left[ \begin{matrix} Q_0 & 0 & \cdots & \cdots & 0 \\ 0 & Q_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 0 & 0 & 0 & \cdots & Q_M \end{matrix} \right] \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right] \\ \Rightarrow \min \left[ \begin{matrix} d_0 \\ d_1 \\ \vdots \\ d_M \end{matrix} \right]^T \left[ \begin{matrix} A_0 & 0 & \cdots & \cdots & 0 \\ 0 & A_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 0 & 0 & 0 & \cdots & A_M \end{matrix} \right]^{-T} \left[ \begin{matrix} Q_0 & 0 & \cdots & \cdots & 0 \\ 0 & Q_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 0 & 0 & 0 & \cdots & Q_M \end{matrix} \right] \left[ \begin{matrix} A_0 & 0 & \cdots & \cdots & 0 \\ 0 & A_1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 0 & 0 & 0 & \cdots & A_M \end{matrix} \right]^{-1} \left[ \begin{matrix} P_0 \\ P_1 \\ \vdots \\ P_M \end{matrix} \right] \\ \Rightarrow \min d^TA^{-T}QA^{-1}d min⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤T⎣⎢⎢⎢⎡Q00⋮00Q1⋮0⋯0⋮0⋯⋯⋱⋯00⋮QM⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤⇒min⎣⎢⎢⎢⎡d0d1⋮dM⎦⎥⎥⎥⎤T⎣⎢⎢⎢⎡A00⋮00A1⋮0⋯0⋮0⋯⋯⋱⋯00⋮AM⎦⎥⎥⎥⎤−T⎣⎢⎢⎢⎡Q00⋮00Q1⋮0⋯0⋮0⋯⋯⋱⋯00⋮QM⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡A00⋮00A1⋮0⋯0⋮0⋯⋯⋱⋯00⋮AM⎦⎥⎥⎥⎤−1⎣⎢⎢⎢⎡P0P1⋮PM⎦⎥⎥⎥⎤⇒mindTA−TQA−1d
上式实际上替换，实际上已经把（10）式中微分约束给带入到了代价函数中，但是连续性约束还没有引入到代价函数中。

在论文【1】中，作者引入一个permutation matrix (置换矩阵)CCC来将连续性约束进行引入到代价函数中。CCC矩阵本身是一个只包含0和1的矩阵，它的作用就是给所有[d0d1⋯dM]\left[\begin{matrix}d_0 & d_1 & \cdots &d_M\end{matrix}\right][d0d1⋯dM]进行一个组合，将固定的变量放在头部，将需要就行优化的变量方位尾部，并且对于连续性约束的变量，由于他们的值是相等的，因此选择其中一个来表达连续性约束。说起来有一些抽象，用数学表达就会清晰很多。

根据之前的介绍，那么紫色框中变量的值都是我们事先已经设置好的，浅蓝色框中的变量就是代价函数在优化时分配的最优值，也就是我们需要优化的变量。

咱们以最小化Jerk为例，假设有如下一条轨迹，它有两个分段轨迹，如上图，那么它的CCC矩阵如下：
[d0,0(0)d0,0(1)d0,0(2)dT,0(0)dT,0(1)dT,0(2)d0,1(0)d0,1(1)d0,1(2)dT,1(0)dT,1(1)dT,1(2)]=[100000000010000000001000000000100000000000010000000001000100000000000010000000001000010000000001000000000100][d0,0(0)d0,0(1)d0,0(2)dT,0(0)dT,1(0)dT,1(1)dT,1(2)dT,0(1)dT,0(2)]⇒dm=CT[dmFdmP]\left[ \begin{matrix} d_{0,0}^{(0)} \\ d_{0,0}^{(1)}\\d_{0,0}^{(2)}\\d_{T,0}^{(0)}\\d_{T,0}^{(1)}\\d_{T,0}^{(2)}\\d_{0,1}^{(0)}\\d_{0,1}^{(1)}\\d_{0,1}^{(2)}\\d_{T,1}^{(0)}\\d_{T,1}^{(1)}\\d_{T,1}^{(2)} \end{matrix} \right] = \left[ \begin{matrix} 1 & 0&0&0&0&0&0&0&0 \\ 0 & 1&0&0&0&0&0&0&0 \\ 0 & 0&1&0&0&0&0&0&0 \\ 0 & 0&0&1&0&0&0&0&0 \\ 0 & 0&0&0&0&0&0&1&0 \\ 0 & 0&0&0&0&0&0&0&1 \\ 0 & 0&0&1&0&0&0&0&0 \\ 0 & 0&0&0&0&0&0&1&0 \\ 0 & 0&0&0&0&0&0&0&1 \\ 0 & 0&0&0&1&0&0&0&0 \\ 0 & 0&0&0&0&1&0&0&0 \\ 0 & 0&0&0&0&0&1&0&0 \\ \end{matrix} \right] \left[ \begin{matrix} d_{0,0}^{(0)} \\ d_{0,0}^{(1)}\\d_{0,0}^{(2)}\\d_{T,0}^{(0)}\\d_{T,1}^{(0)}\\d_{T,1}^{(1)}\\d_{T,1}^{(2)}\\d_{T,0}^{(1)} \\d_{T,0}^{(2)} \end{matrix} \right] \\ \Rightarrow d_m = C^T \left[ \begin{matrix} d_{m_F} \\d_{m_P} \end{matrix} \right] ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡d0,0(0)d0,0(1)d0,0(2)dT,0(0)dT,0(1)dT,0(2)d0,1(0)d0,1(1)d0,1(2)dT,1(0)dT,1(1)dT,1(2)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡100000000000010000000000001000000000000100100000000000000100000000000010000000000001000010010000000001001000⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡d0,0(0)d0,0(1)d0,0(2)dT,0(0)dT,1(0)dT,1(1)dT,1(2)dT,0(1)dT,0(2)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⇒dm=CT[dmFdmP]

其中dmFd_{mF}dmF表示第m段轨迹中的已经确定的变量(对应紫色框中的变量)，dmPd_{mP}dmP表示可以自由变化的变量(浅蓝色框中的变量)。

类似的你可以写出更高阶情况下的置换矩阵CCC.

不知道你是否发现了这其中奥妙，由于使用了置换矩阵CCC，我们可以将连续性约束给引入到代价函数中，因为用了CCC矩阵之后，就把相等的变量给省略了，它们的值相等，所以只需要保留一个进行优化即可。

跳转到生成置换矩阵C的代码

将上式带入代价函数，则最终得到：
J=min⁡[dFdP]TCA−TQA−1CT[dFdP]J = \min \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right]^T CA^{-T}QA^{-1}C^T \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right] J=min[dFdP]TCA−TQA−1CT[dFdP]
最终我们将带约束的QP优化问题，转换成了一个无约束的QP问题，对于上式因为它是一个有全局最小值的凸函数，所以对其求导取极值点就是它的全局最优值：

令
CA−TQA−1CT=RCA^{-T}QA^{-1}C^T = R CA−TQA−1CT=R
并对RRR矩阵根据dFd_FdF和dPd_PdP的尺寸进行分块，则得到如下变换：
J=min⁡[dFdP]TCA−TQA−1CT[dFdP]=J=min⁡[dFdP]TR[dFdP]=min⁡[dFdP]T[RFFRFPRPFRPP][dFdP]=dFTRFFdF+dFTRFPdP+dPTRPFdF+dPTRPPdPJ = \min \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right]^T CA^{-T}QA^{-1}C^T \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right] = J = \min \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right]^T R \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right] \\= \min \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right]^T \left[ \begin{matrix} R_{FF} & R_{FP} \\ R_{PF} & R_{PP} \\ \end{matrix} \right] \left[ \begin{matrix} d_F \\ d_P \end{matrix} \right] \\=d_F^TR_{FF}d_F+d_F^TR_{FP}d_P+d_P^TR_{PF}d_F+d_P^TR_{PP}d_P J=min[dFdP]TCA−TQA−1CT[dFdP]=J=min[dFdP]TR[dFdP]=min[dFdP]T[RFFRPFRFPRPP][dFdP]=dFTRFFdF+dFTRFPdP+dPTRPFdF+dPTRPPdP
对于上式JJJ是一个标量，并且QQQ矩阵是对称矩阵，而CA−TQA−1CT=(CA−TQA−1CT)TCA^{-T}QA^{-1}C^T = (CA^{-T}QA^{-1}C^T)^TCA−TQA−1CT=(CA−TQA−1CT)T，所以RRR也是对称矩阵，则RPF=RFPTR_{PF} = R_{FP}^TRPF=RFPT

上式中：
J=dFTRFFdF+dFTRFPdP+dPTRPFdF+dPTRPPdP=dFTRFFdF+2dFTRFPdP+dPTRPPdPJ=d_F^TR_{FF}d_F+d_F^TR_{FP}d_P+d_P^TR_{PF}d_F+d_P^TR_{PP}d_P \\ = d_F^TR_{FF}d_F+2d_F^TR_{FP}d_P+d_P^TR_{PP}d_P J=dFTRFFdF+dFTRFPdP+dPTRPFdF+dPTRPPdP=dFTRFFdF+2dFTRFPdP+dPTRPPdP
对于上式我们要求的变量是dPd_PdP，因此将其对dPd_PdP求导，并零导数为0时取极值点的值对应的就是最优的dP∗d_P^*dP∗：
∂J∂dP=2dFTRFP+2RPPdP=0⇒dP∗=−RPP−1RFPTdF\frac{\partial J}{\partial d_P} = 2d_F^TR_{FP}+2R_{PP}d_P = 0 \\ \Rightarrow d_P^* = -R_{PP}^{-1}R_{FP}^Td_F ∂dP∂J=2dFTRFP+2RPPdP=0⇒dP∗=−RPP−1RFPTdF
求解出来了dP∗d_P^*dP∗之后，问题就很简单了，我们只需要将它带入CCC和AAA就能求出最终的各个多项式的系数：

P=A−1CT[dFdP∗]P = A^{-1}C^T \left[ \begin{matrix} d_F \\ d_P^* \end{matrix} \right] P=A−1CT[dFdP∗]

6. 从多项式中解算出轨迹

有了多项式的系数之后，我们只要通过离散化时间ttt，就可以得到每一段轨迹。原理比较简单，在此不做深究，直接看代码：跳转到结算轨迹的代码

7.每段多项式的时间分配问题

还记得上面我在推导Minimum Snap时说假设事先已经给定了每段轨迹所用的时间，实际上这并不是一个非常好的做法，因为实际上时间对于最终生成轨迹的好坏有很大影响，下图是来自论文中的结果：

可以看出来，不同的时间轨迹相差还是比较明显的。

在【1】和【2】的论文中都提出了对每段时间也加入到优化的代价函数中，这样固然是好，但是势必会增加优化的时间，并且也可能会导致无法再闭式的求解这个问题。

实际上在工程中，有一个简化的方法去分配每一段轨迹的时间，也就是先假设机器人以最大的加速度加速到最大速度，然后运行一段时间，再以最大的减速度到终点，到达是速度为0。经过我的多次实验，这样的方法分配出的时间基本都能优化出一条非常好的轨迹。在我的代码中就采用了这种思路！跳转到时间分配的代码

8. 如何解决生成的轨迹与障碍物碰撞

这个实际已经超过了本篇博客要介绍的内容，但是为了完整性，我在这里简单地介绍一下在论文【2】中提供的思路：如果发生碰撞，就在两个航迹点的中间再选一个航迹点，然后重新优化生成轨迹，重复这样的操作，直到生成安全的轨迹。如下图：

9. 实验效果

在该项目的仓库中，我提供了两个关于Minimum Snap的应用场景，一个是直接通过rviz中的插件点击空间位置，然后自动生成轨迹。另一个则是编写了一个完整的从前端轨迹搜索(Astar)到后端轨迹优化的全过程的代码。

先来看第一个(GIF图片过大，需要跳转到我的github中观看)：

实验效果GIF
对应的运行文件

第二个：
该实验，我实现了一个自动随机生成栅格地图，然后使用A star进行路径搜索，搜索到路径之后，使用Minimum Snap进行轨迹拟合。

对应的运行文件

为了方便，我省略了从A star搜索的路径中选择关键位置的步骤，所以A star搜索到的所有节点，都会被输入Minimum Snap进行拟合！

跳转到代码运行方法的文档

另外，我实现的Minimum Snap算法其计算用时比论文中低，当然作者的电脑是接近10年之前的，哈哈，不过实现的效果以及效率基本和论文中是一致的！

在写这篇博客时浙江大学湖州实验室的高飞老师给了我很多帮助，同时也借鉴了他的一些代码和课件，在此对他进行感谢！

【1】: Mellinger D, Kumar V. Minimum snap trajectory generation and control for quadrotors[C]//Robotics and Automation (ICRA), 2011 IEEE International Conference on. IEEE, 2011: 2520-2525.
【2】：Polynomial Trajectory Planning for Aggressive Quadrotor Flight in Dense Indoor Environments, Charles Richter, Adam Bry, and Nicholas Roy