Model Predictive Control Toolbox

Designing Model Predictive Controller Equivalent to Infinite-Horizon LQR

This example shows how to design an infinite-horizon model predictive controller by setting the weights on the terminal predicted states.

Setup MPC Controller

Define the sampling time, prediction model, and input and output weights to be used in the MPC controller setup:

Ts = 0.1;                                   % Sampling time
A = [0.8 Ts;0 0.9];
B = [0;Ts];
C = [1 0];
sysd = ss(A,B,C,0,Ts);                      % Discrete-time prediction model
dcg = dcgain(sysd);                         % DC-gain of prediction model

Transform Infinite-Horizon LQR Problem Into Finite-Horizon MPC Problem

Compute the Riccati matrix associated with the LQR problem with output weight Qy and input weight Qu:

Qy = 10;                                    % Output weight: y'*Qy*y
Qu = 0.1;                                   % Input weight:  u'*Qu*y
[K,P] = lqry(sysd,Qy,Qu);                   % LQR gain and Riccati matrix

Weight the terminal state x'(t+p)*P*x(t+p), where p is the prediction horizon of the MPC controller. Compute the Cholesky factor chol(P) of the Riccati matrix P, so that the terminal penalty becomes:

x'(t+p)*P*x(t+p) = [chol(P)*x(t+p)]'*[chol(P)*x(t+p)] = yc'(t+p)*yc(t+p)

where the new output yc(t+p) = chol(P)*x(t+p) and the state x(t) is assumed to be fully measurable.

cholP = chol(P);

Change and augment the output vector to include the full state x and yc:

set(sysd,'C',[eye(2);cholP],'D',zeros(4,1));% Output = state vector x + output yc such that yc'*yc = x'*P*x

Label the new additional output signal yc(t) as unmeasured:

sysd = setmpcsignals(sysd,'UO',[3 4]);      % Cholesky factor is not measured
-->Assuming unspecified output signals are measured outputs.

Define prediction horizons and input saturation constraints:

p = 3;                                      % Prediction horizon (for any p>=1, unconstrained MPC = LQR)
m = p;                                      % Control horizon = prediction horizon
mpc1 = mpc(sysd,Ts,p,m);                    % MPC object
mpc1.MV = struct('Min',-3,'Max',3);         % Input saturation constraints
-->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.
   for output(s) y1 and zero weight for output(s) y2 y3 y4 

Define input and output weights at each step of the prediction horizon (terminal weights are set later). Very small weights on input increments are included to make the QP problem associated with the MPC controller positive definite:

mpc1.Weights.OV = [sqrt(Qy) 0 0 0];         % Output weights (only on original output)
mpc1.Weights.MV = sqrt(Qu);                 % Input weight
mpc1.Weights.MVRate = 1e-5;                 % Very small weight on command input increments

Define set-points for the output and input signals:

ry = 1;                                     % Output set point
mpc1.MV.Target = ry/dcg;                    % Set-point for manipulated variable

Impose the terminal penalty x'(t+p)*P*x(t+p) by specifying a unit weight on yc(t+p) = chol(P)*x(t+p). The terminal weight on u(t+p-1) remains the same, that is sqrt(Qu):

Y = struct('Weight',[0 0 1 1]);             % Weight on y(t+p)
U = struct('Weight',sqrt(Qu));              % Weight on u(t+p-1)
setterminal(mpc1,Y,U);                      % Set terminal weight y'*y = x'*P*x

Since the measured output is the entire state, remove any additional output disturbance integrator inserted by the MPC controller:

setoutdist(mpc1,'remove');                  % Remove additional disturbance integrators
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Remove the state estimator by defining the following measurement update equation:

x[n|n] = x[n|n-1] + I * (x[n]-x[n|n-1]) = x[n]

Because setterminal function resets state estimator to default value, we use setestim function to change state estimator after setterminal is called.

setestim(mpc1,eye(2));                      % State estimate = state measurement

Compare MPC and LQR Controllers

Compute the gain of the MPC controller when constraints are inactive, and compare it to the LQR gain:

mpcgain = dcgain(ss(mpc1));
fprintf('\n(unconstrained) MPC: u(k)=[%8.8g,%8.8g]*x(k)',mpcgain(1),mpcgain(2));
fprintf('\n                LQR: u(k)=[%8.8g,%8.8g]*x(k)\n\n',-K(1),-K(2));
(unconstrained) MPC: u(k)=[-2.8363891,-2.1454028]*x(k)
                LQR: u(k)=[-2.8363891,-2.1454028]*x(k)

The state feedback gains are exactly the same.

Simulate the Model in Simulink®

Define a set-point for the new extended output vector containing x and chol(P)*x:

rx = (eye(2)-A)\B/dcg;
r = [rx;cholP*rx];                          % Set point for extended prediction model

Define initial condition and total simulation time:

x0 = [0;0];                                 % Initial state
Tstop = 5;                                  % Simulation time

Open and simulate the Simulink® model:

if ~mpcchecktoolboxinstalled('simulink')
    disp('Simulink(R) is required to run this example.')
    return
end
open_system('mpc_infinite');                % Open Simulink(R) Model
sim('mpc_infinite',Tstop);                  % Start Simulation

bdclose('mpc_infinite');