address @roryyorke review comments

murrayrm · murrayrm · commit 17ed781721b3 · 2021-02-06T13:17:56.000-08:00
diff --git a/control/descfcn.py b/control/descfcn.py
@@ -3,11 +3,11 @@
 # RMM, 23 Jan 2021
 #
 # This module adds functions for carrying out analysis of systems with
-# static nonlinear feedback functions using describing functions.
+# memoryless nonlinear feedback functions using describing functions.
 #
 
 """The :mod:~control.descfcn` module contains function for performing
-closed loop analysis of systems with static nonlinearities using
+closed loop analysis of systems with memoryless nonlinearities using
 describing function analysis.
 
 """
@@ -26,21 +26,21 @@
 
 # Class for nonlinearities with a built-in describing function
 class DescribingFunctionNonlinearity():
-    """Base class for nonlinear functions with a describing function
+    """Base class for nonlinear systems with a describing function
 
     This class is intended to be used as a base class for nonlinear functions
-    that have a analytically defined describing function (accessed via the
-    :meth:`describing_function` method).  Objects using this class should also
-    implement a `call` method that evaluates the nonlinearity at a given point
-    and an `_isstatic` method that is `True` if the nonlinearity has no
-    internal state.
+    that have an analytically defined describing function.  Subclasses should
+    override the `__call__` and `describing_function` methods and (optionally)
+    the `_isstatic` method (should be `False` if `__call__` updates the
+    instance state).
 
     """
     def __init__(self):
-        """Initailize a describing function nonlinearity"""
+        """Initailize a describing function nonlinearity (optional)"""
         pass
 
     def __call__(self, A):
+        """Evaluate the nonlinearity at a (scalar) input value"""
         raise NotImplementedError(
             "__call__() not implemented for this function (internal error)")
 
@@ -56,9 +56,15 @@ def describing_function(self, A):
             "describing function not implemented for this function")
 
     def _isstatic(self):
-        """Return True if the function has not internal state"""
-        raise NotImplementedError(
-            "_isstatic() not implemented for this function (internal error)")
+        """Return True if the function has no internal state (memoryless)
+
+        This internal function is used to optimize numerical computation of
+        the describing function.  It can be set to `True` if the instance
+        maintains no internal memory of the instance state.  Assumed False by
+        default.
+
+        """
+        return False
 
     # Utility function used to compute common describing functions
     def _f(self, x):
@@ -70,10 +76,10 @@ def describing_function(
         F, A, num_points=100, zero_check=True, try_method=True):
     """Numerical compute the describing function of a nonlinear function
 
-    The describing function of a static nonlinear function is given by
-    magnitude and phase of the first harmonic of the function when evaluated
-    along a sinusoidal input :math:`a \\sin \\omega t`.  This function returns
-    the magnitude and phase of the describing function at amplitude :math:`A`.
+    The describing function of a nonlinearity is given by magnitude and phase
+    of the first harmonic of the function when evaluated along a sinusoidal
+    input :math:`A \\sin \\omega t`.  This function returns the magnitude and
+    phase of the describing function at amplitude :math:`A`.
 
     Parameters
     ----------
@@ -119,11 +125,15 @@ def describing_function(
     """
     # If there is an analytical solution, trying using that first
     if try_method and hasattr(F, 'describing_function'):
-        # Go through all of the amplitudes we were given
-        df = []
-        for a in np.atleast_1d(A):
-            df.append(F.describing_function(a))
-        return np.array(df).reshape(np.shape(A))
+        try:
+            # Go through all of the amplitudes we were given
+            df = []
+            for a in np.atleast_1d(A):
+                df.append(F.describing_function(a))
+            return np.array(df).reshape(np.shape(A))
+        except NotImplementedError:
+            # Drop through and do the numerical computation
+            pass
 
     #
     # The describing function of a nonlinear function F() can be computed by
@@ -136,24 +146,24 @@ def describing_function(
     #
     #    N(A) = M_1(A) e^{j \phi_1(A)} / A
     #
-    # To compute this, we compute F(\theta) for \theta between 0 and 2 \pi,
-    # use the identities
+    # To compute this, we compute F(A \sin\theta) for \theta between 0 and 2
+    # \pi, use the identities
     #
     #   \sin(\theta + \phi) = \sin\theta \cos\phi + \cos\theta \sin\phi
     #   \int_0^{2\pi} \sin^2 \theta d\theta = \pi
     #   \int_0^{2\pi} \cos^2 \theta d\theta = \pi
     #
     # and then integrate the product against \sin\theta and \cos\theta to obtain
     #
-    #   \int_0^{2\pi} F(a\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi
-    #   \int_0^{2\pi} F(a\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi
+    #   \int_0^{2\pi} F(A\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi
+    #   \int_0^{2\pi} F(A\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi
     #
     # From these we can compute M1 and \phi.
     #
 
-    # Evaluate over a full range of angles
-    theta = np.linspace(0, 2*np.pi, num_points)
-    dtheta = theta[1] - theta[0]
+    # Evaluate over a full range of angles (leave off endpoint a la DFT)
+    theta, dtheta = np.linspace(
+        0, 2*np.pi, num_points, endpoint=False, retstep=True)
     sin_theta = np.sin(theta)
     cos_theta = np.cos(theta)
 
@@ -175,10 +185,10 @@ def describing_function(
         elif a < 0:
             raise ValueError("cannot evaluate describing function for A < 0")
 
-        # Save the scaling factor for to make the formulas simpler
+        # Save the scaling factor to make the formulas simpler
         scale = dtheta / np.pi / a
 
-        # Evaluate the function (twice) along a sinusoid (for internal state)
+        # Evaluate the function along a sinusoid
         F_eval = np.array([F(x) for x in a*sin_theta]).squeeze()
 
         # Compute the prjections onto sine and cosine
@@ -353,7 +363,7 @@ def __init__(self, ub=1, lb=None):
         self.ub = ub
 
     def __call__(self, x):
-        return np.maximum(self.lb, np.minimum(x, self.ub))
+        return np.clip(x, self.lb, self.ub)
 
     def _isstatic(self):
         return True
@@ -453,18 +463,12 @@ def __init__(self, b):
         self.center = 0         # current center position
 
     def __call__(self, x):
-        # Convert input to an array
-        x_array = np.array(x)
-
-        y = []
-        for x in np.atleast_1d(x_array):
-            # If we are outside the backlash, move and shift the center
-            if x - self.center > self.b/2:
-                self.center = x - self.b/2
-            elif x - self.center < -self.b/2:
-                self.center = x + self.b/2
-            y.append(self.center)
-        return(np.array(y).reshape(x_array.shape))
+        # If we are outside the backlash, move and shift the center
+        if x - self.center > self.b/2:
+            self.center = x - self.b/2
+        elif x - self.center < -self.b/2:
+            self.center = x + self.b/2
+        return self.center
 
     def _isstatic(self):
         return False
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -522,10 +522,8 @@ def gen_zero_centered_series(val_min, val_max, period):
 
 def nyquist_plot(
         syslist, omega=None, plot=True, omega_limits=None, omega_num=None,
-        label_freq=0, arrowhead_length=0.1, arrowhead_width=0.1,
-        mirror='--', color=None, *args, **kwargs):
-    """
-    Nyquist plot for a system
+        label_freq=0, color=None, mirror='--', arrowhead_length=0.1,
+        arrowhead_width=0.1, *args, **kwargs): """Nyquist plot for a system
 
     Plots a Nyquist plot for the system over a (optional) frequency range.
 
diff --git a/control/iosys.py b/control/iosys.py
@@ -325,6 +325,10 @@ def __neg__(sys):
         # Return the newly created system
         return newsys
 
+    def _isstatic(self):
+        """Check to see if a system is a static system (no states)"""
+        return self.nstates == 0
+
     # Utility function to parse a list of signals
     def _process_signal_list(self, signals, prefix='s'):
         if signals is None:
@@ -450,10 +454,6 @@ def issiso(self):
         """Check to see if a system is single input, single output"""
         return self.ninputs == 1 and self.noutputs == 1
 
-    def isstatic(self):
-        """Check to see if a system is a static system (no states)"""
-        return self.nstates == 0
-
     def feedback(self, other=1, sign=-1, params={}):
         """Feedback interconnection between two input/output systems
 
@@ -812,9 +812,28 @@ def __init__(self, updfcn, outfcn=None, inputs=None, outputs=None,
         self._current_params = params.copy()
 
     # Return the value of a static nonlinear system
-    def __call__(sys, u, squeeze=None, params=None):
+    def __call__(sys, u, params=None, squeeze=None):
+        """Evaluate a (static) nonlinearity at a given input value
+
+        If a nonlinear I/O system has not internal state, then evaluating the
+        system at an input `u` gives the output `y = F(u)`, determined by the
+        output function.
+
+        Parameters
+        ----------
+        params : dict, optional
+            Parameter values for the system. Passed to the evaluation function
+            for the system as default values, overriding internal defaults.
+        squeeze : bool, optional
+            If True and if the system has a single output, return the system
+            output as a 1D array rather than a 2D array.  If False, return the
+            system output as a 2D array even if the system is SISO.  Default
+            value set by config.defaults['control.squeeze_time_response'].
+
+        """
+
         # Make sure the call makes sense
-        if not sys.isstatic():
+        if not sys._isstatic():
             raise TypeError(
                 "function evaluation is only supported for static "
                 "input/output systems")
diff --git a/control/tests/descfcn_test.py b/control/tests/descfcn_test.py
@@ -17,10 +17,10 @@
 
 
 # Static function via a class
-class saturation_class():
+class saturation_class:
     # Static nonlinear saturation function
     def __call__(self, x, lb=-1, ub=1):
-        return np.maximum(lb, np.minimum(x, ub))
+        return np.clip(x, lb, ub)
 
     # Describing function for a saturation function
     def describing_function(self, a):
@@ -33,7 +33,7 @@ def describing_function(self, a):
 
 # Static function without a class
 def saturation(x):
-    return np.maximum(-1, np.minimum(x, 1))
+    return np.clip(x, -1, 1)
 
 
 # Static nonlinear system implementing saturation
@@ -47,7 +47,7 @@ def _satfcn(t, x, u, params):
 
 def test_static_nonlinear_call(satsys):
     # Make sure that the saturation system is a static nonlinearity
-    assert satsys.isstatic()
+    assert satsys._isstatic()
 
     # Make sure the saturation function is doing the right computation
     input = [-2, -1, -0.5, 0, 0.5, 1, 2]
@@ -61,7 +61,7 @@ def test_static_nonlinear_call(satsys):
     np.testing.assert_array_equal(satsys([0.]), [0.])
 
     # Test SIMO nonlinearity
-    def _simofcn(t, x, u, params={}):
+    def _simofcn(t, x, u, params):
         return np.array([np.cos(u), np.sin(u)])
     simo_sys = ct.NonlinearIOSystem(None, outfcn=_simofcn, input=1, output=2)
     np.testing.assert_array_equal(simo_sys([0.]), [1, 0])
@@ -72,7 +72,6 @@ def _misofcn(t, x, u, params={}):
         return np.array([np.sin(u[0]) * np.cos(u[1])])
     miso_sys = ct.NonlinearIOSystem(None, outfcn=_misofcn, input=2, output=1)
     np.testing.assert_array_equal(miso_sys([0, 0]), [0])
-    np.testing.assert_array_equal(miso_sys([0, 0]), [0])
     np.testing.assert_array_equal(miso_sys([0, 0], squeeze=True), [0])
 
 
@@ -85,6 +84,10 @@ def test_saturation_describing_function(satsys):
     df_anal = [satfcn.describing_function(a) for a in amprange]
 
     # Compute describing function for a static function
+    df_fcn = [ct.describing_function(saturation, a) for a in amprange]
+    np.testing.assert_almost_equal(df_fcn, df_anal, decimal=3)
+
+    # Compute describing function for a describing function nonlinearity
     df_fcn = [ct.describing_function(satfcn, a) for a in amprange]
     np.testing.assert_almost_equal(df_fcn, df_anal, decimal=3)
 
@@ -100,6 +103,15 @@ def test_saturation_describing_function(satsys):
     with pytest.raises(ValueError, match="cannot evaluate"):
         ct.describing_function(saturation, -1)
 
+    # Create describing function nonlinearity w/out describing_function method
+    # and make sure it drops through to the underlying computation
+    class my_saturation(ct.DescribingFunctionNonlinearity):
+        def __call__(self, x):
+            return saturation(x)
+    satfcn_nometh = my_saturation()
+    df_nometh = [ct.describing_function(satfcn_nometh, a) for a in amprange]
+    np.testing.assert_almost_equal(df_nometh, df_anal, decimal=3)
+
 
 @pytest.mark.parametrize("fcn, amin, amax", [
     [saturation_nonlinearity(1), 0, 10],
@@ -118,7 +130,7 @@ def test_describing_function(fcn, amin, amax):
 
     # Make sure the describing function method also works
     df_meth = ct.describing_function(fcn, amprange, zero_check=False)
-    np.testing.assert_almost_equal(df_meth, df_anal, decimal=1)
+    np.testing.assert_almost_equal(df_meth, df_anal)
 
     # Make sure that evaluation at negative amplitude generates an exception
     with pytest.raises(ValueError, match="cannot evaluate"):
diff --git a/doc/descfcn.rst b/doc/descfcn.rst
@@ -8,8 +8,8 @@ For nonlinear systems consisting of a feedback connection between a
 linear system and a static nonlinearity, it is possible to obtain a
 generalization of Nyquist's stability criterion based on the idea of
 describing functions.  The basic concept involves approximating the
-response of a static nonlinearity to an input :math:`u = A e^{\omega
-t}` as an output :math:`y = N(A) (A e^{\omega t})`, where :math:`N(A)
+response of a static nonlinearity to an input :math:`u = A e^{j \omega
+t}` as an output :math:`y = N(A) (A e^{j \omega t})`, where :math:`N(A)
 \in \mathbb{C}` represents the (amplitude-dependent) gain and phase
 associated with the nonlinearity.
 
diff --git a/examples/describing_functions.ipynb b/examples/describing_functions.ipynb
@@ -122,7 +122,7 @@
     "backlash = ct.backlash_nonlinearity(0.5)\n",
     "theta = np.linspace(0, 2*np.pi, 50)\n",
     "x = np.sin(theta)\n",
-    "plt.plot(x, backlash(x))\n",
+    "plt.plot(x, [backlash(z) for z in x])\n",
     "plt.xlabel(\"Input, x\")\n",
     "plt.ylabel(\"Output, y = backlash(x)\")\n",
     "plt.title(\"Input/output map for a backlash nonlinearity\");"
