diff --git a/control/freqplot.py b/control/freqplot.py index 881ec93dd..546972298 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,21 +522,23 @@ 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.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, *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. 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). @@ -550,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 @@ -592,6 +596,11 @@ 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`). + 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 the imaginary axis. @@ -679,11 +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) + 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: @@ -692,13 +704,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 @@ -717,46 +729,64 @@ 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 - # TODO: smarter indent radius that depends on dcgain of system - # and timebase of discrete system. + # indent the contour around any poles on/near the imaginary axis 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 - 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") - - # change contour to z-plane if necessary + if omega_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(): contour = splane_contour else: @@ -765,6 +795,15 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None, # Compute the primary curve resp = sys(contour) + 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)) count = int(np.round(np.sum(np.diff(phase)) / np.pi, 0)) @@ -792,6 +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 + 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) @@ -799,12 +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) + 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) + 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+') @@ -838,6 +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") + 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: @@ -873,7 +927,9 @@ 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)] + if len(x) in (0, 1): + return [] arrow_kw = { "arrowstyle": arrowstyle, } diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py index 4667c6219..00e31eb62 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 @@ -137,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", [ @@ -193,13 +208,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 +220,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(); @@ -226,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)