Stability margins for a cascaded system with getLoopTransfer functions

Hello all, I have tuned the control system above and it's stable. I'm trying to find the (classical) stability margins at the signals highlighted in red, and I'm still a bit inexperienced with this. How would I go about this?

I have calculated those at e as follows:

Le = getLoopTransfer(CL,'e',-1); % open-loop function at e
[margins] = allmargin(Le);

which is I think the correct way. I am more unsure however of q and , since there are inner-loops with positive feedback and there is the added complexity of breaking the loops of the outer loop(s). I found that breaking the loop at both and and only at gives the same transfer function when evaluating q. Apologies for inconsistency, in the code below and

%% Verification of margins at plant input 
L_P_input = getLoopTransfer(CL,'u_c',-1,{'y2'});
zpk(L_P_input)
eig(L_P_input)
figure; nyquist(L_P_input)
[margins_input] = allmargin(L_P_input);   % Margins at the plant input
%%

%% Verification of margins at q
L_P_out_q =  getLoopTransfer(CL,'y1',-1,{'y2'});
zpk(L_P_out_q)
eig(L_P_out_q)
[margins_q] = allmargin(L_P_out_q); % Margins at the output of q
figure; nyquist(L_P_out_q)

Futhermore, I also saw that the OL transfer functions for and q are exactly the same (when I break the outer loop at ). Anyone could clarify this?

Lastly, is it enough to find the stability margins at the highlighted locations? I reckoned it wouldnt be necessary to also do it at since the error is already on this signal.

I am wondering if my approach is correct and if not what should be changed.

Thanks in advance for any help!

Answer 1

[Paul]

Hi Daam,

To get the stability margins at 'e', I believe this is correct:

Le = getLoopTransfer(CL,'e',-1); % open-loop function at e
[margins] = allmargin(Le);

To get the stability margins at 'u_c', I believe this is correct.

L_P_input = getLoopTransfer(CL,'u_c',-1);
[margins_input] = allmargin(L_P_input);   % Margins at the plant input

Note that I removed the loop break at 'y2' which is 'a_n'. You could compute the margins at u_c with the loop broken at 'a_n', but the real system has the loop closed at a_n.

Similary, at 'q'

L_P_out_q =  getLoopTransfer(CL,'y1',-1);
[margins_q] = allmargin(L_P_out_q); % Margins at the output of q

If in doubt, it's always easy to verify the gain margins (and not so difficult to verify the phase margins). For example, rebuild your system with the nominal value of Kq multiplied by a gain margin from margins_q.GainMargin. The closed loop system should have poles (or one pole at the origin) on the imaginary axis at the corresponding margins_q.GMFrequency.

"I found that breaking the loop at both a_n and Kp*a_n and only at a_n gives the same transfer function when evaluating q."

Breaking the loop at a_n eliminates the feedback loops through Kp and Ki. So additionally breaking the loop at the output of Kp doesn't change anything. In either case, the loop transfer at q

L_q = getLoopTransfer(CL,'y1',-1,{'a_n'});

should be L_q(s) = -Kq*G21(s).

"Futhermore, I also saw that the OL transfer functions for u_c and q are exactly the same (when I break the outer loop at a_n)."

When you open the loop at a_n, the two outer loops are eliminated when computing the open loop transfer at either u_c or q. The lone remaining loop contains both q and u_c, so the open loop tranfer fucntion at either should be L(s) = -Kq*G21(s).

If you save G, Kp, Ki, Kq in a .mat file and upload that file using the paperclip icon on the Insert menu, then we can verify all of the above, which at this point is based only on my view of the diagram.

"Lastly, is it enough to find the stability margins at the highlighted locations? I reckoned it wouldnt be necessary to also do it at a_n since the error is already on this signal."

As discussed above, opening the loop at a_n is actually breaking two feedback loops and the resulting open loop transfer fucntion is not the same as it would be at e. One apprach is to find the stability margins at the control inputs and feedback outputs, i.e., u_c, q, and a_n. However, you can, of course, compute stabilty margins at any point in the system should you feel there is a need to do so.

[Daam Schoon]

Hi Paul, thanks a lot for helping!

[Q1:] I've verified the gain margin at e by applying the gain around gain margin and checking stability. However, for the other signals the gain margin was Inf, so I couldn't really verify it by gain, I believe (although they could handle infinite gain), and I didn't manage to add pure phase to the control system and still have slLinearizer manage to linearize it. I tried adding a transport delay block, but this didn't seem to change the transfer functions at all, and I tried a Matlab Fcn block which adds the following as a block

exp(1j*deg2rad(Phase));

but this gave an error when linearizing it due to the complex component. What would you suggest to do to add pure phase to the signals to test phase margins?

The latter three parts of your answer make sense to me, thank you for that!

[Q2:] For my own understanding, why do we need to specify -1 in the getLoopTransfer for the inner-loop q and u_c?. Aren't they positive feedback loops?

[Q3:] When closing the outer loops for u_c and q, there is a pole in the origin for u_c leading to no gain margin. I was wondering if you comment on this, what it's caused by (maybe the integrator after e?) and what the implications of this are.

I've attached the required variables. Again, thank you for helping!

Daam

[Daam Schoon]

Here is also the simplified simulink file to test margins to save you some time!

[Paul]

Create the model

s = load('G,Kp,Ki,Kq.mat');
G = s.Gswitched;
G.InputName = 'u_c';
G.OutputName = {'a_n' 'q'};
Ki = ss(s.Ki,'InputName','epsilon','OutputName','ui');
Kp = ss(s.Kp,'InputName','a_n','OutputName','ua');
Kq = ss(s.Kq,'InputName','q','OutputName','uq');
CL = connect(...
    sumblk('e = r - a_n'), ...
    tf(1,[1 0],'InputName','e','OutputName','epsilon'), ...
    Ki, ...
    Kp, ...
    Kq, ...
    sumblk('u_c = ui + ua + uq'), ...
    G, ...
    {'r'}, ...
    {'a_n'}, ...
    {'a_n','q','e','u_c'});

Check the step response

step(CL)

Verify some claims made above.

Open loop TF at q, with loop also broken at a_n.

L_q = getLoopTransfer(CL,'q',-1,{'a_n'});
tf(minreal(L_q))

ans =
 
  From input "q" to output "q":
          38.89 s + 40.02
  -------------------------------
  s^3 + 22.3 s^2 + 42.7 s + 7.142
 
Continuous-time transfer function.

Show that L_q = -Kq*G(2,1).

tf(-Kq*G(2,1))

ans =
 
  From input "u_c" to output "uq":
          38.89 s + 40.02
  -------------------------------
  s^3 + 22.3 s^2 + 42.7 s + 7.142
 
Continuous-time transfer function.

which should also be the same as the open loop transfer function at u_c with the loops also opened at a_n

L_u = getLoopTransfer(CL,'u_c',-1,{'a_n'});
tf(minreal(L_u))

ans =
 
  From input "u_c" to output "u_c":
          38.89 s + 40.02
  -------------------------------
  s^3 + 22.3 s^2 + 42.7 s + 7.142
 
Continuous-time transfer function.

Q1: Verifying stablity margins, done here at e but applicable to any other point in the system.

format short e
margins_e = allmargin(getLoopTransfer(CL,'e',-1))

margins_e = struct with fields:
     GainMargin: [9.8313e+00 7.5615e+01 Inf]
    GMFrequency: [3.2637e+00 6.5644e+00 Inf]
    PhaseMargin: 6.6146e+01
    PMFrequency: 7.4907e-01
    DelayMargin: 1.5412e+00
    DMFrequency: 7.4907e-01
         Stable: 1

Verify the gain margin by putting a gain of 9.8313 on the integrator and check the eigenvalues of the close loop system

CLgain = connect(...
    sumblk('e = r - a_n'), ...
    tf(margins_e.GainMargin(1),[1 0],'InputName','e','OutputName','epsilon'), ...
    Ki, ...
    Kp, ...
    Kq, ...
    sumblk('u_c = ui + ua + uq'), ...
    G, ...
    {'r'}, ...
    {'a_n'}, ...
    {'a_n','q','e','u_c'});
format short e
eig(CLgain)

ans = 
  -1.8516e+01 + 0.0000e+00i
   8.3606e-05 + 3.2636e+00i
   8.3606e-05 - 3.2636e+00i
  -5.3552e+00 + 0.0000e+00i

As expected, closed loop eigenvalues on the imaginary axis at the expected frequency

margins_e.GMFrequency(1)

ans = 
   3.2637e+00

For the phase margin, we have to determine the time delay that corresponds to -66.146 deg of phase at the gain crossover frequency (margins_e.PMFrequency) of 0.74907 rad/sec. In other words, solve for T such that

angle(exp(-1j*w*T)) = -66.146*pi/180 for w = 0.74907

T = pi/180*margins_e.PhaseMargin(1)/margins_e.PMFrequency(1)

T = 
   1.5412e+00

Note that is also the DelayMargin.

Now we add that time delay to the model by adding an InputDelay to the integrator, which would be accomplished by using a Transport Delay block in Simulink (make sure to select an appropriate value for the 'Pade order' parameter if going to use Simulink linearization tools).

CLphase = connect(...
    sumblk('e = r - a_n'), ...
    tf(1,[1 0],'InputName','e','OutputName','epsilon','InputDelay',T), ...
    Ki, ...
    Kp, ...
    Kq, ...
    sumblk('u_c = ui + ua + uq'), ...
    G, ...
    {'r'}, ...
    {'a_n'}, ...
    {'a_n','q','e','u_c'});

Use a Pade approximant of order 3 for the delay to compute the closed loop poles

format short e
pole(pade(CLphase,3))

ans = 
  -1.9760e+01 + 0.0000e+00i
  -4.4986e+00 + 2.7228e+00i
  -4.4986e+00 - 2.7228e+00i
   3.2028e-05 + 7.4910e-01i
   3.2028e-05 - 7.4910e-01i
  -1.4499e+00 + 3.0631e+00i
  -1.4499e+00 - 3.0631e+00i

As predicted, poles on the imaginary axis at the expected frequency

margins_e.PMFrequency(1)

ans = 
   7.4907e-01

Q2: getLoopTransfer, without the sign input (or sign = +1) computes the open loop transfer function, LT(s), including all of the signs at the sum junctions. By Mason's Rule, the closed loop characteristic equation is:

D = 1 - LT(s)

But the Nyquist criterion is based on a characteristic equation of the form

D = 1 + L(s)

or L(s) = -LT(s).

Therefore, to apply the Nyquist criterion to determine stability, and then compute GM and PM if the closed loop system is stable, we need to use the sign = -1 to get the L(s) that we need, regardless of the polarity of the feedback loop(s).

For example, find the open loop system at q, use sign = -1 even though that loop has positive feedback

L_q = getLoopTransfer(CL,'q',-1);

Get the stability margins

margin_q = allmargin(L_q)

margin_q = struct with fields:
     GainMargin: [4.5981e+17 Inf]
    GMFrequency: [0 Inf]
    PhaseMargin: [-1.2092e+02 1.2161e+02]
    PMFrequency: [1.3617e+00 2.7485e+00]
    DelayMargin: [3.0644e+00 7.7223e-01]
    DMFrequency: [1.3617e+00 2.7485e+00]
         Stable: 1

Time delay corresponding to the high-frequency decreasing phase margin at 2.7485 rad/sec

T = pi/180*margin_q.PhaseMargin(2)/margin_q.PMFrequency(2)

T = 
   7.7223e-01

which is the DelayMargin.

Add that delay to Kq, build the system, and get the closed loop poles

KqWithDelay = Kq;
KqWithDelay.InputDelay = T;
CLphaseAtq = connect(...
    sumblk('e = r - a_n'), ...
    tf(1,[1 0],'InputName','e','OutputName','epsilon'), ...
    Ki, ...
    Kp, ...
    KqWithDelay, ...
    sumblk('u_c = ui + ua + uq'), ...
    G, ...
    {'r'}, ...
    {'a_n'}, ...
    {'a_n','q','e','u_c'});

pole(pade(CLphaseAtq,3))

ans = 
  -2.6495e+01 + 0.0000e+00i
  -5.4464e+00 + 9.1957e+00i
  -5.4464e+00 - 9.1957e+00i
  -5.6000e-04 + 2.7495e+00i
  -5.6000e-04 - 2.7495e+00i
  -1.0106e+00 + 4.4982e-01i
  -1.0106e+00 - 4.4982e-01i

As expected, two poles on the imaginary axis 2.7495 rad/sec.

Q3: Yes, the open loop transfer function at u_c breaks all of the loops, so the integrator results in an open loop pole at the origin.

L_u = getLoopTransfer(CL,'u_c',-1);
zpk(minreal(L_u))

ans =
 
  From input "u_c" to output "u_c":
  1.5745 (s+29.2) (s^2 + 3.033s + 2.337)
  --------------------------------------
      s (s+0.185) (s+1.911) (s+20.2)
 
Continuous-time zero/pole/gain model.

The Bode plot shows that there is no phase crossover (i.e., no gain margin), but I'm not sure I'd say that pole at the origin is "leading to no gain margin." That pole at the origin adds 90 deg of phase lag across all frequencies; it's the other compensation paths that shape the remainder of the phase response and keep it away from -180 deg. Not sure what else you might be looking for regarding "implications."

bode(L_u)

[Daam Schoon]

This has been an eye-opener for me. Why to choose this sign of the getLoopTransfer makes sense to me now, as well as the method of verifying the gain and phase margins.

Thanks a lot for your extensive help, Paul!

[Paul]

You're very welcome.