add warning if encirclements is not near an integer

murrayrm · murrayrm · commit 36b0d0ca6efd · 2022-04-13T22:56:55.000-07:00
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -522,6 +522,7 @@ def gen_zero_centered_series(val_min, val_max, period):
     'nyquist.mirror_style': ['--', ':'],        # style for mirror curve
     'nyquist.arrows': 2,                        # number of arrors 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
@@ -610,6 +611,11 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
     arrow_style : matplotlib.patches.ArrowStyle, optional
         Define style used for Nyquist curve arrows (overrides `arrow_size`).
 
+    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
+        be set using config.defaults['nyquist.encirclement_threshold'].
+
     indent_direction : str, optional
         For poles on the imaginary axis, set the direction of indentation to
         be 'right' (default), 'left', or 'none'.
@@ -716,6 +722,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
     arrow_style = config._get_param('nyquist', 'arrow_style', kwargs, None)
     indent_radius = config._get_param(
         'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True)
+    encirclement_threshold = config._get_param(
+        'nyquist', 'encirclement_threshold', kwargs,
+        _nyquist_defaults, pop=True)
     indent_direction = config._get_param(
         'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True)
     indent_points = config._get_param(
@@ -749,8 +758,11 @@ def _parse_linestyle(style_name, allow_false=False):
     if not hasattr(syslist, '__iter__'):
         syslist = (syslist,)
 
+    # Determine the range of frequencies to use, based on args/features
     omega, omega_range_given = _determine_omega_vector(
         syslist, omega, omega_limits, omega_num)
+
+    # If omega was not specified explicitly, start at omega = 0
     if not omega_range_given:
         if omega_num_given:
             # Just reset the starting point
@@ -851,6 +863,7 @@ def _parse_linestyle(style_name, allow_false=False):
                     first_point = 0
                     start_freq = 0
 
+                # Find the frequencies after the pole frequency
                 above_points = np.argwhere(
                     splane_contour.imag - abs(p.imag) > indent_radius)
                 last_point = above_points[0].item()
@@ -899,7 +912,15 @@ def _parse_linestyle(style_name, allow_false=False):
 
         # 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))
+        encirclements = np.sum(np.diff(phase)) / np.pi
+        count = int(np.round(encirclements, 0))
+
+        # Let the user know if the count might not make sense
+        if abs(encirclements - count) > encirclement_threshold:
+            warnings.warn(
+                "number of encirclements was a non-integer value; this can"
+                " happen is contour is not closed, possibly based on a"
+                " frequency range that does not include zero.")
 
         #
         # Make sure that the enciriclements match the Nyquist criterion
@@ -1533,6 +1554,7 @@ def _determine_omega_vector(syslist, omega_in, omega_limits, omega_num,
                                     num=omega_num, endpoint=True)
     else:
         omega_out = np.copy(omega_in)
+
     return omega_out, omega_range_given
 
 
diff --git a/control/tests/descfcn_test.py b/control/tests/descfcn_test.py
@@ -140,7 +140,7 @@ def test_describing_function(fcn, amin, amax):
 def test_describing_function_plot():
     # Simple linear system with at most 1 intersection
     H_simple = ct.tf([1], [1, 2, 2, 1])
-    omega = np.logspace(-1, 2, 100)
+    omega = np.logspace(-2, 2, 100)
 
     # Saturation nonlinearity
     F_saturation = ct.descfcn.saturation_nonlinearity(1)
@@ -160,7 +160,7 @@ def test_describing_function_plot():
 
     # Multiple intersections
     H_multiple = H_simple * ct.tf(*ct.pade(5, 4)) * 4
-    omega = np.logspace(-1, 3, 50)
+    omega = np.logspace(-2, 3, 50)
     F_backlash = ct.descfcn.friction_backlash_nonlinearity(1)
     amp = np.linspace(0.6, 5, 50)
     xsects = ct.describing_function_plot(H_multiple, F_backlash, amp, omega)
diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py
@@ -81,8 +81,9 @@ def test_nyquist_basic(ss_siso):
         tf_siso, plot=False, return_contour=True, omega_num=20)
     assert len(contour) == 20
 
-    count, contour = nyquist_plot(
-        tf_siso, plot=False, omega_limits=(1, 100), return_contour=True)
+    with pytest.warns(UserWarning, match="encirclements was a non-integer"):
+        count, contour = nyquist_plot(
+            tf_siso, plot=False, omega_limits=(1, 100), return_contour=True)
     assert_allclose(contour[0], 1j)
     assert_allclose(contour[-1], 100j)
 
diff --git a/control/tests/nyquist_test.py b/control/tests/nyquist_test.py
@@ -219,6 +219,17 @@ def test_nyquist_encirclements():
     plt.title("Pole at the origin; encirclements = %d" % count)
     assert _Z(sys) == count + _P(sys)
 
+    # Non-integer number of encirclements
+    plt.figure();
+    sys = 1 / (s**2 + s + 1)
+    with pytest.warns(UserWarning, match="encirclements was a non-integer"):
+        count = ct.nyquist_plot(sys, omega_limits=[0.5, 1e3])
+    with pytest.warns(None) as records:
+        count = ct.nyquist_plot(
+            sys, omega_limits=[0.5, 1e3], encirclement_threshold=0.2)
+    assert len(records) == 0
+    plt.title("Non-integer number of encirclements [%g]" % count)
+
 
 @pytest.fixture
 def indentsys():
