MATLAB Examples

MPC in Simulink Tutorial

```% 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') ```