Clear Filters
Clear Filters

allmargin gives incorrect answers for some discrete time systems

13 views (last 30 days)
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
  3 Comments
Carlos Felipe Rengifo
Carlos Felipe Rengifo on 15 Aug 2024
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".

Sign in to comment.

Accepted Answer

Paul
Paul on 15 Aug 2024
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.
  1 Comment
Paul
Paul on 16 Aug 2024
Edited: Paul on 16 Aug 2024
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.

Sign in to comment.

More Answers (1)

Carlos Felipe Rengifo
Carlos Felipe Rengifo on 16 Aug 2024
Hi Paul
Thank you for your valuable comments. They helped me to understand the problem. Yesterday, I reported the bug to Mathworks.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!