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.l≤Ax≤u
其中,
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元素个数为矩阵列数加1,P_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=0∑n−2[wssi2+(wv+pi)s˙i2−Δt22wjerks¨is¨i+1]+i=1∑n−2(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)sn−12+(wv+pn−1+wvend)s˙n−12+(wa+Δt2wjerk)s¨02+(wa+Δt2wjerk+waend)s¨n−12
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=0∑n−1[wssrefisi+(wv+pi)vrefis˙i]
注意:对于终点处的
s
n
−
1
s_{n-1}
sn−1额外增加了惩罚系数
w
s
e
n
d
w_{s_{end}}
wsend,对于终点处的
s
˙
n
−
1
\dot s_{n-1}
s˙n−1额外增加了惩罚系数
w
v
e
n
d
w_{v_{end}}
wvend,对于终点处的
s
¨
n
−
1
\ddot s_{n-1}
s¨n−1额外增加了惩罚系数
w
a
e
n
d
w_{a_{end}}
waend。
预备知识:
- 矩阵与列向量相乘,等于矩阵中的每一列的线性组合,组合系数为列向量中的每个元素:
[
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
s0s1⋮sn−1s˙0s˙1⋮s˙n−1s¨0s¨1⋮s¨n−1
TP
2An∗n0002Bn∗n0002Cn∗n
s0s1⋮sn−1s˙0s˙1⋮s˙n−1s¨0s¨1⋮s¨n−1
+q
−2wssref0−2wssref1⋮−2wssrefn−1−2(wv+p0)vref0−2(wv+p1)vref1⋮−2(wv+pn−1)vrefn−1
T
s0s1⋮sn−1s˙0s˙1⋮s˙n−1
接下来我们写出矩阵
A
n
∗
n
A_{n*n}
An∗n、
B
n
∗
n
B_{n*n}
Bn∗n和
C
n
∗
n
C_{n*n}
Cn∗n的具体表达式(空余部分元素全部为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]
An∗n=
wsws⋱wsws+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] Bn∗n= wv+p0wv+p1⋱wv+pn−2wv+pn−1+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] Cn∗n= wa+Δt2wjerk−Δt22wjerkwa+Δt22wjerk−Δt22wjerk⋱⋱wa+Δt22wjerk−Δt22wjerkwa+Δt22wjerk−Δt22wjerkwa+Δt2wjerk+waend
注意:上述矩阵
C
n
∗
n
C_{n*n}
Cn∗n不是对称的,按照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]
Cn∗n=
wa+Δt2wjerk−Δt2wjerk−Δt2wjerkwa+Δt22wjerk⋱−Δt2wjerk⋱−Δt2wjerk⋱wa+Δ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]
Cn∗n=
wa+Δt2wjerk−Δt2wjerkwa+Δt22wjerk⋱−Δt2wjerk⋱⋱wa+Δt22wjerk−Δt2wjerkwa+Δt2wjerk+waend
内核矩阵
P
P
P中的系数就是PiecewiseJerkSpeedProblem
(在文件modules/planning/math/piecewise_jerk/piecewise_jerk_speed_problem.cc
中)的CalculateKernel
函数中的系数,偏移向量q
中的系数就是PiecewiseJerkSpeedProblem
的CalculateOffset
函数中的系数(注意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]≤si≤sub[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˙i≤s˙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¨i≤s¨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Δt≤s¨i+1−s¨i≤jerkubΔ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 l≤Ax≤u形式。
- 初始条件约束
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[n−1]s˙lb[0]⋮s˙lb[n−1]s¨lb[0]⋮s¨lb[n−1]jerklbΔt⋮jerklbΔt0.0⋮0.0sinitvinitainit
6n∗1≤A6n∗3n
s0⋮sn−1s˙0⋮s˙n−1s¨0s¨1⋮s¨n−1
3n∗1≤
sub[0]⋮sub[n−1]s˙ub[0]⋮s˙ub[n−1]s¨ub[0]⋮s¨ub[n−1]jerkubΔt⋮jerkubΔt0.0⋮0.0sinitvinitainit
6n∗1
**仿射矩阵
A
6
n
∗
3
n
A_{6n*3n}
A6n∗3n**的表达式如下所示(空白部分的元素为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=
1−1111⋱0⋱⋱−1⋯1n∗n1(n−1)∗n01∗n1−1−Δt111⋱⋱0⋱⋱−1−Δt(n−1)∗n⋯1n∗n1(n−1)∗n01∗n1−1−2Δt−3Δt2111⋱−2Δt⋱−6Δt2⋱0⋱⋱−1⋱−2Δt⋱−3Δt2⋯1n∗n1(n−1)∗n−2Δt(n−1)∗n−6Δt2(n−1)∗n01∗n
6n∗3n
注意:上式中的下标
1
n
∗
n
1_{n*n}
1n∗n表示某块矩阵的行列数是
n
∗
n
n*n
n∗n。仿射矩阵
A
A
A中的系数也就是PiecewiseJerkProblem
的CalculateAffineConstraint
函数(在文件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();