0.友情链接

• B站链接1-北京理工大学无人驾驶技术课程
• B站链接2-MATLAB实现模型预测控制
• B站链接3-北京理工大学的学生讲解无人驾驶车辆的模型预测控制
• B站链接4-英国谢菲尔德大学模型预测控制网课（全英文）
• IEEE论文参考：Model Predictive Control for Visual Servo Steering of Nonholonomic Mobile Robots（非完整移动机器人视觉伺服转向的模型预测控制）
• 书籍参考：国之重器出版社《无人驾驶车辆模型预测控制（第二版）》
• 书籍配套代码下载链接
• 本文完整配套代码链接
• 阿克曼转向原理介绍

1.引言

\qquadMATLAB自身是带预测控制模块的，但是看了一些MATLAB的预测模块使用教程视频，仍然不如自己书写代码原理清晰，而且MATLAB的预测控制是基于状态方程的滑模控制，并不包含所有预测控制的种类，使用范围有限。因此我想按照预测控制的原理，逐步为大家讲解预测控制算法及MATLAB代码的书写步骤。
\qquadSimulink自身也携带一个关于无人驾驶的模块Car_sim，但是Car_sim使用的模型比一般的模型更复杂，因此本人参考《无人驾驶车辆的模型预测控制》中的
\qquad预测控制都包含三个主要步骤——预测模型、滚动优化和反馈矫正。按照模型的不同，可以分为动态矩阵算法、广义极点配置、基于状态方程的滑模控制。现代控制理论中，以状态方程为主要模型，本文也以其为例建立模型，其他的模型原理与其类似，只是求解方法有所差异。
\qquad为了方便大家理解，本人不直接列出所有代码，而是按照讲解和代码对照的方式，逐步展现每个环节的代码，最终展现全部代码和效果。那我们就从第一步开始吧。

2.预测模型

\qquad建立预测模型的目的在于获取输入量（控制量）和输出量（以及其他干扰项）的映射关系。在本例（汽车变道）示例中，控制量是前轮角度$\varphi$和汽车速度大小$v$，而输出量则是汽车的质心与待变道直线的距离。控制的目的是让汽车的质心尽快接近待变道直线，并维持其在该直线上运动。
\qquad我们的思路已经很清晰了，但是别着急，我们似乎忽略了什么。例如汽车的速度不能是无穷大，前轮的角度也是有限制的。我们取车辆的后轴中心速度大小限制为$[v_{min},v_{max}]$，前轮角度绝对值限幅为${\varphi }_m$。前轮转动的这个角度，又称为车辆的横摆角。变道的情况如下图所示：

\qquad汽车是前轮驱动还是后轮驱动直接决定了控制量是前轮的速度还是后轮速度，当然前后轮速度是可以通过后轴中心速度和车辆横摆角转换得到的。因此设小汽车后轴速度为$v_d$，车辆横摆角为$\varphi$，同时设车长为$l$。以上图左下角为原点建立坐标系。（不知道阿克曼转向原理的可以看一下友情链接的介绍，这里用的是简化的两轮模型）

\qquad黑色虚线代表了汽车朝向，而红色虚线则是前轮的朝向。二者夹角为横摆角$\varphi$，以逆时针为正方向。汽车可以看成是二维平面的刚体，因此只需要三个参数就可以确定其状态。我们选取它质心的横纵坐标$x,y$加汽车朝向与x轴逆时针方向的夹角$\theta$作为状态变量。（在一篇IEEE的文章上看到选取极坐标作为状态变量的方法，也是可行的，但对于变道这种单移线运动来说，直角坐标更加方便）。采用$S$表示状态向量，$U$表示控制向量。
\qquad汽车的后轮速度大小$v_d$和汽车的朝向$\theta$决定了一小段时间内汽车的位移$(\delta x,\delta y)$，而前轮的转向角$\varphi$则决定了一小段时间汽车朝向的变化量$\delta \theta$。状态方程列写如下：
$$S^{\prime }=\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ {\theta }^{\prime }\end{array}\right]=\left[\begin{array}{c}{v}_{d}cos\theta \\ v_{d}sin\theta \\ v_{d}tan\varphi /l\end{array}\right]=f(S,U)$$
\qquad模拟阿克曼转向的车辆的MATLAB代码如下：

function [X,Y,Th] = car_run(X,Y,Th,L,dt,vd,phi)X = X + vd*cos(Th)*dt;Y = Y + vd*sin(Th)*dt;Th = Th + vd*tan(phi)/L*dt;
end

X和Y为车辆后轴中心的位置信息，Th为车辆的朝向，L为车长，dt为仿真采样间隔，vd为车辆后轴速度，phi为车辆横摆角。

3.滚动优化

3.1.线性化

\qquad滚动优化的目的是根据预测模型和控制目的，建立目标函数和约束条件，并在约束条件内，求得目标函数的最优解（即期望控制量），从而达到预测控制的目的。
\qquad求目标函数极小值的方法有很多，智能优化算法依赖于蒙特卡洛模拟，数值优化算法依赖于函数的导数。前者比后者的实用性更大，但是速度更慢。汽车行驶中的在线规划对运算速度的要求高，因此普遍采用后者进行求函数极小值的操作。
\qquad二次规划是数值优化算法中一类特殊的优化问题，即特指目标函数为二次型$\frac{1}{2}{x}^{T}Hx+c^{T}x$，约束条件为线性约束的（线性不等式约束＋线性灯等式约束）的有约束优化问题。将预测模型的目标函数转换为二次规划问题的模型预测控制（MPC）又称为线性预测控制（LPC），其运算速度是MPC中最快的。本文将以无人驾驶领域为例，介绍模型的线性化过程。
\qquad汽车在变道时，给定期望的轨迹$(x_r,y_r)$，并通过控制变量使汽车在时域上经过的路径逼近期望轨迹。当然，我们期望汽车的朝向$\theta$最终逼近于0（即直线行驶），期望的状态量使用$S_r$表示，则我们的目的是使$|S-S_r|$即$(S-S_r{)}^2/2$最小。当期望状态向量与实际状态向量相差较小且期望控制量与实际控制量相差较小时，可采用在$S_r$处的导数$dy/dx$来代替变化率$\Delta y/\Delta x$：
$$\begin{array}{l}{S}^{\prime }-S_r^{\prime }=f(S,U)-f(S_r,U_r)\\ \approx \frac{\partial f}{\partial S}{|}_{S_r}\cdot (S-S_r)+\frac{\partial f}{\partial U}{|}_{U_r}\cdot (U-U_r)\\ =A\cdot (S-S_r)+B\cdot (U-U_r)\end{array}$$
其中
$$\begin{gathered} A=\frac{\partial f}{\partial S}{|}_{S_r}=\left[\begin{array}{ccc}\frac{\partial f_1}{\partial S_1}& \frac{\partial f_1}{\partial S_2}& \frac{\partial f_1}{\partial S_3}\\ \frac{\partial f_2}{\partial S_1}& \frac{\partial f_2}{\partial S_2}& \frac{\partial f_2}{\partial S_3}\\ \frac{\partial f_3}{\partial S_1}& \frac{\partial f_3}{\partial S_2}& \frac{\partial f_3}{\partial S_3}\end{array}\right]=\left[\begin{array}{ccc}1& 0& -v_{dr}\ast sin({\theta }_r)\\ 0& 1& v_{dr}\ast cos({\theta }_r)\\ 0& 0& 1\end{array}\right] \\ B=\frac{\partial f}{\partial U}{|}_{U_r}=\left[\begin{array}{cc}\frac{\partial f_1}{\partial U_1}& \frac{\partial f_1}{\partial U_2}\\ \frac{\partial f_2}{\partial U_1}& \frac{\partial f_2}{\partial U_2}\\ \frac{\partial f_3}{\partial U_1}& \frac{\partial f_3}{\partial U_2}\end{array}\right]=\left[\begin{array}{cc}cos({\theta }_r)& 0\\ sin({\theta }_r)& 0\\ tan({\varphi }_r)/l& v_{dr}\ast se{c}^2({\varphi }_r)/l\end{array}\right] \end{gathered}$$
记$\widetilde{S^{\prime }}=S^{\prime }-S_r^{\prime },\tilde{S}=S-S_r,\tilde{U}=U-U_r$
则线性化后的模型为
$$\widetilde{S^{\prime }}=A\tilde{S}+B\tilde{U}$$
即我们所熟知的线性状态方程。我们的目标是使得汽车按照既定的轨迹形式，因此希望实际状态逼近于期望状态，即$||\tilde{S}||\to 0$，因此输出方程即为
$$Y=\tilde{S}$$

3.2.$U_r$的求取

\qquad期望状态向量$S_r$是容易理解的，无人驾驶车辆需要提前规划形式的路径$(x_r,y_r)$，有时还对车辆的朝向有所要求（${\theta }_r$），那么期望控制量$U_r$又是怎么得出的呢？其实很简单，在没有约束条件的情况下，根据状态方程即可得到$U_r$的取值。根据状态方程我们可以得知：
$$\{\begin{array}{l}{x}_r^{\prime }=v_{dr}cos({\theta }_r)\\ y_r^{\prime }=v_{dr}sin({\theta }_r)\\ {\theta }_r^{\prime }=v_{dr}tan{\varphi }_r/l\end{array}\Rightarrow \{\begin{array}{l}{v}_{dr}=\sqrt{(x_r^{\prime }{)}^2+(y_r^{\prime }{)}^2}\\ {\varphi }_r=atan(l\cdot {\theta }_r^{\prime }/v_{dr})\end{array}$$
\qquad在实际系统中，若$S=(x_r,y_r,{\theta }_r)$是离散取值的，可以借助三次样条插值的方法补全需要的值，并用差分替代微分求取状态量的导数。样条差值的方法在此不再赘述，仿真时建议使用MATLAB的pchip函数。
本部分对应的MATLAB代码如下

X_r = [0,1,2,3,4,5,6]*3;  % X的期望坐标
Y_r = [3,3,2,1,0,0,0];  % Y的期望坐标
L=2.5; % 车身长度
v0 = 10; % 期望速度
T = 20e-3; % 控制周期20ms
dt = 1e-3; % 仿真周期1ms
Xr = linspace(X_r(1),X_r(numel(X_r)),(X_r(numel(X_r))-X_r(1))/(T*v0));  % 横坐标期望值的插值序列
Yr = pchip(X_r,Y_r,Xr);  % 通过三次样条插值算出纵坐标坐标插值序列
n = numel(Yr);  % 预测模型中所有规划点
FYr = Yr(2:n);
dYr = zeros(1,n);
dYr(1:n-1) = FYr - Yr(1:n-1);
dYr(n) = dYr(n-1);%保证导数连续
FXr = Xr(2:n);
dXr = zeros(1,n);
dXr(1:n-1) = FXr - Xr(1:n-1);
dXr(n) = dXr(n-1);
Thr = atan(dYr./dXr);  % 车朝向角的期望值
FThr = Thr(2:n);
dThr = zeros(1,n);
dThr(1:n-1) = FThr-Thr(1:n-1);
dThr(n)=dThr(n-1);
vdr = sqrt(dXr.^2+dYr.^2)/T;  % 车后轮的期望速度
Phr = atan(L./vdr.*dThr/T);  % 车前轮的偏转角

这里以变道为例，绘制期望轨迹$(X_r,Y_r)$如下图

\qquad可以发现，车辆从上道（y=3）逐步变化至下道(y=0)。样条插值后的轨迹并不是折线，但由于不一定满足车辆的物理约束，因此不能用期望控制量作为实际控制量进行控制。

3.3.离散化与序列化

\qquad预测控制依赖于计算机，计算机的控制信号是离散的，因此我们还需要把线性化后的状态方程离散化。假设控制周期为$T$，则离散化后的状态方程为：
$$\begin{array}{l}\frac{\tilde{S}(k+1)-\tilde{S}(k)}{T}=A_k\tilde{S}(k)+B_k\tilde{U}(k)\\ \Rightarrow \tilde{S}(k+1)=(A_{k}T+I)\tilde{S}(k)+B_{k}T\tilde{U}(k)\\ \Rightarrow \tilde{S}(k+1)=A_k^{\ast }\tilde{S}(k)+B_k^{\ast }\tilde{U}(k)\end{array}$$
\qquad接着我们需要选择接下来需要预测的未来若干个采样时刻的状态量的个数和控制量个数，即预测时域$N$和控制时域$M$（$M\le N$）。我们采取最简单的情形$M=N$讲解。
\qquad首先列写未来$N$个时刻的离散状态方程
$$\begin{gathered} \{\begin{array}{l}\tilde{S}(k+1)=A_k^{\ast }\tilde{S}(k)+B_k^{\ast }\tilde{U}(k)\\ \tilde{S}(k+2)=A_{k+1}^{\ast }\tilde{S}(k+1)+B_{k+1}^{\ast }\tilde{U}(k+1)\\ \cdots \\ \tilde{S}(k+N)=A_{k+N-1}^{\ast }\tilde{S}(k+N-1)+B_{k+N-1}^{\ast }\tilde{U}(k+N-1)\end{array} \\ \Rightarrow \tilde{S}(k+i)=A{A}_i\tilde{S}(k)+A{B}_1\tilde{U}(k)+A{B}_2\tilde{U}(k+1)+...+A{B}_i\tilde{U}(k+i-1) \end{gathered}$$
其中
$$\{\begin{array}{l}A{A}_i=A_{k+i-1}^{\ast }{A}_{k+i-2}^{\ast }...A_k^{\ast }\\ A{B}_1^{(i)}=A_{k+i-1}^{\ast }{A}_{k+i-2}^{\ast }...A_{k+1}^{\ast }{B}_k^{\ast }\\ A{B}_2^{(i)}=A_{k+i-1}^{\ast }{A}_{k+i-2}^{\ast }...A_{k+2}^{\ast }{B}_{k+1}^{\ast }\\ \cdots \\ A{B}_{i-1}^{(i)}=A_{k+i-1}^{\ast }{B}_{k+i-2}^{\ast }\\ A{B}_i^{(i)}=B_{k+i-1}^{\ast }\end{array}$$
\qquad已知量为$\tilde{S}(k),A_k^{\ast },B_k^{\ast }$，后两者是因为期望状态量和期望控制量已知因此也已知，但由于期望状态量是一个时间序列，因此矩阵$A_k^{\ast },B_k^{\ast }$也是一个时间序列。输入为$\tilde{U}(k),\tilde{U}(k+1),...,\tilde{U}(k+N)$，输出为$\tilde{S}(k+1),\tilde{S}(k+2),....,\tilde{S}(k+N)$，记
$$\begin{gathered} \hat{\tilde{S}}=[\tilde{S}(k+1),\tilde{S}(k+2),....,\tilde{S}(k+N){]}^T \\ \hat{\tilde{U}}=[\tilde{U}(k),\tilde{U}(k+1),...,\tilde{U}(k+N-1){]}^T \\ \hat{U}=[U(k),U(k+1),...,U(k+N-1){]}^T \end{gathered}$$
根据上面的推导，可以推出序列化后的状态方程：
$$\begin{gathered} \hat{\tilde{S}}=AA\cdot \tilde{S}(k)+BB\cdot \hat{\tilde{U}} \\ s.t.\{\begin{array}{l}{v}_{min}\le v_d\le v_{max}\\ -{\varphi }_m\le \varphi \le {\varphi }_m\end{array} \end{gathered}$$
其中
$$\begin{gathered} AA=\left[\begin{array}{c}A{A}_1\\ A{A}_2\\ ⋮\\ A{A}_N\end{array}\right] \\ BB=\left[\begin{array}{cccc}A{B}_1^{(1)}& 0& \cdots & 0\\ A{B}_1^{(2)}& A{B}_2^{(2)}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ A{B}_1^{(N)}& A{B}_2^{(N)}& \cdots & A{B}_N^{(N)}\end{array}\right] \end{gathered}$$
优化问题为
$$\hat{U}=arg\underset{\hat{U}}{\min }\frac{1}{2}{\hat{\tilde{S}}}^{T}Q\hat{\tilde{S}}+\frac{1}{2}{\hat{U}}^{T}P\hat{U}$$

3.4.实现增量控制

\qquad除了保证车辆基本的约束外，为了保证车辆的舒适性，还应防止控制量的突变，即车辆的加速度和打方向盘的速度不至于过快。控制量的加速度不宜过大，需要在约束条件中增加关于控制量变化率的硬约束和在目标函数增加关于控制量变化率的软约束。
\qquad硬约束很好理解，即对控制增量进行限幅
$$\{\begin{array}{l}{a}_{min}\le v_d^{\prime }\le a_{max}\\ -d{\varphi }_m\le {\varphi }^{\prime }\le d{\varphi }_m\end{array}$$
\qquad软约束则是通过在目标函数中加入关于控制增量的二次型进行约束，目的是使得目标函数在求取极小值的过程中兼顾到控制增量的最小化。因此上述3.3.节的目标函数则需要修改为
$$\Delta \hat{U}=arg\underset{\Delta \hat{U}}{\min }\frac{1}{2}{\hat{\tilde{S}}}^{T}Q\hat{\tilde{S}}+\frac{1}{2}(\Delta \hat{U}{)}^{T}P(\Delta \hat{U})$$
\qquad仅此而已是不够的，我们还需要重写线性状态方程，记
$$\{\begin{array}{l}{G}_u=[u(k-1)-u_r(k),u(k-1)-u_r(k+1),...,u(k-1)-u_r(k+N-1){]}^T\\ E_N=tril(ones(1,N))\text{(元素均为1的下三角阵)}\\ \hat{\tilde{U}}=G_u+T\ast E_N\cdot \Delta \hat{U}\end{array}$$
因此目标函数为
$$\begin{array}{rl}J=& \frac{1}{2}{\hat{\tilde{S}}}^{T}Q\hat{\tilde{S}}+\frac{1}{2}(\Delta \hat{U}{)}^{T}P(\Delta \hat{U})\\ =& \frac{1}{2}(AA\cdot \tilde{S}(k)+BB\cdot \hat{\tilde{U}}{)}^{T}Q(AA\cdot \tilde{S}(k)+BB\cdot \hat{\tilde{U}})+\frac{1}{2}(\Delta \hat{U}{)}^{T}P(\Delta \hat{U})\\ =& \frac{1}{2}(AA\cdot \tilde{S}(k)+BB(G_u+T\ast E_N\Delta \hat{U}){)}^{T}Q(AA\cdot \tilde{S}(k)+BB(G_u+T\ast E_N\Delta \hat{U}))+\frac{1}{2}(\Delta \hat{U}{)}^{T}P(\Delta \hat{U})\end{array}$$
令
$$\begin{gathered} G_{AB}=AA\cdot \tilde{S}(k)+BB\cdot G_u \\ {E}_{NB}=T\ast BB\cdot E_N \end{gathered}$$
则
$$\begin{array}{rl}J=& \frac{1}{2}(G_{AB}+E_{NB}\cdot \Delta \hat{U}{)}^{T}Q(G_{AB}+E_{NB}\cdot \Delta \hat{U})+\frac{1}{2}(\Delta \hat{U}{)}^{T}P(\Delta \hat{U})\\ =& \frac{1}{2}{G}_{AB}^{T}Q{G}_{AB}+G_{AB}^{T}Q{E}_{NB}\cdot \Delta \hat{U}+\frac{1}{2}(\Delta \hat{U}{)}^T(E_{NB}^{T}Q{E}_{NB}+P)(\Delta \hat{U})\\ ⟹& \Delta \hat{U}=arg\underset{\Delta \hat{U}}{\min }J(\Delta \hat{U})\\ & s.t.\{\begin{array}{l}{a}_{min}\le v_d^{\prime }\le a_{max}\\ -d{\varphi }_m\le {\varphi }^{\prime }\le d{\varphi }_m\end{array}\end{array}$$
\qquad此处参考了友情链接中IEEE的论文的写法，与原书中的写法不一样。这种方法推导出的关于控制增量的目标函数更加直接，也不需要再修改状态向量了。所有MPC的原理讲解部分就到此为止了。
\qquad完整的MPC函数的MATLAB代码如下：

    function u = MPC(Thr_list,ur_list,pSk,u0,L,N,T)u_max = [20;pi/6];u_min = [5;-pi/6];du_max = [2/T;25*pi/180/T];du_min = -du_max;A = cell(N,1);B = cell(N,1);for i=1:NThr = Thr_list(i);   ur = ur_list(:,i);vdr = ur(1);Phr = ur(2);A{i} = [1 0 -vdr*sin(Thr)*T;0 1 vdr*cos(Thr)*T;0 0     1];B{i} = [cos(Thr)*T  0;sin(Thr)*T  0;tan(Phr)*T/L vdr*T*sec(Phr)^2/L];endAA = zeros(3*N,3);  % 序列化之后的ABB = zeros(3*N,2*N);  % 序列化之后的B% S = AA*Sk + BB*u for i = 1:NAA(3*(i-1)+1:3*i,:) = eye(3);for j = 1:iAA(3*(i-1)+1:3*i,:) = A{j}*AA(3*(i-1)+1:3*i,:);endendfor i = 1:N  % 列for j = i:N  % 行BB(3*(j-1)+1:3*j,2*(i-1)+1:2*i) = B{i}; % 3×2for t = i+1:jBB(3*(j-1)+1:3*j,2*(i-1)+1:2*i) = A{t}*BB(3*(j-1)+1:3*j,2*(i-1)+1:2*i);endendendGu = cell(N,1);U0 = cell(N,1);for i = 1:NGu{i} = u0 - ur_list(:,i);U0{i} = u0;endGu = cell2mat(Gu);  %2N*2U0 = cell2mat(U0);  %2N*1EN = zeros(2*N,2*N);for i = 1:Nfor j = 1:iEN(2*(i-1)+1,2*(j-1)+1)=T;EN(2*i,2*j)=T;endendENB = BB*EN;GAB = AA*pSk + BB*(Gu); % (3N×3)×(3×1)+(3N×2N)×(2N×1)→3N×1       Q = zeros(3*N,3*N);P = zeros(2*N,2*N);Qu = [3,3,1]; % 各状态量的加权Pu = [1,0.2];  % 各控制量的加权for i = 1:NQ(3*(i-1)+1:3*i,3*(i-1)+1:3*i) = diag(Qu*1e5*0.9^(i-1));P(2*(i-1)+1:2*i,2*(i-1)+1:2*i) = diag(Pu*0.9^(i-1));endAb = [EN;-EN];  %4N*2Numax = cell(N,1);umax(:)={u_max};umax = cell2mat(umax);  % 2N×1umin = cell(N,1);umin(:)={u_min};umin = cell2mat(umin);  % 2N×1dumax = cell(N,1);dumax(:)={du_max};dumax =cell2mat(dumax);  % 2N×1dumin = cell(N,1);dumin(:)={du_min};dumin =cell2mat(dumin);  % 2N×1b = cat(1,umax-U0,-umin+U0);  % 4N×1H = ENB'*Q*ENB + P;  % 2N×2NH = 1/2*(H+H');  % 使其严格对称c = (GAB'*Q*ENB)'; % 2N×1opt = optimoptions('quadprog','Display','off');du0 = zeros(2*N,1);du0(1:2) = (ur_list(:,1)-u0)/T;for i = 1:N-1du0(2*i+1:2*i+2) = (ur_list(:,i+1)-ur_list(:,i))/T;enddu = quadprog(H,c,Ab,b,[],[],dumin,dumax,du0,opt);u = u0 +du(1:2)*T;  end

Thr_list是车辆朝向角$\theta$的期望值序列，ur_list是车辆的控制量期望值序列，pSk是$S(k)-S_r(k)$，u0是$u(k-1)$，N是预测时域，L是车长，T是控制周期。

4.仿真示例

车辆变道的轨迹图、后轴中心速度图及横摆角图如下

\qquad可以发现车辆的变道轨迹与期望轨迹基本一致，横摆角成功限制在了-30°~30°之间，速度限幅也符合要求（这里限制上限为20，下限为5）。
\qquad希望本文对您有帮助，谢谢阅读。本文的完整代码请见第0节友情链接处。