Bug fix and improvements to Nyquist plots #722
Conversation
* use plot_curve_magnitude to scale response at large magnitudes * add new line styles of scaled points on the curve * add points to contour near imaginary poles to avoid missing encirclements * add offsets to primary/mirror curve at large amplitude to avoid overlap * (still needs some unit tests, documention updates, and comparisons)
f3dfe04 to
5d9c7f0
Compare
* changed return type for poles() and zeros() to complex * updated (missing) discrete-time tests for stability warnings * changed processing of near imaginary poles to pure imaginary poles
| - (s - p).real | ||
|
|
||
| # Figure out which way to offset the contour point | ||
| if p.real < 0 or (p.real == 0 and |
There was a problem hiding this comment.
I think I was the one who changed this code from equality to zero to is close to zero for the reason that sometimes when working with discrete systems you start with a continuous time system, with poles on imag axis, and then discretize it. When these poles are mapped back to the s plane in this code, their real part doesn’t end up to be perfectly zero anymore.
I Probably should have written a unit test!
There was a problem hiding this comment.
The current code should OK because the Nyquist criterion is "robust" with respect to whether you indent to the right or the left. You just have to be consistent in how you count open and closed loop poles.
Here's an example (from tests/nyquist_test.py):
sys_c = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1])
sys = ct.c2d(sys_c, 0.01)
np.abs(sys.poles()) - 1
array([ 9.58553237e-11, 9.58553237e-11, -9.95016656e-03, -1.98013266e-02])
The closed loop system has poles at
sys.feedback().poles()
array([1.00116706+0.01054462j, 1.00116706-0.01054462j,
0.98390719+0.00189796j, 0.98390719-0.00189796j])
so it is unstable (as is the continuous time system, which has a pair of poles with real part 0.122).
If you run the current code, it will map the z-plane poles to
np.log(sys.poles()) / sys.dt
array([ 9.58553380e-09+1.j, 9.58553380e-09-1.j, -1.00000003e+00+0.j,
-1.99999999e+00+0.j])
Note that what used to be a purely imaginary pole (at +/- 1j) is now slightly unstable (as you described). But the Nyquist criterion is still OK: we indent to the left (since they are unstable) and we end up with zero net encirclements:
count = ct.nyquist_plot(sys)
Since Z = N + P = 0 + 2 = 2, we conclude that the closed loop system has two poles in the right half-plane (correct!).
If it had turned out that the poles had mapped to exactly on the unit circle or to just inside the unit circle, then we would get a plot that looks like the continuous time plot for this example:
Note that in this case there are 2 encirclements and we get Z = N + P = 2 + 0 = 2, which is the same as before.
This PR fixes a bug in the the way that the
nyquist_plotcommand was implemented and introduces a bunch of enhancements, including some of the ideas in PR #671.Bug fix: because of the new default parameter settings introduced in #534, it was possible for the indentation around an open-loop pole that was near the imaginary axis to exclude a region in the right half-plane that contained a closed-loop pole (so, for example, if you had a stable open-loop pole at -0.05 and an unstable closed-loop pole at 0.05, the default indentation of 0.1 would cause you to miss the encirclement). This has now been fixed by (1) generating a warning when the indented contour bypasses a closed loop pole and (2) implementing new functionality that allows the default indentation to be reduced (to 1e-4) while still keeping the magnitude of the Nyquist plot within reasonable limits.
Enhancements: in fixing the bug above, it was necessary to find a better balance between the size of the indentation and the magnitude of the Nyquist curve. This is also the goal of PR #671. The following changes have been done to implement this:
max_curve_magnitudefunction that limits the size of the plotted Nyquist curve. The default value for this parameter is 20, which is small enough that that -1 point is still visible in the plot. All points on the curve that are greater than the maximum magnitude specified are scaled so that the angle of the curve is preserved, but the magnitude is limited. The value for this parameter can be passed as a keyword argument or set usingct.defaults.primary_linestyleandmirror_linestylekeywords and default to['-', '-.']and['--', ':'].curve_offsetkeyword can be used so that overlapping sections of the Nyquist curve (typically at the limit values set bymax_curve_magnitudeare now separated out by a small amount (by adjusting the primary curve outward and the mirror curve inward by a small amount that is smooth at the end of the scaled sections).start_markerkeyword and the size of the marker is controlled with thestart_marker_sizekeyword (default = a filled circle of size 4).warn_encirclements=False.use_legacy_defaults('0.9.1').The new features can be seen in the following example
which produces the following plot:
Features to note:
Comparisons between old and new Nyquist plots: