Code covered by the BSD License

# MPC Tutorial III: MPC in Simulink V2

by

### Yi Cao (view profile)

27 Jan 2009 (Updated )

A tutorial on using MPC in Simulink.

## Contents

% Based on the State Space MPC code in MATLAB Central, a Simulink block of
% the State Space MPC is developed using sfunction. Two examples to use
% this MPC block, one with a continous time state space plant and another
% with a discrete space state plant are included in the tutorial. These two
% plant can be replaced by a corresponding nonlinear model in
% continous-time or discrete-time.

State Space MPC code.

## Discrete-time plant model

The Simulink model for discrete-time state-space plant is as follows.

% The plant model
A=[ 0.1555  -13.7665   -0.0604         0         0         0
0.0010    1.0008    0.0068         0         0         0
0    0.0374    0.9232         0         0         0
0.0015   -0.1024   -0.0003    0.1587  -13.6705   -0.0506
0    0.0061         0    0.0006    0.9929    0.0057
0    0.0001         0         0    0.0366    0.9398];
B=[0.0001       0
0       0
-0.0036       0
0  0.0001
0       0
0 -0.0028];
Bd=[      0         0
0         0
0.0013         0
0         0
0         0
0    0.0008];
C=[0 362.995 0 0 0 0
0 0 0 0 362.995 0];
D=zeros(2,2);
% MPC internal model
Am=A;
Bm=B;
Cm=C;
Dm=D;
% Prediction horizon and moving horizon
P=10;
M=3;
% Performance wights
Q=1.5*eye(2*P);
R=eye(2*M);
Ts=1;
x0=zeros(6,1);
d0=[0.5;0.5];
dv=[2;2];
ds=[0;1001];
N=1500;
% Predefined reference
yr=zeros(N,2);
yr(10:N,1)=1;
yr(351:N,1)=3;
yr(600:N,1)=5;
yr(1100:N,1)=3;
yr(210:N,2)=2;
yr(551:N,2)=1;
yr(800:N,2)=4;
yr(1300:N,2)=2;
d=(rand(N,2)-0.5)*2;
t=(0:N-1)';
ny=2;
nd=2;

[T,X,YU]=sim('mpcdss',[0 N-1],[],[t yr d]);
Y=YU(:,1:2);
U=YU(:,3:4);
t=t*0.1;
subplot(211)
plot(t,Y,t,yr,'--','linewidth',2)
title('output and setpoint')
ylabel('temp, C^\circ')
legend('T_1','T_2','T_1 Ref','T_2 Ref', 'location','southeast')
subplot(212)
stairs(t,U,'linewidth',2)
legend('u_1','u_2','location','southeast')
title('input')
ylabel('flow rate, m^3/s')
xlabel('time, s')

## Continuous-time plant model

The Simulink model for continous-time plant model is as follows.

% The continous-time state space model
A=[-17.9527 -295.5969    0.0034   -0.0000   -0.0000   -0.0000
0.0215    0.1960    0.0706   -0.0000    0.0000   -0.0000
-0.0005    0.3832   -0.8005    0.0000    0.0000    0.0000
0.0913   -0.0804   -0.0044  -18.0146 -295.9326    0.0016
-0.0001    0.0601   -0.0002    0.0130    0.0424    0.0589
0.0000   -0.0001    0.0000   -0.0003    0.3755   -0.6220];
B=[    0.0028   -0.0000
0.0001   -0.0000
-0.0375    0.0000
-0.0000    0.0025
-0.0000    0.0001
0.0000   -0.0289];
C=[      0  362.9950         0         0         0         0
0         0         0         0  362.9950         0];
D=zeros(2);
Bd=[-0.0002    0.0000
-0.0000    0.0000
0.0135   -0.0000
0.0000   -0.0001
0.0000   -0.0000
-0.0000    0.0083];

% Discretized model for MPC
Ts=0.1;
sysd = c2d(ss(A,B,C,D),Ts);
Am=sysd.A;
Bm=sysd.B;
Cm=sysd.C;
Dm=sysd.D;

% Prediction horizon and moving horizon
P=10;
M=3;
% Performance wights, change these to see the effect
Q=1.5*eye(2*P);
R=eye(2*M);
Ts=0.2;
x0=zeros(6,1);
d0=[0.5;0.5];
dv=[2;2];
ds=[0;1001];
N=1500;
% Predefined reference
yr=zeros(N,2);
yr(10:N,1)=1;
yr(351:N,1)=3;
yr(600:N,1)=5;
yr(1100:N,1)=3;
yr(210:N,2)=2;
yr(551:N,2)=1;
yr(800:N,2)=4;
yr(1300:N,2)=2;
d=(rand(N,2)-0.5)*2;
t=Ts*(0:N-1)';
ny=2;
nd=2;

[t,X,YU]=sim('mpcss',t,[],[t yr d]);
Y=YU(:,1:2);
U=YU(:,3:4);
subplot(211)
plot(t,Y,t,yr,'--','linewidth',2)
axis([0 t(end) -6 6])
title('output and setpoint')
ylabel('temp, C^\circ')
legend('T_1','T_2','T_2 Ref','T_2 Ref','location','southeast')
subplot(212)
stairs(t,U,'linewidth',2)
legend('u_1','u_2','location','southeast')
title('input')
ylabel('flow rate, m^3/s')
xlabel('time, s')