Apollo 7.0 PiecewiseJerkSpeedOptimizer(分段加加速度优化器)代码解读

Apollo 7.0 PiecewiseJerkSpeedOptimizer(分段加加速度优化器)代码解读

贺志国
2022-7-30

PiecewiseJerkSpeedOptimizer (对应PIECEWISE_JERK_SPEED_OPTIMIZER任务)使用QP来求解速度曲线。Apollo 3.0中的QpSplineStSpeedOptimizer使用多段五次多项式连接的样条曲线(每2s使用一个五次多项式拟合,8s时长的轨迹共有8/2=4段五次多项式。实际求解过程中,为了降低计算量,实际求解中仅使用了3段多项式来求解[0, 6]s区间内的速度曲线。这些五次多项式隐含的控制点和节点向量满足样条曲线的迭代表达式就形成了五次多项式样条曲线)来拟合速度,使用二次规划来优化五次多项式的系数,并使用qpOASES库来求解,之后再使用求解出的多项式系数来确定自变量为时间t,因变量为s的五次多项式样条曲线,再通过样条曲线来计算出不同时刻的位移s,最后借助运动学公式计算出指定位移s处的速度v和加速度a。Apollo 7.0不使用样条曲线而是使用位移s关于时间t的三阶偏导(称为jerk)为分段常量作为优化项来建模二次规划问题,并使用OSQP库来求解,并且求解的结果直接就是多个采样点处的s, v, a(采样时间间隔为0.1 s),采样点以外的s, v, a未给出,直接由控制模块自行插值获取。

一、数学模型

1.1 凸二次规划介绍
1.1.1 凸二次规划的一般表达式

m i n 1 2 x T P x + q T x s . t . l ≤ A x ≤ u min \quad \frac{1}{2}x^T Px + q^Tx \qquad s.t. \quad l \leq Ax \leq u min21xTPx+qTxs.t.lAxu
其中, P P P称为内核矩阵,要求是半正定矩阵。正定矩阵在复数域上是共轭对称矩阵(Hermite矩阵),在实数域上是对称矩阵。因为我们在实数域内求解, P P P也是一个对称矩阵。 q q q称为偏移向量, A A A称为仿射矩阵。

凸二次规划的求解算法主要有如下几种:

  • 有效集法:Active Set Method
  • 内点法:Interior Point Method
  • 一阶方法:First Order Method
  • 交替方向乘子法:Alternating Direction Method of Multipliers(ADMM)

PiecewiseJerkSpeedOptimizer使用的求解库OSQP利用交替方向乘子法(ADMM)。Apollo 3.0中QpSplineStSpeedOptimizer使用的求解库qpOASES则利用有效集法。从网上给出的测评报告看,OSQP库的求解效率比qpOASES库要高很多。

1.1.2 凸二次规划的示例

为了形象地说明凸二次规划的求解过程,下面给出一个简单的示例:
m i n 1 2 x T [ 4 1 1 2 ] x + [ 1 1 ] T x min \quad \frac{1}{2}x^T\left[\begin{matrix} 4 & 1\\ 1 & 2\end{matrix}\right]x + \left[\begin{matrix} 1 \\ 1\end{matrix}\right]^Tx min21xT[4112]x+[11]Tx
s.t. (subject to):
[ 1 0 0 ] ≤ [ 1 1 1 0 0 1 ] x ≤ [ 1 0.7 0.7 ] \left[\begin{matrix}1 \\ 0 \\ 0 \end{matrix}\right] \leq \left[\begin{matrix} 1 & 1 \\ 1 & 0 \\ 0 & 1 \end{matrix}\right]x \leq \left[\begin{matrix}1 \\ 0.7 \\ 0.7 \end{matrix}\right] 100 110101 x 10.70.7
本示例中,内核矩阵 P = [ 4 1 1 2 ] P=\left[\begin{matrix} 4 & 1\\1 & 2\end{matrix}\right] P=[4112]的元素全部为实数,因此它一定是对称矩阵。对于对称矩阵,只需要存储上三角矩阵即可,这也是OSQP存储稀疏矩阵的常见做法,当然也是该库提高计算效率的有效措施。

具体而言,OSQP库使用csc_matrix(indices, indptr, data, csc = Compressed Sparse Column)的方式来存储一个矩阵。csc_matrix存储格式定义如下:

  • indices is array of row indices
  • data is array of corresponding nonzero values
  • indptr points to column starts in indices and data

下面仍然以内核矩阵 P = [ 4 1 1 2 ] P=\left[\begin{matrix} 4 & 1\\1 & 2\end{matrix}\right] P=[4112]为例进行说明。首先,将内核矩阵简化为上三角矩阵:
P = [ 4 1 0 2 ] P=\left[\begin{matrix} 4 & 1\\0 & 2\end{matrix}\right] P=[4012]
P P P的非零元素个数为3,记为:P_nnz = 3 (nnz: number of non-zero elements)

P P P中非零元素为:4, 1, 2,记为:P_x = {4, 1, 2};

P P P中非零元素行索引为:0, 0, 1,这是按照先第一列,再第二列。。。,最后第N列的顺序标记。第一列的非零元素是4,行索引为0;第二列的非零元素为1, 2,行索引分别为:0, 1,于是记为:P_i = {0, 0, 1};

P P P中第一列非零元素是1个,第二列非零元素是2个,记为:P_p = {0, 1, 3}。这个记法有些反人类,也很难形象理解,容我再详细解释一番。P_p元素个数为矩阵列数加1P_p的第一个元素为0,P_p中前后两个元素的差值表示矩阵P相应列的非零元素个数。例如:P_p = {0, 1, 3}表示:1 - 0 = 1 代表矩阵P第一列元素非零个数为1个,3 - 1 = 2 代表矩阵P第二列元素非零个数为2个。助记方法:P_i (i代表indices), P_p(p代表indptr)。

根据上述介绍,下面给出该示例的求解代码:

#include "osqp.h"
int main(int argc, char **argv) {
    // P矩阵是对称矩阵,只存储上三角部分,下三角部分视为零值
    c_float P_x[3] = {4.0, 1.0, 2.0, }; 
    // P矩阵非零元素个数为3
    // nnz = number of non-zero elements.
    c_int P_nnz = 3; 
    // P矩阵非零元素行索引,按照先第一列,再第二列。。。
    // 的顺序排列。下例表示非零元素的行序号为0、0、1
    c_int P_i[3] = {0, 0, 1, }; 
    // 1-0=1表示第0列非零元素个数
    // 3-1=2表示第1列非零元素个数
    c_int P_p[3] = {0, 1, 3, }; 
    c_float q[2] = {1.0, 1.0, };
    // A矩阵非零元素
    c_float A_x[4] = {1.0, 1.0, 1.0, 1.0, };
    // A矩阵非零元素个数为4
    c_int A_nnz = 4;
    // A矩阵非零元素的行序号为0、1、0、2(按列排序)
    c_int A_i[4] = {0, 1, 0, 2, };
    // 2-0=2表示第0列2个元素非零
    // 4-2=2表示第1列2个元素非零
    c_int A_p[3] = {0, 2, 4, };
    c_float l[3] = {1.0, 0.0, 0.0, };
    c_float u[3] = {1.0, 0.7, 0.7, };
    c_int n = 2;
    c_int m = 3;

    // Exitflag
    c_int exitflag = 0;
    // Workspace structures
    OSQPWorkspace *work;
    OSQPSettings  *settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings));
    OSQPData      *data     = (OSQPData *)c_malloc(sizeof(OSQPData));

    // Populate data
    if (data) {
        data->n = n;
        data->m = m;
        data->P = csc_matrix(data->n, data->n, P_nnz, P_x, P_i, P_p);
        data->q = q;
        data->A = csc_matrix(data->m, data->n, A_nnz, A_x, A_i, A_p);
        data->l = l;
        data->u = u;
    }

    // Define solver settings as default
    if (settings) {
        osqp_set_default_settings(settings);
        settings->alpha = 1.0; // Change alpha parameter
    }

    // Setup workspace
    exitflag = osqp_setup(&work, data, settings);

    // Solve Problem
    osqp_solve(work);

    // Cleanup
    if (data) {
        if (data->A) c_free(data->A);
        if (data->P) c_free(data->P);
        c_free(data);
    }
    if (settings) c_free(settings);

    return exitflag;
}
1.2 速度规划待优化的状态变量

在这里插入图片描述

1.3 速度规划目标函数

我们希望速度曲线尽可能贴近启发式搜索给出的位移和限速,同时希望加速度平滑,加加速度为常量,得出目标函数如下:

在这里插入图片描述

c o s t = m i n { c o s t 1 + c o s t 2 + c o s t 3 } cost = min \left\{ cost1 + cost2 + cost3 \right\} cost=min{cost1+cost2+cost3}

其中:

c o s t 1 = ∑ i = 0 n − 2 [ w s s i 2 + ( w v + p i ) s ˙ i 2 − 2 w j e r k Δ t 2 s ¨ i s ¨ i + 1 ] + ∑ i = 1 n − 2 ( w a + 2 w j e r k Δ t 2 ) s ¨ i 2 cost1 = \sum_{i=0}^{n-2} \left[w_s s_i^2 + (w_v+p_i)\dot s_i^2 - \frac{2w_{jerk}}{\Delta t^2}\ddot s_i \ddot s_{i+1}\right] + \sum_{i=1}^{n-2}(w_a + \frac{2w_{jerk}}{\Delta t^2})\ddot s_i^2 cost1=i=0n2[wssi2+(wv+pi)s˙i2Δt22wjerks¨is¨i+1]+i=1n2(wa+Δt22wjerk)s¨i2
c o s t 2 = ( w s + w s e n d ) s n − 1 2 + ( w v + p n − 1 + w v e n d ) s ˙ n − 1 2 + ( w a + w j e r k Δ t 2 ) s ¨ 0 2 + ( w a + w j e r k Δ t 2 + w a e n d ) s ¨ n − 1 2 cost2 = (w_s + w_{s_{end}})s_{n-1}^2 + (w_v + p_{n-1} + w_{v_{end}})\dot s_{n-1}^2+(w_a + \frac{w_{jerk}}{\Delta t^2})\ddot s_0^2 + (w_a + \frac{w_{jerk}}{\Delta t^2} + w_{a_{end}})\ddot s_{n-1}^2 cost2=(ws+wsend)sn12+(wv+pn1+wvend)s˙n12+(wa+Δt2wjerk)s¨02+(wa+Δt2wjerk+waend)s¨n12
c o s t 3 = − 2 ∑ i = 0 n − 1 [ w s s r e f i s i + ( w v + p i ) v r e f i s ˙ i ] cost3 = -2\sum_{i=0}^{n-1} \left[ w_s s_{ref_i} s_i + (w_v + p_i) v_{ref_i} \dot s_i\right] cost3=2i=0n1[wssrefisi+(wv+pi)vrefis˙i]
注意:对于终点处的 s n − 1 s_{n-1} sn1额外增加了惩罚系数 w s e n d w_{s_{end}} wsend,对于终点处的 s ˙ n − 1 \dot s_{n-1} s˙n1额外增加了惩罚系数 w v e n d w_{v_{end}} wvend,对于终点处的 s ¨ n − 1 \ddot s_{n-1} s¨n1额外增加了惩罚系数 w a e n d w_{a_{end}} waend

预备知识

  1. 矩阵与列向量相乘,等于矩阵中的每一列的线性组合,组合系数为列向量中的每个元素:

[ a b c d e f ] [ m n o ] = m [ a d ] + n [ b e ] + o [ c f ] \left[\begin{matrix} a & b & c \\ d & e & f\end{matrix}\right] \left[\begin{matrix} m \\ n \\ o \end{matrix}\right] = m \left[\begin{matrix} a \\ d \end{matrix}\right] + n \left[\begin{matrix} b \\ e \end{matrix}\right] + o \left[\begin{matrix} c \\ f \end{matrix}\right] [adbecf] mno =m[ad]+n[be]+o[cf]
2. 行向量与矩阵相乘,等于矩阵中的每一行的线性组合,组合系数为行向量中的每个元素:
[ m n ] [ a b c d e f ] = m [ a b c ] + n [ d e f ] \left[\begin{matrix} m & n \end{matrix}\right] \left[\begin{matrix} a & b & c \\ d & e & f\end{matrix}\right] = m \left[\begin{matrix} a & b & c \end{matrix}\right] + n \left[\begin{matrix} d & e & f \end{matrix}\right] [mn][adbecf]=m[abc]+n[def]
基于上述预备知识,将上述表达式写成矩阵形式:
m i n 1 2 [ s 0 s 1 ⋮ s n − 1 s ˙ 0 s ˙ 1 ⋮ s ˙ n − 1 s ¨ 0 s ¨ 1 ⋮ s ¨ n − 1 ] T [ 2 A n ∗ n 0 0 0 2 B n ∗ n 0 0 0 2 C n ∗ n ] ⏟ P [ s 0 s 1 ⋮ s n − 1 s ˙ 0 s ˙ 1 ⋮ s ˙ n − 1 s ¨ 0 s ¨ 1 ⋮ s ¨ n − 1 ] + [ − 2 w s s r e f 0 − 2 w s s r e f 1 ⋮ − 2 w s s r e f n − 1 − 2 ( w v + p 0 ) v r e f 0 − 2 ( w v + p 1 ) v r e f 1 ⋮ − 2 ( w v + p n − 1 ) v r e f n − 1 ] T ⏟ q [ s 0 s 1 ⋮ s n − 1 s ˙ 0 s ˙ 1 ⋮ s ˙ n − 1 ] min \quad \frac{1}{2}\left[\begin{matrix} s_0 \\ s_1 \\ \vdots \\ s_{n-1} \\ \dot s_0 \\ \dot s_1 \\ \vdots \\ \dot s_{n-1} \\ \ddot s_0 \\ \ddot s_1 \\ \vdots \\ \ddot s_{n-1} \end{matrix}\right]^T \underbrace{\left[\begin{matrix} 2A_{n*n} & 0 & 0 & \\ 0 & 2B_{n*n} & 0 \\ 0 & 0 & 2C_{n*n} \end{matrix}\right]}_{\mathbf{P}}\left[\begin{matrix} s_0 \\ s_1 \\ \vdots \\ s_{n-1} \\ \dot s_0 \\ \dot s_1 \\ \vdots \\ \dot s_{n-1} \\ \ddot s_0 \\ \ddot s_1 \\ \vdots \\ \ddot s_{n-1} \end{matrix}\right] + \underbrace{\left[\begin{matrix} -2w_s s_{ref_0} \\ -2w_s s_{ref_1} \\ \vdots \\ -2w_s s_{ref_{n-1}} \\ -2(w_v+p_0) v_{ref_0} \\ -2(w_v+p_1) v_{ref_1} \\ \vdots \\ -2(w_v+p_{n-1}) v_{ref_{n-1}} \end{matrix}\right]^T}_{\mathbf{q}}\left[\begin{matrix} s_0 \\ s_1 \\ \vdots \\ s_{n-1} \\ \dot s_0 \\ \dot s_1 \\ \vdots \\ \dot s_{n-1}\end{matrix}\right] min21 s0s1sn1s˙0s˙1s˙n1s¨0s¨1s¨n1 TP 2Ann0002Bnn0002Cnn s0s1sn1s˙0s˙1s˙n1s¨0s¨1s¨n1 +q 2wssref02wssref12wssrefn12(wv+p0)vref02(wv+p1)vref12(wv+pn1)vrefn1 T s0s1sn1s˙0s˙1s˙n1

接下来我们写出矩阵 A n ∗ n A_{n*n} Ann B n ∗ n B_{n*n} Bnn C n ∗ n C_{n*n} Cnn的具体表达式(空余部分元素全部为0):
A n ∗ n = [ w s w s ⋱ w s w s + w s e n d ] A_{n*n}=\left[\begin{matrix}w_s\\&w_s\\&&\ddots\\&&&w_s\\&&&&w_s+w_{s_{end}}\end{matrix}\right] Ann= wswswsws+wsend

B n ∗ n = [ w v + p 0 w v + p 1 ⋱ w v + p n − 2 w v + p n − 1 + w v e n d ] B_{n*n}=\left[\begin{matrix}w_v+p_0\\&w_v+p_1\\&&\ddots\\&&&w_v+p_{n-2}\\&&&&w_v+p_{n-1}+w_{v_{end}}\end{matrix}\right] Bnn= wv+p0wv+p1wv+pn2wv+pn1+wvend


C n ∗ n = [ w a + w j e r k Δ t 2 − 2 w j e r k Δ t 2 w a + 2 w j e r k Δ t 2 − 2 w j e r k Δ t 2 ⋱ ⋱ w a + 2 w j e r k Δ t 2 − 2 w j e r k Δ t 2 w a + 2 w j e r k Δ t 2 − 2 w j e r k Δ t 2 w a + w j e r k Δ t 2 + w a e n d ] C_{n*n}=\left[\begin{matrix}w_a+\frac{w_{jerk}}{\Delta t^2}&-\frac{2w_{jerk}}{\Delta t^2}\\&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{2w_{jerk}}{\Delta t^2}\\&&\ddots&\ddots\\&&&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{2w_{jerk}}{\Delta t^2}\\&&&&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{2w_{jerk}}{\Delta t^2}\\&&&&&w_a+\frac{w_{jerk}}{\Delta t^2} + w_{a_{end}}\end{matrix}\right] Cnn= wa+Δt2wjerkΔt22wjerkwa+Δt22wjerkΔt22wjerkwa+Δt22wjerkΔt22wjerkwa+Δt22wjerkΔt22wjerkwa+Δt2wjerk+waend

注意:上述矩阵 C n ∗ n C_{n*n} Cnn不是对称的,按照OSQP要求转换为对称矩阵:
C n ∗ n = [ w a + w j e r k Δ t 2 − w j e r k Δ t 2 − w j e r k Δ t 2 w a + 2 w j e r k Δ t 2 − w j e r k Δ t 2 ⋱ ⋱ ⋱ − w j e r k Δ t 2 w a + 2 w j e r k Δ t 2 − w j e r k Δ t 2 − w j e r k Δ t 2 w a + w j e r k Δ t 2 + w a e n d ] C_{n*n}=\left[\begin{matrix}w_a+\frac{w_{jerk}}{\Delta t^2}&-\frac{w_{jerk}}{\Delta t^2}\\-\frac{w_{jerk}}{\Delta t^2}&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{w_{jerk}}{\Delta t^2}\\&\ddots&\ddots&\ddots\\&&-\frac{w_{jerk}}{\Delta t^2}&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{w_{jerk}}{\Delta t^2}\\&&&-\frac{w_{jerk}}{\Delta t^2} &w_a+\frac{w_{jerk}}{\Delta t^2} + w_{a_{end}}\end{matrix}\right] Cnn= wa+Δt2wjerkΔt2wjerkΔt2wjerkwa+Δt22wjerkΔt2wjerkΔt2wjerkwa+Δt22wjerkΔt2wjerkΔt2wjerkwa+Δt2wjerk+waend
根据OSQP对于内核矩阵的存储要求,保留上三角矩阵,将下三角矩阵的非零元素全部抹去,结果为(Apollo转换对称矩阵的做法是错误的):
C n ∗ n = [ w a + w j e r k Δ t 2 − w j e r k Δ t 2 w a + 2 w j e r k Δ t 2 − w j e r k Δ t 2 ⋱ ⋱ ⋱ w a + 2 w j e r k Δ t 2 − w j e r k Δ t 2 w a + w j e r k Δ t 2 + w a e n d ] C_{n*n}=\left[\begin{matrix}w_a+\frac{w_{jerk}}{\Delta t^2}&-\frac{w_{jerk}}{\Delta t^2}\\&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{w_{jerk}}{\Delta t^2}\\&\ddots&\ddots&\ddots\\&&&w_a+\frac{2w_{jerk}}{\Delta t^2}&-\frac{w_{jerk}}{\Delta t^2}\\&&&&w_a+\frac{w_{jerk}}{\Delta t^2} + w_{a_{end}}\end{matrix}\right] Cnn= wa+Δt2wjerkΔt2wjerkwa+Δt22wjerkΔt2wjerkwa+Δt22wjerkΔt2wjerkwa+Δt2wjerk+waend
内核矩阵 P P P中的系数就是PiecewiseJerkSpeedProblem(在文件modules/planning/math/piecewise_jerk/piecewise_jerk_speed_problem.cc中)的CalculateKernel函数中的系数,偏移向量q中的系数就是PiecewiseJerkSpeedProblemCalculateOffset函数中的系数(注意Apollo对于偏移向量的实现忽略了曲率惩罚系数 p i p_i pi,我认为这是不正确的)。

1.4 速度规划约束条件

在这里插入图片描述

s l b [ i ] ≤ s i ≤ s u b [ i ] s_{lb}[i] \leq s_i \leq s_{ub}[i] slb[i]sisub[i]
s ˙ l b [ i ] ≤ s ˙ i ≤ s ˙ u b [ i ] \dot s_{lb}[i] \leq \dot s_i \leq \dot s_{ub}[i] s˙lb[i]s˙is˙ub[i]
s ¨ l b [ i ] ≤ s ¨ i ≤ s ¨ l b [ i ] \ddot s_{lb}[i] \leq \ddot s_i \leq \ddot s_{lb}[i] s¨lb[i]s¨is¨lb[i]
j e r k l b Δ t ≤ s ¨ i + 1 − s ¨ i ≤ j e r k u b Δ t jerk_{lb} \Delta t \leq \ddot s_{i+1} - \ddot s_{i} \leq jerk_{ub} \Delta t jerklbΔts¨i+1s¨ijerkubΔt

其中,位移下、上限 s l b [ i ] s_{lb}[i] slb[i] s u b [ i ] s_{ub}[i] sub[i]通过动态规划算法在ST图上搜索得到,速度下、上限 s ˙ l b [ i ] \dot s_{lb}[i] s˙lb[i] s ˙ l b [ i ] \dot s_{lb}[i] s˙lb[i]通过限速决策得到,加速度下、上限 s ¨ l b [ i ] \ddot s_{lb}[i] s¨lb[i] s ¨ u b [ i ] \ddot s_{ub}[i] s¨ub[i]通过车辆运动学模型参数给出,加加速度下、下限 j e r k l b jerk_{lb} jerklb j e r k u b jerk_{ub} jerkub通过调参经验得到。

  • 连续性约束
    在这里插入图片描述

只要令上、下限为0,上述条件就可转换为OSQP要求的 l ≤ A x ≤ u l \leq Ax \leq u lAxu形式。

  • 初始条件约束

s 0 = s i n i t s_0 = s_{init} s0=sinit
s ˙ 0 = v i n i t \dot s_0 = v_{init} s˙0=vinit
s ¨ 0 = a i n i t \ddot s_0 = a_{init} s¨0=ainit

将上述表达式写成矩阵为:
[ s l b [ 0 ] ⋮ s l b [ n − 1 ] s ˙ l b [ 0 ] ⋮ s ˙ l b [ n − 1 ] s ¨ l b [ 0 ] ⋮ s ¨ l b [ n − 1 ] j e r k l b Δ t ⋮ j e r k l b Δ t 0.0 ⋮ 0.0 s i n i t v i n i t a i n i t ] 6 n ∗ 1 ≤ A 6 n ∗ 3 n [ s 0 ⋮ s n − 1 s ˙ 0 ⋮ s ˙ n − 1 s ¨ 0 s ¨ 1 ⋮ s ¨ n − 1 ] 3 n ∗ 1 ≤ [ s u b [ 0 ] ⋮ s u b [ n − 1 ] s ˙ u b [ 0 ] ⋮ s ˙ u b [ n − 1 ] s ¨ u b [ 0 ] ⋮ s ¨ u b [ n − 1 ] j e r k u b Δ t ⋮ j e r k u b Δ t 0.0 ⋮ 0.0 s i n i t v i n i t a i n i t ] 6 n ∗ 1 \left[\begin{matrix} s_{lb}[0] \\ \vdots \\ s_{lb}[n-1] \\ \dot s_{lb}[0] \\ \vdots \\ \dot s_{lb}[n-1] \\ \ddot s_{lb}[0] \\ \vdots \\ \ddot s_{lb}[n-1] \\ jerk_{lb} \Delta t \\ \vdots \\ jerk_{lb} \Delta t \\ 0.0 \\ \vdots \\ 0.0 \\ s_{init} \\ v_{init} \\ a_{init} \end{matrix}\right]_{6n*1} \leq A_{6n*3n} \left[\begin{matrix} s_0 \\\vdots \\ s_{n-1} \\ \dot s_0 \\ \vdots \\ \dot s_{n-1} \\ \ddot s_0 \\ \ddot s_1 \\ \vdots \\ \ddot s_{n-1} \end{matrix}\right]_{3n*1} \leq \left[\begin{matrix} s_{ub}[0] \\ \vdots \\ s_{ub}[n-1] \\ \dot s_{ub}[0] \\ \vdots \\ \dot s_{ub}[n-1] \\ \ddot s_{ub}[0] \\ \vdots \\ \ddot s_{ub}[n-1] \\ jerk_{ub} \Delta t \\ \vdots \\ jerk_{ub} \Delta t \\ 0.0 \\ \vdots \\ 0.0 \\ s_{init} \\ v_{init} \\ a_{init} \end{matrix}\right]_{6n*1} slb[0]slb[n1]s˙lb[0]s˙lb[n1]s¨lb[0]s¨lb[n1]jerklbΔtjerklbΔt0.00.0sinitvinitainit 6n1A6n3n s0sn1s˙0s˙n1s¨0s¨1s¨n1 3n1 sub[0]sub[n1]s˙ub[0]s˙ub[n1]s¨ub[0]s¨ub[n1]jerkubΔtjerkubΔt0.00.0sinitvinitainit 6n1

**仿射矩阵 A 6 n ∗ 3 n A_{6n*3n} A6n3n**的表达式如下所示(空白部分的元素为0):
A = [ 1 1 ⋱ 1 n ∗ n 1 1 ⋱ 1 n ∗ n 1 1 ⋱ 1 n ∗ n − 1 1 ⋱ ⋱ − 1 1 ( n − 1 ) ∗ n − 1 1 − Δ t 2 − Δ t 2 ⋱ ⋱ ⋱ ⋱ − 1 1 ( n − 1 ) ∗ n − Δ t 2 − Δ t 2 ( n − 1 ) ∗ n − 1 1 − Δ t − Δ t 2 3 − Δ t 2 6 ⋱ ⋱ ⋱ ⋱ ⋱ − 1 1 ( n − 1 ) ∗ n − Δ t ( n − 1 ) ∗ n − Δ t 2 3 − Δ t 2 6 ( n − 1 ) ∗ n 1 0 ⋯ 0 1 ∗ n 1 0 ⋯ 0 1 ∗ n 1 0 ⋯ 0 1 ∗ n ] 6 n ∗ 3 n A = \left[\begin{matrix}1\\ &1 \\ &&\ddots\\ &&&1_{n*n} \\&&&&1 \\&&&&&1 \\ &&&&&&\ddots \\ &&&&&&&1_{n*n} \\ &&&&&&&&1 \\ &&&&&&&&&1 \\ &&&&&&&&&&\ddots \\ &&&&&&&&&&&1_{n*n} \\ &&&&&&&&-1 &1 \\ &&&&&&&&&\ddots &\ddots \\ &&&&&&&&&&-1 &1_{(n-1)*n} \\ &&&&-1 &1 &&&-\frac{\Delta{t}}{2} &-\frac{\Delta{t}}{2} \\ &&&&&\ddots &\ddots &&&\ddots &\ddots \\ &&&&&&-1 &1_{(n-1)*n} &&&-\frac{\Delta{t}}{2} &-\frac{\Delta{t}}{2}_{(n-1)*n} \\ -1 &1 &&&-\Delta{t} &&&& -\frac{\Delta{t}^2}{3} &-\frac{\Delta{t}^2}{6} \\ &\ddots &\ddots &&&\ddots &&&&\ddots &\ddots \\ && -1 & 1_{(n-1)*n} &&&-\Delta{t}_{(n-1)*n} &&&& -\frac{\Delta{t}^2}{3} &-\frac{\Delta{t}^2}{6}_{(n-1)*n} \\ 1 & 0 &\cdots &0_{1*n} \\ &&&&1 & 0 &\cdots &0_{1*n}\\ &&&&&&&&1 & 0 &\cdots &0_{1*n} \end{matrix}\right]_{6n*3n} A= 11111011nn1(n1)n01n11Δt11101Δt(n1)n1nn1(n1)n01n112Δt3Δt21112Δt6Δt2012Δt3Δt21nn1(n1)n2Δt(n1)n6Δt2(n1)n01n 6n3n

注意:上式中的下标 1 n ∗ n 1_{n*n} 1nn表示某块矩阵的行列数是 n ∗ n n*n nn。仿射矩阵 A A A中的系数也就是PiecewiseJerkProblemCalculateAffineConstraint函数(在文件modules/planning/math/piecewise_jerk/piecewise_jerk_problem.cc中)的系数。仿射矩阵不是对称矩阵,OSQP不要求使用上三角矩阵的形式存储。

优化求解采用OSQP库来实现。

二、流程图

流程图

三、主要代码解读

  • 调用入口
  Status PiecewiseJerkSpeedOptimizer::Process(
           const PathData& path_data, 
           const TrajectoryPoint& init_point, 
           SpeedData* const speed_data)
  • 如果抵达终点, 直接返回OK
  if (reference_line_info_->ReachedDestination()) {
    return Status::OK();
  }
  • 如果无路径点,返回错误
  if (path_data.discretized_path().empty()) {
    const std::string msg = "Empty path data";
    AERROR << msg;
    return Status(ErrorCode::PLANNING_ERROR, msg);
  }
  • 获取ST图数据和车辆参数
  StGraphData& st_graph_data = *reference_line_info_->mutable_st_graph_data();

  const auto& veh_param =
      common::VehicleConfigHelper::GetConfig().vehicle_param();
  • 设置状态变量的初始值
  std::array<double, 3> init_s = {0.0, st_graph_data.init_point().v(),
                                  st_graph_data.init_point().a()};
  • 设置时间采样间隔、路径总长度、总时长、状态变量的个数
  double delta_t = 0.1;
  double total_length = st_graph_data.path_length();
  double total_time = st_graph_data.total_time_by_conf();
  int num_of_knots = static_cast<int>(total_time / delta_t) + 1;
  • 创建PiecewiseJerkSpeedProblem对象
 PiecewiseJerkSpeedProblem piecewise_jerk_problem(num_of_knots, delta_t,
                                                   init_s);
  • 设置权重系数、位移、速度、加速度、加加速度的范围(即上、下限)、巡航速度
  const auto& config = config_.piecewise_jerk_speed_optimizer_config();
  piecewise_jerk_problem.set_weight_ddx(config.acc_weight());
  piecewise_jerk_problem.set_weight_dddx(config.jerk_weight());

  piecewise_jerk_problem.set_x_bounds(0.0, total_length);
  piecewise_jerk_problem.set_dx_bounds(
      0.0, std::fmax(FLAGS_planning_upper_speed_limit,
                     st_graph_data.init_point().v()));
  piecewise_jerk_problem.set_ddx_bounds(veh_param.max_deceleration(),
                                        veh_param.max_acceleration());
  piecewise_jerk_problem.set_dddx_bound(FLAGS_longitudinal_jerk_lower_bound,
                                        FLAGS_longitudinal_jerk_upper_bound);

  piecewise_jerk_problem.set_dx_ref(config.ref_v_weight(),
                                    reference_line_info_->GetCruiseSpeed());
  • 根据前面动态规划获取到的初始位移点,更新位移范围(即上、下限)。
   // Update STBoundary
  std::vector<std::pair<double, double>> s_bounds;
  for (int i = 0; i < num_of_knots; ++i) {
    double curr_t = i * delta_t;
    double s_lower_bound = 0.0;
    double s_upper_bound = total_length;
    for (const STBoundary* boundary : st_graph_data.st_boundaries()) {
      double s_lower = 0.0;
      double s_upper = 0.0;
      if (!boundary->GetUnblockSRange(curr_t, &s_upper, &s_lower)) {
        continue;
      }
      switch (boundary->boundary_type()) {
        case STBoundary::BoundaryType::STOP:
        case STBoundary::BoundaryType::YIELD:
          s_upper_bound = std::fmin(s_upper_bound, s_upper);
          break;
        case STBoundary::BoundaryType::FOLLOW:
          // TODO(Hongyi): unify follow buffer on decision side
          s_upper_bound = std::fmin(s_upper_bound, s_upper - 8.0);
          break;
        case STBoundary::BoundaryType::OVERTAKE:
          s_lower_bound = std::fmax(s_lower_bound, s_lower);
          break;
        default:
          break;
      }
    }
    if (s_lower_bound > s_upper_bound) {
      const std::string msg =
          "s_lower_bound larger than s_upper_bound on STGraph";
      AERROR << msg;
      speed_data->clear();
      return Status(ErrorCode::PLANNING_ERROR, msg);
    }
    s_bounds.emplace_back(s_lower_bound, s_upper_bound);
  }
  piecewise_jerk_problem.set_x_bounds(std::move(s_bounds));
  • 根据ST图中的限速决策数据,更新速度范围(即上、下限)
 // Update SpeedBoundary and ref_s
  std::vector<double> x_ref;
  std::vector<double> penalty_dx;
  std::vector<std::pair<double, double>> s_dot_bounds;
  const SpeedLimit& speed_limit = st_graph_data.speed_limit();
  for (int i = 0; i < num_of_knots; ++i) {
    double curr_t = i * delta_t;
    // get path_s
    SpeedPoint sp;
    reference_speed_data.EvaluateByTime(curr_t, &sp);
    const double path_s = sp.s();
    x_ref.emplace_back(path_s);
    // get curvature
    PathPoint path_point = path_data.GetPathPointWithPathS(path_s);
    penalty_dx.push_back(std::fabs(path_point.kappa()) *
                         config.kappa_penalty_weight());
    // get v_upper_bound
    const double v_lower_bound = 0.0;
    double v_upper_bound = FLAGS_planning_upper_speed_limit;
    v_upper_bound = speed_limit.GetSpeedLimitByS(path_s);
    s_dot_bounds.emplace_back(v_lower_bound, std::fmax(v_upper_bound, 0.0));
  }
  piecewise_jerk_problem.set_x_ref(config.ref_s_weight(), std::move(x_ref));
  piecewise_jerk_problem.set_penalty_dx(penalty_dx);
  piecewise_jerk_problem.set_dx_bounds(std::move(s_dot_bounds));
  • 求解QP问题
      if (!piecewise_jerk_problem.Optimize()) {
    const std::string msg = "Piecewise jerk speed optimizer failed!";
    AERROR << msg;
    speed_data->clear();
    return Status(ErrorCode::PLANNING_ERROR, msg);
  }
  • 输出结果,如果速度曲线的总时长小于FLAGS_fallback_total_time = 3s,则使用SpeedProfileGenerator::FillEnoughSpeedPoints(speed_data)在后面补齐速度全部为0的点。
  // Extract output
  const std::vector<double>& s = piecewise_jerk_problem.opt_x();
  const std::vector<double>& ds = piecewise_jerk_problem.opt_dx();
  const std::vector<double>& dds = piecewise_jerk_problem.opt_ddx();
  for (int i = 0; i < num_of_knots; ++i) {
    ADEBUG << "For t[" << i * delta_t << "], s = " << s[i] << ", v = " << ds[i]
           << ", a = " << dds[i];
  }
  speed_data->clear();
  speed_data->AppendSpeedPoint(s[0], 0.0, ds[0], dds[0], 0.0);
  for (int i = 1; i < num_of_knots; ++i) {
    // Avoid the very last points when already stopped
    if (ds[i] <= 0.0) {
      break;
    }
    speed_data->AppendSpeedPoint(s[i], delta_t * i, ds[i], dds[i],
                                 (dds[i] - dds[i - 1]) / delta_t);
  }
  SpeedProfileGenerator::FillEnoughSpeedPoints(speed_data);
  RecordDebugInfo(*speed_data, st_graph_data.mutable_st_graph_debug());
  return Status::OK();
  • 3
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值