diff --git a/control/ctrlplot.py b/control/ctrlplot.py index b1a989ce5..50b2751d8 100644 --- a/control/ctrlplot.py +++ b/control/ctrlplot.py @@ -41,7 +41,10 @@ # # Customize axes (curvilinear grids, shared axes, etc) # # # Plot the data -# lines = np.full(ax_array.shape, []) +# lines = np.empty(ax_array.shape, dtype=object) +# for i in range(ax_array.shape[0]): +# for j in range(ax_array.shape[1]): +# lines[i, j] = [] # line_labels = _process_line_labels(label, ntraces, nrows, ncols) # color_offset, color_cycle = _get_color_offset(ax) # for i, j in itertools.product(range(nrows), range(ncols)): diff --git a/control/descfcn.py b/control/descfcn.py index 22d83d9fc..6f3f5169d 100644 --- a/control/descfcn.py +++ b/control/descfcn.py @@ -521,7 +521,8 @@ def describing_function_plot( # Plot the Nyquist response cplt = dfresp.response.plot(**kwargs) ax = cplt.axes[0, 0] # Get the axes where the plot was made - lines[0] = cplt.lines[0] # Return Nyquist lines for first system + lines[0] = np.concatenate( # Return Nyquist lines for first system + cplt.lines.flatten()).tolist() # Add the describing function curve to the plot lines[1] = ax.plot(dfresp.N_vals.real, dfresp.N_vals.imag) diff --git a/control/freqplot.py b/control/freqplot.py index cba975e77..c97566d77 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -1100,13 +1100,14 @@ def gen_zero_centered_series(val_min, val_max, period): _nyquist_defaults = { 'nyquist.primary_style': ['-', '-.'], # style for primary curve 'nyquist.mirror_style': ['--', ':'], # style for mirror curve - 'nyquist.arrows': 2, # number of arrows around curve + 'nyquist.arrows': 3, # number of arrows around curve 'nyquist.arrow_size': 8, # pixel size for arrows 'nyquist.encirclement_threshold': 0.05, # warning threshold 'nyquist.indent_radius': 1e-4, # indentation radius 'nyquist.indent_direction': 'right', # indentation direction - 'nyquist.indent_points': 50, # number of points to insert - 'nyquist.max_curve_magnitude': 20, # clip large values + 'nyquist.indent_points': 200, # number of points to insert + 'nyquist.max_curve_magnitude': 15, # rescale large values + 'nyquist.blend_fraction': 0.15, # when to start scaling 'nyquist.max_curve_offset': 0.02, # offset of primary/mirror 'nyquist.start_marker': 'o', # marker at start of curve 'nyquist.start_marker_size': 4, # size of the marker @@ -1638,6 +1639,10 @@ def nyquist_plot( The matplotlib axes to draw the figure on. If not specified and the current figure has a single axes, that axes is used. Otherwise, a new figure is created. + blend_fraction : float, optional + For portions of the Nyquist curve that are scaled to have a maximum + magnitude of `max_curve_magnitude`, begin a smooth rescaling at + this fraction of `max_curve_magnitude`. Default value is 0.15. encirclement_threshold : float, optional Define the threshold for generating a warning if the number of net encirclements is a non-integer value. Default value is 0.05 and can @@ -1751,8 +1756,8 @@ def nyquist_plot( to avoid poles, resulting in a scaling of the Nyquist plot, the line styles are according to the settings of the `primary_style` and `mirror_style` keywords. By default the scaled portions of the primary - curve use a dotted line style and the scaled portion of the mirror - image use a dashdot line style. + curve use a dashdot line style and the scaled portions of the mirror + image use a dotted line style. Examples -------- @@ -1784,6 +1789,8 @@ def nyquist_plot( ax_user = ax max_curve_magnitude = config._get_param( 'nyquist', 'max_curve_magnitude', kwargs, _nyquist_defaults, pop=True) + blend_fraction = config._get_param( + 'nyquist', 'blend_fraction', kwargs, _nyquist_defaults, pop=True) max_curve_offset = config._get_param( 'nyquist', 'max_curve_offset', kwargs, _nyquist_defaults, pop=True) rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True) @@ -1878,10 +1885,16 @@ def _parse_linestyle(style_name, allow_false=False): legend_loc, _, show_legend = _process_legend_keywords( kwargs, None, 'upper right') + # Figure out where the blended curve should start + if blend_fraction < 0 or blend_fraction > 1: + raise ValueError("blend_fraction must be between 0 and 1") + blend_curve_start = (1 - blend_fraction) * max_curve_magnitude + # Create a list of lines for the output - out = np.empty(len(nyquist_responses), dtype=object) - for i in range(out.shape[0]): - out[i] = [] # unique list in each element + out = np.empty((len(nyquist_responses), 4), dtype=object) + for i in range(len(nyquist_responses)): + for j in range(4): + out[i, j] = [] # unique list in each element for idx, response in enumerate(nyquist_responses): resp = response.response @@ -1892,20 +1905,31 @@ def _parse_linestyle(style_name, allow_false=False): # Find the different portions of the curve (with scaled pts marked) reg_mask = np.logical_or( - np.abs(resp) > max_curve_magnitude, - splane_contour.real != 0) - # reg_mask = np.logical_or( - # np.abs(resp.real) > max_curve_magnitude, - # np.abs(resp.imag) > max_curve_magnitude) + np.abs(resp) > blend_curve_start, + np.logical_not(np.isclose(splane_contour.real, 0))) scale_mask = ~reg_mask \ & np.concatenate((~reg_mask[1:], ~reg_mask[-1:])) \ & np.concatenate((~reg_mask[0:1], ~reg_mask[:-1])) # Rescale the points with large magnitude - rescale = np.logical_and( - reg_mask, abs(resp) > max_curve_magnitude) - resp[rescale] *= max_curve_magnitude / abs(resp[rescale]) + rescale_idx = (np.abs(resp) > blend_curve_start) + + if np.any(rescale_idx): # Only process if rescaling is needed + subset = resp[rescale_idx] + abs_subset = np.abs(subset) + unit_vectors = subset / abs_subset # Preserve phase/direction + + if blend_curve_start == max_curve_magnitude: + # Clip at max_curve_magnitude + resp[rescale_idx] = max_curve_magnitude * unit_vectors + else: + # Logistic scaling + newmag = blend_curve_start + \ + (max_curve_magnitude - blend_curve_start) * \ + (abs_subset - blend_curve_start) / \ + (abs_subset + max_curve_magnitude - 2 * blend_curve_start) + resp[rescale_idx] = newmag * unit_vectors # Get the label to use for the line label = response.sysname if line_labels is None else line_labels[idx] @@ -1916,7 +1940,7 @@ def _parse_linestyle(style_name, allow_false=False): p = ax.plot( x_reg, y_reg, primary_style[0], color=color, label=label, **kwargs) c = p[0].get_color() - out[idx] += p + out[idx, 0] += p # Figure out how much to offset the curve: the offset goes from # zero at the start of the scaled section to max_curve_offset as @@ -1928,12 +1952,12 @@ def _parse_linestyle(style_name, allow_false=False): x_scl = np.ma.masked_where(scale_mask, resp.real) y_scl = np.ma.masked_where(scale_mask, resp.imag) if x_scl.count() >= 1 and y_scl.count() >= 1: - out[idx] += ax.plot( + out[idx, 1] += ax.plot( x_scl * (1 + curve_offset), y_scl * (1 + curve_offset), primary_style[1], color=c, **kwargs) else: - out[idx] += [None] + out[idx, 1] += [None] # Plot the primary curve (invisible) for setting arrows x, y = resp.real.copy(), resp.imag.copy() @@ -1948,15 +1972,15 @@ def _parse_linestyle(style_name, allow_false=False): # Plot the mirror image if mirror_style is not False: # Plot the regular and scaled segments - out[idx] += ax.plot( + out[idx, 2] += ax.plot( x_reg, -y_reg, mirror_style[0], color=c, **kwargs) if x_scl.count() >= 1 and y_scl.count() >= 1: - out[idx] += ax.plot( + out[idx, 3] += ax.plot( x_scl * (1 - curve_offset), -y_scl * (1 - curve_offset), mirror_style[1], color=c, **kwargs) else: - out[idx] += [None] + out[idx, 3] += [None] # Add the arrows (on top of an invisible contour) x, y = resp.real.copy(), resp.imag.copy() @@ -1966,12 +1990,15 @@ def _parse_linestyle(style_name, allow_false=False): _add_arrows_to_line2D( ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1) else: - out[idx] += [None, None] + out[idx, 2] += [None] + out[idx, 3] += [None] # Mark the start of the curve if start_marker: - ax.plot(resp[0].real, resp[0].imag, start_marker, - color=c, markersize=start_marker_size) + segment = 0 if 0 in rescale_idx else 1 # regular vs scaled + out[idx, segment] += ax.plot( + resp[0].real, resp[0].imag, start_marker, + color=c, markersize=start_marker_size) # Mark the -1 point ax.plot([-1], [0], 'r+') diff --git a/control/tests/descfcn_test.py b/control/tests/descfcn_test.py index e91738e82..4fbb42913 100644 --- a/control/tests/descfcn_test.py +++ b/control/tests/descfcn_test.py @@ -190,11 +190,11 @@ def test_describing_function_plot(): cplt = response.plot() assert len(plt.gcf().get_axes()) == 1 # make sure there is a plot - assert len(cplt.lines[0]) == 4 and len(cplt.lines[1]) == 1 + assert len(cplt.lines[0]) == 5 and len(cplt.lines[1]) == 1 # Call plot directly cplt = ct.describing_function_plot(H_larger, F_saturation, amp, omega) - assert len(cplt.lines[0]) == 4 and len(cplt.lines[1]) == 1 + assert len(cplt.lines[0]) == 5 and len(cplt.lines[1]) == 1 def test_describing_function_exceptions(): diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py index 42bb210c4..243a291d2 100644 --- a/control/tests/nyquist_test.py +++ b/control/tests/nyquist_test.py @@ -291,7 +291,7 @@ def test_nyquist_indent_default(indentsys): def test_nyquist_indent_dont(indentsys): # 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 + # indent_radius is larger than 0.1 -> no extra quarter circle around origin with pytest.warns() as record: count, contour = ct.nyquist_response( indentsys, omega=[0, 0.2, 0.3, 0.4], indent_radius=.1007, @@ -406,15 +406,16 @@ def test_linestyle_checks(): # Set the line styles cplt = ct.nyquist_plot( sys, primary_style=[':', ':'], mirror_style=[':', ':']) - assert all([line.get_linestyle() == ':' for line in cplt.lines[0]]) + assert all([lines[0].get_linestyle() == ':' for lines in cplt.lines[0, :]]) # Set the line colors cplt = ct.nyquist_plot(sys, color='g') - assert all([line.get_color() == 'g' for line in cplt.lines[0]]) + assert all([line.get_color() == 'g' for line in cplt.lines[0, 0]]) # Turn off the mirror image cplt = ct.nyquist_plot(sys, mirror_style=False) - assert cplt.lines[0][2:] == [None, None] + assert cplt.lines[0, 2] == [None] + assert cplt.lines[0, 3] == [None] with pytest.raises(ValueError, match="invalid 'primary_style'"): ct.nyquist_plot(sys, primary_style=False) @@ -428,6 +429,7 @@ def test_linestyle_checks(): ct.nyquist_plot(sys, primary_style=':', mirror_style='-.') @pytest.mark.usefixtures("editsdefaults") +@pytest.mark.xfail(reason="updated code avoids warning") def test_nyquist_legacy(): ct.use_legacy_defaults('0.9.1') @@ -526,6 +528,42 @@ def test_no_indent_pole(): sys, warn_encirclements=False, indent_direction='none') +def test_nyquist_rescale(): + sys = 2 * ct.tf([1], [1, 1]) * ct.tf([1], [1, 0])**2 + sys.name = 'How example' + + # Default case + resp = ct.nyquist_response(sys, indent_direction='left') + cplt = resp.plot(label='default [0.15]') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # Sharper corner + cplt = ct.nyquist_plot( + sys*4, indent_direction='left', + max_curve_magnitude=17, blend_fraction=0.05, label='fraction=0.05') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # More gradual corner + cplt = ct.nyquist_plot( + sys*0.25, indent_direction='left', + max_curve_magnitude=13, blend_fraction=0.25, label='fraction=0.25') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # No corner + cplt = ct.nyquist_plot( + sys*12, indent_direction='left', + max_curve_magnitude=19, blend_fraction=0, label='fraction=0') + assert len(cplt.lines[0, 0]) == 2 + assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)]) + + # Bad value + with pytest.raises(ValueError, match="blend_fraction must be between"): + ct.nyquist_plot(sys, indent_direction='left', blend_fraction=1.2) + + if __name__ == "__main__": # # Interactive mode: generate plots for manual viewing @@ -566,8 +604,8 @@ def test_no_indent_pole(): sys = 3 * (s+6)**2 / (s * (s**2 + 1e-4 * s + 1)) plt.figure() ct.nyquist_plot(sys) - ct.nyquist_plot(sys, max_curve_magnitude=15) - ct.nyquist_plot(sys, indent_radius=1e-6, max_curve_magnitude=25) + ct.nyquist_plot(sys, max_curve_magnitude=10) + ct.nyquist_plot(sys, indent_radius=1e-6, max_curve_magnitude=20) print("Unusual Nyquist plot") sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1]) @@ -595,3 +633,29 @@ def test_no_indent_pole(): plt.figure() cplt = ct.nyquist_plot([sys, sys1, sys2]) cplt.set_plot_title("Mixed FRD, tf data") + + plt.figure() + print("Jon How example") + test_nyquist_rescale() + + # + # Save the figures in a PDF file for later comparisons + # + import subprocess + from matplotlib.backends.backend_pdf import PdfPages + from datetime import date + + # Create the file to store figures + try: + git_info = subprocess.check_output( + ['git', 'describe'], text=True).strip() + except subprocess.CalledProcessError: + git_info = 'UNKNOWN-REPO-INFO' + pdf = PdfPages( + f'nyquist_gallery-{git_info}-{date.today().isoformat()}.pdf') + + # Go through each figure and save it + for fignum in plt.get_fignums(): + pdf.savefig(plt.figure(fignum)) + + pdf.close() diff --git a/doc/figures/freqplot-nyquist-custom.png b/doc/figures/freqplot-nyquist-custom.png index 5cd2c19d0..29e7e476d 100644 Binary files a/doc/figures/freqplot-nyquist-custom.png and b/doc/figures/freqplot-nyquist-custom.png differ diff --git a/doc/figures/freqplot-nyquist-default.png b/doc/figures/freqplot-nyquist-default.png index c511509fa..968987d2c 100644 Binary files a/doc/figures/freqplot-nyquist-default.png and b/doc/figures/freqplot-nyquist-default.png differ diff --git a/doc/response.rst b/doc/response.rst index 0058a500d..8ccdccba8 100644 --- a/doc/response.rst +++ b/doc/response.rst @@ -547,7 +547,7 @@ the computation of the Nyquist curve and the way the data are plotted: sys = ct.tf([1, 0.2], [1, 0, 1]) * ct.tf([1], [1, 0]) nyqresp = ct.nyquist_response(sys) nyqresp.plot( - max_curve_magnitude=6, max_curve_offset=1, + max_curve_magnitude=6, max_curve_offset=1, blend_fraction=0.05, arrows=[0, 0.15, 0.3, 0.6, 0.7, 0.925], title='Custom Nyquist plot') print("Encirclements =", nyqresp.count) diff --git a/examples/cds110-L8b_pvtol-complete-limits.ipynb b/examples/cds110-L8b_pvtol-complete-limits.ipynb index 0b482c865..b0ae80d76 100644 --- a/examples/cds110-L8b_pvtol-complete-limits.ipynb +++ b/examples/cds110-L8b_pvtol-complete-limits.ipynb @@ -830,13 +830,13 @@ "\n", "# Gain margin for Lx\n", "neg1overgm_x = -0.67 # vary this manually to find intersection with curve\n", - "color = cplt.lines[0][0].get_color()\n", + "color = cplt.lines[0, 0][0].get_color()\n", "plt.plot(neg1overgm_x, 0, color=color, marker='o', fillstyle='none')\n", "gm_x = -1/neg1overgm_x\n", "\n", "# Gain margin for Ly\n", "neg1overgm_y = -0.32 # vary this manually to find intersection with curve\n", - "color = cplt.lines[1][0].get_color()\n", + "color = cplt.lines[1, 0][0].get_color()\n", "plt.plot(neg1overgm_y, 0, color=color, marker='o', fillstyle='none')\n", "gm_y = -1/neg1overgm_y\n", "\n", @@ -885,13 +885,13 @@ "# Phase margin of Lx:\n", "th_pm_x = 0.14*np.pi\n", "th_plt_x = np.pi + th_pm_x\n", - "color = cplt.lines[0][0].get_color()\n", + "color = cplt.lines[0, 0][0].get_color()\n", "plt.plot(np.cos(th_plt_x), np.sin(th_plt_x), color=color, marker='o')\n", "\n", "# Phase margin of Ly\n", "th_pm_y = 0.19*np.pi\n", "th_plt_y = np.pi + th_pm_y\n", - "color = cplt.lines[1][0].get_color()\n", + "color = cplt.lines[1, 0][0].get_color()\n", "plt.plot(np.cos(th_plt_y), np.sin(th_plt_y), color=color, marker='o')\n", "\n", "print('Margins obtained visually:')\n", @@ -936,12 +936,12 @@ "\n", "# Stability margin:\n", "sm_x = 0.3 # vary this manually to find min which intersects\n", - "color = cplt.lines[0][0].get_color()\n", + "color = cplt.lines[0, 0][0].get_color()\n", "sm_circle = plt.Circle((-1, 0), sm_x, color=color, fill=False, ls=':')\n", "cplt.axes[0, 0].add_patch(sm_circle)\n", "\n", "sm_y = 0.5 # vary this manually to find min which intersects\n", - "color = cplt.lines[1][0].get_color()\n", + "color = cplt.lines[1, 0][0].get_color()\n", "sm_circle = plt.Circle((-1, 0), sm_y, color=color, fill=False, ls=':')\n", "cplt.axes[0, 0].add_patch(sm_circle)\n", "\n", @@ -1023,8 +1023,22 @@ } ], "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, "language_info": { - "name": "python" + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.1" } }, "nbformat": 4,