allmargin gives incorrect answers for some discrete time systems

Hi,

Consider the following code:

s = tf("s");
gs = 1/(s*(s+1));
gz = c2d(gs, 4, "zoh");
allmargin(gz)
gclz = zpk(feedback(gz,1,-1));
abs(gclz.P{1})

In R2024a, the stable flag given by allmargin is equal to one, indicating that the closed loop system is stable. On the other hand, the last sentence of the above code shows that the closed-loop system has a pole outside the unit circle, and therefore it is unstable. I do not have this issue in R2022b. Why the difference?

Sincerely

Carlos

[Les Beckham]

Your gs has a pole at -1, but when you discretize it, you specify a sample period of 4 seconds which can't possibly capture the dynamics of that pole. I know this isn't really an answer to your question, I just wanted to point out that it doesn't make sense to use a sampling period that slow when discretizing this transfer function.

[Carlos Felipe Rengifo]

Dear friend, no matter whether the sampling time is well chosen or not, a MATLAB function should not assess incorreclty the stability of a dynamic system. Imagine a Toolbox that gives correct answers only when sampling time is "well chosen".

[Paul]

Bug confirmed: allmargin can misdiagnose closed-loop stability for some discrete-time systems (3408625)

Answer 1

[Paul]

Hi Carlos,

In 2024a, allmargin introduced a new "Focus" input. Check the doc page for more information.

As far as I can tell, what's happening is as follows:

For a discrete time system, the default Focuse is Focus = [0 pi/Ts], that is the upper bound is the Nyquist frequency.

Then, allmargin uses an undocumented call to damp to compute the "natural frequency" of the closed loop poles. In your example, we have

s = tf("s");
gs = 1/(s*(s+1));
gz = c2d(gs, 4, "zoh");
clp = pole(feedback(gz,1))

clp = 2x1

-1.2707 -0.7293

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

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

wn = damp(clp,4)  % Ts = 4;

wn = 2x1

0.7877 0.7894

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

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

I'm not exactly sure yet what wn is supposed to represent in the context of this problem.

Nevertheless, both elements of wn are greater than the upper bound of Focus, which is

pi/4

ans = 0.7854

Hence, allmargin excludes both closed loop poles from the unit circle test, and since there are no poles left it returns Stable = 1, as you've observed.

I'd have to do some more digging to get a better idea of what they mean with that wn computation, or anyone can as damp is an ordinary m-file.

Also, at present I don't see a way to use the explicit Focus name/value pair to override the default. As best I can tell, if Focus(2) on input is greater than pi/Ts, then it's set to pi/Ts internally.

I'd say that this situation looks very suspicious, but I'd need to get a better understanding of what's happening with damp before going so far as to say that there is a bug, though it certainly looks like that's the case.

[Paul]

The bottom of the doc page for damp show how it computes the natural frequency of the poles of a discrete-time sytem:

s = log(z)/T

wn = abs(s)

where z is the discrete time pole. Clearly, the idea is use the inverse of z = exp(s*T) and then get the natural frequency of the "equivalent" continuous time pole.

We can use the symbolic toolbox to see exactly what's happening:

syms Ts positive
z = sym(-1.2);  % one of the closed loop poles in question
s = log(z)/Ts

s = 

We see that the complex part of s is pi/Ts, independent of the actual value of z (for z real and less than 0). Hence, the abs(s) will ALWAYS be greater than pi/Ts, and the allmargin closed-loop stability test will not include any real, negative poles.

Of course, the above is also true for poles that are inside the unit circle, so those won't be included in the stability test either, though I guess that's not important.

I think MathWorks screwed this up and it's a SERIOUS bug that you should report. If you do so, please post back as a comment to your Question with a summary of their response.

Answer 2

[Carlos Felipe Rengifo]

Hi Paul

Thank you for your valuable comments. They helped me to understand the problem. Yesterday, I reported the bug to Mathworks.

[Paul]

Please let us know what their response is. I am very interested.

[Paul]

@Carlos Felipe Rengifo

What was MathWorks response to your bug report, if you don't mind me asking?