【轨迹规划论文整理(1)】UAV轨迹规划的开山之作Minimum Snap Trajectory_minimum snap trajectory generation and control for

【轨迹规划论文整理(1)】UAV轨迹规划的开山之作Minimum Snap Trajectory Generation and Control for Quadrotors

本系列主要是对精读的一些关于无人机、无人车的轨迹搜索论文的整理，包括了论文所拓展的其他一些算法的改进思路。
这是本系列的第一篇文章，也是UAV轨迹规划的开山之作，是所有学习无人机方向的需要精读的第一篇文章，两个作者来自于宾夕法尼亚大学的GRASP Lab，是全球顶尖的机器人团队，也是UAV方向研究的开山鼻祖之一。

Mellinger D , Kumar V .Minimum snap trajectory generation and control for quadrotors[J].IEEE, 2011.DOI:10.1109/ICRA.2011.5980409.

1. Introduction

1.1 轨迹规划的定义

在机器人导航过程中，如何控制机器人从A点移动到B点，通常称之为运动规划。运动规划一般分为两步：

• 路径搜索：在地图(包括给栅格、八叉树、点云地图等等) 搜索一条从A到B点的路径，这条路径通常是由一系列离散的空间点(waypoint)组成。该部分称之为前端。
• 轨迹规划：由于路径点可能比较稀疏、而且不平滑，为了更好地控制机器人运动，需要将稀疏的路径点变成平滑的曲线或稠密的轨迹点，即轨迹。这部分称之为后端。

本文的主要工作就是一种无人机轨迹生成的方法，额外还包括了建模和控制器的设计。

2. Model

图中是本文定义的坐标关系，利用Z-X-Y欧拉角定义横滚角、俯仰角、偏航角(roll,pitch,yaw)作为本地坐标系统。从B至W的旋转矩阵通过$R_B^W=R_C^W{R}_B^C$来表示，$R_C^W$代表偏航角旋转至中间坐标系(Z轴方向与世界坐标Z轴方向一致)。UAV机体在世界坐标系的角速度计算如下，p、q、r分量在记载坐标系中的分量为：

$${\omega }_B^W=p{X}_B+q{Y}_B+r{Z}_B$$

每个电机有角速度${\omega }_i$，产生一个力$F_i$和$M_i$，关系为:

$$\begin{gathered} F_i=k_F{\omega }_i^2 \\ {M}_i=k_M{\omega }_i^2 \end{gathered}$$

控制输入记为u，$u_1$为净升力，$u_2、u_3、u_4$是机体力矩可以通过电机转速表示为(L是电机旋转中心到电机中心的距离)

u = [ k F k F k F k F 0 k F L 0 − k F L − k F L 0 k F L 0 k M − k M k M − k M ] [ ω 1 2 ω 2 2 ω 3 2 ω 4 2 ] \mathbf{u}=

$$\begin{bmatrix}k_F&k_F&k_F&k_F\\0&k_FL&0&-k_FL\\-k_FL&0&k_FL&0\\k_M&-k_M&k_M&-k_M\end{bmatrix}$$ $$\begin{bmatrix}\omega_1^2\\\omega_2^2\\\omega_3^2\\\omega_4^2\end{bmatrix}$$ u= ​kF​0−kF​LkM​​kF​kF​L0−kM​​kF​0kF​LkM​​kF​−kF​L0−kM​​ ​ ​ω12​ω22​ω32​ω42​​ ​

机体的线加速度，其中系统受力为$-z_W$方向上的重力和四个电机的合力为$z_B$方向上的$u_1$。

m r ¨ = − m g z W + u 1 z B

$$\begin{aligned}m\ddot{\mathbf{r}}=-mg\mathbf{z}_W+u_1\mathbf{z}_B\end{aligned}$$ mr¨=−mgzW​+u1​zB​​

用欧拉公式计算角速度如下，其中I是惯性矩矩阵：

ω ˙ B W = I − 1 [ − ω B W × I ω B W + [ u 2 u 3 u 4 ] ] \left.\dot{\omega}_{\mathcal{BW}}=\mathcal{I}^{-1}\left[-\omega_{\mathcal{BW}}\times\mathcal{I}\omega_{\mathcal{BW}}+\left[

$$\begin{array}{c}u_2\\u_3\\u_4\end{array}$$ \right.\right]\right] ω˙BW​=I−1 ​−ωBW​×IωBW​+ ​u2​u3​u4​​ ​ ​

UAV系统的状态向量由质心的位置和速度、方向和角速度给出：

$$X=[x,y,z,\varphi ,\theta ,\psi ,\dot{x},\dot{y},\dot{z},p,q,r]$$

3. 微分平坦性

UAV的平坦输出可以写作：

$$\sigma =[x,y,z,\psi {]}^T$$

其中$r=[x,y,z{]}^T$是世界坐标系中质心的坐标，而ψ是偏航角，这里定义轨迹$\sigma (t)$为平坦输出空间中的平滑曲线：

$$\sigma (t):[t_0,t_m]\to {\mathbb{R}}^3\times SO(2)$$

根据微分平坦性，系统的状态和控制输入可以用$\sigma$及其导数表示，具体推导参考博客。

4. 控制

本文设计了一个控制器去跟踪特定轨迹${\sigma }_T(t)=[r_T(t{)}^T,{\psi }_T(t){]}^T$。首先定义位置和速度误差为

$$e_p=r-r_T,e_v=\dot{r}-{\dot{r}}_T$$

计算控制器所需要的力矢量：

$$F_{des}=-K_p{e}_p-K_v{e}_v+mg{z}_W+m{\ddot{r}}_T$$

将所需的力矢量投影到实际的身体坐标系z轴上，从而计算第一输入的所需力：

$$u_1=F_{des}\cdot z_B$$

理想的$Z_{B,des}$一定沿着理想推力方向：

$$Z_{B,des}=\frac{F_{des}}{||F_{des}||}$$

此时，理想的旋转$R_B^W$为：

$$R_{des}{e}_3=Z_{B,des}$$

其他理想状态量：

$$\begin{gathered} {\mathbf{x}}_{C,des}={\left[\cos {\psi }_T,\sin {\psi }_T,0\right]}^T, \\ \text{and} \\ {\mathbf{y}}_{B,des}=\frac{{\mathbf{z}}_{B,des}\times {\mathbf{x}}_{C,des}}{\Vert {\mathbf{z}}_{B,des}\times {\mathbf{x}}_{C,des}\Vert },{\mathbf{x}}_{B,des}={\mathbf{y}}_{B,des}\times {\mathbf{z}}_{B,des}, \end{gathered}$$

$Z_B,X_B,Y_B$构成了旋转矩阵，从而获得了理想的旋转矩阵

定义旋转误差：

$${\mathbf{e}}_R=\frac{1}{2}(R_{des}^T{R}_B^W-R_B^W{}^T{R}_{des}{)}^{\vee }$$

角速度误差（机身系）：

$$e_{\omega }=\omega -{\omega }_T$$

另外三个输入的计算如下，两个K都是对角增益矩阵：

$$[u_2,u_3,u_4{]}^T=-K_R{e}_R-K_{\omega }{e}_{\omega }$$

5. 轨迹生成(本文最重要部分)

生成的轨迹可以写成m个时间间隔上的n阶分段多项式函数：

σ T ( t ) = { ∑ i = 0 n σ T i 1 t i t 0 ≤ t < t 1 ∑ i = 0 n σ T i 2 t i t 1 ≤ t < t 2 ⋮ ∑ i = 0 n σ T i m t i t m − 1 ≤ t ≤ t m \sigma_T(t)=\left\{

$$\begin{array}{cc}\sum_{i=0}^n\sigma_{Ti1}t^i&t_0\leq t<t_1\\\sum_{i=0}^n\sigma_{Ti2}t^i&t_1\leq t<t_2\\\vdots\\\sum_{i=0}^n\sigma_{Tim}t^i&t_{m-1}\leq t\leq t_m\end{array}$$ \right. σT​(t)=⎩ ⎨ ⎧​∑i=0n​σTi1​ti∑i=0n​σTi2​ti⋮∑i=0n​σTim​ti​t0​≤t<t1​t1​≤t<t2​tm−1​≤t≤tm​​

根据上一节的微分平坦性，这里的${\sigma }_T=[x_T,y_T,z_T,{\psi }_T{]}^T$，${\sigma }_i=[x_i,y_i,z_i,{\psi }_i{]}^T$。上述中要求相邻两段轨迹相同的时间节点是$r_T$微分连续。

建立二次规划问题，c是一个4nm*1的向量，包含了${\sigma }_{T_{ij}}=[x_{T_{ij}},y_{T_{ij}},z_{T_{ij}},{\psi }_{T_{ij}}{]}^T$及它们的导数：

$$\begin{gathered} \min c^{T}Hc+f^{T}c \\ s.t.Ac\le b \end{gathered}$$

接下来的解析我觉得文章本身不是很清晰，但是结合深蓝学院的课程和其他博客的思路，重新解读一下。

5.1 一段多项式轨迹

N+1个已知参数点可以拟合N次函数，对于首尾位置、速度、加速度已知的条件(相当于6个已知参数)，可以唯一确定一个五次多项式。

5.2 多段多项式轨迹

本文提出了一个经过多个中间航迹点的轨迹，分别使用多个轨迹段进行拟合。但是相对于一段轨迹我们指定了节点的PVA，对于多段轨迹的中间点我们很难给定UAV经过这些航点的VAJ等信息，此时我们将中间速度、加速度等信息看作一个可以自由变化的状态，对于中间状态不确定的问题，可以构建一个优化问题，使得中间状态取得某一个值时，制定的代价函数达到最小值。

Minimum Snap的核心思路就是多端轨迹的中间航点的变量如何取值可以使机器人能量耗费最少。

5.3 Minimum Snap

UAV的轨迹再各个坐标轴上是独立的，因此可以对其分别进行轨迹拟合。即xyz轴上路径生成可以直接将三个轴合成就得到一个完整的空间轨迹。下面的推导均在一维上进行。

第m段轨迹表示如下：

P m ( t ) = p m , 0 t 0 + p m , 1 t 1 + ⋯ + p m , N t N = [ p m , 0 p m , 1 ⋯ p m , N ] [ t 0 t 1 ⋮ t N ] \mathrm{P_m\left(t\right)=p_{m,0}t^0+p_{m,1}t^1+\cdots+p_{m,N}t^N}\\=

$$\begin{bmatrix}\mathrm{p_{m,0}}&\mathrm{p_{m,1}}&\cdots&\mathrm{p_{m,N}}\end{bmatrix}$$ $$\begin{bmatrix}\mathrm{t^0}\\\mathrm{t^1}\\\vdots\\\mathrm{t^N}\end{bmatrix}$$ Pm​(t)=pm,0​t0+pm,1​t1+⋯+pm,N​tN=[pm,0​​pm,1​​⋯​pm,N​​] ​t0t1⋮tN​ ​

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

对其分别求一阶、二阶、三阶导，获得它的速度、加速度、Jerk函数函数如下：

d P ( t ) d t = [ p m , 0 p m , 1 ⋯ p m , N ] [ 0 1 ⋮ N t N − 1 ]   d 2 P ( t ) d t 2 = [ p m , 0 p m , 1 ⋯ p m , N ] [ 0 1 ∗ 0 2 ∗ 1 ⋮ N ( N − 1 ) t N − 2 ]

$$\begin{aligned}\frac{\mathrm{dP}\left(\mathrm{t}\right)}{\mathrm{dt}}&=\begin{bmatrix}\mathrm{p}_{\mathrm{m},0}&\mathrm{p}_{\mathrm{m},1}&\cdots&\mathrm{p}_{\mathrm{m},\mathrm{N}}\end{bmatrix}\begin{bmatrix}0\\1\\\vdots\\\mathrm{N} \mathrm{t}^{\mathrm{N}-1}\end{bmatrix}\\\\\ \frac{\mathrm{d}^2\mathrm{P}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}^2}&=\begin{bmatrix}\mathrm{p}_{\mathrm{m},0}&\mathrm{p}_{\mathrm{m},1}&\cdots&\mathrm{p}_{\mathrm{m},\mathrm{N}}\end{bmatrix}\begin{bmatrix}0\\1*0\\2*1\\\vdots\\\mathrm{N}\left(\mathrm{N}-1\right)\mathrm{t}^{\mathrm{N}-2}\end{bmatrix}\end{aligned}$$ dtdP(t)​ dt2d2P(t)​​=[pm,0​​pm,1​​⋯​pm,N​​] ​01⋮NtN−1​ ​=[pm,0​​pm,1​​⋯​pm,N​​] ​01∗02∗1⋮N(N−1)tN−2​ ​​

d 3 P ( t ) d t 3 = [ p m , 0 p m , 1 ⋯ p m , N ] [ 0 ∗ 0 ∗ 0 1 ∗ 0 ∗ 0 2 ∗ 1 ∗ 0 ⋮ N ( N − 1 ) ( N − 2 ) t N − 3 ] d 4 P ( t ) d t 4   =   [ p m , 0 p m , 1 ⋯ p m , N ]   [ 0 ∗ 0 ∗ 0 ∗ 0 1 ∗ 0 ∗ 0 ∗ 0 2 ∗ 1 ∗ 0 ∗ 0 3 ∗ 2 ∗ 1 ∗ 0 ⋮ N ( N − 1 ) ( N − 2 ) ( N − 3 ) t N − 4 ] \frac{\mathrm{d}^3\mathrm{P}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}^3}=

$$\begin{bmatrix}\mathrm{p}_{\mathrm{m},0}&\mathrm{p}_{\mathrm{m},1}&\cdots&\mathrm{p}_{\mathrm{m},\mathrm{N}}\end{bmatrix}$$ $$\begin{bmatrix}0*0*0\\1*0*0\\2*1*0\\\vdots\\\mathrm{N}(\mathrm{N}-1)(\mathrm{N}-2)\mathrm{t}^{\mathrm{N}-3}\end{bmatrix}$$ \\\mathrm{\frac{d^4P(t)}{dt^4}}\:=\: $$\begin{bmatrix}\mathrm{p_{m,0}}&\mathrm{p_{m,1}}&\cdots&\mathrm{p_{m,N}}\end{bmatrix}$$ \: $$\begin{bmatrix}0*0*0*0\\1*0*0*0\\2*1*0*0\\3*2*1*0\\\vdots\\\mathrm{N(N-1)(N-2)(N-3)t^{N-4}}\end{bmatrix}$$ dt3d3P(t)​=[pm,0​​pm,1​​⋯​pm,N​​] ​0∗0∗01∗0∗02∗1∗0⋮N(N−1)(N−2)tN−3​ ​dt4d4P(t)​=[pm,0​​pm,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​ ​

在论文中提出了一个代价函数$J_m$，对应的是一段多项式轨迹的四阶导数(snap)平方的积分，具体公式如下：

J m = ∫ T m − 1 T m ( d 4 P ( t ) d t 4 ) 2 = ∫ T m − 1 T m [ p N p N − 1 ⋮ p 0 ] T [ N ( N − 1 ) ( N − 2 ) ( N − 3 ) t N − 4 ⋮ 3 ∗ 2 ∗ 1 ∗ 0 2 ∗ 1 ∗ 0 ∗ 0 1 ∗ 0 ∗ 0 ∗ 0 0 ∗ 0 ∗ 0 ∗ 0 ] [ N ( N − 1 ) ( N − 2 ) ( N − 3 ) t N − 4 ⋮ 3 ∗ 2 ∗ 1 ∗ 0 2 ∗ 1 ∗ 0 ∗ 0 1 ∗ 0 ∗ 0 ∗ 0 0 ∗ 0 ∗ 0 ∗ 0 ] T [ p N ⋮ p N − 1 ⋮ p 0 ] = P m T Q m P m \mathrm{J_m=\int_{T_{m-1}}^{T_m}(\frac{d^4P\left(t\right)}{dt^4})^2}\\=\int_{\mathrm{T_{m-1}}}^{\mathrm{T_m}}

$$\begin{bmatrix}\mathrm{p_N}\\\mathrm{p_{N-1}}\\\vdots\\\mathrm{p_0}\end{bmatrix}$$ ^{\mathrm{T}} $$\begin{bmatrix}\mathrm{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{bmatrix}$$ $$\begin{bmatrix}\mathrm{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{bmatrix}$$ ^{\mathrm{T}} $$\begin{bmatrix}\mathrm{p_N}\\\vdots\\\mathrm{p_{N-1}}\\\vdots\\\mathrm{p_0}\end{bmatrix}$$ \\=P_m^TQ_mP_m Jm​=∫Tm−1​Tm​​(dt4d4P(t)​)2=∫Tm−1​Tm​​ ​pN​pN−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 ​pN​⋮pN−1​⋮p0​​ ​=PmT​Qm​Pm​

其中$Q_m$如下：

[ 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 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 ⋯ ⋮ ⋮ ⋮ ⋯ i ( i − 1 ) ( i − 2 ) ( i − 3 ) ( i − 1 ) ( i − 2 ) ( T m − T m − 1 ) i + 1 − 7 i + 1 − 7 ⋯ ⋮ ⋮ ⋮ 0 0 ⋯ ]

$$\begin{bmatrix}\frac{\mathrm{N}\left(\mathrm{N}-1\right)\left(\mathrm{N}-2\right)\left(\mathrm{N}-3\right)\mathrm{N}\left(\mathrm{N}-1\right)\left(\mathrm{N}-2\right)\left(\mathrm{N}-3\right)\left(\mathrm{T}_{ {\mathrm{m}}}-\mathrm{T}_{ {\mathrm{m}-1}}\right)^{ {\mathrm{N}+\mathrm{N}-7}}}{\mathrm{N}+\mathrm{N}-7}&\frac{\mathrm{N}\left(\mathrm{N}-1\right)\left(\mathrm{N}-2\right)\left(\mathrm{N}-3\right)\left(\mathrm{N}-1\right)\left(\mathrm{N}-2\right)\left(\mathrm{N}-3\right)\left(\mathrm{N}-4\right)\left(\mathrm{T}_{ {\mathrm{m}}}-\mathrm{T}_{ {\mathrm{m}-1}}\right)^{ {\mathrm{N}+\mathrm{N}-1-7}}}{\mathrm{N}+\mathrm{N}-1-7}&\cdots\\\varvdots&\varvdots&\varvdots\\\cdots&\frac{\mathrm{i}\left(\mathrm{i}-1\right)\left(\mathrm{i}-2\right)\left(\mathrm{i}-3\right)\left(\mathrm{i}-1\right)\left(\mathrm{i}-2\right)\left(\mathrm{T}_{ {\mathrm{m}}}-\mathrm{T}_{ {\mathrm{m}-1}}\right)^{ {\mathrm{i}+1-7}}}{\mathrm{i}+1-7}&\cdots\\\varvdots&\varvdots&\varvdots\\0&0&\cdots\end{bmatrix}$$ ​N+N−7N(N−1)(N−2)(N−3)N(N−1)(N−2)(N−3)(Tm​−Tm−1​)N+N−7​⋮⋯⋮0​N+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+1−7i(i−1)(i−2)(i−3)(i−1)(i−2)(Tm​−Tm−1​)i+1−7​⋮0​⋯⋮⋯⋮⋯​ ​

上式中$T_m-T_{m-1}$实际就是该段轨迹运动所用的时长，为了不让每一段轨迹的时间差过大，这里采用相对时间来表示，对于每一段轨迹认为它们的时间都是$(0,T_m)$。

总的代价函数

假设现在有M段轨迹，那么每段轨迹都是一个高阶的多项式函数，写出代价函数如下：

J = J 0 + J 1 + ⋯ + J M = ∑ m = 0 M P m ⃗ T Q m P m ⃗ = [ P 0 P 1 ⋮ P M ] T [ Q 0 0 ⋯ ⋯ 0 0 Q 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋯ ⋮ 0 0 0 ⋯ Q M ] [ P 0 P 1 ⋮ P M ]

$$\begin{gathered}\mathrm{J}=\mathrm{J}_{0}+\mathrm{J}_{1}+\cdots+\mathrm{J}_{\mathrm{M}} \\=\sum_{\mathrm{m=0}}^\mathrm{M}\vec{\mathrm{P_m}}^\mathrm{T}\mathrm{Q_m}\vec{\mathrm{P_m}} \\\left.=\begin{bmatrix}\mathrm{P}_0\\\mathrm{P}_1\\\vdots\\\mathrm{P}_\mathrm{M}\end{bmatrix}\right.^\mathrm{T}\begin{bmatrix}\mathrm{Q}_0&0&\cdots&\cdots&0\\0&\mathrm{Q}_1&0&\cdots&0\\\vdots&\vdots&\vdots&\cdots&\vdots\\0&0&0&\cdots&\mathrm{Q}_\mathrm{M}\end{bmatrix}\begin{bmatrix}\mathrm{P}_0\\\mathrm{P}_1\\\vdots\\\mathrm{P}_\mathrm{M}\end{bmatrix} \end{gathered}$$ J=J0​+J1​+⋯+JM​=m=0∑M​Pm​ ​TQm​Pm​ ​= ​P0​P1​⋮PM​​ ​T ​Q0​0⋮0​0Q1​⋮0​⋯0⋮0​⋯⋯⋯⋯​00⋮QM​​ ​ ​P0​P1​⋮PM​​ ​​

代价函数完成了，但是还缺少约束，即每个中间航点在前后两段轨迹的速度、加速度或者jerk是连续的，如此运动轨迹才是平滑的。

5.4 约束

连续性约束

限制中间的航迹点在两段轨迹的接合处是左右两边状态连续的，写成数学形式如下：

$$P_m^{(k)}(T_m)=P_{m+1}^{(k)}(T_m)$$

即轨迹段$P_m$和$P_{m+1}$在$T_m$时刻具有相同的k阶导数。

P m ( k ) ( T m ) = P m + 1 ( k ) ( T m ) ⇒ ∑ i ≥ k i ! ( i − k ) ! T m i − k p m , i − ∑ l ≥ k l ! ( l − k ) ! T m l − k p m + 1 , l = 0 ⇒ [ A m − A m + 1 ] [ P m P m + 1 ] = 0

$$\begin{aligned}\mathrm{P_{m}^{(k)}(T_{m})=P_{m+1}^{(k)}(T_{m})}\\\Rightarrow\sum_{ {\mathrm{i\geq k}}}\frac{\mathrm{i!}}{(\mathrm{i-k})!}\mathrm{T_{m}^{i-k}p_{m,i}-\sum_{ {\mathrm{l\geq k}}}\frac{l!}{(l-k)!}T_{m}^{l-k}p_{m+1,l}=0}\\\Rightarrow\begin{bmatrix}\mathrm{A_{m}}&-\mathrm{A_{m+1}}\end{bmatrix}\begin{bmatrix}\mathrm{P_{m}}\\\mathrm{P_{m+1}}\end{bmatrix}=0\end{aligned}$$ Pm(k)​(Tm​)=Pm+1(k)​(Tm​)⇒i≥k∑​(i−k)!i!​Tmi−k​pm,i​−l≥k∑​(l−k)!l!​Tml−k​pm+1,l​=0⇒[Am​​−Am+1​​][Pm​Pm+1​​]=0​

微分约束

启动和终点的PVA使用确定值限制，即$P_0(0),P_0^{(1)}(0),P_0^{(2)}(0),P_M(T_M),P_M^{(1)}(T_M),P_M^{(2)}(T_M)$已知。

中间航迹点的位置确定，即$P_m(0)和P_M(T_m)$已知。

融合上述的约束得到：

[ A 0 0 0 0 ⋯ 0 A 0 − A 1 0 0 ⋯ 0 0 A 1 0 0 ⋯ 0 0 A 1 − A 2 0 ⋯ 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ 0 0 0 0 0 A M ] [ P 0 P 1 P 2 P 3 ⋮ P M ] = [ d 0 0 d 1 0 ⋮ d M ]

$$\begin{bmatrix}\mathrm{A}_0&0&0&0&\cdots&0\\\mathrm{A}_0&-\mathrm{A}_1&0&0&\cdots&0\\0&\mathrm{A}_1&0&0&\cdots&0\\0&\mathrm{A}_1&-\mathrm{A}_2&0&\cdots&0\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\0&0&0&0&0&\mathrm{A}_\mathrm{M}\end{bmatrix}$$ $$\begin{bmatrix}\mathrm{P}_0\\\mathrm{P}_1\\\mathrm{P}_2\\\mathrm{P}_3\\\vdots\\\mathrm{P}_\mathrm{M}\end{bmatrix}$$ = $$\begin{bmatrix}\mathrm{d}_0\\0\\\mathrm{d}_1\\0\\\vdots\\\mathrm{d}_\mathrm{M}\end{bmatrix}$$ ​A0​A0​00⋮0​0−A1​A1​A1​⋮0​000−A2​⋮0​0000⋮0​⋯⋯⋯⋯⋮0​0000⋮AM​​ ​ ​P0​P1​P2​P3​⋮PM​​ ​= ​d0​0d1​0⋮dM​​ ​

到此，整个QP问题构建完毕。

5.5 求解QP

多项式轨迹中的系数$p$没有具体的物理意义，通过无约束方程求解的p可能会出现一些特别小的值，受计算精度影响，导致数值不稳定。

Polynomial Trajectory Planning for Aggressive Quadrotor Flight in Dense Indoor Environments, Charles Richter, Adam Bry, and Nicholas Roy.

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

以最小化Jerk为例，此时是多项式是一个5次多项式轨迹，那么它有6个未知的系数，需要提供轨迹两端的位置、速度、加速度约束，对轨迹进行过次求导：

[ P m ( t ) P m ( 1 ) ( t ) P m ( 2 ) ( t ) ] = [ 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 5 p m , 5 t 4 4 p m , 4 t 3 3 p m , 3 t 2 2 p m , 2 t 1 p m , 1 t 0 0 5 ∗ 4 p m , 5 t 3 4 ∗ 3 p m , 4 t 2 3 ∗ 2 p m , 3 t 1 2 ∗ 1 p m , 2 t 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 ] [ p m , 5 p m , 4 p m , 3 p m , 2 p m , 1 p m , 0 ]

$$\begin{aligned}\begin{bmatrix}\mathrm{P}_{\mathrm{m}}\left(\mathrm{t}\right)\\\mathrm{P}_{\mathrm{m}}^{\left(1\right)}\left(\mathrm{t}\right)\\\mathrm{P}_{\mathrm{m}}^{\left(2\right)}\left(\mathrm{t}\right)\end{bmatrix}&=\begin{bmatrix}\mathrm{p}_{\mathrm{m},5}\mathrm{t}^5&\mathrm{p}_{\mathrm{m},4}\mathrm{t}^4&\mathrm{p}_{\mathrm{m},3}\mathrm{t}^3&\mathrm{p}_{\mathrm{m},2}\mathrm{t}^2&\mathrm{p}_{\mathrm{m},1}\mathrm{t}^1&\mathrm{p}_{\mathrm{m},0}\mathrm{t}^0\\5\mathrm{p}_{\mathrm{m},5}\mathrm{t}^4&4\mathrm{p}_{\mathrm{m},4}\mathrm{t}^3&3\mathrm{p}_{\mathrm{m},3}\mathrm{t}^2&2\mathrm{p}_{\mathrm{m},2}\mathrm{t}^1&\mathrm{p}_{\mathrm{m},1}\mathrm{t}^0&0\\5*4\mathrm{p}_{\mathrm{m},5}\mathrm{t}^3&4*3\mathrm{p}_{\mathrm{m},4}\mathrm{t}^2&3*2\mathrm{p}_{\mathrm{m},3}\mathrm{t}^1&2*1\mathrm{p}_{\mathrm{m},2}\mathrm{t}^0&0&0\end{bmatrix}\\&=\begin{bmatrix}\mathrm{t}^5&\mathrm{t}^4&\mathrm{t}^3&\mathrm{t}^2&\mathrm{t}^1&\mathrm{t}^0\\5\mathrm{t}^4&4\mathrm{t}^3&3\mathrm{t}^2&2\mathrm{t}^1&\mathrm{t}^0&0\\5*4\mathrm{t}^3&4*3\mathrm{t}^2&3*2\mathrm{t}^1&2*1\mathrm{t}^0&0&0\end{bmatrix}\begin{bmatrix}\mathrm{p}_{\mathrm{m},5}\\\mathrm{p}_{\mathrm{m},4}\\\mathrm{p}_{\mathrm{m},3}\\\mathrm{p}_{\mathrm{m},2}\\\mathrm{p}_{\mathrm{m},1}\\\mathrm{p}_{\mathrm{m},0}\end{bmatrix}\end{aligned}$$ ​Pm​(t)Pm(1)​(t)Pm(2)​(t)​ ​​= ​pm,5​t55pm,5​t45∗4pm,5​t3​pm,4​t44pm,4​t34∗3pm,4​t2​pm,3​t33pm,3​t23∗2pm,3​t1​pm,2​t22pm,2​t12∗1pm,2​t0​pm,1​t1pm,1​t00​pm,0​t000​ ​= ​t55t45∗4t3​t44t34∗3t2​t33t23∗2t1​t22t12∗1t0​t1t00​t000​ ​ ​pm,5​pm,4​pm,3​pm,2​pm,1​pm,0​​ ​​

带入两端的时间，得到多项式系数到运动微分的映射关系：

[ P ( 0 ) P ( 1 ) ( 0 ) P ( 2 ) ( 0 ) P ( T ) P ( 1 ) ( T ) P ( 2 ) ( T ) ] = [ 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 ] [ p 5 p 4 p 3 p 2 p 1 p 0 ] ⇒ d m = A m P m ⇒ P m = A m − 1 d m

$$\begin{aligned}\begin{bmatrix}\mathrm{P}\left(0\right)\\\mathrm{P}^{\left(1\right)}\left(0\right)\\\mathrm{P}^{\left(2\right)}\left(0\right)\\\mathrm{P}\left(\mathrm{T}\right)\\\mathrm{P}^{\left(1\right)}\left(\mathrm{T}\right)\\\mathrm{P}^{\left(2\right)}\left(\mathrm{T}\right)\end{bmatrix}&=\begin{bmatrix}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\\\mathrm{T}^5&\mathrm{T}^4&\mathrm{T}^3&\mathrm{T}^2&\mathrm{T}^1&\mathrm{T}^0\\5*\mathrm{T}^4&4*\mathrm{T}^3&3*\mathrm{T}^2&2*\mathrm{T}^1&\mathrm{T}^0&0\\5*4*\mathrm{T}^3&4*3*\mathrm{T}^2&3*2*\mathrm{T}^1&2*1*\mathrm{T}^0&0&0\end{bmatrix}\begin{bmatrix}\mathrm{p}_5\\\mathrm{p}_4\\\mathrm{p}_3\\\mathrm{p}_2\\\mathrm{p}_1\\\mathrm{p}_0\end{bmatrix}\\&\Rightarrow\mathrm{d}_{\mathrm{m}}=\mathrm{A}_{\mathrm{m}}\mathrm{P}_{\mathrm{m}}\\&\Rightarrow\mathrm{P}_{\mathrm{m}}=\mathrm{A}_{\mathrm{m}}^{-1}\mathrm{d}_{\mathrm{m}}\end{aligned}$$ ​P(0)P(1)(0)P(2)(0)P(T)P(1)(T)P(2)(T)​ ​​= ​055∗045∗4∗03T55∗T45∗4∗T3​044∗034∗3∗02T44∗T34∗3∗T2​033∗023∗2∗01T33∗T23∗2∗T1​022∗012∗1∗00T22∗T12∗1∗T0​01000T1T00​0000T000​ ​ ​p5​p4​p3​p2​p1​p0​​ ​⇒dm​=Am​Pm​⇒Pm​=Am−1​dm​​

得到映射矩阵：

A m = [ 0 N 0 N − 1 ⋯ 0 2 0 1 0 0 N ∗ 0 N − 1 ( N − 1 ) ∗ 0 N − 2 ⋯ 2 ∗ 0 1 0 0 0 N ( N − 1 ) ∗ 0 N − 2 ( N − 1 ) ( N − 2 ) ∗ 0 N − 3 ⋯ 2 ∗ 1 ∗ 0 0 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ T N T N − 1 ⋯ T 2 T 1 T 0 N ∗ T N − 1 ( N − 1 ) ∗ T N − 2 ⋯ 2 ∗ T 1 T 0 0 N ( N − 1 ) ∗ T N − 2 ( N − 1 ) ( N − 2 ) ∗ T N − 3 ⋯ 2 ∗ 1 ∗ T 0 0 0 ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ] \mathrm{A_m=

$$\begin{bmatrix}0^N&0^{N-1}&\cdots&0^2&0^1&0^0\\\mathrm{N*0^{N-1}}&(\mathrm{N-1})*0^{N-2}&\cdots&2*0^1&0^0&0\\\mathrm{N(N-1)*0^{N-2}}&(\mathrm{N-1})(\mathrm{N-2})*0^{N-3}&\cdots&2*1*0^0&0&0\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\mathrm{T^N}&\mathrm{T^{N-1}}&\cdots&\mathrm{T^2}&\mathrm{T^1}&\mathrm{T^0}\\\mathrm{N*T^{N-1}}&(\mathrm{N-1})*\mathrm{T^{N-2}}&\cdots&2*\mathrm{T^1}&\mathrm{T^0}&0\\\mathrm{N(N-1)*T^{N-2}}&(\mathrm{N-1})(\mathrm{N-2})*\mathrm{T^{N-3}}&\cdots&2*1*\mathrm{T^0}&0&0\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\end{bmatrix}$$ } 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_m$进行替换，得到下式：

min ⁡ [ P 0 P 1 ⋮ P M ] T [ Q 0 0 ⋯ ⋯ ⋯ 0 0 Q 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ Q M ] [ P 0 P 1 ⋮ P M ] min ⁡ [ d 1 d 1 ⋮ d M ] T [ A 0 0 ⋯ ⋯ 0 0 A 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ A M ] T [ Q 0 0 ⋯ ⋯ 0 0 Q 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ Q M ] [ A 0 0 ⋯ ⋯ 0 0 A 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ A M ] − 1 [ P 0 P 1 ⋮ P M ] ⇒ min ⁡ d T A − T Q A − 1 d \min\left[

$$\begin{array}{c}\mathrm{P_{0}}\\\mathrm{P_{1}}\\\vdots\\\mathrm{P_{M}}\end{array}$$ \right]^{\mathrm{T}} $$\begin{bmatrix}\mathrm{Q_{0}}&0&\cdots&\cdots&\cdots&0\\0&\mathrm{Q_{1}}&0&\cdots&0\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&\cdots&\mathrm{Q_{M}}\end{bmatrix}$$ $$\begin{bmatrix}\mathrm{P_{0}}\\\mathrm{P_{1}}\\\vdots\\\mathrm{P_{M}}\end{bmatrix}$$ \\\min\left[ $$\begin{array}{c}\mathrm{d_{1}}\\\mathrm{d_{1}}\\\vdots\\\mathrm{d_{M}}\end{array}$$ \right]^{\mathrm{T}} $$\begin{bmatrix}\mathrm{A_{0}}&0&\cdots&\cdots&0\\0&\mathrm{A_{1}}&0&\cdots&0\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\mathrm{A_{M}}\end{bmatrix}$$ ^{\mathrm{T}} $$\begin{bmatrix}\mathrm{Q_{0}}&0&\cdots&\cdots&0\\0&\mathrm{Q_{1}}&0&\cdots&0\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&\cdots&\mathrm{Q_{M}}\end{bmatrix}$$ $$\begin{bmatrix}\mathrm{A_{0}}&0&\cdots&\cdots&0\\0&\mathrm{A_{1}}&0&\cdots&0\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&\cdots&\mathrm{A_{M}}\end{bmatrix}$$ ^{-1} $$\begin{bmatrix}\mathrm{P_{0}}\\\mathrm{P_{1}}\\\vdots\\\mathrm{P_{M}}\end{bmatrix}$$ \\\Rightarrow\min\mathrm{d^{T}A^{-T}QA^{-1}d} min ​P0​P1​⋮PM​​ ​T ​Q0​0⋮0​0Q1​⋮0​⋯0⋮0​⋯⋯⋱⋯​⋯0⋮QM​​0 ​ ​P0​P1​⋮PM​​ ​min ​d1​d1​⋮dM​​ ​T ​A0​0⋮0​0A1​⋮0​⋯0⋮⋯​⋯⋯⋱AM​​00⋮​ ​T ​Q0​0⋮0​0Q1​⋮0​⋯0⋮0​⋯⋯⋱⋯​00⋮QM​​ ​ ​A0​0⋮0​0A1​⋮0​⋯0⋮0​⋯⋯⋱⋯​00⋮AM​​ ​−1 ​P0​P1​⋮PM​​ ​⇒mindTA−TQA−1d

如此一来，微分约束已经融入到了PQ问题当中，但是连续性约束还没有引入到代价函数当中。

本文中引入了一个permutation Matrix(置换矩阵)C来将连续性约束引入到代价函数当中。C矩阵本身是一个只包含0和1的矩阵，它的作用是将向量$[d_0,d_1,...,d_M]$进行一个组合，将固定的变量放在头部，将需要优化的变量放在尾部，并且对于连续性约束的变量值相等，因此只选择其中一个来表达连续性变量。

换言之，首尾状态和中间点航点是已知的设为$d_{mF}$，中间点除位置信息外其他状态是代价函数再优化时需要分配的最优值，即需要优化的变量设为$d_{mP}$。

将置换矩阵带入到代价函数当中：

J = min ⁡ [ d F d P ] T C A − T Q A − 1 C T [ d F d P ] \mathrm{J}=\min

$$\begin{bmatrix}\mathrm{d}_{\mathrm{F}}\\\mathrm{d}_{\mathrm{P}}\end{bmatrix}$$ ^{\mathrm{T}}\mathrm{C}\mathrm{A}^{{-\mathrm{T}}}\mathrm{Q}\mathrm{A}^{-1}\mathrm{C}^{\mathrm{T}} $$\begin{bmatrix}\mathrm{d}_{\mathrm{F}}\\\mathrm{d}_{\mathrm{P}}\end{bmatrix}$$ J=min[dF​dP​​]TCA−TQA−1CT[dF​dP​​]

令$R=C{A}^{-T}Q{A}^{-1}{C}^T$。对R矩阵根据$d_F$和$d_P$的尺寸进行分块，得到如下变换：

J = min ⁡ [ d F d P ] T C A − T Q A − 1 C T [ d F d P ] = J = min ⁡ [ d F d P ] T R [ d F d P ] = min ⁡ [ d F d P ] T [ R F F R F P R P F R P P ] [ d F d P ] = d F T R F F   d F   + d F T R F P   d P   + d P T R P F   d F   + d P T R P P   d P

$$\begin{gathered}\mathrm J=\operatorname*{min}\left[\begin{matrix}{\mathrm d_{\mathrm F}}\\{\mathrm d_{\mathrm P}}\\\end{matrix}\right]^{\mathrm T}\mathrm C\mathrm A^{-\mathrm T}\mathrm Q\mathrm A^{-1}\mathrm C^{\mathrm T}\left[\begin{matrix}{\mathrm d_{\mathrm F}}\\{\mathrm d_{\mathrm P}}\\\end{matrix}\right]=\mathrm J=\operatorname*{min}\left[\begin{matrix}{\mathrm d_{\mathrm F}}\\{\mathrm d_{\mathrm P}}\\\end{matrix}\right]^{\mathrm T}\mathrm R\left[\begin{matrix}{\mathrm d_{\mathrm F}}\\{\mathrm d_{\mathrm P}}\\\end{matrix}\right] \\=\min\begin{bmatrix}\mathrm{d_F}\\\mathrm{d_P}\end{bmatrix}^\mathrm{T}\begin{bmatrix}\mathrm{R_FF}&\mathrm{R_FP}\\\mathrm{R_PF}&\mathrm{R_PP}\end{bmatrix}\begin{bmatrix}\mathrm{d_F}\\\mathrm{d_P}\end{bmatrix} \\=\mathrm{d_F^TR_{FF}~d_F~+d_F^TR_{FP}~d_P~+d_P^TR_{PF}~d_F~+d_P^TR_{PP}~d_P} \end{gathered}$$ J=min[dF​dP​​]TCA−TQA−1CT[dF​dP​​]=J=min[dF​dP​​]TR[dF​dP​​]=min[dF​dP​​]T[RF​FRP​F​RF​PRP​P​][dF​dP​​]=dFT​RFF​ dF​ +dFT​RFP​ dP​ +dPT​RPF​ dF​ +dPT​RPP​ dP​​

对于上式J是一个标量，并且Q矩阵是对称阵，所以R矩阵也是对称阵，则$R_{PF}=R_{FP}^T$。对于上式我们要求的变量是$d_p$，对其求导得到

∂ J ∂ d P = 2 d F T R F P + 2 R P P   d P = 0 ⇒ d P ∗ = − R P P − 1 R F P T   d F

$$\begin{aligned}\frac{\partial\mathcal{J}}{\partial\mathcal{d}_{ {\mathcal{P}}}}&=2\mathrm{d}_{ {\mathcal{F}}}^{ {\mathcal{T}}}\mathcal{R}_{ {\mathcal{F}\mathcal{P}}}+2\mathrm{R}_{ {\mathcal{P}\mathcal{P}}}\mathrm{~d}_{ {\mathcal{P}}}=0\\&\Rightarrow\mathrm{d}_{ {\mathcal{P}}}^{*}=-\mathrm{R}_{ {\mathcal{P}\mathcal{P}}}^{-1}\mathrm{R}_{ {\mathcal{F}\mathcal{P}}}^{ {\mathcal{T}}}\mathrm{~d}_{ {\mathcal{F}}}\end{aligned}$$ ∂dP​∂J​​=2dFT​RFP​+2RPP​ dP​=0⇒dP∗​=−RPP−1​RFPT​ dF​​

然后反求得到多项式系数：

P   = A − 1 C T   [ d F d P ∗ ] \mathrm{P~=A^{-1}C^T~

$$\begin{bmatrix}d_F\\d_P^*\end{bmatrix}$$ } P =A−1CT [dF​dP∗​​]