refactoring of describing_function_plot in response/plot + updated un…

…it tests
murrayrm · Sep 16, 2023 · Jul 1, 2023 · Jul 4, 2023 · Jul 4, 2023 · Jul 7, 2023
commit 83c9e4e8e80fe09d5b5797dc25cd247ff0fa9f86
diff --git a/control/descfcn.py b/control/descfcn.py
@@ -19,10 +19,12 @@
 from warnings import warn
 
 from .freqplot import nyquist_response
+from . import config
 
 __all__ = ['describing_function', 'describing_function_plot',
-           'DescribingFunctionNonlinearity', 'friction_backlash_nonlinearity',
-           'relay_hysteresis_nonlinearity', 'saturation_nonlinearity']
+           'describing_function_response', 'DescribingFunctionNonlinearity',
+           'friction_backlash_nonlinearity', 'relay_hysteresis_nonlinearity',
+           'saturation_nonlinearity']
 
 # Class for nonlinearities with a built-in describing function
 class DescribingFunctionNonlinearity():
@@ -205,14 +207,41 @@ def describing_function(
     # Return the values in the same shape as they were requested
     return retdf
 
+#
+# Describing function response/plot
+#
 
-def describing_function_plot(
-        H, F, A, omega=None, refine=True, label="%5.2g @ %-5.2g",
-        warn=None, **kwargs):
-    """Plot a Nyquist plot with a describing function for a nonlinear system.
+# Simple class to store the describing function response
+class DescribingFunctionResponse:
+    def __init__(self, response, N_vals, positions, intersections):
+        self.response = response
+        self.N_vals = N_vals
+        self.positions = positions
+        self.intersections = intersections
+
+    def plot(self, **kwargs):
+        return describing_function_plot(self, **kwargs)
+
+    # Implement iter, getitem, len to allow recovering the intersections
+    def __iter__(self):
+        return iter(self.intersections)
+
+    def __getitem__(self, index):
+        return list(self.__iter__())[index]
+
+    def __len__(self):
+        return len(self.intersections)
 
-    This function generates a Nyquist plot for a closed loop system consisting
-    of a linear system with a static nonlinear function in the feedback path.
+
+# Compute the describing function response + intersections
+def describing_function_response(
+        H, F, A, omega=None, refine=True, warn_nyquist=None,
+        plot=False, check_kwargs=True, **kwargs):
+    """Compute the describing function response of a system.
+
+    This function uses describing function analysis to analyze a closed
+    loop system consisting of a linear system with a static nonlinear
+    function in the feedback path.
 
     Parameters
     ----------
@@ -226,10 +255,7 @@ def describing_function_plot(
         List of amplitudes to be used for the describing function plot.
     omega : list, optional
         List of frequencies to be used for the linear system Nyquist curve.
-    label : str, optional
-        Formatting string used to label intersection points on the Nyquist
-        plot.  Defaults to "%5.2g @ %-5.2g".  Set to `None` to omit labels.
-    warn : bool, optional
+    warn_nyquist : bool, optional
         Set to True to turn on warnings generated by `nyquist_plot` or False
         to turn off warnings.  If not set (or set to None), warnings are
         turned off if omega is specified, otherwise they are turned on.
@@ -249,31 +275,27 @@ def describing_function_plot(
     >>> H_simple = ct.tf([8], [1, 2, 2, 1])
     >>> F_saturation = ct.saturation_nonlinearity(1)
     >>> amp = np.linspace(1, 4, 10)
-    >>> ct.describing_function_plot(H_simple, F_saturation, amp)  # doctest: +SKIP
+    >>> ct.describing_function_response(H_simple, F_saturation, amp)  # doctest: +SKIP
     [(3.343844998258643, 1.4142293090899216)]
 
     """
     # Decide whether to turn on warnings or not
-    if warn is None:
+    if warn_nyquist is None:
         # Turn warnings on unless omega was specified
-        warn = omega is None
+        warn_nyquist = omega is None
 
     # Start by drawing a Nyquist curve
     response = nyquist_response(
-        H, omega, warn_encirclements=warn, warn_nyquist=warn,
-        check_kwargs=False, **kwargs)
-    response.plot(**kwargs)
+        H, omega, warn_encirclements=warn_nyquist, warn_nyquist=warn_nyquist,
+        check_kwargs=check_kwargs, **kwargs)
     H_omega, H_vals = response.contour.imag, H(response.contour)
 
     # Compute the describing function
     df = describing_function(F, A)
     N_vals = -1/df
 
-    # Now add the describing function curve to the plot
-    plt.plot(N_vals.real, N_vals.imag)
-
     # Look for intersection points
-    intersections = []
+    positions, intersections = [], []
     for i in range(N_vals.size - 1):
         for j in range(H_vals.size - 1):
             intersect = _find_intersection(
@@ -306,17 +328,99 @@ def _cost(x):
                 else:
                     a_final, omega_final = res.x[0], res.x[1]
 
-            # Add labels to the intersection points
-            if isinstance(label, str):
-                pos = H(1j * omega_final)
-                plt.text(pos.real, pos.imag, label % (a_final, omega_final))
-            elif label is not None or label is not False:
-                raise ValueError("label must be formatting string or None")
+            pos = H(1j * omega_final)
 
             # Save the final estimate
+            positions.append(pos)
             intersections.append((a_final, omega_final))
 
-    return intersections
+    return DescribingFunctionResponse(
+        response, N_vals, positions, intersections)
+
+
+def describing_function_plot(
+        *sysdata, label="%5.2g @ %-5.2g", **kwargs):
+    """Plot a Nyquist plot with a describing function for a nonlinear system.
+
+    This function generates a Nyquist plot for a closed loop system
+    consisting of a linear system with a static nonlinear function in the
+    feedback path.
+
+    Parameters
+    ----------
+    H : LTI system
+        Linear time-invariant (LTI) system (state space, transfer function, or
+        FRD)
+    F : static nonlinear function
+        A static nonlinearity, either a scalar function or a single-input,
+        single-output, static input/output system.
+    A : list
+        List of amplitudes to be used for the describing function plot.
+    omega : list, optional
+        List of frequencies to be used for the linear system Nyquist curve.
+    refine : bool, optional
+        If True (default), refine the location of the intersection of the
+        Nyquist curve for the linear system and the describing function to
+        determine the intersection point
+    label : str, optional
+        Formatting string used to label intersection points on the Nyquist
+        plot.  Defaults to "%5.2g @ %-5.2g".  Set to `None` to omit labels.
+
+    Returns
+    -------
+    intersections : 1D array of 2-tuples or None
+        A list of all amplitudes and frequencies in which :math:`H(j\\omega)
+        N(a) = -1`, where :math:`N(a)` is the describing function associated
+        with `F`, or `None` if there are no such points.  Each pair represents
+        a potential limit cycle for the closed loop system with amplitude
+        given by the first value of the tuple and frequency given by the
+        second value.
+
+    Examples
+    --------
+    >>> H_simple = ct.tf([8], [1, 2, 2, 1])
+    >>> F_saturation = ct.saturation_nonlinearity(1)
+    >>> amp = np.linspace(1, 4, 10)
+    >>> ct.describing_function_plot(H_simple, F_saturation, amp)  # doctest: +SKIP
+    [(3.343844998258643, 1.4142293090899216)]
+
+    """
+    # Process keywords
+    warn_nyquist = config._process_legacy_keyword(
+        kwargs, 'warn', 'warn_nyquist', kwargs.pop('warn_nyquist', None))
+
+    if label not in (False, None) and not isinstance(label, str):
+        raise ValueError("label must be formatting string, False, or None")
+
+    # Get the describing function response
+    if len(sysdata) == 3:
+        sysdata = sysdata + (None, )    # set omega to default value
+    if len(sysdata) == 4:
+        dfresp = describing_function_response(
+            *sysdata, refine=kwargs.pop('refine', True),
+            warn_nyquist=warn_nyquist)
+    elif len(sysdata) == 1:
+        dfresp = sysdata[0]
+    else:
+        raise TypeError("1, 3, or 4 position arguments required")
+
+    # Create a list of lines for the output
+    out = np.empty(2, dtype=object)
+
+    # Plot the Nyquist response
+    out[0] = dfresp.response.plot(**kwargs)[0]
+
+    # Add the describing function curve to the plot
+    lines = plt.plot(dfresp.N_vals.real, dfresp.N_vals.imag)
+    out[1] = lines
+
+    # Label the intersection points
+    if label:
+        for pos, (a, omega) in zip(dfresp.positions, dfresp.intersections):
+            # Add labels to the intersection points
+            plt.text(pos.real, pos.imag, label % (a, omega))
+
+    return out
 
 
 # Utility function to figure out whether two line segments intersection

diff --git a/control/tests/descfcn_test.py b/control/tests/descfcn_test.py
@@ -12,6 +12,7 @@
 import numpy as np
 import control as ct
 import math
+import matplotlib.pyplot as plt
 from control.descfcn import saturation_nonlinearity, \
     friction_backlash_nonlinearity, relay_hysteresis_nonlinearity
 
@@ -137,7 +138,7 @@ def test_describing_function(fcn, amin, amax):
         ct.describing_function(fcn, -1)
 
 
-def test_describing_function_plot():
+def test_describing_function_response():
     # Simple linear system with at most 1 intersection
     H_simple = ct.tf([1], [1, 2, 2, 1])
     omega = np.logspace(-1, 2, 100)
@@ -147,12 +148,12 @@ def test_describing_function_plot():
     amp = np.linspace(1, 4, 10)
 
     # No intersection
-    xsects = ct.describing_function_plot(H_simple, F_saturation, amp, omega)
-    assert xsects == []
+    xsects = ct.describing_function_response(H_simple, F_saturation, amp, omega)
+    assert len(xsects) == 0
 
     # One intersection
     H_larger = H_simple * 8
-    xsects = ct.describing_function_plot(H_larger, F_saturation, amp, omega)
+    xsects = ct.describing_function_response(H_larger, F_saturation, amp, omega)
     for a, w in xsects:
         np.testing.assert_almost_equal(
             H_larger(1j*w),
@@ -163,12 +164,38 @@ def test_describing_function_plot():
     omega = np.logspace(-1, 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)
+    xsects = ct.describing_function_response(H_multiple, F_backlash, amp, omega)
     for a, w in xsects:
         np.testing.assert_almost_equal(
             -1/ct.describing_function(F_backlash, a),
             H_multiple(1j*w), decimal=5)
 
+
+def test_describing_function_plot():
+    # Simple linear system with at most 1 intersection
+    H_larger = ct.tf([1], [1, 2, 2, 1]) * 8
+    omega = np.logspace(-1, 2, 100)
+
+    # Saturation nonlinearity
+    F_saturation = ct.descfcn.saturation_nonlinearity(1)
+    amp = np.linspace(1, 4, 10)
+
+    # Plot via response
+    plt.clf()                                   # clear axes
+    response = ct.describing_function_response(
+        H_larger, F_saturation, amp, omega)
+    assert len(response.intersections) == 1
+    assert len(plt.gcf().get_axes()) == 0       # make sure there is no plot
+
+    out = response.plot()
+    assert len(plt.gcf().get_axes()) == 1       # make sure there is a plot
+    assert len(out[0]) == 4 and len(out[1]) == 1
+
+    # Call plot directly
+    out = ct.describing_function_plot(H_larger, F_saturation, amp, omega)
+    assert len(out[0]) == 4 and len(out[1]) == 1
+
+
 def test_describing_function_exceptions():
     # Describing function with non-zero bias
     with pytest.warns(UserWarning, match="asymmetric"):
@@ -194,3 +221,8 @@ def test_describing_function_exceptions():
     amp = np.linspace(1, 4, 10)
     with pytest.raises(ValueError, match="formatting string"):
         ct.describing_function_plot(H_simple, F_saturation, amp, label=1)
+
+    # Unrecognized keyword
+    with pytest.raises(TypeError, match="unrecognized keyword"):
+        ct.describing_function_response(
+            H_simple, F_saturation, amp, None, unknown=None)
diff --git a/control/tests/kwargs_test.py b/control/tests/kwargs_test.py
@@ -27,6 +27,7 @@
 import control.tests.stochsys_test as stochsys_test
 import control.tests.trdata_test as trdata_test
 import control.tests.timeplot_test as timeplot_test
+import control.tests.descfcn_test as descfcn_test
 
 @pytest.mark.parametrize("module, prefix", [
     (control, ""), (control.flatsys, "flatsys."), (control.optimal, "optimal.")
@@ -210,6 +211,8 @@ def test_response_plot_kwargs(data_fcn, plot_fcn, mimo):
     'create_estimator_iosystem': stochsys_test.test_estimator_errors,
     'create_statefbk_iosystem': statefbk_test.TestStatefbk.test_statefbk_errors,
     'describing_function_plot': test_matplotlib_kwargs,
+    'describing_function_response':
+        descfcn_test.test_describing_function_exceptions,
     'dlqe': test_unrecognized_kwargs,
     'dlqr': test_unrecognized_kwargs,
     'drss': test_unrecognized_kwargs,