Clear Filters
Clear Filters

LQI(): integrator response not as expected

39 views (last 30 days)
John
John on 1 Mar 2023
Commented: Paul on 25 Mar 2023
I have a simple 2nd-order spring/mass/damper system in closed-loop configuration, that works as expected under LQR control; see first image.
I'm trying to use LQI() to add integral control -- and also my own discrete method -- but the integrator response, integral(e), is not as expected.
The 2nd image shows a ref input; I'm monitoring the integral state to see if it's working as expected, before augmented state-space with a disurbance input.
The 3rd image includes an augmente disturbance input + ref input, via [B F]. However, the integrator doesn't seem to accumulate when a disturbance is present at t = ~0.37. The integratorit acts as though the r-y error didn't exists.
See the attached .m file.
Mathworks' lqi() page -- the documentation is very sparse -- says the state vector x = [pos; vel; xi)], with xi = integral(error).
So, xi_dot = -x1+r --> xi = integral(r-y). This seems correct.
Also, xi seems correctly added to x2 (ie the input velocity_dot level) via a non-zero A(2,3) term, so the integral action should ultimately impact y as desired via x1_dot_dot = ...+A(2,3)*xi.
3 questions:
(1) Y isn't tracking to ref with a disturbance. Am I misunderstanding what this lqi() structure should do?
(2) x3 seems very small when ref steps to 1. Instantaneous r-y = 1, so since errors accumulate, x3 should be >1. Instead, it's 2.5e-3. What SS structure is lqi() assuming? Is the integrator set up correctly?
(3) Why is the computed LQI using lqi() gain so small at -10? To make a step that's about the same risetime as LQR (no integrator), i needed to multiply what lqi() comes up with, by ~10,000. Otherwise if using k_lqr(3) that lqi() came up with, the risetime was 130 seconds instead of 10 ms.
The latter is also confusing -- the way integral addition is included in state-space, it shouldn't counter-act the non-integral control response; it's only an addition. When starting from 0 state, wouldn't the integrator only be an additive error to get to ref faster? After that, the accumulator (int(e)) needs to wind down, so lag seems expected.
Either way, this makes me think that I don't understand what the lqi() gains actually are. For example, maybe the integrator is not added in a standard way, and ref signal must pass through the integrator before getting to the feedback K gains, so the plant response to a step will always be inherently slower than a non-int version with a ref step.

Accepted Answer

Paul
Paul on 2 Mar 2023
Before getting to the LQI part of this, I'm curious about the approach for making an LQR regulator into a tracking system.
Here is the original code, with some changes I made to illustrate some points.
%% Reference and Disturbance vectors
tf_s = 0.4;
ts_s = 1e-6;
t = (0 : ts_s : tf_s)'; % changed tout to t
u = ones(1,numel(t));
d = zeros(1,numel(t));
d(t > 0.25 * tf_s) = -1;
refStartTime = 0.1;
distPosTime = 0.25;
distVelTime = 0.5;
ref = zeros(numel(t), 1);
ref(t > refStartTime) = 1;
ref(t > 2*refStartTime) = 0;
ref(t > 3*refStartTime) = -1;
dPos = zeros(numel(t), 1);
% Force ramp for velocity disturbance
dVel = -10 * (t-distVelTime) .* (t > distVelTime);
% Both pos and vel dist (unused
dPosAndVel = dPos + dVel;
%% (1) Make 2nd-0 system
zeta = 0.3;
wn_rPs = 250;
ampl = 1;
m = 1e-5;
b = zeta * 2 * wn_rPs * m;
k = wn_rPs^2 * m;
Add the velocity state to the outputs.
A = [0 1; -k/m -b/m];
B = [0; 1/m];
C = eye(2); % output both states
D = [0];
Form the second order model. Don't normalize by the dcgain to illustrate a point later.
massSys = ss(A, B, C, D);
% massSys = massSys/dcgain(massSys); % got rid of dc gain adjustment
massSys.StateName = {'p', 'v'};
massSys.InputName = {'u'}; % changed to u
massSys.OutputName = {'pos','vel'};
%figure
%step(massSys);
%hold on;
Define the LQR weights. How were the values in Q and R selected?
%% (2) Make CL LQR
Q = [1000 0 ; 0 0.001];
R = [ 0.01 ];
[K_lqr, S1, P1] = lqr(massSys, Q, R);
% Stability check
A_CL_lqr = (massSys.A - massSys.B * K_lqr);
[T, Dd] = eig(A_CL_lqr);
Here is the original approach to forming the closed loop system. This approach assumes a feedback control law of the form:
u = -Kx + r, where r is the reference command.
sysClLqr = ss(A_CL_lqr, massSys.B, massSys.C, massSys.D);
However, for tracking problems it's probably better to form an error signal, in which case the control is
u = K1*(r - pos) - K2*vel = K1*([1;0]*r - x) = K1*[1;0]*r - K*x.
With this control we have
sysClLqr1= ss(A_CL_lqr, massSys.B*K_lqr*[1;0], massSys.C, massSys.D);
Compare the transfer functions of both systems from r to pos
zpk(sysClLqr(1,1))
ans = 1e+05 --------------------- (s+1001) (s+3.161e04) Continuous-time zero/pole/gain model.
zpk(sysClLqr1(1,1))
ans = 3.156e+07 --------------------- (s+1001) (s+3.161e04) Continuous-time zero/pole/gain model.
We see that the poles are the same, as must be the case, but the dc gain is different. One pole is at s = -1000 and the other is way out in the LHP, suggesting there is a very high bandwidth loop in the system, which might not be desirable. Are these poles consistent with design objectives? These pole locations must be related to the selection of Q and R, however that selection was made.
Compare dc gains
dcgain(sysClLqr(1,1))
ans = 0.0032
dcgain(sysClLqr1(1,1))
ans = 0.9980
The original approach has a dcgain on the order of 1e-3 and the latter is neary unity, as I think would be desired.
Another way to form the (new) closed-loop system would be
K = ss([K_lqr(1) -K_lqr(2)],'InputName',{'e' 'vel'},'OutputName',{'u'});
s = sumblk('e = r - pos');
sysClLqr2 = connect(massSys,K,s,'r','pos');
Verify its the same as above.
zpk(sysClLqr2)
ans = From input "r" to output "pos": 3.156e+07 --------------------- (s+1001) (s+3.161e04) Continuous-time zero/pole/gain model.
% sysClLqrDc = ss(A_CL_lqr, massSys.B/dcgain(sysClLqr), massSys.C, massSys.D, ...
% 'StateName', {'p', 'v'}, ...
% 'InputName', {'r(ref)'}, ...
% 'OutputName', {'pos'});
Now plot the response to the reference command
%% (3) Plot to verify
x0 = [0;0];
[y, tout, x] = lsim(sysClLqr2, ref, t, x0);
legStr = {'LQR', 'Pos Ref'};
subplot(3,1,1)
plot(tout, y(:,1), 'LineWidth', 1.5);
hold on;
plot(tout, ref, '--', 'LineWidth', 1.5)
grid on;
legend(legStr);
title('y');
subplot(3,1,2)
plot(tout, x(:,1), 'LineWidth', 1.5);
hold on;
title('pos');
ylabel('x1');
grid on;
subplot(3,1,3)
plot(tout, x(:,2), 'LineWidth', 1.5);
title('vel');
ylabel('x2');
grid on;
sgtitle('CL LQR: Ref response', 'FontSize', 14)
xlabel('Time [s]', 'FontSize', 14)
The response looks pretty good.
I guess the next step might be too look at the response of this system to the disturbance.
Now we could try LQI. But what is the motivation for trying LQI with this system? I'm not saying there isn't one, but I'm curious. Also probably need to get a better understanding how the Q and R matrices are selected for either LQR or LQI relative to actual design objectives.
  30 Comments
John
John on 24 Mar 2023
That's very helpful, thanks @Paul for the explanation. (Also I didn't know matrices could be created symbolically like that)
1) If rlocus gain is 1/R, then physically, why is the dominant CL pole (CL BW) becoming slower, as R gets lower (cheaper), after the gain causes CL poles to hit the Real axis? As you've mentioned before, one mode moves far left to -inf (fast), but the dominant pole slows down and moves right, towards the Im axis.
I'm not sure why easier control would slow things down -- unless it causes instability of some kind, but I thought LQI sets an optimum gainset that guarantees gain/phase margins (ie stability).
2) "We see that the zeros of G(s) are the roots of
H2*s^2 + H1*s - H3 = 0
Using the standard second order form
s^2 + 2*zeta*wn*s + wn^2
we see that we must have
H2 = 1, H1 = 2*zeta*wn, H3 = -wn^2."
Thanks. This is the clearest explanation i've seen of this (though i read a bit online).
How do you know to pick the standard 2nd-O form? Ie how do you know these CL zeros should correspond to a physical 2nd-O system's form, vs anything else? (as far as I know not every 2nd-O CL TF is a 2nd-O physical representation)
Paul
Paul on 25 Mar 2023
1) Are you referring to the rlocus plot in this comment? If so, the reason the dominant closed loop pole migrates to the right is becasuse it has to migrate toward that zero at s = -1000. By the rules of the root locus, the two branches have to meet on the real axis to the left of that zero, and then one root moves off to the left towards -inf and the other to the right towards the zero at s = -1000. That's the only direction it can go. You seem to be expecting that as R gets smaller both closed loop roots should migrate further to the left, but that's not how the math works when Q = H'H (H a row vector), in which case m roots will migrate to the m zeros of G(s) and the remaining n-m roots will migrate to infinity in a Butterworth pattern IIRC (n is the number of states in the plant). Depending on the selection of H, m could be anywhere from 0 to n-1.
If you're not referring to that comment, please clarify what you are referring to.
2. I picked the standard 2nd order form because that seemed to be the kind of response I thought you were trying to achieve, i.e., a well-damped second order response with a time constant faster than that of the plant. The zeros of G(s) aren't closed-loop zeros, unless it happens to be the case that H = C.

Sign in to comment.

More Answers (1)

Sam Chak
Sam Chak on 1 Mar 2023
If I'm not mistaken, the result is due to the disturbance d injected at the e_dot level, which I think it should be injected at the x_dot level.
Try if you get the expected result with this one:
F = [0; 1; 0];
  4 Comments
Sam Chak
Sam Chak on 3 Mar 2023
Hi @John, don't mention it. I also learn new things about control systems from the answers and comments posted by @Paul.

Sign in to comment.

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!