Update Nyquist rescaling + other improvements

murrayrm · murrayrm · commit ca6f3dafd20c · 2025-06-21T01:34:58.000-04:00
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
@@ -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)
@@ -1891,21 +1898,36 @@ def _parse_linestyle(style_name, allow_false=False):
             splane_contour = np.log(response.contour) / response.dt
 
         # Find the different portions of the curve (with scaled pts marked)
+        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
         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 is None or \
+               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]
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,
@@ -428,6 +428,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 +527,34 @@ 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
+    cplt = ct.nyquist_plot(sys, indent_direction='left', label='default [0.15]')
+
+    # Sharper corner
+    cplt = ct.nyquist_plot(
+        sys*4, indent_direction='left',
+        max_curve_magnitude=17, blend_fraction=0.05, label='fraction=0.05')
+
+    # 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')
+
+    # No corner
+    cplt = ct.nyquist_plot(
+        sys*12, indent_direction='left',
+        max_curve_magnitude=19, blend_fraction=0, label='fraction=0')
+
+    # Bad value
+    with pytest.raises(ValueError, match="blend_fraction must be between"):
+        cplt = ct.nyquist_plot(
+            sys, indent_direction='left', blend_fraction=1.2)
+
+
 if __name__ == "__main__":
     #
     # Interactive mode: generate plots for manual viewing
@@ -566,8 +595,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 +624,25 @@ 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
+    git_info = subprocess.check_output(['git', 'describe'], text=True).strip()
+    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()
