Main Content

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

Define a plant model. For this example use continuous time transfer functions from each input to the output. Display the plant model at the command line.

```
plantTF = tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}) % display transfer functions
```

plantTF = From input 1 to output: 1 --------------- s^2 + 0.5 s + 1 From input 2 to output: 1 ----- s + 1 From input 3 to output: 1 ------------------- 0.7 s^2 + 0.5 s + 1 Continuous-time transfer function.

For this example, explicitly convert the plant to a discrete-time state space form before passing in to the MPC controller creation function.

The controller creation function can accept either continuous-time or discrete-time plants. When the optimization problem to find the optimal value of the manipulated variable is set up, the MPC controller automatically converts a continuous-time plant (in any format) to a discrete-time state space model, using Zero Order Hold.

Converting the plant to discrete-time is useful when you need the discrete-time system matrices for analysis or simulation (as in this example) or want the controller to use a discrete-time plant converted with a method other than ZOH.

plantCSS = ss(plantTF); % convert plant from transfer function to continuous-time state space Ts = 0.2; % specify a sample time of 0.2 seconds plantDSS = c2d(plantCSS,Ts); % convert plant to discrete-time state space, using Zero Order Hold

By default, all the plant input signals are assumed to be manipulated variables. Use `setmpcsignal`

to specify which signals are measured and unmeasured disturbances. In this example, the first input signal is a manipulated variable, the second is a measured disturbance, the third one is an unmeasured disturbance. This information is stored in the plant model `plantDSS`

and later used by the mpc controller.

plantDSS = setmpcsignals(plantDSS,'MV',1,'MD',2,'UD',3); % specify signal types

Create the controller object, specifying the sampling time, as well as the prediction and control horizons (10 and 3 steps respectively).

mpcobj = mpc(plantDSS,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.

Because you have not specified the weights of the quadratic cost function to be minimized by the controller, their value is assumed to be the default one (0 for the manipulated variables, 0.1 for the manipulated variable rates, 1 for the output variables). Also, at this point the MPC problem is still unconstrained as you have not specified any constraint yet.

Define hard constraints on the manipulated variable.

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

The input and output disturbance models specify the dynamic characteristics of the unmeasured disturbances on the input and output, respectively, so they can be better rejected. By default, these disturbance models are assumed to be integrators unless you specify them otherwise. The mpc object also has a noise model that specifies the dynamics of the noise on the measured output variable. By default this is assumed to be a unity static gain, which is equivalent to assume that the noise is white. In this example, there is no measured output disturbance, so there is no need to specify an output disturbance model, and the noise on the measured output signal is assumed to be white.

Specify the disturbance model for the unmeasured input as an integrator driven by white noise with variance = `1000`

.

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

Display the MPC controller object `mpcobj`

to review its properties.

mpcobj

MPC object (created on 23-Feb-2021 13:06:50): --------------------------------------------- Sampling time: 0.2 (seconds) Prediction Horizon: 10 Control Horizon: 3 Plant Model: -------------- 1 manipulated variable(s) -->| 5 states | | |--> 1 measured output(s) 1 measured disturbance(s) -->| 3 inputs | | |--> 0 unmeasured output(s) 1 unmeasured disturbance(s) -->| 1 outputs | -------------- Indices: (input vector) Manipulated variables: [1 ] Measured disturbances: [2 ] Unmeasured disturbances: [3 ] (output vector) Measured outputs: [1 ] Disturbance and Noise Models: Output disturbance model: default (type "getoutdist(mpcobj)" for details) Input disturbance model: user specified (type "getindist(mpcobj)" for more details) Measurement noise model: default (unity gain after scaling) Weights: ManipulatedVariables: 0 ManipulatedVariablesRate: 0.1000 OutputVariables: 1 ECR: 100000 State Estimation: Default Kalman Filter (type "getEstimator(mpcobj)" for details) Constraints: 0 <= MV1 <= 1, -10 <= MV1/rate <= 10, MO1 is unconstrained

Display the input disturbance model. As expected is the specified integrator converted to discrete time.

getindist(mpcobj)

ans = A = x1 x1 1 B = Noise#1 x1 0.8 C = x1 UD1 7.906 D = Noise#1 UD1 0 Sample time: 0.2 seconds Discrete-time state-space model.

Display the output disturbance model. As expected it is empty.

getoutdist(mpcobj)

Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel. ans = Empty state-space model.

To examine whether the MPC controller will be able to reject constant output disturbances and track constant setpoint with zero offsets in steady-state, compute the closed loop DC gain from output disturbances to controlled outputs using the `cloffset`

command. This is also known as the steady state output sensitivity of the closed loop.

```
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 (=5.77316e-15)

A zero gain (which is typically the result of the controller having integrators as input or output disturbance models) means that the measured plant output will track the desired output reference setpoint.

`sim`

CommandThe `sim`

command provides a quick way to simulate an MPC controller in a closed loop with a linear time-invariant plant when constraints and weights stay constants and the disturbance and reference signals can be easily and completely specified ahead of time.

First, specify simulation time and the reference and disturbance signals

Tstop = 30; % simulation time Nf = round(Tstop/Ts); % number of simulation steps r = ones(Nf,1); % output reference signal v = [zeros(Nf/3,1);ones(2*Nf/3,1)]; % measured input disturbance signal

Run the closed-loop simulation and plot results. The plant specified in `mpcobj.Model.Plant`

is used both as a plant model in the closed loop and as the internal plant model used by the controller to predict the response over the prediction horizon. Use `sim`

to simulate the closed loop system for Nf steps with reference r and measured input disturbance v.

```
sim(mpcobj,Nf,r,v) % simulate plant and controller in closed loop
```

The manipulated variable hits the upper bound initially, and brings the plant output to the reference value within a few seconds. It then settles in at its maximum allowed value, 1. After 10 seconds, the measured disturbance signal rises from 0 to 1, which causes the plant output to exceed its reference value by about 30%. The manipulated variable hits the lower bound in an effort to reject the disturbance. The controller is able to bring the plant output back to the reference value after a few seconds, and the manipulated variable settles in at its minimum value. The unmeasured disturbance signal is always zero, because no unmeasured disturbance has been specified yet.

You can use a simulation options objects to specify additional simulation options as well as noise and disturbance signals that feed into the plant but are unknown to the controller. For this example, use a simulation option object to specify the unmeasured input disturbance signal as well as additive noise on both manipulated variables and measured outputs. You will also later use this options object to specify a plant to be used in simulation, which differs from the one used internally by the controller. Create a simulation option object with default options and no signal.

```
SimOptions = mpcsimopt; % create object
```

Create a disturbance signal and specify it in the simulation options object

d = [zeros(2*Nf/3,1);-0.5*ones(Nf/3,1)]; % define a step disturbance signal SimOptions.UnmeasuredDisturbance = d; % specify unmeasured input disturbance signal

Specify noise signals in the simulation options object. At simulation time, the simulation function directly adds the specified output noise to the measured output before feeding it to the controller. It also directly adds the specified input noise to the manipulated variable (not to any disturbance signals) before feeding it to the plant.

SimOptions.OutputNoise=.001*(rand(Nf,1)-.5); % specify output measurement noise SimOptions.InputNoise=.05*(rand(Nf,1)-.5); % specify noise on manipulated variables

Note that you can also use the `OutputNoise`

field of the simulation option object to specify a more general additive output disturbance signal (such as a step, for example) on the measured plant output.

Use `sim`

with the additional `SimOptions`

argument to simulate the closed loop system and save the results to the workspace variables `y`

,|t|,|u|, and `xp`

. This allows you to selectively plot signals in a new figure window and in any given color and order.

```
[y,t,u,xp] = sim(mpcobj,Nf,r,v,SimOptions); % simulate closed loop
```

Plot results.

figure % create new figure subplot(2,1,1) % create upper subplot plot(0:Nf-1,y,0:Nf-1,r) % plot plant output and reference title('Output') % add title so upper subplot ylabel('MO1') % add a label to the upper y axis grid % add a grid to upper subplot subplot(2,1,2) % create lower subplot plot(0:Nf-1,u) % plot manipulated variable title('Input'); % add title so lower subplot xlabel('Simulation steps') % add a label to the lower x axis ylabel('MV1') % add a label to the lower y axis grid % add a grid to lower subplot

Despite the added noise (which is especially visible on the manipulated variable plot) and despite the measured and unmeasured disturbances kicking in after 50 and 100 steps respectively, the controller is able to achieve and maintain good tracking. The manipulated variable settles in at about 1 after the initial part of the simulation (steps from 20 to 50), at about 0 to reject the measured disturbance (steps from 70 to 100), and at about 0.5 to reject both disturbances (steps from 120 to 150).

Test the robustness of the MPC controller against a model mismatch. Specify the true plant to be used in simulation as `truePlantCSS`

.

truePlantTF = tf({1,1,1},{[1 .8 1],[1 2],[.6 .6 1]}) % specify and display transfer functions truePlantCSS = ss(truePlantTF); % convert to continuous state space truePlantCSS = setmpcsignals(truePlantCSS,'MV',1,'MD',2,'UD',3); % specify signal types

truePlantTF = From input 1 to output: 1 --------------- s^2 + 0.8 s + 1 From input 2 to output: 1 ----- s + 2 From input 3 to output: 1 ------------------- 0.6 s^2 + 0.6 s + 1 Continuous-time transfer function.

Update the simulation option object by specifying `SimOptions.Model`

as a structure with two fields, `Plant`

(containing the true plant model) and `Nominal`

(containing the operating point values for the true plant). For this example, the nominal values for the state derivatives and the inputs are not specified, so they are assumed to be zero, resulting in `y = SimOptions.Model.Nominal.Y + C*(x-SimOptions.Model.Nominal.X)`

where `x`

and `y`

are the state and measured output of the plant, respectively.

SimOptions.Model = struct('Plant',truePlantCSS); % create the structure and assign the 'Plant' field SimOptions.Model.Nominal.Y = 0.1; % create and assign the 'Nominal.Y' field SimOptions.Model.Nominal.X = -.1*[1 1 1 1 1]; % create and assign the 'Nominal.X' field SimOptions.PlantInitialState = [0.1 0 -0.1 0 .05]; % specify the initial state of the true plant

remove any signal previously added to the measured output and to the manipulated variable

SimOptions.OutputNoise = []; % remove output measurement noise SimOptions.InputNoise = []; % remove noise on manipulated variable

Run the closed-loop simulation and plot the results. Since `SimOptions.Model`

is not empty, `SimOptions.Model.Plant`

is converted to discrete time (using zero order hold) and used as the plant model in the closed loop simulation, while the plant in `mpcobj.Model.Plant`

is only used by the controller to predict the response over the prediction horizon.

```
sim(mpcobj,Nf,r,v,SimOptions) % simulate the closed loop
```

-->Converting model to discrete time.

As a result of the model mismatch, some degradation in the response is visible, notably the controller needs a little more time to achieve tracking and the manipulated variable now settles at about 0.5 to reject the measured disturbance (see values from 5 to 10 seconds) and settles at about 0.9 to reject both input disturbances (from 25 to 30 seconds). Despite this, the controller is still able to track the output reference.

Every constraint is associated to a dimensionless parameter called ECR value. A constraint with larger ECR value is allowed to be violated more than a constraint with smaller ECR value. By default all constraints on the manipulated variables have an ECR value of zero, making them hard. You can specify a nonzero ECR value for a constraint to make it soft.

Note that is possible to refer to the `ManipulatedVariables`

field also as `MV`

. Relax the constraints on manipulated variables from hard to soft.

mpcobj.MV.MinECR = 1; % ECR for the lower bound on the manipulated variable mpcobj.MV.MaxECR = 1; % ECR for the upped bound on the manipulated variable

Define an output constraint. By default all constraints on output variables (measured outputs) have an ECR value of one, making them soft. You can reduce the ECR value for an output constraint to make it harder, however it is always advised to keep output constraints soft. This is mainly because the plant outputs depend on plant states as well as on the measured disturbances, therefore, if a large disturbance occurs, the plant outputs can violate the constraints independently on any control action taken by the mpc controller, especially when the manipulated variables are themselves hard bounded. Such an unavoidable violation of an hard constraint results in an infeasible MPC problem, for which no MV can be calculated.

mpcobj.OV.Max = 1.1; % define the (soft) output constraint % Note that is possible to refer to the |OutputVariables| field also as |OV|.

Run a new closed-loop simulation, without including the simulation option object, and therefore without any model mismatch, unmeasured disturbance, or noise added to the manipulated variable or measured output.

```
sim(mpcobj,Nf,r,v) % simulate the closed loop
```

Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

In an effort to reject the measured disturbance, achieve tracking, and prevent the output to rises above its 1.1 soft constraint, the controller slightly violates the soft constraint on the manipulated variable, which reaches small negative values from seconds 10 to 11 (you can zoom in the picture to see this violation more clearly). The constraint on the measured output is violated more than the one on the manipulated variable.

Penalize more the output constraint violations and rerun the simulation.

mpcobj.OV.MaxECR = 0.001; % the closer to zero, the harder the constraint sim(mpcobj,Nf,r,v) % run a new closed-loop simulation.

Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Now the controller violates the output constraint only very slighly. As expected, this output performance improvement comes at the cost of violating the manipulated variable constraint a lot more (the manipulated variable reaches -3 for a couple of steps).

At each time step, the MPC controller obtains the manipulated variable by multiplying the discrete-time plant state (if available) by a gain matrix (computed by solving a constrained quadratic optimization problem). Since the plant state is normally not available, by default, the controller uses a linear Kalman filter as an observer to estimate the state of the discrete-time augmented plant (that is the plant comprehensive of the disturbance and noise models). Therefore, the states of the controller are the states of this Kalman filter, which are in turn the estimates of the states of the augmented discrete-time plant.

Run a closed loop simulation with model mismatch and unmeasured disturbance, using the default estimator, and return the controller states structure `xc`

in the workspace, so you can later compare the state sequences.

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

-->Converting model to discrete time.

display the component of the controller states structure

xc

xc = struct with fields: Plant: [150x5 double] Disturbance: [150x1 double] Noise: [150x0 double] LastMove: [150x1 double]

Plot the plant output response as well as the plant states estimated by the default observer

figure; % create figure subplot(2,1,1) % create upper subplot axis plot(t,y) % plot y versus time title('Plant output'); % add title to upper plot ylabel('y') % add a label to the upper y axis grid % add grid to upper plot subplot(2,1,2) % create lower subplot axis plot(t,xc.Plant) % plot xc.Plant versus time title('Estimated plant states'); % add title to lower plot xlabel('Time (sec)') % add a label to the lower x axis ylabel('xc') % add a label to the lower y axis grid % add grid to lower plot

As expected, there are sudden changes at 10 and 20 seconds due to the measured and unmeasured disturbance signals kicking in, respectively.

You can change the gains of the Kalman filter used as observer. To do so, first, retrieve the default Kalman gains and state-space matrices.

```
[L,M,A1,Cm1] = getEstimator(mpcobj); % retrieve observer matrices
```

Calculate and display the poles of the default observer. Note that, as expected, they are all inside the unit circle, though a few of them seem relatively close to the border. Note that there are six states, five belonging to the plant model, the sixth belonging to the input disturbance model.

```
e = eig(A1-A1*M*Cm1) % eigenvalues of observer state matrix
```

e = 0.5708 + 0.4144i 0.5708 - 0.4144i 0.4967 + 0.0000i 0.9334 + 0.1831i 0.9334 - 0.1831i 0.8189 + 0.0000i

Design a new state estimator by pole-placement.

poles = [.8 .75 .7 .85 .6 .81]; % specify desired positions for the new poles L = place(A1',Cm1',poles)'; % calculate Kalman gain for time update M = A1\L; % calculate Kalman gain for measurement update

Set the new matrix gains in the mpc controller object

```
setEstimator(mpcobj,L,M); % set the new estimation gains
```

re-run the closed loop simulation

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

Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel. -->Converting model to discrete time.

Plot the plant output response as well as the plant states estimated by the new observer

figure; % create figure subplot(2,1,1) % create upper subplot axis plot(t,y) % plot y versus time title('Plant output'); % add title to upper plot ylabel('y') % add a label to the upper y axis grid % add grid to upper plot subplot(2,1,2) % create lower subplot axis plot(t,xc.Plant) % plot xc.Plant versus time title('Estimated plant states'); % add title to lower plot xlabel('Time (sec)') % add a label to the lower x axis ylabel('xc') % add a label to the lower y axis grid % add grid to lower plot

As expected the controller states are different from the ones previously plotted, and the overall closed loop response is somewhat slower.

Test the behavior of the plant and controller in open-loop, using the `sim`

command. Set the `OpenLoop`

flag to `on`

, and provide the sequence of manipulated variables to excite the system.

SimOptions.OpenLoop = 'on'; % set open loop option SimOptions.MVSignal = sin((0:Nf-1)'/10); % define the manipulated variable signal

Simulate the open loop system containing the true plant (previously specified in `SimOptions.Model`

) followed by the controller. As the reference signal will be ignored, you can avoid specifying it.

```
sim(mpcobj,Nf,[],v,SimOptions) % simulate the open loop system
```

-->Converting model to discrete time.

`mpcmove`

In the more general case of a closed loop in which the controller is applied to a nonlinear or time varying plant, constraints or weights can change at run time, or the disturbance and reference signals are hard to specify completely ahead of time, you can use the `mpcmove`

function at each step in a for loop to simulate the MPC controller in closed loop. If your plant is continuous-time, you will need to convert it to discrete-time (using any method) to simulate it with a for loop.

First, obtain the discrete-time state-space matrices of the plant, and define simulation time and initial states for plant and controller.

[A,B,C,D] = ssdata(plantDSS); % discrete-time plant plant ss matrices Tstop = 30; % 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

Display the initial state of the controller. Note that this is an `mpcstate`

object (not a simple structure) and contains the controller states only at the current time (not the whole history up to the current time, as the structure returned by sim). It also contains the estimator covariance matrix. At each simulation step, you need to update `xmpc`

with the new controller states.

```
xmpc % controller states
```

MPCSTATE object with fields Plant: [0 0 0 0 0] Disturbance: 0 Noise: [1x0 double] LastMove: 0 Covariance: []

Define workspace arrays `YY`

and `UU`

to store the closed-loop trajectories so that they can be plotted after the simulation.

YY=[]; UU=[];

Run the simulation loop

for k=0:round(Tstop/Ts)-1 % Define measured disturbance signal v(k). You might specify a more % complex dependence on time or previous states here. v = 0; if k*Ts>=10 % raising to 1 after 10 seconds v = 1; end % Define unmeasured disturbance signal d(k). You might specify a more % complex dependence on time or previous states here. d = 0; if k*Ts>=20 % falling to -0.5 after 20 seconds d = -0.5; end % Plant equations: current output % If you have a more complex plant output behavior (including, for example, % model mismatch or nonlinearities) you might want to simulate it here. % Note that there cannot be any direct feedthrough between u and y. y = C*x + D(:,2)*v + D(:,3)*d; % calculate current output (D(:,1)=0) YY = [YY,y]; % store current output % Compute the MPC action (u) and update the controller states (xmpc) using mpcmove % Note that xmpc is updated by mpcmove even though it is an input argument! u = mpcmove(mpcobj,xmpc,y,r,v); % xmpc,y,r,v are values at current step k % Plant equations: state update % You might want to simulate a more complex plant state behavior here. x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d; % calculate next state and update current state UU = [UU,u]; % store current input end

Plot the results.

figure % create figure subplot(2,1,1) % create upper subplot axis plot(0:Ts:Tstop-Ts,YY) % plot YY versus time ylabel('y') % add a label to the upper y axis grid % add grid to upper plot title('Output') % add title to upper plot subplot(2,1,2) % create lower subplot axis plot(0:Ts:Tstop-Ts,UU) % plot UU versus time ylabel('y') % add a label to the lower y axis xlabel('Time (sec)') % add a label to the lower x axis grid % add grid to lower plot title('Input') % add title to lower plot

If at any time during the simulation you want to check the optimal predicted trajectories from that point, they can be returned by `mpcmove`

as well. Assume you start from the current state (`x`

and `xmpc`

), while the reference set-point changes to 0.5 and that the disturbance is gone (and that both signals remain constant during the prediction horizon).

r = 0.5; % reference v = 0; % disturbance [~,info] = mpcmove(mpcobj,xmpc,y,r,v); % solve over prediction horizon

Display the info variable.

info

info = struct with fields: Uopt: [11x1 double] Yopt: [11x1 double] Xopt: [11x6 double] Topt: [11x1 double] Slack: 0 Iterations: 1 QPCode: 'feasible' Cost: 0.1399

`info`

is a structure containing the predicted optimal sequences of manipulated variables, plant states, and outputs over the prediction horizon. This sequence is calculated, together with the optimal move, by solving a quadratic optimization problem to minimize the cost function. The plant states and outputs in `info`

result from applying the optimal manipulated variable sequence directly to `mpcobj.Model.Plant`

, in an open-loop fashion. Therefore, this open-loop optimization process is therefore not equivalent to simulating the closed loop consisting of plant, estimator and controller using either the `sim`

command or `mpcmove`

iteratively in a for loop.

Extract the predicted optimal trajectories.

topt = info.Topt; % time yopt = info.Yopt; % predicted optimal plant model outputs uopt = info.Uopt; % predicted optimal mv sequence

Plot the trajectories using stairstep plots. The stairstep plots are used because a normal plot would simply join each calculated point with the following one using a direct straight line, while the optimal discrete-time signals jump instantly at each step and stay constant until the next step. This difference is especially visible since there are only 10 intervals to be plotted for the plant output (and only 3 for the manipulated variable).

figure % create new figure subplot(2,1,1) % create upper subplot stairs(topt,yopt) % plot yopt in a stairstep graph title('Optimal sequence of predicted outputs') % add title to upper subplot grid % add grid to upper subplot subplot(2,1,2) % create lower subplot stairs(topt,uopt) % plot uopt in a stairstep graph title('Optimal sequence of manipulated variables') % add title to upper subplot grid % add grid to upper subplot

When the constraints are not active, the MPC controller behaves like a linear controller. Note that for a finite-time unconstrained Linear Quadratic Regulator problem with a finite non-receding horizon, the value function is time-dependent, which causes the optimal feedback gain to be time varying. In contrast, in MPC the horizon has a constant lenght, because it is always receding, resulting in a time-invariant value function, and consequently in a time-invariant optimal feedback gain.

You can get the state-space form of the MPC controller.

LTI = ss(mpcobj,'rv'); % get state space representation

Get the state-space matrices to simulate the linearized controller.

```
[AL,BL,CL,DL] = ssdata(LTI); % get state space matrices
```

Initialize variables for a closed loop simulation of both the original MPC controller without constraints and the linearized controller

mpcobj.MV = []; % remove input constraints mpcobj.OV = []; % remove output constraints Tstop = 5; % set simulation time y = 0; % set nitial measured output r = 1; % set output reference set-point (constant) u = 0; % set previous (initial) input command x = [0 0 0 0 0]'; % set initial state of plant xmpc = mpcstate(mpcobj); % set initial state of unconstrained MPC controller xL = zeros(size(BL,1),1); % set initial state of linearized MPC controller YY = []; % define workspace array to store plant outputs

Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Simulate both controllers in a closed loop with the same plant model

for k = 0:round(Tstop/Ts)-1 YY = [YY,y]; % store current output for plotting purposes % Define measured disturbance signal v(k). v = 0; if k*Ts>=10 v = 1; % raising to 1 after 10 seconds end % Define unmeasured disturbance signal d(k). d = 0; if k*Ts>=20 d = -0.5; % falling to -0.5 after 20 seconds end % Compute the control actions of both (unconstrained) MPC and linearized MPC uMPC = mpcmove(mpcobj,xmpc,y,r,v); % unconstrained MPC (this also updates xmpc) u = CL*xL+DL*[y;r;v]; % unconstrained linearized MPC % Compare the two control actions dispStr(k+1) = {sprintf('t=%5.2f, u=%7.4f (provided by LTI), u=%7.4f (provided by MPCOBJ)',k*Ts,u,uMPC)}; %#ok<*SAGROW> % Update state of unconstrained linearized MPC controller xL = AL*xL + BL*[y;r;v]; % Update plant state x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d; % Calculate plant output y = C*x + D(:,1)*u + D(:,2)*v + D(:,3)*d; % D(:,1)=0 end

Display the character arrays containing the control actions

for k=0:round(Tstop/Ts)-1 disp(dispStr{k+1}); % display each string as k increases end

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)

Plot the results.

figure % create figure plot(0:Ts:Tstop-Ts,YY) % plot plant outputs grid % add grid title('Unconstrained MPC control: Plant output') % add title xlabel('Time (sec)') % add label to x axis ylabel('y') % add label to y axis

Running a closed-loop simulation in which all controller constraints are turned off is easier using `sim`

, as you just need to specify 'off' in the `Constraint`

field of the related `mpcsimopt`

simulation option object.

SimOptions = mpcsimopt; % create simulation options object SimOptions.Constraints = 'off'; % remove all MPC constraints SimOptions.UnmeasuredDisturbance = d; % unmeasured input disturbance sim(mpcobj,Nf,r,v,SimOptions); % run closed loop simulation

Simulink is a graphical block diagram environment for multidomain system simulation. You can connect blocks representing dynamical systems (in this case the plant and the MPC controller) and simulate the closed loop. Besides the fact that the system and its interconnections are represented visually, one key difference is that Simulink can use a wider range of integration algorithms to solve the underlying system differential equations. Simulink can also automatically generate C (or PLC) code and deploy it for real-time applications.

% To run this part of the example, Simulink(R) is required. Check that % Simulink is installed, otherwise display a message and return if ~mpcchecktoolboxinstalled('simulink') % if Simulink not installed disp('Simulink(R) is required to run this part of the example.') return % end execution here end

Recreate the original mpc object with the original constraint and the original default estimator, so you can easily compare results.

mpcobj = mpc(plantDSS,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.

Extract the state-space matrices of the (continuous-time) plant to be controlled, since they are needed by the Simulink model to be opened.

```
[A,B,C,D] = ssdata(plantCSS); % get state-space realization
```

Open the pre-existing Simulink model for the closed-loop simulation. The plant model is implemented with a continuous state space block, and the ode45 integration algorithm is used to calculate its continuous time response.

%The block inputs are u(t), v(t) and d(t) representing the manipulated % variable, measured input disturbance, and unmeasured input disturbance, % respectively, while y(t) is the measured output. The block parameters % are the matrices forming the state space realization of the continuous-time % plant, and the initial conditions for the 5 states. % The MPC controller is implemented with an MPC controller block, which % has the workspace mpc object |mpcobj| as parameter, the manipulated % variable as output, and then the measured plant output, reference signal, % and measured plant input disturbance as inputs, respectively. % The four Scopes blocks plot the five loop signals, which are also saved % (except the reference signal) by four To-Workspace blocks. open_system('mpc_miso')

Simulate the system using the simulink `sim`

command. Note that this command (which simulates a Simulink model, and is equivalent to clicking the "Run" button in the model) is different from the `sim`

command provided by the mpc toolbox (which instead simulates an MPC controller in a loop with an LTI plant).

```
sim('mpc_miso')
```

The plots in the `Scopes`

blocks are equivalent to the ones in figures 1,2 and 3, with minor differences due to the fact that in the first simulation (Figures 1 and 2) the unmeasured disturbance signal was zero, and that is the second simulation (Figure 3) noise was added to the plant input and output. Also note that that, while the MPC `sim`

command internally discretizes any continuous plant model using ZOH, Simulink typically uses an integration algorithm (in this example ode45) to simulate the closed loop when a continuous-time block is present.

Assume output measurements are affected by a sinusoidal disturbance (a single tone sensor noise) of frequency 0.1 Hz on the measured output.

```
omega = 2*pi/10; % disturbance radial frequency
```

Open the `mpc_misonoise`

Simulink model. It's similar to the 'mpc_miso' model except for the sinusoidal disturbance injected on the measured output. What is also different is that the simulation time is now set to 150 seconds, and that unmeasured and measured input disturbances start to act at after 60 and 120 seconds, respectively. Note that, differently from the previous simulations, the unmeasured disturbance start first. This allows the response to unmeasured disturbance, taking place from 60 to 120 seconds, to be plotted and analyzed more clearly.

open_system('mpc_misonoise') % open new Simulink model

Since this noise is expected, specifying it in a measurement noise model will help the state estimator to ignore it, thereby improving the quality of the state estimates and allowing the whole controller to better reject disturbances while tracking the output reference.

```
mpcobj.Model.Noise = 0.5*tf(omega^2,[1 0 omega^2]); % measurement noise model
```

Revise the MPC design by specifing a disturbance model on the unmeasured input as a white Gaussian noise with zero mean and variance `0.1`

. This allows a better rejection of generic non constant disturbances (at the expense of a worse rejection of input constant disturbances).

```
mpcobj.Model.Disturbance = tf(0.1); % static gain
```

Note that when you specify a static gain as a model for the input disturbance, an output disturbance model consisting in a discretized integrator is automatically added to the controller. This helps rejecting constant (and slow varying) output disturbances.

getoutdist(mpcobj)

-->Assuming output disturbance added to measured output channel #1 is integrated white noise. -->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design. ans = A = x1 x1 1 B = u1 x1 0.2 C = x1 MO1 1 D = u1 MO1 0 Sample time: 0.2 seconds Discrete-time state-space model.

Decrease the weight on the tracking of the output variable, (thereby placing more emphasis on minimizing the rate of change of the manipulated variables).

mpcobj.weights = struct('MV',0,'MVRate',0.1,'OV',0.005); % new weights

By themselves, the new weights would penalize the manipulated variable change too much, resulting in a very low bandwidth (slow) controller. This would result in the output variables lagging behind the reference for many steps. Increasing the predicion horizon to 40 allows the controller to fully take into account the cost of the output error for 40 steps instead of just 10, thereby placing back some emphasis on tracking.

```
mpcobj.predictionhorizon = 40; % new prediction horizon
```

Run the simulation for 145 seconds

sim('mpc_misonoise',145) % the second argument is the simulation duration

-->Assuming output disturbance added to measured output channel #1 is integrated white noise. -->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design.

The kalman filter successfully learns to ignore the measurement noise after 50 seconds. The unmeasured and measured disturbances get rejected in a 10 to 20 second timespan. Finally, as expected, the manipulated variable stays in the interval between 0 and 1.

bdclose('all') % close all open Simulink models without saving any change close all % close all open figures

`mpc`

| MPC
Controller | MPC
Designer