MATLAB Examples

Design an MPC Controller Equivalent to LQR Controller

This example shows shows how to design an unconstrained MPC controller that provides performance equivalent to an LQR controller.

Contents

Define the Plant Model

The plant is a double integrator, represented as a state-space model in discrete time with a sample time of 0.1 seconds. In this case the two plant states are measurable at the plant outputs.

A = [1 0;0.1 1];
B = [0.1;0.005];
C = eye(2);
D = zeros(2,1);
Ts = 0.1;
Plant = ss(A,B,C,D,Ts);
Plant.InputName = {'u'};
Plant.OutputName = {'x_1','x_2'};

Design an LQR Controller

Design an LQR controller with output feedback for the plant.

Q = eye(2);
R = 1;
[K,Qp] = lqry(Plant,Q,R);

Q and R are output and input weight matrices, respectively. Qp is the Ricatti matrix.

Design an Equivalent MPC Controller

To implement the MPC cost function, compute $L$, the Cholesky decomposition of $Q_p$, such that ${L^T}L = {Q_p}$. Then define auxilliary unmeasured output variables ${y_a}\left( k \right) = Lx\left( k\right)$, such that $y_a^T{y_a} = {x^T}{Q_p}x$.

NewPlant = Plant;
L = chol(Qp);
set(NewPlant,'C',[C;L],'D',[D;zeros(2,1)],...
    'OutputName',{'x_1','x_2','Cx_1','Cx_2'})
NewPlant.InputGroup.MV = 1;
NewPlant.OutputGroup.MO = [1 2];
NewPlant.OutputGroup.UO = [3 4];

Create an MPC controller with equal prediction and control horizons.

P = 3;
M = 3;
MPCobj = mpc(NewPlant,Ts,P,M);
-->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 

Specify weights for the manipulated variable and output variables for the first $p-1$ prediction horizon steps. Use the square roots of the diagonal elements of the $Q$ and $R$ weight matrices from the LQR controller design. The standard $Q$ weight matrix values apply to $y$, and $y_a$ has a zero penalty.

ywt = sqrt(diag(Q))';
uwt = sqrt(diag(R))';
MPCobj.Weights.OV = [ywt 0 0];
MPCobj.Weights.MV = uwt;

Specify terminal weights for the final prediction horizon step. On step $p$, the original $y$ has a zero penalty, and $y_a$ has a unit penalty. The input weight remains the same for the terminal step.

U = struct('Weight',uwt);
Y = struct('Weight',[0 0 1 1]);
setterminal(MPCobj,Y,U)

Remove the default state estimator. Since the model states are measured directly, the default state estimator is unnecessary.

setoutdist(MPCobj,'model',tf(zeros(4,1)))
setEstimator(MPCobj,[],C)

Compare Controller Performance

Compare the performance of the LQR controller, the MPC controller with terminal weights, and a standard MPC controller.

Compute the closed-loop response for the LQR controller.

clsys = feedback(Plant,K);
Tstop = 6;
x0 = [0.2;0.2];
[yLQR,tLQR] = initial(clsys,x0,Tstop);

Compute the closed-loop response for the MPC controller with terminal weights.

SimOptions = mpcsimopt(MPCobj);
SimOptions.PlantInitialState = x0;
r = zeros(1,4);
[y,t,u] = sim(MPCobj,ceil(Tstop/Ts),r,SimOptions);
Cost = sum(sum(y(:,1:2)*diag(ywt).*y(:,1:2))) + sum(u*diag(uwt).*u);
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Compute the closed-loop response for a standard MPC controller with default prediction and control horizons ($p=10$, $m=3$). To match the other controllers, remove the default state estimator from the standard MPC controller.

MPCobjSTD = mpc(Plant,Ts);
MPCobjSTD.Weights.MV = uwt;
MPCobjSTD.Weights.OV = ywt;
setoutdist(MPCobjSTD,'model',tf(zeros(2,1)))
setEstimator(MPCobjSTD,[],C)
SimOptions = mpcsimopt(MPCobjSTD);
SimOptions.PlantInitialState = x0;
r = zeros(1,2);
[ySTD,tSTD,uSTD] = sim(MPCobjSTD,ceil(Tstop/Ts),r,SimOptions);
CostSTD = sum(sum(ySTD*diag(ywt).*ySTD)) + sum(uSTD*uwt.*uSTD);
-->The "PredictionHorizon" property of "mpc" object is empty. Trying PredictionHorizon = 10.
-->The "ControlHorizon" property of the "mpc" object is empty. Assuming 2.
-->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 
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Compare the controller responses.

figure
h1 = line(tSTD,ySTD,'color','r');
h2 = line(t,y(:,1:2),'color','b');
h3 = line(tLQR,yLQR,'color','m','marker','o','linestyle','none');
xlabel('Time')
ylabel('Plant Outputs')
legend([h1(1) h2(1) h3(1)],'Standard MPC','MPC with Terminal Weights','LQR','Location','NorthEast')

The MPC controller with terminal weights has a faster settling time compared to the standard MPC controller. The LQR controller and the MPC controller with terminal weights perform identically.

As reported in docid:mpc_ug.bthh9z6-1, the computed Cost value of 2.23 for the MPC controller with terminal weights is identical to the LQR controller cost. The cost for the standard MPC controller, CostSTD, is 4.82, more than double the value of Cost.

You can improve the standard MPC controller performance by adjusting the horizons. For example, if you increase the prediction and control horizons ($p=20$, $m=5$), the standard MPC controller performs almost identically to the MPC controller with terminal weights.

This example shows that using terminal penalty weights can eliminate the need to tune the prediction and control horizons for the unconstrained MPC case. If your application includes constraints, using a terminal weight is insufficient to guarantee nominal stability. You must also choose appropriate horizons and possibly add terminal constraints. For more information, see Rawlings and Mayne docid:mpc_ug.bu0qu79-1.