自动驾驶-控制算法(三)
本文主要记录模型预测控制(MPC)的原理,以及基于MPC实现自动驾驶车辆路径跟踪。
前言
本文主要记录模型预测控制(MPC)的原理,以及基于MPC实现自动驾驶车辆路径跟踪。
一、模型预测控制
1.1 基本概念
模型预测控制(MPC)是一种控制方案,它使用一个模型来预测系统在有限时间窗口(视界)内的未来行为。根据这些预测和系统当前的测量/估计状态,计算出与确定的控制目标相关的最佳控制输入,并遵守系统约束条件。经过一定的时间间隔后,测量、估计和计算过程将以移动的视界重复进行。
MPC的主要优势:
- 主动控制行动: 控制器预测未来的干扰、设定点等;
- 非线性控制: MPC 可明确考虑非线性系统,而无需线性化;
- 任意控制目标: 传统的设定点跟踪和调节或经济型 MPC;
- 受约束的表述: 明确考虑物理、安全或运行系统约束。
1.2 整体流程
1 预测区间与控制区间
对于一般的离散系统,在 k k k时刻,可以测量出系统当前状态 g ( k ) g(k) g(k),通过计算可得到 u ( k ) , u ( k + 1 ) , u ( k + 2 ) , . . . , u ( k + j ) u(k), u(k+1), u(k+2),..., u(k+j) u(k),u(k+1),u(k+2),...,u(k+j),并得到系统未来状态的估计值 y ( k ) , y ( k + 1 ) , y ( k + 2 ) , . . . , y ( k + j ) y(k), y(k+1), y(k+2),..., y(k+j) y(k),y(k+1),y(k+2),...,y(k+j)。
预测区间为一次优化后预测未来输出的时间步个数,如下图中的 [ k , k + P ] [k, k+P] [k,k+P]区间;
控制区间为进行控制估计的部分,如下图中 [ k , k + M − 1 ] [k, k+M-1] [k,k+M−1]区间。
过小的控制区间,可能无法做到较好的控制;而较大的控制区间,会导致只有控制范围的前一部分才会有较好的效果,后一部分的效果甚微,且带来大量的计算开销。
2 约束
约束包括Hard约束和Soft约束,Hard约束是不可违背必须遵守的,如车轮转角;Soft约束是可以违反的(会带来一定的惩罚代价),如道路限速。
在控制系统中,输入输出都可能会有约束,但是在设计时不建议将输入输出都设为Hard约束,因为这两部分的约束有可能存在重叠,导致优化器产生不可行解。
建议输出采用Soft约束,输入中的输入参数、输入参数变化率建议设置一个Hard约束和一个Soft约束。
1.3 MPC设计
当模型时线性的时候,MPC的设计求解一般使用二次规划方法。
设线性模型为以下形式:
x k + 1 = A x k + B u k + C \begin{align} x_{k+1} = Ax_k + Bu_k +C \end{align} xk+1=Axk+Buk+C
假定未来 m m m步的控制输入已知,为 u ( k ) , u ( k + 1 ) , u ( k + 2 ) , . . . , u ( k + m ) u(k), u(k+1), u(k+2),..., u(k+m) u(k),u(k+1),u(k+2),...,u(k+m),根据以上模型,可以计算未来 m m m步的状态:
x k + 1 = A x k + B u k + C x_{k+1} = Ax_k + Bu_k +C xk+1=Axk+Buk+C
x k + 2 = A x k + 1 + B u k + 1 + C = A ( A x k + B u k + C ) + B u k + 1 + C = A 2 x k + A B u k + B u k + 1 + A C + C x_{k+2} = Ax_{k+1} + Bu_{k+1} +C=A(Ax_k + Bu_k +C) + Bu_{k+1} +C = A^2x_k+ABu_k+Bu_{k+1}+AC+C xk+2=Axk+1+Buk+1+C=A(Axk+Buk+C)+Buk+1+C=A2xk+ABuk+Buk+1+AC+C
x k + 3 = A x k + 2 + B u k + 2 + C = A 3 x k + A 2 B u k + A B u k + 1 + B u k + 2 + A 2 C + A C + C x_{k+3} = Ax_{k+2} + Bu_{k+2} +C=A^3x_k+A^2Bu_k +ABu_{k+1} +Bu_{k+2}+A^2C+AC+C xk+3=Axk+2+Buk+2+C=A3xk+A2Buk+ABuk+1+Buk+2+A2C+AC+C
. . . ... ...
x k + m = A m x k + A m − 1 B u k + . . . A m − i B u k + i − 1 + . . . + B u k + m − 1 + A m − 1 C + A m − 2 C + . . . + C x_{k+m}=A^mx_k+A^{m-1}Bu_k+...A^{m-i}Bu_{k+i-1}+...+Bu_{k+m-1}+A^{m-1}C+A^{m-2}C+...+C xk+m=Amxk+Am−1Buk+...Am−iBuk+i−1+...+Buk+m−1+Am−1C+Am−2C+...+C
将上述方程组写成矩阵向量形式可得:
X = A x k + B u + C \begin{align} X=\mathcal{A}x_k + \mathcal{B} \mathbf{u} +\mathcal{C} \end{align} X=Axk+Bu+C
其中:
X = [ x k + 1 x k + 2 x k + 3 ⋮ x k + m ] ; u = [ u k u k + 1 u k + 2 ⋮ u k + m − 1 ] ; A = [ A A 2 A 3 ⋮ A m ] X = \begin{bmatrix} x_{k+1} \\ x_{k+2} \\ x_{k+3} \\ \vdots \\ x_{k+m} \end{bmatrix};\ \ \mathbf{u} = \begin{bmatrix} u_{k} \\ u_{k+1} \\ u_{k+2} \\ \vdots \\ u_{k+m-1} \end{bmatrix}; \ \ \mathcal{A} = \begin{bmatrix} A \\ A^2 \\ A^3 \\ \vdots \\ A^m \end{bmatrix} X=
xk+1xk+2xk+3⋮xk+m
; u=
ukuk+1uk+2⋮uk+m−1
; A=
AA2A3⋮Am
B = [ B 0 0 . . . 0 A B B 0 . . . 0 A 2 B A B B . . . 0 ⋮ ⋮ ⋮ ⋱ ⋮ A m − 1 B A m − 2 B A m − 3 B . . . B ] \mathcal{B} = \begin{bmatrix} B & 0 & 0 & ... & 0 \\ AB & B & 0 & ... & 0 \\ A^2B & AB & B & ... & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ A^{m-1}B & A^{m-2}B & A^{m-3}B & ... & B \end{bmatrix} B= BABA2B⋮Am−1B0BAB⋮Am−2B00B⋮Am−3B.........⋱...000⋮B
C = [ C A C + C A 2 C + A C + C . . . A m − 1 C + . . . + C ] \mathcal{C} = \begin{bmatrix} C \\ AC+C \\ A^2C+AC+C \\ ... \\ A^{m-1}C+...+ C \end{bmatrix} C= CAC+CA2C+AC+C...Am−1C+...+C
上式 B \mathcal{B} B中的下三角形式,直接反应了系统在时间上的因果关系,即 k + 1 k+1 k+1时刻的输入对 k k k时刻的输出没有影响, k + 2 k+2 k+2时刻的输入对 k k k和 k + 1 k+1 k+1时刻没有影响。
假定参考轨迹为 X ˉ = [ x ˉ k + 1 , x ˉ k + 2 , . . . , x ˉ k + m ] T \bar{X}=[\bar{x}_{k+1}, \bar{x}_{k+2},...,\bar{x}_{k+m}]^T Xˉ=[xˉk+1,xˉk+2,...,xˉk+m]T,则MPC的一个简单的目标代价函数如下:
min J = E T Q E + u T R u s . t . u m i n ≤ u ≤ u m a x \begin{align} \min \mathcal{J} = \mathcal{E} ^TQ \mathcal{E} +\mathbf{u}^TR\mathbf{u} \\ s.t. \ u_{min} \le \mathbf{u} \le u_{max} \end{align} minJ=ETQE+uTRus.t. umin≤u≤umax
其中, E = X − X ˉ = [ x k + 1 − x ˉ k + 1 , x k + 2 − x ˉ k + 2 , . . . , x k + m − x ˉ k + m ] T \mathcal{E} = X - \bar{X}=[x_{k+1} - \bar{x}_{k+1}, x_{k+2} -\bar{x}_{k+2},...,x_{k+m} -\bar{x}_{k+m}]^T E=X−Xˉ=[xk+1−xˉk+1,xk+2−xˉk+2,...,xk+m−xˉk+m]T
以上最优化问题可通过二次规划求解,得到满足目标代价函数最小的最优控制序列 u = [ u k u k + 1 u k + 2 . . . u k + m − 1 ] T \mathbf{u} = \begin{bmatrix} u_{k} & u_{k+1} & u_{k+2} & ... & u_{k+m-1} \end{bmatrix} ^T u=[ukuk+1uk+2...uk+m−1]T
二、基于MPC的自动驾驶车辆轨迹跟踪
基于运动学模型的离散状态空间方程如下:
X ~ ( k + 1 ) = A ~ X ~ ( k ) + B ~ u ~ ( k ) \begin{align} \tilde{X}(k+1) = \tilde{A}\tilde{X}(k)+\tilde{B}\tilde{\mathrm{u}}(k) \end{align} X~(k+1)=A~X~(k)+B~u~(k)
其中:
X ~ ( k ) = [ x ( k ) − x r y ( k ) − y r ψ ( k ) − ψ r ] \tilde{X}(k)=\begin{bmatrix} x(k) - x_r \\ y(k) - y_r \\ \psi(k) - \psi_r \end{bmatrix} X~(k)=
x(k)−xry(k)−yrψ(k)−ψr
u ~ ( k ) = [ v ( k ) − v r δ ( k ) − δ r ] \tilde{\mathrm{u}}(k)=\begin{bmatrix} v(k) - v_r \\ \delta(k) - \delta_r \end{bmatrix} u~(k)=[v(k)−vrδ(k)−δr]
A ~ = A ⋅ T + I = [ 1 0 − v r ⋅ T ⋅ sin ψ r 0 1 v r ⋅ T ⋅ cos ψ r 0 0 1 ] \tilde{A}=A\cdot T+I=\begin{bmatrix} 1& 0& -v_r \cdot T \cdot \sin \psi _r \\ 0& 1& v_r \cdot T \cdot \cos \psi _r \\ 0& 0& 1 \end{bmatrix} A~=A⋅T+I=
100010−vr⋅T⋅sinψrvr⋅T⋅cosψr1
B ~ = B ⋅ T = [ T cos ψ r 0 T sin ψ r 0 T ⋅ tan δ r L v r ⋅ T L ⋅ cos 2 δ r ] \tilde{B}=B\cdot T=\begin{bmatrix} T\cos \psi _r & 0 \\ T\sin \psi _r & 0 \\ \frac{T \cdot \tan \delta _r}{L} & \frac{v_r \cdot T}{L\cdot \cos ^2 \delta _r} \end{bmatrix} B~=B⋅T= TcosψrTsinψrLT⋅tanδr00L⋅cos2δrvr⋅T
式中, X r = [ x r , y r , ψ r ] X_r=[x_r, y_r, \psi_r] Xr=[xr,yr,ψr]为参考点处的状态, u r = [ v r , δ r ] \mathrm{u}_r=[v_r, \delta_r] ur=[vr,δr]为参考点处的控制量; T T T为采样步长, I I I为单位矩阵,维度与矩阵 A A A一致。
将上述状态方程进行改写可得:
X ( k + 1 ) = A ~ X ( k ) + B ~ u ( k ) + X r − A ~ X r − B ~ u r \begin{align} X(k+1) = \tilde{A}X(k)+\tilde{B}\mathrm{u}(k)+X_r-\tilde{A}X_r-\tilde{B}\mathrm{u}_r \end{align} X(k+1)=A~X(k)+B~u(k)+Xr−A~Xr−B~ur
其中 X ( k ) = [ x ( k ) , y ( k ) , ψ ( k ) ] X(k)=[x(k), y(k), \psi(k)] X(k)=[x(k),y(k),ψ(k)], u ( k ) = [ v ( k ) , δ ( k ) ] \mathrm{u}(k)=[v(k), \delta(k)] u(k)=[v(k),δ(k)]表示 k k k点处的状态量和控制量。
令 C ~ = X r − A ~ X r − B ~ u r \tilde{C}=X_r-\tilde{A}X_r-\tilde{B}\mathrm{u}_r C~=Xr−A~Xr−B~ur,可得:
X ( k + 1 ) = A ~ X ( k ) + B ~ u ( k ) + C ~ \begin{align} X(k+1) = \tilde{A}X(k)+\tilde{B}\mathrm{u}(k)+\tilde{C} \end{align} X(k+1)=A~X(k)+B~u(k)+C~
MPC控制的代价函数定义如下:
J = ( X ( N ) − X ˉ ( N ) ) T Q f ( X ( N ) − X ˉ ( N ) ) + ∑ k = 0 N − 1 ( X ( k ) − X ˉ ( k ) ) T Q ( X ( k ) − X ˉ ( k ) ) + u ( k ) T R u ( k ) \begin{align} J = \left (X(N) - \bar{X}(N) \right)^TQ_f \left (X(N) - \bar{X}(N) \right) + \sum_{k=0}^{N-1} \left (X(k) - \bar{X}(k) \right) ^T Q \left(X(k) - \bar{X}(k)\right)+\mathrm{u}(k)^TR\mathrm{u}(k) \end{align} J=(X(N)−Xˉ(N))TQf(X(N)−Xˉ(N))+k=0∑N−1(X(k)−Xˉ(k))TQ(X(k)−Xˉ(k))+u(k)TRu(k)
s u b j e c t t o : X ( k + 1 ) = A ~ X ( k ) + B ~ u ( k ) + C ~ u m i n ≤ u ( k ) ≤ u m a x X ( 0 ) = X 0 subject \ to: \ X(k+1) = \tilde{A}X(k)+\tilde{B}\mathrm{u}(k)+\tilde{C}\\ \mathrm{u}_{min} \le \mathrm{u}(k) \le \mathrm{u}_{max} \\ X(0)=X_0 subject to: X(k+1)=A~X(k)+B~u(k)+C~umin≤u(k)≤umaxX(0)=X0
其中,矩阵 Q Q Q, Q f Q_f Qf, R R R为正定对称矩阵, X ( k ) , u ( k ) X(k), \mathrm{u}(k) X(k),u(k)分别表示状态量和控制量,其中控制量受 [ u m i n , u m a x ] [\mathrm{u}_{min}, \mathrm{u}_{max}] [umin,umax]约束。 X 0 X_0 X0表示初始状态, X ( N ) X(N) X(N)表示终端状态。 X ˉ ( k ) \bar{X}(k) Xˉ(k)表示参考路径在 k k k时刻的状态。
求解上述MPC问题需要将其改写为二次规划的形式,然后用求解器进行求解。
二次规划的标准形式为:
m i n i m i z e 1 2 x T P x + q T x s u b j e c t t o l ≤ A c x ≤ u \begin{align} \mathrm{minimize} \ \frac{1}{2} x^TPx +q^Tx \\ subject \ to \ l\le A_c x \le u \end{align} minimize 21xTPx+qTxsubject to l≤Acx≤u
其中hessian矩阵 P P P为:
P = d i a g ( Q , Q , . . . , Q f , R , . . . , R ) \begin{align} P=diag(Q, Q, ..., Q_f, R, ..., R) \end{align} P=diag(Q,Q,...,Qf,R,...,R)
矩阵维度为: P ∈ R ( 3 ∗ ( N + 1 ) + 2 ∗ N ) × ( 3 ∗ ( N + 1 ) + 2 ∗ N ) P \in R^{(3*(N+1)+2*N)\times (3*(N+1)+2*N)} P∈R(3∗(N+1)+2∗N)×(3∗(N+1)+2∗N)。(状态变量个数为3,控制变量个数为2)
梯度向量 q q q为:
q = [ − Q X ˉ ( 0 ) − Q X ˉ ( 1 ) ⋮ − Q X ˉ ( N ) 0 0 ⋮ 0 ] \begin{align} q=\begin{bmatrix} -Q\bar{X}(0) \\ -Q\bar{X}(1) \\ \vdots \\ -Q\bar{X}(N) \\ 0\\ 0\\ \vdots \\ 0 \end{bmatrix} \end{align} q=
−QXˉ(0)−QXˉ(1)⋮−QXˉ(N)00⋮0
矩阵维度为: Q ∈ R ( 3 ∗ ( N + 1 ) + 2 ∗ N ) × 1 Q \in R^{(3*(N+1)+2*N) \times 1} Q∈R(3∗(N+1)+2∗N)×1
优化变量 x x x为:
x = [ X ( 0 ) X ( 1 ) ⋮ X ( N ) u ( 0 ) u ( 1 ) ⋮ u ( N − 1 ) ] \begin{align} x=\begin{bmatrix} X(0) \\ X(1) \\ \vdots \\ X(N) \\ \mathrm{u}(0) \\ \mathrm{u}(1) \\ \vdots \\ \mathrm{u}(N-1) \end{bmatrix} \end{align} x=
X(0)X(1)⋮X(N)u(0)u(1)⋮u(N−1)
矩阵维度为: x ∈ R ( 3 ∗ ( N + 1 ) + 2 ∗ N ) × 1 x \in R^{(3*(N+1)+2*N) \times 1} x∈R(3∗(N+1)+2∗N)×1
线性约束矩阵 A c A_c Ac为:
A c = [ − I 0 0 ⋯ 0 0 0 ⋯ 0 A − I 0 ⋯ 0 B 0 ⋯ 0 0 A − I ⋯ 0 0 B ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ − I 0 0 ⋯ B 0 0 0 ⋯ 0 I 0 ⋯ 0 0 0 0 ⋯ 0 0 I ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ 0 0 0 ⋯ I ] \begin{align} A_c=\begin{bmatrix} -I & 0 & 0 & \cdots & 0& 0& 0& \cdots& 0 \\ A & -I & 0 & \cdots & 0& B& 0& \cdots& 0 \\ 0 & A & -I & \cdots & 0 & 0& B& \cdots& 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots& \vdots& \vdots& \ddots& \vdots \\ 0 & 0 & 0 & \cdots & -I & 0& 0& \cdots& B \\ 0 & 0 & 0 & \cdots & 0 & I& 0& \cdots& 0 \\ 0 & 0 & 0 & \cdots & 0 & 0& I& \cdots& 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots& \vdots& \vdots& \ddots& \vdots \\ 0 & 0 & 0 & \cdots & 0 & 0& 0& \cdots& I \\ \end{bmatrix} \end{align} Ac=
−IA0⋮000⋮00−IA⋮000⋮000−I⋮000⋮0⋯⋯⋯⋱⋯⋯⋯⋱⋯000⋮−I00⋮00B0⋮0I0⋮000B⋮00I⋮0⋯⋯⋯⋱⋯⋯⋯⋱⋯000⋮B00⋮I
矩阵维度为: A c ∈ R ( 3 ∗ ( N + 1 ) + 2 ∗ N ) × ( 3 ∗ ( N + 1 ) + 2 ∗ N ) A_c \in R^{(3*(N+1)+2*N) \times (3*(N+1)+2*N)} Ac∈R(3∗(N+1)+2∗N)×(3∗(N+1)+2∗N)
约束条件的上界 u u u 和下界 l l l 分别为:
l = [ − X ( 0 ) − C ~ ⋮ − C ~ u m i n u m i n ⋮ u m i n ] u = [ − X ( 0 ) − C ~ ⋮ − C ~ u m a x u m a x ⋮ u m a x ] \begin{align} l=\begin{bmatrix} -X(0) \\ -\tilde{C} \\ \vdots \\ -\tilde{C} \\ \mathrm{u}_{min} \\ \mathrm{u}_{min} \\ \vdots \\ \mathrm{u}_{min} \end{bmatrix} \qquad u=\begin{bmatrix} -X(0) \\ -\tilde{C} \\ \vdots \\ -\tilde{C} \\ \mathrm{u}_{max} \\ \mathrm{u}_{max} \\ \vdots \\ \mathrm{u}_{max} \end{bmatrix} \end{align} l=
−X(0)−C~⋮−C~uminumin⋮umin
u=
−X(0)−C~⋮−C~umaxumax⋮umax
矩阵维度为: l ∈ R ( 3 ∗ ( N + 1 ) + 2 ∗ N ) × 1 l \in R^{(3*(N+1)+2*N) \times 1} l∈R(3∗(N+1)+2∗N)×1, u ∈ R ( 3 ∗ ( N + 1 ) + 2 ∗ N ) × 1 u \in R^{(3*(N+1)+2*N) \times 1} u∈R(3∗(N+1)+2∗N)×1
参考文献
1、Basics of model predictive control
2、An Introduction to Model-based PredictiveControl (MPC)
3、【自动驾驶】模型预测控制(MPC)实现轨迹跟踪
4、Using osqp-eigen in MPC fashion
更多推荐




所有评论(0)