Very small u control values in LQR initial response

I have a system described by xdot=Adx+Bdu and y=Cdx+Ddu, where dx and du are the changes in the states and control variables, based on the matrices A, B, C, and D which contain the linearized dynamics of the system. The lqr controller is set up the following way:

The A, B, C, and D matrices:

A =

   -0.0204    0.0994   -0.1976         0  -50.0000   12.0000         0    7.6460    0.6088
   -0.0994    0.0006    0.3023   50.0000         0  -12.0000   -7.5940   -1.7931    9.5517
    0.1927   -0.3085   -0.0359  -12.0000   12.0000         0    1.3239   -5.7965   -2.1517
   -0.0159   -0.0054   -0.0201         0         0         0         0         0         0
    0.2922    0.6400    2.7503         0         0         0         0         0         0
    0.0170    0.0215    0.0906         0         0         0         0         0         0
         0         0         0    1.0000    0.3724    1.2039    0.2035    0.4002         0
         0         0         0         0    0.9553   -0.2955   -0.1546         0         0
         0         0         0         0    0.4754    1.5369    0.2598    0.3135         0


B =

    0.0186   -0.0473    0.0202    0.0100
         0         0   -0.0316   -0.0817
    0.0903    0.0792    0.0843    0.0416
    0.0544   -1.1601    0.5681   -1.4700
   -0.0907   -1.9335    1.5596   -0.7709
   -1.6589    1.1871   -0.9469   -2.4499
         0         0         0         0
         0         0         0         0
         0         0         0         0

C =

   1.0e+03 *

   -0.0159   -0.0054   -0.0201         0         0         0         0         0         0
    0.2922    0.6400    2.7503         0         0         0         0         0         0
    0.0170    0.0215    0.0906         0         0         0         0         0         0


D =

   1.0e+03 *

    0.0544   -1.1601    0.5681   -1.4700
   -0.0907   -1.9335    1.5596   -0.7709
   -1.6589    1.1871   -0.9469   -2.4499

The lqr controller:

Q=1e-3*eye(9,9);
% Q=rand(10,10);
R=1e10*eye(4,4);
% N=eye(10,13)+1e-5;
% N=zeros(10,13);
% R=rand(13,13);
K=lqr(A,B,Q,R);
x0=ones(9,1);
sys=ss((A-B*K),B,C,D);
t=0:1/10:30;
[y,t,dx]=initial(sys,x0,t);

I then retrieve the values of du and apply it to u the following way:

u0 =

   -0.1014
    0.2692
   -0.6107
   -0.2415

   
for ii=1:length(t)
    du(:,ii)=pinv(D)*(y(ii,:)'-C*(dx(ii,:))');
end

u=u0+cumtrapz(t,du(1:4,:)')';

This then gives me the attached results. The controller returns to error, and the y values reach steady state of 0 after 15 or so seconds. The issue however Is that the control values seem incredibly small. In fact, the total change for the control values are in the order of 1e-14 units of control. I have tried also setting the R matrix to be larger. I understand that it may just be that "it is what it is", but I am wondering if I have defined the controls incorrectly or there is a output option for the initial function that gives the control?

[Alessandro Maria Laspina]

I've noticed I may have been calculating the control variables incorrectly.

For an lqr system, the u variables (or du in this case) can be found by:

x=-K*x (mentioned in this video: What Is Gain Scheduling? | Control Systems in Practice, Part 2 - YouTube)

I am not entirely sure, so a confirmation would be appreciated.

Answer 1

[Paul]

For this dicussion, I'm going to just use x and u as the symbols for the state and control variables of the linearized model, instead of dx and du, but with the understand that x and u represent pertrubations to the state and control around which the system is linearized.

We start with the linear model:

xdot = A*x + B*u

y = C*x + D*u

A =[
   -0.0204    0.0994   -0.1976         0  -50.0000   12.0000         0    7.6460    0.6088
   -0.0994    0.0006    0.3023   50.0000         0  -12.0000   -7.5940   -1.7931    9.5517
    0.1927   -0.3085   -0.0359  -12.0000   12.0000         0    1.3239   -5.7965   -2.1517
   -0.0159   -0.0054   -0.0201         0         0         0         0         0         0
    0.2922    0.6400    2.7503         0         0         0         0         0         0
    0.0170    0.0215    0.0906         0         0         0         0         0         0
         0         0         0    1.0000    0.3724    1.2039    0.2035    0.4002         0
         0         0         0         0    0.9553   -0.2955   -0.1546         0         0
         0         0         0         0    0.4754    1.5369    0.2598    0.3135         0];
B =[
    0.0186   -0.0473    0.0202    0.0100
         0         0   -0.0316   -0.0817
    0.0903    0.0792    0.0843    0.0416
    0.0544   -1.1601    0.5681   -1.4700
   -0.0907   -1.9335    1.5596   -0.7709
   -1.6589    1.1871   -0.9469   -2.4499
         0         0         0         0
         0         0         0         0
         0         0         0         0];
C =   1.0e+03 * [
   -0.0159   -0.0054   -0.0201         0         0         0         0         0         0
    0.2922    0.6400    2.7503         0         0         0         0         0         0
    0.0170    0.0215    0.0906         0         0         0         0         0         0];
D =   1.0e+03 * [
    0.0544   -1.1601    0.5681   -1.4700
   -0.0907   -1.9335    1.5596   -0.7709
   -1.6589    1.1871   -0.9469   -2.4499];

Now we introduce the state feedback control

u = -K*x

Q=1e-3*eye(9,9);
R=1e10*eye(4,4);
K=lqr(A,B,Q,R);

and substitute:

xdot = A*x + B*(-K*x) = (A - B*K)*x

y = C*x + D*(-K*x) = (C - D*K)*x

So the resulting closed loop system is:

xdot = (A - B*K)*x

y = (C - D*K)*x

Note that the closed loop system no longer has a "B" or "D" matrix because the system is not being driven by an external input, like a reference command (which can be introduced if needed, but it's a different input than u, which is the control input that actually drives the plant). So the closed loop model in Matlab should be:

sys = ss(A-B*K,zeros(size(A,1)),C-D*K,0); % contrast with the code in the Question, especially the C-argument

Note that I just made the "B" and "D" matrices 0 because we have to put something in the ss model. Note also that ss automatically expands the scalar 0 in the D-argument to appropriate dimensions.

Now to simulate the IC response

x0=ones(9,1);
t=0:1/10:30;  % might want to consider a smaller dt?
[y,t,x]=initial(sys,x0,t);

To recover the control input that actually was driving the plant dynamics, the approach in your comment is correct:

u = (-K*x.').';
subplot(211);plot(t,y),grid
subplot(212);plot(t,u),grid

[Alessandro Maria Laspina]

Thank you for the response. One thing I don't understand is what why the system is formed as:

sys = ss(A-B*K,zeros(size(A,1)),C-D*K,0);

and not as the one I wrote before.

[Paul]

I had a typo in my answer that I just corrected that might clarify why sys is defined that way. Is the portion of the answer between "and substitute" and "Now to simulate" still unclear?

[Alessandro Maria Laspina]

Yes, sorry, it should already have been clear I did not devote enough attention to your answer. Thank you very much for the help.

I may open up a new question to ask about this, but in any case, the controller objective is to reduce the values of y to zero, correct?

If I wanted it to set the values of y to something that was non-zero, which are constants, would the most straightforward way be to introduce the non-zero objectives as noise? What types of controllers can be used in MATLAB to build that?

[Paul]

Actually, the LQR controller objective is to drive the states to zero in a way that minimizes the LQR objective function, which essentially trades off the manner in which the states are regulated to zero with the control effort it takes to do so. Of course, if the states go to zero, then so will the control input, and so will the output.

Having said that, if the plant has D = 0 (your's does not), then you can select

Q = C.'*C % assuming C is such that Q is positive definite

in which case

x.' * Q * x = y.' * y

i.e., the LQR cost function is actually a specification on the output response.

If you want the output, y, to track some reference signal r, one approach (there may be others) would be to augment the plant model with integrators, and make each element of y an input to an integrator. Then design the LQR controller for the augmented plant. Then introduce the reference input r so that the input to each integrator becomes the error in y. effectively, the control scheme becomes:

u = -K1*x + K2/s*(r - y)

where K1 and K2 are partitions of the LQR feeback gain matrix K and assumes that full state vector x and the output vector y are measured for feedback. It's easier to explain in a picture than in words. It may be easier to understand for the case where the ouput, y, is one or more state variables, but that's not the case in the problem in the Question.

[Alessandro Maria Laspina]

How is K partitioned into K1 and K2? and what is s?

[Alessandro Maria Laspina]

Also, how could I add control contraints and control rate constraints?

[Paul]

The augmented plant model will have 12 states. It can be structured so that its states are in the following order

xa = [x; xi]

where xi are the integrator outputs.

The LQR feedback gain is partitioned as

K = [K1 K2]

where K1 has 9 columns and K2 has 3 columns. The s in the equation is as used in the Laplace transform, so I should have written that the control input expressed in the frequency domain is

U(s) = -K1*X(s) + K2/s*(R(s) - Y(s))

Recall that in the frequency domain 1/s represents integration. So it's basically taking the integrators from the augmented plant and implemeting them in the controller.

Anyway, I'm sure an internet search on something like "LQR with integral control" will turn up some good references. Or you might want to start with:

doc lqi

I'm not sure what yout're trying to do with control and control rate constraints so can't comment on that.

[Alessandro Maria Laspina]

My states are perturbed states, so I would need some form of integral constraints on the on the control, i.e. my state-space is:

xdot=Adx+Bdu;

y=Cdx+Bdu;

I don't have any constraints on the states (dx) or the integral of the states (x), or the derivative (rates, which in this case they would be accelerations)

What are my integrator outputs, and why would they be 3?

I am simply trying to limit the values of the integral of my control du. How can I specify the integral of the control for an lqi controller?

[Paul]

There are three additional integrators because the plant has three outputs and I assumed you wanted each to track its own reference signal.

As for the other questions, they've become somewhat far afield of Matlab itself. Again, suggest you undertake a search for good references on LQ control as it applies to the problem(s) that you are trying to solve.