本代碼基于無人駕駛車輛模型預測控制第一版本第五章中的代碼予以糾正，下面代碼已經過測試。
通過函數來生成參考軌迹（雙移線），然後用模型預測控制器去跟蹤。
需要全套資源和課程學習（5套simulink模型+10份代碼+5篇參考論文）的朋友加我微信（Jeossirey-LSJ）擷取。
代碼段：
function [sys,x0,str,ts] = chapter5_2_2(t,x,u,flag)
switch flag,
case 0
[sys,x0,str,ts] = mdlInitializeSizes; % Initialization
case 2
sys = mdlUpdates(t,x,u); % Update discrete states
case 3
sys = mdlOutputs(t,x,u); % Calculate outputs
% case 4
% sys = mdlGetTimeOfNextVarHit(t,x,u); % Get next sample time
case {1,4,9} % Unused flags
sys = [];
otherwise
error(['unhandled flag = ',num2str(flag)]); % Error handling
end
% End of dsfunc.
%==============================================================
% Initialization
%==============================================================
function [sys,x0,str,ts] = mdlInitializeSizes
% Call simsizes for a sizes structure, fill it in, and convert it
% to a sizes array.
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 6;
sizes.NumOutputs = 1;
sizes.NumInputs = 8;
sizes.DirFeedthrough = 1; % Matrix D is non-empty.
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 =[0.001;0.0001;0.0001;0.00001;0.00001;0.00001];
global U;%U為我們的控制量
U=[0];%控制量初始化,這裡面加了一個期望軌迹的輸出，如果去掉，U為一維的
% Initialize the discrete states.
str = []; % Set str to an empty matrix.
ts = [0.05 0]; % sample time: [period, offset]，采樣時間影響最大。
% ts = [0.1 0]; % sample time: [period, offset]，采樣時間影響最大。
%End of mdlInitializeSizes
%==============================================================
% Update the discrete states
%==============================================================
function sys = mdlUpdates(t,x,u)
sys = x;
%End of mdlUpdate.
%==============================================================
% Calculate outputs
%==============================================================
function sys = mdlOutputs(t,x,u)
global a b;
global U;
tic
Nx=6;
Nu=1;
Ny=2;
Np =20;
Nc=5;
Row=1000;%松弛因子權重
fprintf('Update start, t=%6.3f\n',t)
%以下為我們的狀态量誤差矩陣，由于我們的初始值給的很小，這裡予以忽略
y_dot=u(1)/3.6; %橫向速度化為m/s
x_dot=u(2)/3.6+0.0001;%CarSim輸出的是km/h，轉換為m/s
phi=u(3)*3.141592654/180; %CarSim輸出的為角度，角度轉換為弧度，程式設計時我們都是使用弧度
phi_dot=u(4)*3.141592654/180;
Y=u(5);%機關為m
X=u(6);%機關為米
Y_dot=u(7);
X_dot=u(8);
%% 車輛參數輸入
%syms sf為前輪滑移率，sr為後輪滑移率
Sf=0.2; Sr=0.2;
%syms lf%前輪距離車輛質心的距離，lr為後輪距離車輛質心的距離
lf=1.232;lr=1.468;
%syms C_cf前輪線性橫向側偏剛度； C_cr後輪線性橫向側偏剛度 ；C_lf 前輪縱向側偏剛度； C_lr 後輪縱向側偏剛度
Ccf=66900;Ccr=62700;Clf=66900;Clr=62700;
%syms m g I;%m為車輛品質，g為重力加速度，I為車輛繞Z軸的轉動慣量，車輛固有參數
m=1723;g=9.8;I=4175;
%% 參考軌迹生成
shape=2.4;%參數名稱，用于參考軌迹生成
dx1=25;dx2=21.95;%沒有任何實際意義，隻是參數名稱
dy1=4.05;dy2=5.7;%沒有任何實際意義，隻是參數名稱
Xs1=27.19;Xs2=56.46;%參數名稱，以上參數是為了生成我們的雙移線
X_predict=zeros(Np,1);%用于儲存預測時域内的縱向位置資訊，這是計算期望軌迹的基礎
phi_ref=zeros(Np,1);%用于儲存預測時域内的參考橫擺角資訊
Y_ref=zeros(Np,1);%用于儲存預測時域内的參考橫向位置資訊
% % 參考軌迹圖像
% shape=2.4;%參數名稱，用于參考軌迹生成
% dx1=25;dx2=21.95;%沒有任何實際意義，隻是參數名稱
% dy1=4.05;dy2=5.7;%沒有任何實際意義，隻是參數名稱
% Xs1=27.19;Xs2=56.46;%參數名稱
%畫出參考軌迹的圖像
% X_phi=1:1:220;%這個點的區間是根據縱向速度（x_dot）來定的，如果速度為10m/s則區間=10*0.1=1
% z1=shape/dx1*(X_phi-Xs1)-shape/2;
% z2=shape/dx2*(X_phi-Xs2)-shape/2;
% Y_ref=dy1/2.*(1+tanh(z1))-dy2/2.*(1+tanh(z2));
% figure(1);
% plot(X_phi, Y_ref,'r--','LineWidth',2);
% hold on;
% 以下計算kesi,即狀态量與控制量合在一起
kesi=zeros(Nx+Nu,1);
kesi(1)=y_dot;
kesi(2)=x_dot;
kesi(3)=phi;
kesi(4)=phi_dot;
kesi(5)=Y;
kesi(6)=X;
kesi(7)=U(1);
delta_f=U(1);
fprintf('Update start, u(1)=%4.2f\n',U(1))
T=0.05;%仿真步長
T_all=20;%總的仿真時間，主要功能是防止計算期望軌迹越界
%權重矩陣設定
Q_cell=cell(Np,Np);
for i=1:1:Np;
for j=1:1:Np;
if i==j
%Q_cell{i,j}=[200 0;0 100;];
Q_cell{i,j}=[2000 0;0 10000;];
else
Q_cell{i,j}=zeros(Ny,Ny);
end
end
end
R=5*10^5*eye(Nu*Nc);
%最基本也最重要的矩陣，是控制器的基礎，采用動力學模型，該矩陣與車輛參數密切相關，通過對動力學方程求解雅克比矩陣得到，a為6*6，b為6*1
a=[ 1 - (259200*T)/(1723*x_dot), -T*(phi_dot + (2*((460218*phi_dot)/5 - 62700*y_dot))/(1723*x_dot^2) - (133800*((154*phi_dot)/125 + y_dot))/(1723*x_dot^2)), 0, -T*(x_dot - 96228/(8615*x_dot)), 0, 0
T*(phi_dot - (133800*delta_f)/(1723*x_dot)), (133800*T*delta_f*((154*phi_dot)/125 + y_dot))/(1723*x_dot^2) + 1, 0, T*(y_dot - (824208*delta_f)/(8615*x_dot)), 0, 0
0, 0, 1, T, 0, 0
(33063689036759*T)/(7172595384320*x_dot), T*(((2321344006605451863*phi_dot)/8589934592000 - (6325188028897689*y_dot)/34359738368)/(4175*x_dot^2) + (5663914248162509*((154*phi_dot)/125 + y_dot))/(143451907686400*x_dot^2)), 0, 1 - (813165919007900927*T)/(7172595384320000*x_dot), 0, 0
T*cos(phi), T*sin(phi), T*(x_dot*cos(phi) - y_dot*sin(phi)), 0, 1, 0
-T*sin(phi), T*cos(phi), -T*(y_dot*cos(phi) + x_dot*sin(phi)), 0, 0, 1];
b=[ 133800*T/1723
T*((267600*delta_f)/1723 - (133800*((154*phi_dot)/125 + y_dot))/(1723*x_dot))
0
5663914248162509*T/143451907686400
0
0];
d_k=zeros(Nx,1);
state_k1=zeros(Nx,1);
%以下即為根據離散非線性模型預測下一時刻狀态量
%注意，為避免前後軸距的表達式（a,b）與控制器的a,b矩陣沖突，将前後軸距的表達式改為lf和lr
state_k1(1,1)=y_dot+T*(-x_dot*phi_dot+2*(Ccf*(delta_f-(y_dot+lf*phi_dot)/x_dot)+Ccr*(lr*phi_dot-y_dot)/x_dot)/m);
state_k1(2,1)=x_dot+T*(y_dot*phi_dot+2*(Clf*Sf+Clr*Sr+Ccf*delta_f*(delta_f-(y_dot+phi_dot*lf)/x_dot))/m);
state_k1(3,1)=phi+T*phi_dot;
state_k1(4,1)=phi_dot+T*((2*lf*Ccf*(delta_f-(y_dot+lf*phi_dot)/x_dot)-2*lr*Ccr*(lr*phi_dot-y_dot)/x_dot)/I);
state_k1(5,1)=Y+T*(x_dot*sin(phi)+y_dot*cos(phi));
state_k1(6,1)=X+T*(x_dot*cos(phi)-y_dot*sin(phi));
d_k=state_k1-a*kesi(1:6,1)-b*kesi(7,1);
d_piao_k=zeros(Nx+Nu,1)
d_piao_k(1:6,1)=d_k;
d_piao_k(7,1)=0;
A_cell=cell(2,2);
B_cell=cell(2,1);
A_cell{1,1}=a;
A_cell{1,2}=b;
A_cell{2,1}=zeros(Nu,Nx);
A_cell{2,2}=eye(Nu);
B_cell{1,1}=b;
B_cell{2,1}=eye(Nu);
A=cell2mat(A_cell);
B=cell2mat(B_cell);
C=[0 0 1 0 0 0 0;0 0 0 0 1 0 0;];
PSI_cell=cell(Np,1);
THETA_cell=cell(Np,Nc);
GAMMA_cell=cell(Np,Np);%次元 20*20
PHI_cell=cell(Np,1);%次元 20*1
for p=1:1:Np;
PHI_cell{p,1}=d_piao_k;
for q=1:1:Np;
if q<=p; %下三角矩陣
GAMMA_cell{p,q}=C*A^(p-q);
else
GAMMA_cell{p,q}=zeros(Ny,Nx+Nu);
end
end
end
for j=1:1:Np
PSI_cell{j,1}=C*A^j;
for k=1:1:Nc
if k<=j
THETA_cell{j,k}=C*A^(j-k)*B;
else
THETA_cell{j,k}=zeros(Ny,Nu);
end
end
end
PSI=cell2mat(PSI_cell);
THETA=cell2mat(THETA_cell);
GAMMA=cell2mat(GAMMA_cell);
PHI=cell2mat(PHI_cell);
Q=cell2mat(Q_cell);
H_cell=cell(2,2);
H_cell{1,1}=THETA'*Q*THETA+R;
H_cell{1,2}=zeros(Nu*Nc,1);
H_cell{2,1}=zeros(1,Nu*Nc);
H_cell{2,2}=Row;
H=cell2mat(H_cell);
H=(H+H')/2;
error_1=zeros(Ny*Np,1);%40*7*7*1=40*1
Yita_ref_cell=cell(Np,1);%參考的元胞數組為20*1
for p=1:1:Np
if t+p*T>T_all
X_DOT=x_dot*cos(phi)-y_dot*sin(phi);
X_predict(Np,1)=X+X_DOT*Np*T;
z1=shape/dx1*(X_predict(Np,1)-Xs1)-shape/2;
z2=shape/dx2*(X_predict(Np,1)-Xs2)-shape/2;
Y_ref(p,1)=dy1/2*(1+tanh(z1))-dy2/2*(1+tanh(z2));
phi_ref(p,1)=atan(dy1*(1/cosh(z1))^2*(1.2/dx1)-dy2*(1/cosh(z2))^2*(1.2/dx2));
Yita_ref_cell{p,1}=[phi_ref(p,1);Y_ref(p,1)];
else
X_DOT=x_dot*cos(phi)-y_dot*sin(phi);
X_predict(p,1)=X+X_DOT*p*T;
z1=shape/dx1*(X_predict(p,1)-Xs1)-shape/2;
z2=shape/dx2*(X_predict(p,1)-Xs2)-shape/2;
Y_ref(p,1)=dy1/2*(1+tanh(z1))-dy2/2*(1+tanh(z2));
phi_ref(p,1)=atan(dy1*(1/cosh(z1))^2*(1.2/dx1)-dy2*(1/cosh(z2))^2*(1.2/dx2));
Yita_ref_cell{p,1}=[phi_ref(p,1);Y_ref(p,1)];
end
end
Yita_ref=cell2mat(Yita_ref_cell);%将我們得到的輸出元胞矩陣轉換為矩陣
error_1=Yita_ref-PSI*kesi-GAMMA*PHI; %求偏差
f_cell=cell(1,2);
f_cell{1,1}=2*error_1'*Q*THETA;
f_cell{1,2}=0;
% f=(cell2mat(f_cell))';
f=-cell2mat(f_cell);
%% 以下為限制生成區域
%控制量限制
A_t=zeros(Nc,Nc);%見falcone論文 P181
for p=1:1:Nc
for q=1:1:Nc
if q<=p %下三角矩陣包含對角線
A_t(p,q)=1;
else
A_t(p,q)=0;
end
end
end
A_I=kron(A_t,eye(Nu));%求克羅内克積
Ut=kron(ones(Nc,1),U(1));
umin=-0.1744;%維數與控制變量的個數相同，前輪偏角的上限制
umax=0.1744;%前輪偏角的下限制
delta_umin=-0.0148*0.4;%前輪偏角變化量的下限制
delta_umax=0.0148*0.4;%前輪偏角變化量的上限制
Umin=kron(ones(Nc,1),umin);
Umax=kron(ones(Nc,1),umax);
%輸出量限制
ycmax=[0.21;5]; %橫擺角和縱向位移的限制
ycmin=[-0.3;-3];
Ycmax=kron(ones(Np,1),ycmax);
Ycmin=kron(ones(Np,1),ycmin);
%結合控制量限制和輸出量限制
A_cons_cell={A_I zeros(Nu*Nc,1);-A_I zeros(Nu*Nc,1)};
b_cons_cell={Umax-Ut;-Umin+Ut};
% A_cons_cell={A_I zeros(Nu*Nc,1);-A_I zeros(Nu*Nc,1);THETA zeros(Ny*Np,1);-THETA zeros(Ny*Np,1)};
% b_cons_cell={Umax-Ut;-Umin+Ut;Ycmax-PSI*kesi-GAMMA*PHI;-Ycmin+PSI*kesi+GAMMA*PHI};
A_cons=cell2mat(A_cons_cell);
b_cons=cell2mat(b_cons_cell);
%控制增量限制
M=10;
delta_Umin=kron(ones(Nc,1),delta_umin);
delta_Umax=kron(ones(Nc,1),delta_umax);
lb=[delta_Umin;0];
ub=[delta_Umax;M];
%% 開始求解過程
options = optimset('Algorithm','interior-point-convex');
x_start=zeros(Nc+1,1);%加入一個起始點
[X,fval,exitflag]=quadprog(H,f,A_cons,b_cons,[],[],lb,ub,x_start,options);
fprintf('exitflag=%d\n',exitflag);
fprintf('H=%4.2f\n',H(1,1));
fprintf('f=%4.2f\n',f(1,1));
%% 計算輸出
u_piao=X(1);
U(1)=kesi(7,1)+u_piao;
sys= U;
toc
% End of mdlOutputs.
跟蹤結果：
