LQR for MPC with input rate weight

Referring to MATLAB's documentation https://www.mathworks.com/help/mpc/ug/using-terminal-penalty-to-provide-lqr-performance.html# . Different to the example, my MPC has zero input weight, R (versus 1 in the example) and 1.83 input rate weight, mpcobj.Weights.MVRate (versus 1e-5). I followed the step-by-step instruction, but I couldn't get the value of DC gain of mpcobject and K of LQR to be the same (as I should get if I follow the instruction). This is the code I used:

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);

Q = eye(2);

R = 0;

[K,Qp] = lqry(plant,Q,R);

L = chol(Qp);

newPlant = plant;

set(newPlant,'C',[C;L],'D',[D;zeros(2,1)]);

newPlant = setmpcsignals(newPlant,'MO',[1 2],'UO',[3 4]);

p = 3;

m = p;

mpcobj = mpc(newPlant,Ts,p,m);

ywt = sqrt(diag(Q))';

uwt = sqrt(diag(R))';

mpcobj.Weights.OV = [ywt 0 0];

mpcobj.Weights.MV = uwt;

mpcobj.Weights.MVRate = 1.83;

Y = struct('Weight',[0 0 1 1]);

U = struct('Weight',uwt);

setterminal(mpcobj,Y,U);

setoutdist(mpcobj,'model',ss(zeros(4,1)));

setEstimator(mpcobj,[],eye(2));

mpcgain = dcgain(ss(mpcobj));

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));

I suspect that this is due to the LQR isn't considering the input rate in its cost function calculation. I tried modifying R to be the value of 1.83^2 [(mpcobj.Weights.MVRate)^2] but the gains (K and mpcgain) are still not the same. Please advise the correct method to get the equivalent LQR for my MPC. Thank you.

[Paul]

Are you sure that there actually is an equivalent LQR? If you compare mpcgain and K with full precision with mpcopb.Weights.MVRate = 1e-5 as in the doc, you'll see that they are not equal, presumably because the MVRate is non-zero. As you say, the LQR cannot penalize the change in the control,

If you want to use LQR to penalize the rate of the control, you can augment the plant with an integrator so that the control input for the LQR design is udot, but then you end up with 1x3 feedback gain matrix and feed the compensation includes feedback of the true control input, u.

[aurochs]

Thank you for your quick response Paul.

--Are you sure that there actually is an equivalent LQR?--

Do you mean that there's no equivalent LQR for my MPC?

--If you compare mpcgain and K with full precision with mpcopb.Weights.MVRate = 1e-5 as in the doc, you'll see that they are not equal, presumably because the MVRate is non-zero. As you say, the LQR cannot penalize the change in the control--

Yes, you're correct.

--If you want to use LQR to penalize the rate of the control, you can augment the plant with an integrator so that the control input for the LQR design is udot, but then you end up with 1x3 feedback gain matrix and feed the compensation includes feedback of the true control input, u.--

I'm sorry, I'm a bit new to this..can you kindly help to elaborate more? Does it mean that if my previous x matrix is x=[x1;x2], then my augmented x would be xnew=[x1;x2;u]? And if my previous u is u1, then unew=Δu1? How will I be able to get K of u=-Kx from this augmented plant, since the control now is Δu (which means Δu=-Kx instead of u=-Kx)?

Your kind respsonse is much appreciated.

[Paul]

I don't know enough about MPC to answer that first question.

As for the plant augmentation approach ... yes, conceptually, the state vector of the augmented plant would be [x1 ; x2; u] and the control input to the augmented plant would be udot. However, for a discrete time implementation of an integrator, H(z), the output of H(z) typically is not the internal state, so in discrete time the augented state vector is [x1; x2; z], where z is the state variable of the integrator difference equation, and the the control input is the discrete equivalent of udot, which I'm going to still call udot for simplicity. With the augmented plant, the LQR control is udot = -K * x = - K * [x1 ; x2; z]. But for implementation with the actual plant, z is pulled back into the compensation, so the compensation will be dynamic, rather than just feedback gains. Here is the code to illustrate.

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);
integrator = ss(c2d(tf(1,[1 0]),Ts,'impulse')); 
augmentedplant = series(integrator,plant);
set(augmentedplant,'C',eye(3),'D',zeros(3,1));
Q = blkdiag(eye(2),0);
R = 1.83;
[K,Qp] = lqry(augmentedplant,Q,R);
L = chol(Qp);
newPlant = augmentedplant;
set(newPlant,'C',[augmentedplant.C;L],'D',[augmentedplant.D;zeros(3,1)]);
newPlant = setmpcsignals(newPlant,'MO',[1 2 3],'UO',[4 5 6]);
p = 3;
m = p;
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 y5 y6 

ywt = sqrt(diag(Q))';
uwt = sqrt(diag(R))';
mpcobj.Weights.OV = [ywt 0 0 0];
mpcobj.Weights.MV = uwt;
mpcobj.Weights.MVRate = 1e-5; %1.83;
Y = struct('Weight',[0 0 0 1 1 1]);
U = struct('Weight',uwt);
setterminal(mpcobj,Y,U);
setoutdist(mpcobj,'model',ss(zeros(6,1)));
setEstimator(mpcobj,[],eye(3));
mpcgain = dcgain(ss(mpcobj));

-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

fprintf('\n(unconstrained) MPC: u(k)=[%8.8g,%8.8g,%8.8g]*x(k)',mpcgain(1),mpcgain(2),mpcgain(3));

(unconstrained) MPC: u(k)=[-1.6719771,-0.67161995,-0.69814305]*x(k)

fprintf('\n                LQR: u(k)=[%8.8g,%8.8g,%8.8g]*x(k)\n\n',-K(1),-K(2),-K(3));

                LQR: u(k)=[-1.6719771,-0.67161995,-0.69814305]*x(k)

x0 = [0.2;0.2;0];  % initialize integrator to zero

% closed loop response of augmented plant and LQR gains
clsys = feedback(augmentedplant,K);
Tstop = 6;
[yLQR,tLQR] = initial(clsys,x0,Tstop);

% the actual implemtation would have the integrator as part of the compensation, not the plant
compensator = series(ss(K(1:2)),ss(integrator.a - integrator.b*K(3) , integrator.b , integrator.c - integrator.d*K(3) , integrator.d, Ts)); 
clsysimplemented = feedback(plant,compensator);
[yimp,timp] = initial(clsysimplemented,x0,Tstop);

% close loop repsonse of MPC controller
simOpt = mpcsimopt(mpcobj);
simOpt.PlantInitialState = x0;
r = zeros(1,6);
[ympc,tmpc,umpc] = sim(mpcobj,ceil(Tstop/Ts),r,simOpt);

% compare
plot(tLQR,yLQR,timp,yimp,'x',tmpc,ympc(:,1:3),'o')

[aurochs]

Thanks again Paul for the detailed explanation.

I learned a lot from your post even though I still have some questions unanswered.

Can you help to verify if I understand your answer correctly?

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);
integrator = ss(c2d(tf(1,[1 0]),Ts,'impulse')); 
augmentedplant = series(integrator,plant);
set(augmentedplant,'C',eye(3),'D',zeros(3,1));
Q = blkdiag(eye(2),0);
R = 1.83;
[K,Qp] = lqry(augmentedplant,Q,R);

The control signal would be Δu because of the embedded integrator (Δu=-Kx). Δu cost function weight is 1.83.

L = chol(Qp);
newPlant = augmentedplant;
set(newPlant,'C',[augmentedplant.C;L],'D',[augmentedplant.D;zeros(3,1)]);
newPlant = setmpcsignals(newPlant,'MO',[1 2 3],'UO',[4 5 6]);
p = 3;
m = p;
mpcobj = mpc(newPlant,Ts,p,m);
ywt = sqrt(diag(Q))';
uwt = sqrt(diag(R))';
mpcobj.Weights.OV = [ywt 0 0 0];
mpcobj.Weights.MV = uwt;
mpcobj.Weights.MVRate = 1e-5; %1.83;
Y = struct('Weight',[0 0 0 1 1 1]);
U = struct('Weight',uwt);
setterminal(mpcobj,Y,U);
setoutdist(mpcobj,'model',ss(zeros(6,1)));
setEstimator(mpcobj,[],eye(3));
mpcgain = dcgain(ss(mpcobj));

The control signal of mpcobj is Δu as well, instead of u, because it uses the same prediction model as LQR (original plant plus integrator). Δu cost function weight is 1.83.

% closed loop response of augmented plant and LQR gains
clsys = feedback(augmentedplant,K);
Tstop = 6;
[yLQR,tLQR] = initial(clsys,x0,Tstop);

% the actual implemtation would have the integrator as part of the compensation, not the plant
compensator = series(ss(K(1:2)),ss(integrator.a - integrator.b*K(3) , integrator.b , integrator.c - integrator.d*K(3) , integrator.d, Ts)); 
clsysimplemented = feedback(plant,compensator);
[yimp,timp] = initial(clsysimplemented,x0,Tstop);

Can you help to explain how the compensator implementation different to the LQR implementation, since it looks to me they are both using the same control signal with LQR gains (K versus K(1:2)) to control the same plant (augmentedplant versus integrator in series with the original plant).

Your kind response is very much appreciated.

[Paul]

Yes, you're clear on the first two parts.

As to the third part ... the clsys is the result of wrapping the static feedback gains around the augmented plant. However, the augmented plant contains that integrator, which is not part of the real plant. To control the real plant, which does not include an integrator, we need a compensator that inlcudes the integrator and the feedback of that integrator through K3 in order to realize the same loop. It's just a block diagram manipulation. If still unclear on why the compensator is defined the way it is, I suggest you draw the block diagram with the discrete integrator in series with the plant, and then add K1, K2, and K3 feedbacks. Draw a box around the plant, then show that the remainder of the diagram, which is the compensator, must be defined as I've shown.

Note that the "control" input to the augmented plant is delta-u, but the control input to the real plant is u, which is the output of the compensator.

[aurochs]

Thanks Paul! Really appreciate your effort in explaining the answer.