From 5a5ca77d804edb421e3c4a89b6cfbf9637dbad97 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Sun, 21 Nov 2021 00:11:08 -0800 Subject: [PATCH 1/6] nyquist plot above a certain radius is also plotted at specfied infinity-radius to help visualization --- control/freqplot.py | 42 +++++++++++++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 11 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index 881ec93dd..dc4a1e137 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -326,8 +326,8 @@ def bode_plot(syslist, omega=None, if nyquistfrq_plot: # append data for vertical nyquist freq indicator line. - # if this extra nyquist lime is is plotted in a single plot - # command then line order is preserved when + # by including this in the same plot command, line order + # is preserved when # creating a legend eg. legend(('sys1', 'sys2')) omega_nyq_line = np.array( (np.nan, nyquistfrq_plot, nyquistfrq_plot)) @@ -522,14 +522,16 @@ def gen_zero_centered_series(val_min, val_max, period): 'nyquist.mirror_style': '--', 'nyquist.arrows': 2, 'nyquist.arrow_size': 8, - 'nyquist.indent_radius': 1e-1, + 'nyquist.indent_radius': 1e-6, + 'nyquist.infinity_radius': 10, 'nyquist.indent_direction': 'right', } def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, omega_num=None, label_freq=0, color=None, - return_contour=False, warn_nyquist=True, *args, **kwargs): + return_contour=False, warn_nyquist=True, + *args, **kwargs): """Nyquist plot for a system Plots a Nyquist plot for the system over a (optional) frequency range. @@ -592,6 +594,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, arrow_style : matplotlib.patches.ArrowStyle Define style used for Nyquist curve arrows (overrides `arrow_size`). + infinity_radius : float + + indent_radius : float Amount to indent the Nyquist contour around poles that are at or near the imaginary axis. @@ -679,6 +684,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, 'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True) indent_direction = config._get_param( 'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True) + infinity_radius = config._get_param( + 'nyquist', 'infinity_radius', kwargs, _nyquist_defaults, pop=True) # If argument was a singleton, turn it into a list if not hasattr(syslist, '__iter__'): @@ -692,13 +699,13 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Go through each system and keep track of the results counts, contours = [], [] - for sys in syslist: + for isys, sys in enumerate(syslist): if not sys.issiso(): # TODO: Add MIMO nyquist plots. raise ControlMIMONotImplemented( "Nyquist plot currently only supports SISO systems.") - # Figure out the frequency range + # local copy of freq range that may be truncated above nyquist freq omega_sys = np.asarray(omega) # Determine the contour used to evaluate the Nyquist curve @@ -718,15 +725,14 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, splane_contour = 1j * omega_sys # Bend the contour around any poles on/near the imaginary axis - # TODO: smarter indent radius that depends on dcgain of system - # and timebase of discrete system. if isinstance(sys, (StateSpace, TransferFunction)) \ and indent_direction != 'none': if sys.isctime(): splane_poles = sys.pole() else: - # map z-plane poles to s-plane, ignoring any at the origin - # because we don't need to indent for them + # map z-plane poles to s-plane, ignoring any at z-plane origin + # that don't map to finite s-plane location, because we don't + # need to indent for them zplane_poles = sys.pole() zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)] splane_poles = np.log(zplane_poles)/sys.dt @@ -756,7 +762,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, else: ValueError("unknown value for indent_direction") - # change contour to z-plane if necessary + # map contour to z-plane if necessary if sys.isctime(): contour = splane_contour else: @@ -765,6 +771,14 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Compute the primary curve resp = sys(contour) + # compute infinity radius contour + # fill with nans that don't plot + infinity_contour = np.full_like(contour, np.nan + 1j * np.nan) + # where too big, use angle of resp but specifified mag + too_big = abs(resp) > infinity_radius + infinity_contour[too_big] = \ + infinity_radius *(1+isys/40) * np.exp(1j * np.angle(resp[too_big])) + # Compute CW encirclements of -1 by integrating the (unwrapped) angle phase = -unwrap(np.angle(resp + 1)) count = int(np.round(np.sum(np.diff(phase)) / np.pi, 0)) @@ -792,6 +806,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Save the components of the response x, y = resp.real, resp.imag + x_inf, y_inf = infinity_contour.real, infinity_contour.imag # Plot the primary curve p = plt.plot(x, y, '-', color=color, *args, **kwargs) @@ -799,12 +814,15 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, ax = plt.gca() _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) + # plot infinity contour + plt.plot(x_inf, y_inf, ':', color='red', *args, **kwargs) # Plot the mirror image if mirror_style is not False: p = plt.plot(x, -y, mirror_style, color=c, *args, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) + plt.plot(x_inf, -y_inf, ':', color='red', *args, **kwargs) # Mark the -1 point plt.plot([-1], [0], 'r+') @@ -838,6 +856,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, ax.set_xlabel("Real axis") ax.set_ylabel("Imaginary axis") ax.grid(color="lightgray") + ax.set_xlim(-infinity_radius*1.1, infinity_radius*1.1) + ax.set_ylim(-infinity_radius*1.1, infinity_radius*1.1) # "Squeeze" the results if len(syslist) == 1: From c008083c81f25bff16e14ea1c9f62885f8097729 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Sun, 21 Nov 2021 00:42:04 -0800 Subject: [PATCH 2/6] improve terminology and docstring --- control/freqplot.py | 40 +++++++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index dc4a1e137..82d3231ef 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -522,8 +522,8 @@ def gen_zero_centered_series(val_min, val_max, period): 'nyquist.mirror_style': '--', 'nyquist.arrows': 2, 'nyquist.arrow_size': 8, - 'nyquist.indent_radius': 1e-6, - 'nyquist.infinity_radius': 10, + 'nyquist.indent_radius': 1e-3, + 'nyquist.maximum_magnitude': 5, 'nyquist.indent_direction': 'right', } @@ -594,8 +594,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, arrow_style : matplotlib.patches.ArrowStyle Define style used for Nyquist curve arrows (overrides `arrow_size`). - infinity_radius : float - + maximum_magnitude : float + Magnitude/radius at which to draw contours whose magnitude exceeds + this value. indent_radius : float Amount to indent the Nyquist contour around poles that are at or near @@ -684,8 +685,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, 'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True) indent_direction = config._get_param( 'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True) - infinity_radius = config._get_param( - 'nyquist', 'infinity_radius', kwargs, _nyquist_defaults, pop=True) + maximum_magnitude = config._get_param( + 'nyquist', 'maximum_magnitude', kwargs, _nyquist_defaults, pop=True) # If argument was a singleton, turn it into a list if not hasattr(syslist, '__iter__'): @@ -771,13 +772,13 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Compute the primary curve resp = sys(contour) - # compute infinity radius contour + # compute large radius contour where it exceeds # fill with nans that don't plot - infinity_contour = np.full_like(contour, np.nan + 1j * np.nan) + large_mag_contour = np.full_like(contour, np.nan + 1j * np.nan) # where too big, use angle of resp but specifified mag - too_big = abs(resp) > infinity_radius - infinity_contour[too_big] = \ - infinity_radius *(1+isys/40) * np.exp(1j * np.angle(resp[too_big])) + too_big = abs(resp) > maximum_magnitude + large_mag_contour[too_big] = \ + maximum_magnitude *(1+isys/20) * np.exp(1j * np.angle(resp[too_big])) # Compute CW encirclements of -1 by integrating the (unwrapped) angle phase = -unwrap(np.angle(resp + 1)) @@ -806,7 +807,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Save the components of the response x, y = resp.real, resp.imag - x_inf, y_inf = infinity_contour.real, infinity_contour.imag + x_lg, y_lg = large_mag_contour.real, large_mag_contour.imag # Plot the primary curve p = plt.plot(x, y, '-', color=color, *args, **kwargs) @@ -814,15 +815,19 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, ax = plt.gca() _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) - # plot infinity contour - plt.plot(x_inf, y_inf, ':', color='red', *args, **kwargs) + # plot large magnitude contour + p = plt.plot(x_lg, y_lg, ':', color='red', *args, **kwargs) + _add_arrows_to_line2D( + ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) # Plot the mirror image if mirror_style is not False: p = plt.plot(x, -y, mirror_style, color=c, *args, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) - plt.plot(x_inf, -y_inf, ':', color='red', *args, **kwargs) + p = plt.plot(x_lg, -y_lg, ':', color='red', *args, **kwargs) + _add_arrows_to_line2D( + ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) # Mark the -1 point plt.plot([-1], [0], 'r+') @@ -856,8 +861,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, ax.set_xlabel("Real axis") ax.set_ylabel("Imaginary axis") ax.grid(color="lightgray") - ax.set_xlim(-infinity_radius*1.1, infinity_radius*1.1) - ax.set_ylim(-infinity_radius*1.1, infinity_radius*1.1) + ax.set_xlim(-maximum_magnitude*1.1, maximum_magnitude*1.1) + ax.set_ylim(-maximum_magnitude*1.1, maximum_magnitude*1.1) # "Squeeze" the results if len(syslist) == 1: @@ -893,6 +898,7 @@ def _add_arrows_to_line2D( if not isinstance(line, mpl.lines.Line2D): raise ValueError("expected a matplotlib.lines.Line2D object") x, y = line.get_xdata(), line.get_ydata() + x, y = x[np.isfinite(x)], y[np.isfinite(x)] arrow_kw = { "arrowstyle": arrowstyle, From 5b047f75eab8e4a9077529e388e6f76fa9ab1426 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Sun, 21 Nov 2021 00:57:35 -0800 Subject: [PATCH 3/6] fix error when no large contours --- control/freqplot.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/control/freqplot.py b/control/freqplot.py index 82d3231ef..d3fc9ace3 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -899,7 +899,8 @@ def _add_arrows_to_line2D( raise ValueError("expected a matplotlib.lines.Line2D object") x, y = line.get_xdata(), line.get_ydata() x, y = x[np.isfinite(x)], y[np.isfinite(x)] - + if len(x) == 0: + return [] arrow_kw = { "arrowstyle": arrowstyle, } From 372efb4ac0ea07c3c550ecc09e0e6676b9afe944 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Sun, 21 Nov 2021 01:24:53 -0800 Subject: [PATCH 4/6] bugfixes --- control/freqplot.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index d3fc9ace3..2ff6dd409 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -522,7 +522,7 @@ def gen_zero_centered_series(val_min, val_max, period): 'nyquist.mirror_style': '--', 'nyquist.arrows': 2, 'nyquist.arrow_size': 8, - 'nyquist.indent_radius': 1e-3, + 'nyquist.indent_radius': 1e-2, 'nyquist.maximum_magnitude': 5, 'nyquist.indent_direction': 'right', } @@ -899,7 +899,7 @@ def _add_arrows_to_line2D( raise ValueError("expected a matplotlib.lines.Line2D object") x, y = line.get_xdata(), line.get_ydata() x, y = x[np.isfinite(x)], y[np.isfinite(x)] - if len(x) == 0: + if len(x) in (0, 1): return [] arrow_kw = { "arrowstyle": arrowstyle, From 23b3e8c74cdbfa10176c0d19bf72cde442e9bade Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Mon, 22 Nov 2021 12:58:15 -0800 Subject: [PATCH 5/6] new indent code that adds a consistent number of points at every indent if omega_range_given is false --- control/freqplot.py | 80 ++++++++++++++++++++++------------- control/tests/nyquist_test.py | 17 ++++---- 2 files changed, 60 insertions(+), 37 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index 2ff6dd409..47401e6e1 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -522,7 +522,7 @@ def gen_zero_centered_series(val_min, val_max, period): 'nyquist.mirror_style': '--', 'nyquist.arrows': 2, 'nyquist.arrow_size': 8, - 'nyquist.indent_radius': 1e-2, + 'nyquist.indent_radius': 1e-6, 'nyquist.maximum_magnitude': 5, 'nyquist.indent_direction': 'right', } @@ -552,7 +552,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, If True, plot magnitude omega : array_like - Set of frequencies to be evaluated, in rad/sec. + Set of frequencies to be evaluated, in rad/sec. Remark: if omega is + specified, you may need to specify specify a larger `indent_radius` + than the default to get reasonable contours. omega_limits : array_like of two values Limits to the range of frequencies. Ignored if omega is provided, and @@ -725,7 +727,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # do indentations in s-plane where it is more convenient splane_contour = 1j * omega_sys - # Bend the contour around any poles on/near the imaginary axis + # indent the contour around any poles on/near the imaginary axis if isinstance(sys, (StateSpace, TransferFunction)) \ and indent_direction != 'none': if sys.isctime(): @@ -738,30 +740,49 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)] splane_poles = np.log(zplane_poles)/sys.dt - if splane_contour[1].imag > indent_radius \ - and np.any(np.isclose(abs(splane_poles), 0)) \ - and not omega_range_given: - # add some points for quarter circle around poles at origin - splane_contour = np.concatenate( - (1j * np.linspace(0., indent_radius, 50), - splane_contour[1:])) - for i, s in enumerate(splane_contour): - # Find the nearest pole - p = splane_poles[(np.abs(splane_poles - s)).argmin()] - # See if we need to indent around it - if abs(s - p) < indent_radius: - if p.real < 0 or (np.isclose(p.real, 0) \ - and indent_direction == 'right'): - # Indent to the right - splane_contour[i] += \ - np.sqrt(indent_radius ** 2 - (s-p).imag ** 2) - elif p.real > 0 or (np.isclose(p.real, 0) \ - and indent_direction == 'left'): - # Indent to the left - splane_contour[i] -= \ - np.sqrt(indent_radius ** 2 - (s-p).imag ** 2) - else: - ValueError("unknown value for indent_direction") + if omega_range_given: + # indent given contour + for i, s in enumerate(splane_contour): + # Find the nearest pole + p = splane_poles[(np.abs(splane_poles - s)).argmin()] + # See if we need to indent around it + if abs(s - p) < indent_radius: + if p.real < 0 or (np.isclose(p.real, 0) \ + and indent_direction == 'right'): + # Indent to the right + splane_contour[i] += \ + np.sqrt(indent_radius ** 2 - (s-p).imag ** 2) + elif p.real > 0 or (np.isclose(p.real, 0) \ + and indent_direction == 'left'): + # Indent to the left + splane_contour[i] -= \ + np.sqrt(indent_radius ** 2 - (s-p).imag ** 2) + else: + ValueError("unknown value for indent_direction") + else: + # find poles that are near imag axis and replace contour + # near there with indented contour + if indent_direction == 'right': + indent = indent_radius * \ + np.exp(1j * np.linspace(-np.pi/2, np.pi/2, 51)) + elif indent_direction == 'left': + indent = -indent_radius * \ + np.exp(1j * np.linspace(np.pi/2, -np.pi/2, 51)) + else: + ValueError("unknown value for indent_direction") + for p in splane_poles[np.isclose(splane_poles.real, 0) \ + & (splane_poles.imag >= 0)]: + start = np.searchsorted( + splane_contour.imag, p.imag - indent_radius) + end = np.searchsorted( + splane_contour.imag, p.imag + indent_radius) + # only keep indent parts that are > 0 and within contour + indent_mask = ((indent + p).imag >= 0) \ + & ((indent + p).imag <= splane_contour[-1].imag) + splane_contour = np.concatenate(( + splane_contour[:start], + indent[indent_mask] + p, + splane_contour[end:])) # map contour to z-plane if necessary if sys.isctime(): @@ -816,7 +837,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) # plot large magnitude contour - p = plt.plot(x_lg, y_lg, ':', color='red', *args, **kwargs) + p = plt.plot(x_lg, y_lg, color='red', *args, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) @@ -825,7 +846,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, p = plt.plot(x, -y, mirror_style, color=c, *args, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) - p = plt.plot(x_lg, -y_lg, ':', color='red', *args, **kwargs) + p = plt.plot(x_lg, -y_lg, mirror_style, + color='red', *args, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py index 4667c6219..8afd6958e 100644 --- a/control/tests/nyquist_test.py +++ b/control/tests/nyquist_test.py @@ -69,7 +69,8 @@ def test_nyquist_basic(): # Make sure that we can turn off frequency modification count, contour_indented = ct.nyquist_plot( - sys, np.linspace(1e-4, 1e2, 100), return_contour=True) + sys, np.linspace(1e-4, 1e2, 100), indent_radius=0.01, + return_contour=True) assert not all(contour_indented.real == 0) count, contour = ct.nyquist_plot( sys, np.linspace(1e-4, 1e2, 100), return_contour=True, @@ -87,14 +88,17 @@ def test_nyquist_basic(): assert _Z(sys) == count + _P(sys) # Nyquist plot with poles on imaginary axis, omega specified + # when specifying omega, usually need to specify larger indent radius sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1]) - count = ct.nyquist_plot(sys, np.linspace(1e-3, 1e1, 1000)) + count = ct.nyquist_plot(sys, np.linspace(1e-3, 1e1, 1000), + indent_radius=0.1) assert _Z(sys) == count + _P(sys) # Nyquist plot with poles on imaginary axis, omega specified, with contour sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1]) count, contour = ct.nyquist_plot( - sys, np.linspace(1e-3, 1e1, 1000), return_contour=True) + sys, np.linspace(1e-3, 1e1, 1000), indent_radius=0.1, + return_contour=True) assert _Z(sys) == count + _P(sys) # Nyquist plot with poles on imaginary axis, return contour @@ -193,13 +197,10 @@ def test_nyquist_indent(): plt.title("Pole at origin; indent_radius=default") assert _Z(sys) == count + _P(sys) - # first value of default omega vector was 0.1, replaced by 0. for contour - # indent_radius is larger than 0.1 -> no extra quater circle around origin + # confirm that first element of indent is properly indented count, contour = ct.nyquist_plot(sys, plot=False, indent_radius=.1007, return_contour=True) np.testing.assert_allclose(contour[0], .1007+0.j) - # second value of omega_vector is larger than indent_radius: not indented - assert np.all(contour.real[2:] == 0.) plt.figure(); count, contour = ct.nyquist_plot(sys, indent_radius=0.01, @@ -208,7 +209,7 @@ def test_nyquist_indent(): assert _Z(sys) == count + _P(sys) # indent radius is smaller than the start of the default omega vector # check that a quarter circle around the pole at origin has been added. - np.testing.assert_allclose(contour[:50].real**2 + contour[:50].imag**2, + np.testing.assert_allclose(contour[:25].real**2 + contour[:25].imag**2, 0.01**2) plt.figure(); From 5e9b1791ca10ea7a235696df10f2688bc3d9ff34 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Tue, 23 Nov 2021 10:24:46 -0800 Subject: [PATCH 6/6] rename keyword to plot_maximum_magintude and option to turn it off by setting it to zero; docstring improvement; old indentation method still available if omega is explicitly specified --- control/freqplot.py | 61 +++++++++++++++++++---------------- control/tests/nyquist_test.py | 17 +++++++++- 2 files changed, 50 insertions(+), 28 deletions(-) diff --git a/control/freqplot.py b/control/freqplot.py index 47401e6e1..546972298 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -523,14 +523,14 @@ def gen_zero_centered_series(val_min, val_max, period): 'nyquist.arrows': 2, 'nyquist.arrow_size': 8, 'nyquist.indent_radius': 1e-6, - 'nyquist.maximum_magnitude': 5, + 'nyquist.plot_maximum_magnitude': 5, 'nyquist.indent_direction': 'right', } def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, omega_num=None, label_freq=0, color=None, - return_contour=False, warn_nyquist=True, + return_contour=False, warn_nyquist=True, *args, **kwargs): """Nyquist plot for a system @@ -538,7 +538,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, The curve is computed by evaluating the Nyqist segment along the positive imaginary axis, with a mirror image generated to reflect the negative imaginary axis. Poles on or near the imaginary axis are avoided using a - small indentation. The portion of the Nyquist contour at infinity is not + small indentation. If portions of the Nyquist plot The portion of the Nyquist contour at infinity is not explicitly computed (since it maps to a constant value for any system with a proper transfer function). @@ -596,9 +596,10 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, arrow_style : matplotlib.patches.ArrowStyle Define style used for Nyquist curve arrows (overrides `arrow_size`). - maximum_magnitude : float - Magnitude/radius at which to draw contours whose magnitude exceeds - this value. + plot_maximum_magnitude : float + If this value is not 0, an additional curve showing the path of + the Nyquist plot when it leaves the plot window is drawn at a radius + equal to this value. indent_radius : float Amount to indent the Nyquist contour around poles that are at or near @@ -687,13 +688,14 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, 'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True) indent_direction = config._get_param( 'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True) - maximum_magnitude = config._get_param( - 'nyquist', 'maximum_magnitude', kwargs, _nyquist_defaults, pop=True) + plot_maximum_magnitude = config._get_param( + 'nyquist', 'plot_maximum_magnitude', kwargs, _nyquist_defaults, pop=True) # If argument was a singleton, turn it into a list if not hasattr(syslist, '__iter__'): syslist = (syslist,) + omega_given = True if omega is not None else False omega, omega_range_given = _determine_omega_vector( syslist, omega, omega_limits, omega_num) if not omega_range_given: @@ -740,7 +742,7 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)] splane_poles = np.log(zplane_poles)/sys.dt - if omega_range_given: + if omega_given: # indent given contour for i, s in enumerate(splane_contour): # Find the nearest pole @@ -793,13 +795,14 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Compute the primary curve resp = sys(contour) - # compute large radius contour where it exceeds - # fill with nans that don't plot - large_mag_contour = np.full_like(contour, np.nan + 1j * np.nan) - # where too big, use angle of resp but specifified mag - too_big = abs(resp) > maximum_magnitude - large_mag_contour[too_big] = \ - maximum_magnitude *(1+isys/20) * np.exp(1j * np.angle(resp[too_big])) + if plot_maximum_magnitude != 0: + # compute large radius contour where it exceeds + # fill with nans that don't plot + large_mag_contour = np.full_like(contour, np.nan + 1j * np.nan) + # where too big, use angle of resp but specifified mag + too_big = abs(resp) > plot_maximum_magnitude + large_mag_contour[too_big] = plot_maximum_magnitude *\ + (1+isys/20) * np.exp(1j * np.angle(resp[too_big])) # Compute CW encirclements of -1 by integrating the (unwrapped) angle phase = -unwrap(np.angle(resp + 1)) @@ -828,7 +831,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Save the components of the response x, y = resp.real, resp.imag - x_lg, y_lg = large_mag_contour.real, large_mag_contour.imag + if plot_maximum_magnitude != 0: + x_lg, y_lg = large_mag_contour.real, large_mag_contour.imag # Plot the primary curve p = plt.plot(x, y, '-', color=color, *args, **kwargs) @@ -836,20 +840,22 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, ax = plt.gca() _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) - # plot large magnitude contour - p = plt.plot(x_lg, y_lg, color='red', *args, **kwargs) - _add_arrows_to_line2D( - ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) + if plot_maximum_magnitude != 0: + # plot large magnitude contour + p = plt.plot(x_lg, y_lg, color='red', *args, **kwargs) + _add_arrows_to_line2D( + ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1) # Plot the mirror image if mirror_style is not False: p = plt.plot(x, -y, mirror_style, color=c, *args, **kwargs) _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) - p = plt.plot(x_lg, -y_lg, mirror_style, - color='red', *args, **kwargs) - _add_arrows_to_line2D( - ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) + if plot_maximum_magnitude != 0: + p = plt.plot(x_lg, -y_lg, mirror_style, + color='red', *args, **kwargs) + _add_arrows_to_line2D( + ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) # Mark the -1 point plt.plot([-1], [0], 'r+') @@ -883,8 +889,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, ax.set_xlabel("Real axis") ax.set_ylabel("Imaginary axis") ax.grid(color="lightgray") - ax.set_xlim(-maximum_magnitude*1.1, maximum_magnitude*1.1) - ax.set_ylim(-maximum_magnitude*1.1, maximum_magnitude*1.1) + if plot_maximum_magnitude != 0: + ax.set_xlim(-plot_maximum_magnitude*1.1,plot_maximum_magnitude*1.1) + ax.set_ylim(-plot_maximum_magnitude*1.1,plot_maximum_magnitude*1.1) # "Squeeze" the results if len(syslist) == 1: diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py index 8afd6958e..00e31eb62 100644 --- a/control/tests/nyquist_test.py +++ b/control/tests/nyquist_test.py @@ -141,12 +141,23 @@ def test_nyquist_fbs_examples(): count = ct.nyquist_plot(sys) assert _Z(sys) == count + _P(sys) + # when omega is specified, no points are added at indentations, leading + # to incorrect encirclement count + plt.figure() + plt.title("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2) [zoom]") + count = ct.nyquist_plot(sys, omega=np.linspace(1.5, 1e3, 1000)) + # assert _Z(sys) == count + _P(sys) + assert count == -1 + + # when only the range of the frequency is provided, this permits + # extra points to be added at indentations, leading to correct count plt.figure() plt.title("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2) [zoom]") count = ct.nyquist_plot(sys, omega_limits=[1.5, 1e3]) # Frequency limits for zoom give incorrect encirclement count # assert _Z(sys) == count + _P(sys) - assert count == -1 + assert count == 0 + @pytest.mark.parametrize("arrows", [ @@ -227,6 +238,10 @@ def test_nyquist_indent(): plt.title("Imaginary poles; encirclements = %d" % count) assert _Z(sys) == count + _P(sys) + # confirm we can turn off maximum_magitude plot + plt.figure(); + count = ct.nyquist_plot(sys, plot_maximum_magnitude=0) + # Imaginary poles with indentation to the left plt.figure(); count = ct.nyquist_plot(sys, indent_direction='left', label_freq=300)