python-control · murrayrm · Sep 16, 2023 · Jul 1, 2023 · Jul 4, 2023 · Jul 4, 2023
diff --git a/LICENSE b/LICENSE
@@ -1,4 +1,5 @@
 Copyright (c) 2009-2016 by California Institute of Technology
+Copyright (c) 2016-2023 by python-control developers
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without

diff --git a/control/config.py b/control/config.py
@@ -326,13 +326,13 @@ def use_legacy_defaults(version):
 #
 # Use this function to handle a legacy keyword that has been renamed.  This
 # function pops the old keyword off of the kwargs dictionary and issues a
-# warning.  if both the old and new keyword are present, a ControlArgument
+# warning.  If both the old and new keyword are present, a ControlArgument
 # exception is raised.
 #
 def _process_legacy_keyword(kwargs, oldkey, newkey, newval):
     if kwargs.get(oldkey) is not None:
         warnings.warn(
-            f"keyworld '{oldkey}' is deprecated; use '{newkey}'",
+            f"keyword '{oldkey}' is deprecated; use '{newkey}'",
             DeprecationWarning)
         if newval is not None:
             raise ControlArgument(

diff --git a/control/ctrlutil.py b/control/ctrlutil.py
@@ -86,18 +86,9 @@ def unwrap(angle, period=2*math.pi):
     return angle
 
 def issys(obj):
-    """Return True if an object is a Linear Time Invariant (LTI) system,
-    otherwise False.
+    """Deprecated function to check if an object is an LTI system.
 
-    Examples
-    --------
-    >>> G = ct.tf([1], [1, 1])
-    >>> ct.issys(G)
-    True
-
-    >>> K = np.array([[1, 1]])
-    >>> ct.issys(K)
-    False
+    Use isinstance(obj, ct.LTI)
 
     """
     warnings.warn("issys() is deprecated; use isinstance(obj, ct.LTI)",

diff --git a/control/descfcn.py b/control/descfcn.py
@@ -18,9 +18,11 @@
 import scipy
 from warnings import warn
 
-from .freqplot import nyquist_plot
+from .freqplot import nyquist_response
+from . import config
 
 __all__ = ['describing_function', 'describing_function_plot',
+           'describing_function_response', 'DescribingFunctionResponse',
            'DescribingFunctionNonlinearity', 'friction_backlash_nonlinearity',
            'relay_hysteresis_nonlinearity', 'saturation_nonlinearity']
 
@@ -205,14 +207,74 @@ 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:
+    """Results of describing function analysis.
+
+    Describing functions allow analysis of a linear I/O systems with a
+    static nonlinear feedback function.  The DescribingFunctionResponse
+    class is used by the :func:`~control.describing_function_response`
+    function to return the results of a describing function analysis.  The
+    response object can be used to obtain information about the describing
+    function analysis or generate a Nyquist plot showing the frequency
+    response of the linear systems and the describing function for the
+    nonlinear element.
+
+    Attributes
+    ----------
+    response : :class:`~control.FrequencyResponseData`
+        Frequency response of the linear system component of the system.
+    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.
+    N_vals : complex array
+        Complex value of the describing function.
+    positions : list of complex
+        Location of the intersections in the complex plane.
+
+    """
+    def __init__(self, response, N_vals, positions, intersections):
+        """Create a describing function response data object."""
+        self.response = response
+        self.N_vals = N_vals
+        self.positions = positions
+        self.intersections = intersections
+
+    def plot(self, **kwargs):
+        """Plot the results of a describing function analysis.
+
+        See :func:`~control.describing_function_plot` for details.
+        """
+        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,53 +288,53 @@ 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.
 
     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.
+    response : :class:`~control.DescribingFunctionResponse` object
+        Response object that contains the result of the describing function
+        analysis.  The following information can be retrieved from this
+        object:
+    response.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
+    >>> response = ct.describing_function_response(H_simple, F_saturation, amp)
+    >>> response.intersections  # doctest: +SKIP
     [(3.343844998258643, 1.4142293090899216)]
+    >>> lines = response.plot()
 
     """
     # 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
-    count, contour = nyquist_plot(
-        H, omega, plot=True, return_contour=True,
-        warn_encirclements=warn, warn_nyquist=warn, **kwargs)
-    H_omega, H_vals = contour.imag, H(contour)
+    response = nyquist_response(
+        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(
@@ -305,17 +367,114 @@ 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):
+    """describing_function_plot(data, *args, **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.
+
+    The function may be called in one of two forms:
+
+        describing_function_plot(response[, options])
+
+        describing_function_plot(H, F, A[, omega[, options]])
+
+    In the first form, the response should be generated using the
+    :func:`~control.describing_function_response` function.  In the second
+    form, that function is called internally, with the listed arguments.
+
+    Parameters
+    ----------
+    data : :class:`~control.DescribingFunctionData`
+        A describing function response data object created by
+        :func:`~control.describing_function_response`.
+    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. If not specified (or None), frequencies are computed
+        automatically based on the properties of the linear system.
+    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
+    -------
+    lines : 1D array of Line2D
+        Arrray of Line2D objects for each line in the plot.  The first
+        element of the array is a list of lines (typically only one) for
+        the Nyquist plot of the linear I/O styem.  The second element of
+        the array is a list of lines (typically only one) for the
+        describing function curve.
+
+    Examples
+    --------
+    >>> H_simple = ct.tf([8], [1, 2, 2, 1])
+    >>> F_saturation = ct.saturation_nonlinearity(1)
+    >>> amp = np.linspace(1, 4, 10)
+    >>> lines = ct.describing_function_plot(H_simple, F_saturation, amp)
+
+    """
+    # 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