% Matlab Example
%https://se.mathworks.com/help/mpc/gs/control-of-a-multi-input-single-output-plant.html

% Define plant model
% one measured output, one manipulated variable, one measured disturbance, and one unmeasured disturbance
sys = ss(tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}));
Ts = 0.2;
model = c2d(sys,Ts);
%Design MPC controller
% the first signal is a manipulated variable, the second signal is a measured disturbance, the third one is an unmeasured disturbance.
model = setmpcsignals(model,'MV',1,'MD',2,'UD',3);
P = 10;     % prediction horizon
M = 3;      % control horizon
mpcobj = mpc(model,Ts,P,M); % M is either an integer (<= P) or a vector of blocked moves such that sum(M) <= P.
% Constraints for MV
mpcobj.MV = struct('Min',-2,'Max',2,'RateMin',-10,'RateMax',10);
% Weights for MV
mpcobj.Weights.ManipulatedVariables = [0];

%For unmeasured input disturbances, its model is an integrator driven by white noise with variance = 1000.
mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);
% Simulation time and steps
Tstop = 60;                                 % simulation time
Tf = round(Tstop/Ts);                       % number of simulation steps

%Reference and measured disturbance
r = ones(Tf,1);                           % reference signal vector
v = [zeros(2*Tf/3,1);ones(3*Tf/3,1)];       % measured disturbance signal

% Add unmeasured disturbance
SimOptions = mpcsimopt(mpcobj);
d = [zeros(2*Tf/3,1);-0.5*ones(Tf/3,1)];
SimOptions.Unmeas = d;                          % unmeasured input disturbance signal
SimOptions.OutputNoise=.001*(rand(Tf,1)-.5);    % output measurement noise
SimOptions.InputNoise=.05*(rand(Tf,1)-.5);      % noise on manipulated variables

%Simulate
[y,t,u,xp] = sim(mpcobj,Tf,r,v);             %output, time, MV, states (initial model was a tf, so states have no physical  meaning

%Plot
figure
subplot(2,1,1)
plot(0:Tf-1,y,0:Tf-1,r)
title('Output')
grid
subplot(2,1,2)
plot(0:Tf-1,u)
title('Input')
grid