Merge pull request #1155 from murrayrm/fix_nyquist_rescaling-24Mar2025

murrayrm · web-flow · commit bb47d594eea3 · 2025-07-04T14:18:34.000-07:00
Update Nyquist rescaling + other improvements
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
diff --git a/doc/figures/freqplot-nyquist-default.png b/doc/figures/freqplot-nyquist-default.png
diff --git 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
