Contents

Control of a Multi-Input Single-Output Plant

This example shows how to design model predictive controller with one measured output, one manipulated variable, one measured disturbance, and one unmeasured disturbance in a typical workflow.

Define Plant Model

The discrete time linear open-loop dynamic model is defined below.

sys = ss(tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}));
Ts = 0.2;               % sampling time
model = c2d(sys,Ts);

Design MPC Controller

Define type of input signals: 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);
% Create the controller object with sampling period, prediction and control
% horizons:
mpcobj = mpc(model,Ts,10,3);
-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000.
-->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000.
-->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.

Define constraints on the manipulated variable.

mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10);

For unmeasured input disturbances, its model is an integrator driven by white noise with variance = 1000.

mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);

Simulate Closed-Loop Response Using the SIM Command

Specify simulation parameters.

Tstop = 30;                               % simulation time
Tf = round(Tstop/Ts);                     % number of simulation steps
r = ones(Tf,1);                           % reference signal
v = [zeros(Tf/3,1);ones(2*Tf/3,1)];       % measured disturbance signal
% Run the closed-loop simulation and plot results.
sim(mpcobj,Tf,r,v);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Specify disturbance and noise signals in simulation option object.

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

Run the closed-loop simulation and save the results to workspace.

[y,t,u,xp] = sim(mpcobj,Tf,r,v,SimOptions);

Plot results.

figure
subplot(211)
plot(0:Tf-1,y,0:Tf-1,r)
title('Output');
grid
subplot(212)
plot(0:Tf-1,u)
title('Input');
grid

Simulate Closed-Loop Response with Model Mismatch

Test the robustness of the MPC controller against a model mismatch. Assume the true plant is below:

simModel = ss(tf({1,1,1},{[1 .8 1],[1 2],[.6 .6 1]}));
simModel = setmpcsignals(simModel,'MV',1,'MD',2,'UD',3);
simModel = struct('Plant',simModel);
simModel.Nominal.Y = 0.1; % The nominal value of the output of the true plant is 0.1
simModel.Nominal.X = -.1*[1 1 1 1 1];

Create option object.

SimOptions.Model = simModel;
SimOptions.plantinit = [0.1 0 -0.1 0 .05];  % Initial state of the true plant
SimOptions.OutputNoise = [];                % remove output measurement noise
SimOptions.InputNoise = [];                 % remove noise on manipulated variables
% Run the closed-loop simulation and plot results.
sim(mpcobj,Tf,r,v,SimOptions);
-->Converting model to discrete time.

Soften Constraints

Relax the constraints on manipulated variables from hard to soft.

mpcobj.MV.MinECR = 1;
mpcobj.MV.MaxECR = 1;

Define an output constraint.

mpcobj.OV.Max = 1.1;

Run a new closed-loop simulation.

sim(mpcobj,Tf,r,v);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

MV constraint has been slightly violated while MO constraint has been violated more. Penalize more output constraint and rerun the simulation.

mpcobj.OV.MaxECR = 0.001;  % The closer to zero, the harder the constraint
sim(mpcobj,Tf,r,v);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Now MO constraint has been slightly violated while MV constraint has been violated more as expected.

Change Kalman Gains Used in the Built-In State Estimator

Model Predictive Control Toolbox™ provides a default Kalman filter to estimate the state of plant, disturbance, and noise models. You can change the Kalman gains.

First retrieve the default Kalman gains and state-space matrices.

[L,M,A1,Cm1] = getEstimator(mpcobj);

The default observer poles are:

e = eig(A1-A1*M*Cm1);
fprintf('\nDefault observer poles: [%s]\n',sprintf('%5.4f ',e));
Default observer poles: [0.5708 0.5708 0.9334 0.9334 0.4967 0.8189 ]

Design a new state estimator by pole-placement.

poles = [.8 .75 .7 .85 .6 .81];
L = place(A1',Cm1',poles)';
M = A1\L;
setEstimator(mpcobj,L,M);

Simulate Open-Loop Response

Test the behavior of plant in open-loop using the SIM command. Set the 'OpenLoop' flag on, and provide the sequence of manipulated variables that excite the system.

SimOptions.OpenLoop = 'on';
SimOptions.MVSignal = sin((0:Tf-1)'/10);

As the reference signal will be ignored, we can avoid specifying it.

sim(mpcobj,Tf,[],v,SimOptions);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
-->Converting model to discrete time.

Examine Steady-State Offset

To examine whether the MPC controller will be able to reject constant output disturbances and/or track constant set-point with zero offsets in steady-state, compute the DC gain from output disturbances to controlled outputs using the CLOFFSET command.

DC = cloffset(mpcobj);
fprintf('DC gain from output disturbance to output = %5.8f (=%g) \n',DC,DC);
DC gain from output disturbance to output = -0.00000000 (=-3.33067e-15) 

A zero gain means that the output will track the desired set-point.

Simulate MPC with MPCMOVE

First get the discrete-time state-space matrices of the plant.

[A,B,C,D] = ssdata(model);
Tstop = 5;                  % Simulation time
x = [0 0 0 0 0]';           % Initial state of the plant
xmpc = mpcstate(mpcobj);    % Initial state of the MPC controller
r = 1;                      % Output reference signal

Store the closed-loop MPC trajectories in arrays YY,UU,XX.

YY=[];
UU=[];
XX=[];

Create simulation loop

for t=0:round(Tstop/Ts)-1
    % Store states
    XX = [XX,x]; %#ok<*AGROW>
    % Define measured disturbance signal
    v = 0;
    if t*Ts>=10,
        v = 1;
    end
    % Define unmeasured disturbance signal
    d = 0;
    if t*Ts>=20,
       d = -0.5;
    end
    % Plant equations: output update (no feedthrough from MV to Y)
    y = C*x + D(:,2)*v + D(:,3)*d;
    YY = [YY,y];
    % Compute MPC action
    u = mpcmove(mpcobj,xmpc,y,r,v);
    % Plant equations: state update
    x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d;
    UU = [UU,u];
end

Plot results.

figure
subplot(211)
plot(0:Ts:Tstop-Ts,YY)
grid
title('Output');
subplot(212)
plot(0:Ts:Tstop-Ts,UU)
grid
title('Input');

If at any time during the simulation you want to check the optimal predicted trajectories, it can be returned by MPCMOVE as well. Assume you want to start from the current state and have a set-point change to 0.5, and assume the measured disturbance is gone.

r = 0.5;
v = 0;
[~,Info] = mpcmove(mpcobj,xmpc,y,r,v);

We now extract the optimal predicted trajectories.

topt = Info.Topt;
yopt = Info.Yopt;
uopt = Info.Uopt;
figure
subplot(211)
stairs(topt,yopt);
title('Optimal sequence of predicted outputs')
grid
subplot(212)
stairs(topt,uopt);
title('Optimal sequence of manipulated variables')
grid

Linear Representation of MPC Controller

When the constraints are not active, MPC controller behaves like a linear controller. You can get the state-space form of the MPC controller.

LTI = ss(mpcobj,'rv');

Get state-space matrices of linearized controller.

[AL,BL,CL,DL] = ssdata(LTI);

Simulate linear MPC closed-loop system and compare the linearized MPC controller with the original MPC controller with constraints turned off.

mpcobj.MV = [];           % No input constraints
mpcobj.OV = [];           % No output constraints
Tstop = 5;                %Simulation time
xL = zeros(size(BL,1),1); % Initial state of linearized MPC controller
x = [0 0 0 0 0]';         % Initial state of plant
y = 0;                    % Initial measured output
r = 1;                    % Output reference set-point
u = 0;                    % Previous input command
YY = [];
XX = [];
xmpc = mpcstate(mpcobj);
for t = 0:round(Tstop/Ts)-1
    YY = [YY,y];
    XX = [XX,x];
    v = 0;
    if t*Ts>=10,
        v = 1;
    end
    d = 0;
    if t*Ts>=20,
        d = -0.5;
    end
    uold = u;
    % Compute the linear MPC control action
    u = CL*xL+DL*[y;r;v];
    % Compare the input move with the one provided by MPCMOVE
    uMPC = mpcmove(mpcobj,xmpc,y,r,v);
    dispStr(t+1) = {sprintf('t=%5.2f, u=%7.4f (provided by LTI), u=%7.4f (provided by MPCOBJ)',t*Ts,u,uMPC)}; %#ok<*SAGROW>
    % Update plant equations
    x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d;
    % Update controller equations
    xL = AL*xL + BL*[y;r;v];
    % Update output equations
    y = C*x + D(:,1)*u + D(:,2)*v + D(:,3)*d;
end
% Display results.
for t=0:round(Tstop/Ts)-1,
    disp(dispStr{t+1});
end
% Plot results.
figure
plot(0:Ts:Tstop-Ts,YY)
grid
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
t= 0.00, u= 5.2478 (provided by LTI), u= 5.2478 (provided by MPCOBJ)
t= 0.20, u= 3.0134 (provided by LTI), u= 3.0134 (provided by MPCOBJ)
t= 0.40, u= 0.2281 (provided by LTI), u= 0.2281 (provided by MPCOBJ)
t= 0.60, u=-0.9952 (provided by LTI), u=-0.9952 (provided by MPCOBJ)
t= 0.80, u=-0.8749 (provided by LTI), u=-0.8749 (provided by MPCOBJ)
t= 1.00, u=-0.2022 (provided by LTI), u=-0.2022 (provided by MPCOBJ)
t= 1.20, u= 0.4459 (provided by LTI), u= 0.4459 (provided by MPCOBJ)
t= 1.40, u= 0.8489 (provided by LTI), u= 0.8489 (provided by MPCOBJ)
t= 1.60, u= 1.0192 (provided by LTI), u= 1.0192 (provided by MPCOBJ)
t= 1.80, u= 1.0511 (provided by LTI), u= 1.0511 (provided by MPCOBJ)
t= 2.00, u= 1.0304 (provided by LTI), u= 1.0304 (provided by MPCOBJ)
t= 2.20, u= 1.0053 (provided by LTI), u= 1.0053 (provided by MPCOBJ)
t= 2.40, u= 0.9920 (provided by LTI), u= 0.9920 (provided by MPCOBJ)
t= 2.60, u= 0.9896 (provided by LTI), u= 0.9896 (provided by MPCOBJ)
t= 2.80, u= 0.9925 (provided by LTI), u= 0.9925 (provided by MPCOBJ)
t= 3.00, u= 0.9964 (provided by LTI), u= 0.9964 (provided by MPCOBJ)
t= 3.20, u= 0.9990 (provided by LTI), u= 0.9990 (provided by MPCOBJ)
t= 3.40, u= 1.0002 (provided by LTI), u= 1.0002 (provided by MPCOBJ)
t= 3.60, u= 1.0004 (provided by LTI), u= 1.0004 (provided by MPCOBJ)
t= 3.80, u= 1.0003 (provided by LTI), u= 1.0003 (provided by MPCOBJ)
t= 4.00, u= 1.0001 (provided by LTI), u= 1.0001 (provided by MPCOBJ)
t= 4.20, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ)
t= 4.40, u= 0.9999 (provided by LTI), u= 0.9999 (provided by MPCOBJ)
t= 4.60, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ)
t= 4.80, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ)

Running a closed-loop where all constraints are turned off is easy using SIM. We just specify an option in the SimOptions structure:

SimOptions = mpcsimopt(mpcobj);
SimOptions.Constr = 'off';    % Remove all MPC constraints
SimOptions.Unmeas = d;        % unmeasured input disturbance
sim(mpcobj,Tf,r,v,SimOptions);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Simulate Using Simulink®

To run this example, Simulink® is required.

if ~mpcchecktoolboxinstalled('simulink')
    disp('Simulink(R) is required to run this part of the example.')
    return
end
% Recreate MPC controller.
mpcobj = mpc(model,Ts,10,3);
mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10);
mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);
-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000.
-->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000.
-->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.

The continuous-time plant to be controlled has the following state-space realization:

[A,B,C,D] = ssdata(sys);

Now simulate closed-loop MPC in Simulink.

mdl1 = 'mpc_miso';
open_system(mdl1)
sim(mdl1);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Next, you run a simulation with sinusoidal output noise. Assume output measurements are affected by a sinusoidal measurement noise of frequency 0.1 Hz. You want to inform the MPC object about this so that state estimates can be improved.

omega = 2*pi/10;
mpcobj.Model.Noise = 0.5*tf(omega^2,[1 0 omega^2]);
% Also revise the MPC design.
mpcobj.Model.Disturbance = tf(0.1); % Model for unmeasured disturbance = white Gaussian noise with zero mean and variance 0.1
mpcobj.weights = struct('MV',0,'MVRate',0.1,'OV',.005);
mpcobj.predictionhorizon = 40;
mpcobj.controlhorizon = 3;
%Simulation time
Tstop=150;
mdl2 = 'mpc_misonoise';
open_system(mdl2)   % Open new Simulink(R) Model
sim(mdl2);          % Start Simulation
-->Integrated white noise added on measured output channel #1.
-->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design.

bdclose(mdl1);
bdclose(mdl2);
Was this topic helpful?