Stability Margins for MIMO Feedback Loop

Hello,

I have this NDI-Controller for the lateral motion of an aircraft modelled in Simuink and want to determine its stabiity margins. I am strugging in how to linearize the open feedback loop. I placed an 'looptransfer' analysis point after my Control Allocation, hence the control input. When I now linearize the open-loop with the Model Linearizer, and determine the margins in MATLAB with 'allmargin(-linsys)', it says, that my open-loop is not stable.

For me it seems, like I am using the Matlab Linearizer not correctly, as my outputs of the feedback loop seem stable. The same to all the poles.

Does anyone have an idea, how I should linearize the feedback loop to receive the correct gain & phase margins? Would be so helpful.

Thanks a lot.

Cheers,

Flo

[Paul]

allmargin doesn't tell anything about open-loop stability, only closed loop stability via the Stable field of the output.

I think your plant has 2 inputs, so the result from allmargin should be a two-element structure. The Stable field is zero for both elements of the allmargin output?

You might get more help if you upload linsys in a .mat file using the Paperclip icon on the Insert menu.

[Florian]

Yes, the result of allmargin is a two-element structure, where the stable field is zero for both elements. Attached is the resulting linsys.mat, which I receive if I perform a linearization with the depcited signal marked as 'looptransfer'.

Answer 1

[Sam Chak]

Hi @Florian

If you wish to determine the stability margin of a linearized system, you may consider using the Bode Plot block directly. Please refer to the example provided. The water tank system is nonlinear, and the Bode Plot block can perform linearization at the specified input-output points. However, your aircraft model is somewhat unconventional, as the Plant's output state vector x is fed back to the Reference Model. Should the Linearizer interpret this as an external feedback loop?

https://www.mathworks.com/help/slcontrol/ug/visualize-bode-response-of-simulink-model-during-simulation.html

[Florian]

Thanks for the hint, that makes linearizing easier. But, it does not change my result. Is it correct to mark the output of the controller as 'looptransfer'?

The state vector x is fed into the reference model, as accoring to the current bank angle (which is stored in the state vector), the commanded yaw rate is calculated

[Sam Chak]

Hi @Florian

I am not familiar with using the allmargin() function to analyze the stability of a nonlinear control system. If I were to perform linearization in Simulink, I would choose the I/O signals approach because it allows me to know exactly what to expect. Also, if the commanded yaw rate depends on the bank angle (one of the state variables), it should probably be placed in the Error Controller subsystem. Reference signals should be independent of the states.

Are you required to conduct stability analysis on the nonlinear control system? You may already be aware that the stability indicated by allmargin() derived from the linearized model only provides information about local stability, as the linear approximation is valid only within a limited range around a specific operating point. It is also important to consider: "Does the aircraft practically have only a specific operating point?"

If this is for a journal article, you need to plan and decide what to report regarding stability, as it is generally insufficient to claim that the designed nonlinear control system is stable solely based on the indications from MATLAB's allmargin(). Some people have used the MIMO Nyquist Criterion; however, it can potentially lead to incorrect stability analysis as shown in this article.

[Paul]

Hi Sam,

Is linsys1 at the end of your comment that which was posted by @Florian? It doesn't seem to be.

load linsys
S = allmargin(-linsys1);
S.Stable

ans = logical
   0

ans = logical
   0

[Sam Chak]

Hi @Paul

That was a snapshot from your post. I understand that it is acceptable to use allmargin() to investigate stability. However, instead of merely showing the logical value for the stability of NDI-controlled aircraft in the journal paper, I encourage @Florian to conduct a more comprehensive analysis, either in the form of graphs or through a rigorous mathematical proof.

load linsys
linsys1 = minreal(linsys1);

7 states removed.

S = allmargin(-linsys1);
S.Stable

ans = logical
   1

ans = logical
   1

[Florian]

Alright, nevertheless I don’t use allmargin for the stability analysis. I just wanted to check and compare the margins of different controller, but was wondering, why it said Stable=0, as I didn’t expected the loop to be unstable.

[Paul]

@Sam Chak

Thanks for brining that paper to my attention. I only skimmed it, but I think its premise is highly misleading.

I don't think the problem is with the MIMO NC per se. Rather, it seems to be with how the authors attempted to implement the MIMO NC in Matlab (I'm not aware of any MathWorks-provided function that develops the MIMO Nyquist plot, but would like to learn of one if such exists).

The authors make statements like "MATLAB directly connects the two points at 59.8 Hz and 60.2 Hz." The software only does what it's told, so it really sounds like the issue is with the authors' implementation and subsquent interpretation.

I'd like to dig in more, but the paper doesn't have enough information in it to recreate their results, at least for me.

Having said that, the paper did get me curious as to how nyquist and nyquistplot handle SISO loop transfer functions that have poles on the imaginary axis. After some investigation, I think those functions can, at a minimum, yield misleading results, though I'm not sure if such results would ever lead to an incorrect assessment of closed-loop stability. I'm not saying they wouldn't, I just don't know.

[Sam Chak]

Hmm... It's certainly interesting, @Paul. I seldom use Nyquist plots because I do not have a direct visualization of the frequency response when designing compensators to achieve the desired gain and phase margins.

Meanwhile, in evaluating stability for MIMO systems, I have found two examples: 'MIMO Stability Margins for Spinning Satellite' and 'Robust MIMO Controller for Two-Loop Autopilot,' both of which utilize the diskmargin() command. @Florian might find this info relevant.

s   = tf('s');

%% Case 1: Open-Loop Poles on the Imaginary Axis
G1  = 1/(s*(s^2 + 1));
G1  = zpk(G1)

G1 =
 
       1
  -----------
  s (s^2 + 1)
 
Continuous-time zero/pole/gain model.

P1  = pole(G1)

P1 = 
   0.0000 + 0.0000i
   0.0000 + 1.0000i
   0.0000 - 1.0000i

figure(1)
nyquist(G1), grid on, xlim([-2 2]), ylim([-2 2])
title({'Case 1. $G_{1}(s) = \frac{1}{s \left(s^{2} + 1\right)}$'}, 'interpreter', 'latex')

%% Case 2: One unstable Open-Loop Pole in the right-half plane
G2  = 1/(s*(s^2 - 1));
G2  = zpk(G2)

G2 =
 
        1
  -------------
  s (s+1) (s-1)
 
Continuous-time zero/pole/gain model.

P2  = pole(G2)

P2 = 3×1

0 -1 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure(2)
nyquist(G2), grid on, xlim([-2 2]), ylim([-2 2])
title({'Case 2. $G_{2}(s) = \frac{1}{s \left(s^{2} - 1\right)}$'}, 'interpreter', 'latex')

%% Case 3: Open-Loop Poles on the Real and Imaginary Axes
G3  = 1/((s^2 + 1)*(s^2 + 2*s + 1));
G3  = zpk(G3)

G3 =
 
          1
  -----------------
  (s+1)^2 (s^2 + 1)
 
Continuous-time zero/pole/gain model.

P3  = pole(G3)

P3 = 
   0.0000 + 1.0000i
   0.0000 - 1.0000i
  -1.0000 + 0.0000i
  -1.0000 - 0.0000i

figure(3)
nyquist(G3), grid on, xlim([-2 2]), ylim([-2 2])
title({'Case 3. $G_{3}(s) = \frac{1}{\left(s^{2} + 1\right) \left(s^{2} + 2 s + 1\right)}$'}, 'interpreter', 'latex')

%% Case 4: Non-minimum phase system with Open-Loop Poles on the Re and Im Axes
G4  = (s - 1)/((s^2 + 1)*(s^2 + 2*s + 1));
G4  = zpk(G4)

G4 =
 
        (s-1)
  -----------------
  (s+1)^2 (s^2 + 1)
 
Continuous-time zero/pole/gain model.

P4  = pole(G4)

P4 = 
   0.0000 + 1.0000i
   0.0000 - 1.0000i
  -1.0000 + 0.0000i
  -1.0000 - 0.0000i

figure(4)
nyquist(G4), grid on, xlim([-2 2]), ylim([-2 2])
title({'Case 4. $G_{3}(s) = \frac{\left(s - 1\right)}{\left(s^{2} + 1\right) \left(s^{2} + 2 s + 1\right)}$'}, 'interpreter', 'latex')

[Paul]

Yes, it looks like neither nyquist nor nyquistplot takes the indentation to the right around poles on the imaginary axis and so do not display the resulting arcs at infinity on the Nyquist plot, which can lead to an incorrect conclusion regarding closed-loop stability.

Let's look at your G3 example

s = zpk('s');
G3  = 1/((s^2 + 1)*(s^2 + 2*s + 1));

The open loop system has zero poles in the open RHP

pole(G3)

ans = 
   0.0000 + 1.0000i
   0.0000 - 1.0000i
  -1.0000 + 0.0000i
  -1.0000 - 0.0000i

The closed loop has two unstable poles

pole(feedback(G3,1))

ans = 
   0.1898 + 1.0432i
   0.1898 - 1.0432i
  -1.1898 + 0.6028i
  -1.1898 - 0.6028i

Draw the Nyquist plots using nyquist. The plot on the left is the "raw" plot and the plot on the right is zoomed in.

figure
hax11 = subplot(221);
nyquist(G3);title('nyquist');
hax12 = subplot(222);
nyquist(G3);title('nyquist');
axis([-2,2,-2,2]);

Now construct our own frequency vector that takes the indent to the right around the pole at s = 1j. The indent has a pretty large radius of 1e-2 just to show what the plot should look like.

s = [1j*(0:.01:0.99) , 1j+1e-2*exp(1j*linspace(-1,1,50)*pi/2), 1j*(1.01:0.01:2),1j*(3:100)];

Evaluate the open-loop transfer function.

for ii = 1:numel(s)
    n(ii) = evalfr(G3,s(ii));
end

Make the same plots

subplot(223);
plot(real(n),imag(n),real(n),-imag(n)),title('w/indent')
hold on; plot(-1,0,'r+','MarkerSize',5)
subplot(224);
plot(real(n),imag(n),real(n),-imag(n)),title('w/indent')
hold on; plot(-1,0,'r+','MarkerSize',5)
axis([-2,2,-2,2])

As we can see, the graphs from nyquist do not show any encirclements of the -1 point, which would indicate that the closed-loop system is stable because the open-loop system has no poles in the open RHP. But drawing the plots with the indentation, at full scale, shows the arcs at infinity that results in two CW encirclements, which means that there are two unstable closed-loop poles, as we showed above.

Maybe nyquist should at least warn in this case, something like "poles detected on imaginary axis, results might be suspect" or something like that.

I suspect the authors of that paper were only having a problem because whatever Matlab code they were using didn't take the indent around the imaginary poles and therefore their MIMO NC plots were not showing the arcs at infinity. But that's just a guess on my part.

[Paul]

Upon further reading the cited paper, the authors say that taking the indent around the imaginary poles is one solution to their "problem." But taking that indented path is standard procedure for Nyquist analysis of systems with poles on the imaginary axis (even if apparently not done by Control Systems Toolbox), so I'm still not sure why the authors were having any issues in the first place.

[Sam Chak]

Hi @Paul

Thank you for your detailed analysis and demonstration. I agree that the functions nyquist() and nyquistplot() should be improved to provide warnings when poles are detected on the imaginary axis.

Of course, it may be beneficial to strategically modify the Nyquist path by adding semicircular indentations around these poles, shifting to the left side of the imaginary axis to ensure that the contour does not pass through them.

[Paul]

Here is the full scale plot from nyquist

s = zpk('s');
G3  = 1/((s^2 + 1)*(s^2 + 2*s + 1));
figure
nyquist(G3)

Now zoom in really tight in the x-axis

figure
nyquist(G3)
xlim(0.25+[-1,1]/1e6)

We see that the plot isn't closed, which is at least a clue that something peculiar might be happening at infinity. But that pecularity would be pretty easy to overlook (as I did originally)

To be fair to the developers, it might be tricky to determine what actually constitutes a pole on the imaginary axis and the implications of that on how to select the contour, how to display the results, how the results are returned to the user if output arguments are requested, etc. I'm sure it's not a trivial issue.

Answer 2

[Paul]

Load in the data

load linsys
format short e

The open loop system has a pole just to the right of the real axis. Maybe this is really supposed to be an open loop integrator?

max(real(pole(linsys1)))

ans = 
   2.1349e-09

If we close the loop

H = feedback(-linsys1,eye(2));

the closed-loop system still has (probably) that same pole, which is unstable.

max(real(pole(H)))

ans = 
   2.8572e-08

But linsys1 is non-minimal

linsys1 = minreal(linsys1);

7 states removed.

max(real(pole(linsys1)))

ans = 
  -7.3799e-13

allmargin claims the closed-loop is stable

S = allmargin(-linsys1);
S.Stable

ans = logical
   1

ans = logical
   1

Which is verified by the closed-loop poles.

H = feedback(-linsys1,eye(2));
max(real(pole(H)))

ans = 
  -2.4630e-02

For your design, is a closed-loop pole so close to the origin expected?

[Florian]

Thanks a lot Paul, that makes sense. Yeah, there can be a pole close to the origin. So if I understand it correctly, minreal() removes all uncontrollable poles which don’t influence my I/O behavior

[Paul]

Each pole removed by minreal is one of: uncontrollable and observable, controllable and unobservable, or unobservable and uncontrollable.