rename method _rhs in iosys to dynamics

sawyerbfuller · Feb 16, 2021 · Feb 16, 2021 · Feb 16, 2021 · Feb 17, 2021 · Feb 17, 2021
commit 033b236c4df611987e19de550b433af282708af2
diff --git a/control/iosys.py b/control/iosys.py
@@ -111,7 +111,7 @@ class for a set of subclasses that are used to implement specific
     The :class:`~control.InputOuputSystem` class (and its subclasses) makes
     use of two special methods for implementing much of the work of the class:
 
-    * _rhs(t, x, u): compute the right hand side of the differential or
+    * dynamics(t, x, u): compute the right hand side of the differential or
       difference equation for the system.  This must be specified by the
       subclass for the system.
 
@@ -353,28 +353,69 @@ def _process_signal_list(self, signals, prefix='s'):
     # Find a signal by name
     def _find_signal(self, name, sigdict): return sigdict.get(name, None)
 
-    # Update parameters used for _rhs, _out (used by subclasses)
+    # Update parameters used for dynamics, _out (used by subclasses)
     def _update_params(self, params, warning=False):
         if (warning):
             warn("Parameters passed to InputOutputSystem ignored.")
 
-    def _rhs(self, t, x, u):
-        """Evaluate right hand side of a differential or difference equation.
+    def dynamics(self, t, x, u):
+        """Compute the dynamics of a differential or difference equation.
 
-        Private function used to compute the right hand side of an
-        input/output system model.
+        Given time `t`, input `u` and state `x`, returns the dynamics of the
+        system. If the system is continuous, returns the time derivative
 
+            ``dx/dt = f(t, x, u)``
+
+        where `f` is the system's (possibly nonlinear) dynamics function.
+        If the system is discrete-time, returns the next value of `x`:
+
+            ``x[t+dt] = f(t, x[t], u[t])``
+
+        Where `t` is a scalar.
+
+        The inputs `x` and `u` must be of the correct length.
+
+        Parameters
+        ----------
+        t : float
+            the time at which to evaluate
+        x : array_like
+            current state
+        u : array_like
+            input
+
+        Returns
+        -------
+        `dx/dt` or `x[t+dt]` : ndarray
         """
-        NotImplemented("Evaluation not implemented for system of type ",
+
+        NotImplemented("Dynamics not implemented for system of type ",
                        type(self))
 
     def _out(self, t, x, u, params={}):
-        """Evaluate the output of a system at a given state, input, and time
+        """Compute the output of the system
+
+        Given time `t`, input `u` and state `x`, returns the output of the
+        system:
+
+            ``y = g(t, x, u)``
 
-        Private function used to compute the output of of an input/output
-        system model given the state, input, parameters, and time.
+        The inputs `x` and `u` must be of the correct length.
 
+        Parameters
+        ----------
+        t : float
+            the time at which to evaluate
+        x : array_like
+            current state
+        u : array_like
+            input
+
+        Returns
+        -------
+        y : ndarray
         """
+
         # If no output function was defined in subclass, return state
         return x
 
@@ -533,7 +574,7 @@ def linearize(self, x0, u0, t=0, params={}, eps=1e-6,
         """
         #
         # If the linearization is not defined by the subclass, perform a
-        # numerical linearization use the `_rhs()` and `_out()` member
+        # numerical linearization use the `dynamics()` and `_out()` member
         # functions.
         #
 
@@ -554,7 +595,7 @@ def linearize(self, x0, u0, t=0, params={}, eps=1e-6,
         self._update_params(params)
 
         # Compute the nominal value of the update law and output
-        F0 = self._rhs(t, x0, u0)
+        F0 = self.dynamics(t, x0, u0)
         H0 = self._out(t, x0, u0)
 
         # Create empty matrices that we can fill up with linearizations
@@ -567,14 +608,14 @@ def linearize(self, x0, u0, t=0, params={}, eps=1e-6,
         for i in range(nstates):
             dx = np.zeros((nstates,))
             dx[i] = eps
-            A[:, i] = (self._rhs(t, x0 + dx, u0) - F0) / eps
+            A[:, i] = (self.dynamics(t, x0 + dx, u0) - F0) / eps
             C[:, i] = (self._out(t, x0 + dx, u0) - H0) / eps
 
         # Perturb each of the input variables and compute linearization
         for i in range(ninputs):
             du = np.zeros((ninputs,))
             du[i] = eps
-            B[:, i] = (self._rhs(t, x0, u0 + du) - F0) / eps
+            B[:, i] = (self.dynamics(t, x0, u0 + du) - F0) / eps
             D[:, i] = (self._out(t, x0, u0 + du) - H0) / eps
 
         # Create the state space system
@@ -694,7 +735,7 @@ def _update_params(self, params={}, warning=True):
         if params and warning:
             warn("Parameters passed to LinearIOSystems are ignored.")
 
-    def _rhs(self, t, x, u):
+    def dynamics(self, t, x, u):
         # Convert input to column vector and then change output to 1D array
         xdot = np.dot(self.A, np.reshape(x, (-1, 1))) \
             + np.dot(self.B, np.reshape(u, (-1, 1)))
@@ -863,7 +904,7 @@ def _update_params(self, params, warning=False):
         self._current_params = self.params.copy()
         self._current_params.update(params)
 
-    def _rhs(self, t, x, u):
+    def dynamics(self, t, x, u):
         xdot = self.updfcn(t, x, u, self._current_params) \
             if self.updfcn is not None else []
         return np.array(xdot).reshape((-1,))
@@ -1033,7 +1074,7 @@ def _update_params(self, params, warning=False):
             local.update(params)        # update with locally passed parameters
             sys._update_params(local, warning=warning)
 
-    def _rhs(self, t, x, u):
+    def dynamics(self, t, x, u):
         # Make sure state and input are vectors
         x = np.array(x, ndmin=1)
         u = np.array(u, ndmin=1)
@@ -1047,7 +1088,7 @@ def _rhs(self, t, x, u):
         for sys in self.syslist:
             # Update the right hand side for this subsystem
             if sys.nstates != 0:
-                xdot[state_index:state_index + sys.nstates] = sys._rhs(
+                xdot[state_index:state_index + sys.nstates] = sys.dynamics(
                     t, x[state_index:state_index + sys.nstates],
                     ulist[input_index:input_index + sys.ninputs])
 
@@ -1497,14 +1538,14 @@ def input_output_response(sys, T, U=0., X0=0, params={}, method='RK45',
 
     # Create a lambda function for the right hand side
     u = sp.interpolate.interp1d(T, U, fill_value="extrapolate")
-    def ivp_rhs(t, x): return sys._rhs(t, x, u(t))
+    def ivp_dynamics(t, x): return sys.dynamics(t, x, u(t))
 
     # Perform the simulation
     if isctime(sys):
         if not hasattr(sp.integrate, 'solve_ivp'):
             raise NameError("scipy.integrate.solve_ivp not found; "
                             "use SciPy 1.0 or greater")
-        soln = sp.integrate.solve_ivp(ivp_rhs, (T0, Tf), X0, t_eval=T,
+        soln = sp.integrate.solve_ivp(ivp_dynamics, (T0, Tf), X0, t_eval=T,
                                       method=method, vectorized=False)
 
         # Compute the output associated with the state (and use sys.out to
@@ -1549,7 +1590,7 @@ def ivp_rhs(t, x): return sys._rhs(t, x, u(t))
             y.append(sys._out(T[i], x, u(T[i])))
 
             # Update the state for the next iteration
-            x = sys._rhs(T[i], x, u(T[i]))
+            x = sys.dynamics(T[i], x, u(T[i]))
 
         # Convert output to numpy arrays
         soln.y = np.transpose(np.array(soln.y))
@@ -1670,8 +1711,8 @@ def find_eqpt(sys, x0, u0=[], y0=None, t=0, params={},
         if y0 is None:
             # Take u0 as fixed and minimize over x
             # TODO: update to allow discrete time systems
-            def ode_rhs(z): return sys._rhs(t, z, u0)
-            result = root(ode_rhs, x0, **kw)
+            def ode_dynamics(z): return sys.dynamics(t, z, u0)
+            result = root(ode_dynamics, x0, **kw)
             z = (result.x, u0, sys._out(t, result.x, u0))
         else:
             # Take y0 as fixed and minimize over x and u
@@ -1680,7 +1721,7 @@ def rootfun(z):
                 x, u = np.split(z, [nstates])
                 # TODO: update to allow discrete time systems
                 return np.concatenate(
-                    (sys._rhs(t, x, u), sys._out(t, x, u) - y0), axis=0)
+                    (sys.dynamics(t, x, u), sys._out(t, x, u) - y0), axis=0)
             z0 = np.concatenate((x0, u0), axis=0)   # Put variables together
             result = root(rootfun, z0, **kw)        # Find the eq point
             x, u = np.split(result.x, [nstates])    # Split result back in two
@@ -1782,7 +1823,7 @@ def rootfun(z):
             u[input_vars] = z[nstate_vars:]
 
             # Compute the update and output maps
-            dx = sys._rhs(t, x, u) - dx0
+            dx = sys.dynamics(t, x, u) - dx0
             if dtime:
                 dx -= x           # TODO: check
             dy = sys._out(t, x, u) - y0

diff --git a/control/sisotool.py b/control/sisotool.py
@@ -101,7 +101,7 @@ def sisotool(sys, kvect=None, xlim_rlocus=None, ylim_rlocus=None,
 
     # First time call to setup the bode and step response plots
     _SisotoolUpdate(sys, fig,
-        1 if kvect is None else kvect[0], bode_plot_params)
+        1 if kvect is None else kvect[0], bode_plot_params, tvect=tvect)
 
     # Setup the root-locus plot window
     root_locus(sys, kvect=kvect, xlim=xlim_rlocus,

diff --git a/control/tests/iosys_test.py b/control/tests/iosys_test.py
@@ -56,7 +56,7 @@ def test_linear_iosys(self, tsys):
         # Make sure that the right hand side matches linear system
         for x, u in (([0, 0], 0), ([1, 0], 0), ([0, 1], 0), ([0, 0], 1)):
             np.testing.assert_array_almost_equal(
-                np.reshape(iosys._rhs(0, x, u), (-1, 1)),
+                np.reshape(iosys.dynamics(0, x, u), (-1, 1)),
                 np.dot(linsys.A, np.reshape(x, (-1, 1))) + np.dot(linsys.B, u))
 
         # Make sure that simulations also line up
@@ -687,7 +687,7 @@ def test_find_eqpts(self, tsys):
         assert result.success
         np.testing.assert_array_almost_equal(xeq, [1.64705879, 1.17923874])
         np.testing.assert_array_almost_equal(
-            nlsys._rhs(0, xeq, ueq), np.zeros((2,)))
+            nlsys.dynamics(0, xeq, ueq), np.zeros((2,)))
 
         # Ducted fan dynamics with output = velocity
         nlsys = ios.NonlinearIOSystem(pvtol, lambda t, x, u, params: x[0:2])
@@ -697,15 +697,15 @@ def test_find_eqpts(self, tsys):
             nlsys, [0, 0, 0, 0], [0, 4*9.8], return_result=True)
         assert result.success
         np.testing.assert_array_almost_equal(
-            nlsys._rhs(0, xeq, ueq), np.zeros((4,)))
+            nlsys.dynamics(0, xeq, ueq), np.zeros((4,)))
         np.testing.assert_array_almost_equal(xeq, [0, 0, 0, 0])
 
         # Use a small lateral force to cause motion
         xeq, ueq, result = ios.find_eqpt(
             nlsys, [0, 0, 0, 0], [0.01, 4*9.8], return_result=True)
         assert result.success
         np.testing.assert_array_almost_equal(
-            nlsys._rhs(0, xeq, ueq), np.zeros((4,)), decimal=5)
+            nlsys.dynamics(0, xeq, ueq), np.zeros((4,)), decimal=5)
 
         # Equilibrium point with fixed output
         xeq, ueq, result = ios.find_eqpt(
@@ -715,7 +715,7 @@ def test_find_eqpts(self, tsys):
         np.testing.assert_array_almost_equal(
             nlsys._out(0, xeq, ueq), [0.1, 0.1], decimal=5)
         np.testing.assert_array_almost_equal(
-            nlsys._rhs(0, xeq, ueq), np.zeros((4,)), decimal=5)
+            nlsys.dynamics(0, xeq, ueq), np.zeros((4,)), decimal=5)
 
         # Specify outputs to constrain (replicate previous)
         xeq, ueq, result = ios.find_eqpt(
@@ -725,15 +725,15 @@ def test_find_eqpts(self, tsys):
         np.testing.assert_array_almost_equal(
             nlsys._out(0, xeq, ueq), [0.1, 0.1], decimal=5)
         np.testing.assert_array_almost_equal(
-            nlsys._rhs(0, xeq, ueq), np.zeros((4,)), decimal=5)
+            nlsys.dynamics(0, xeq, ueq), np.zeros((4,)), decimal=5)
 
         # Specify inputs to constrain (replicate previous), w/ no result
         xeq, ueq = ios.find_eqpt(
             nlsys, [0, 0, 0, 0], [0.01, 4*9.8], y0=[0.1, 0.1], iu = [])
         np.testing.assert_array_almost_equal(
             nlsys._out(0, xeq, ueq), [0.1, 0.1], decimal=5)
         np.testing.assert_array_almost_equal(
-            nlsys._rhs(0, xeq, ueq), np.zeros((4,)), decimal=5)
+            nlsys.dynamics(0, xeq, ueq), np.zeros((4,)), decimal=5)
 
         # Now solve the problem with the original PVTOL variables
         # Constrain the output angle and x velocity
@@ -746,7 +746,7 @@ def test_find_eqpts(self, tsys):
         np.testing.assert_array_almost_equal(
             nlsys_full._out(0, xeq, ueq)[[2, 3]], [0.1, 0.1], decimal=5)
         np.testing.assert_array_almost_equal(
-            nlsys_full._rhs(0, xeq, ueq)[-4:], np.zeros((4,)), decimal=5)
+            nlsys_full.dynamics(0, xeq, ueq)[-4:], np.zeros((4,)), decimal=5)
 
         # Fix one input and vary the other
         nlsys_full = ios.NonlinearIOSystem(pvtol_full, None)
@@ -759,7 +759,7 @@ def test_find_eqpts(self, tsys):
         np.testing.assert_array_almost_equal(
             nlsys_full._out(0, xeq, ueq)[[3]], [0.1], decimal=5)
         np.testing.assert_array_almost_equal(
-            nlsys_full._rhs(0, xeq, ueq)[-4:], np.zeros((4,)), decimal=5)
+            nlsys_full.dynamics(0, xeq, ueq)[-4:], np.zeros((4,)), decimal=5)
 
         # PVTOL with output = y velocity
         xeq, ueq, result = ios.find_eqpt(
@@ -771,7 +771,7 @@ def test_find_eqpts(self, tsys):
         np.testing.assert_array_almost_equal(
             nlsys_full._out(0, xeq, ueq)[-3:], [0.1, 0, 0], decimal=5)
         np.testing.assert_array_almost_equal(
-            nlsys_full._rhs(0, xeq, ueq)[-5:], np.zeros((5,)), decimal=5)
+            nlsys_full.dynamics(0, xeq, ueq)[-5:], np.zeros((5,)), decimal=5)
 
         # Unobservable system
         linsys = ct.StateSpace(

diff --git a/control/tests/type_conversion_test.py b/control/tests/type_conversion_test.py
@@ -19,7 +19,7 @@ def sys_dict():
     sdict['frd'] = ct.frd([10+0j, 9 + 1j, 8 + 2j], [1,2,3])
     sdict['lio'] = ct.LinearIOSystem(ct.ss([[-1]], [[5]], [[5]], [[0]]))
     sdict['ios'] = ct.NonlinearIOSystem(
-        sdict['lio']._rhs, sdict['lio']._out, 1, 1, 1)
+        sdict['lio'].dynamics, sdict['lio']._out, 1, 1, 1)
     sdict['arr'] = np.array([[2.0]])
     sdict['flt'] = 3.
     return sdict
@@ -66,7 +66,7 @@ def sys_dict():
     ('add',     'ios', ['xos', 'xos', 'E',   'ios', 'ios', 'xos', 'xos']),
     ('add',     'arr', ['ss',  'tf',  'xrd', 'xio', 'xos', 'arr', 'arr']),
     ('add',     'flt', ['ss',  'tf',  'xrd', 'xio', 'xos', 'arr', 'flt']),
-    
+
     # op        left     ss     tf    frd    lio    ios    arr    flt
     ('sub',     'ss',  ['ss',  'ss',  'xrd', 'ss',  'xos', 'ss',  'ss' ]),
     ('sub',     'tf',  ['tf',  'tf',  'xrd', 'tf',  'xos', 'tf',  'tf' ]),
@@ -75,7 +75,7 @@ def sys_dict():
     ('sub',     'ios', ['xos', 'xio', 'E',   'ios', 'xos'  'xos', 'xos']),
     ('sub',     'arr', ['ss',  'tf',  'xrd', 'xio', 'xos', 'arr', 'arr']),
     ('sub',     'flt', ['ss',  'tf',  'xrd', 'xio', 'xos', 'arr', 'flt']),
-    
+
     # op        left     ss     tf    frd    lio    ios    arr    flt
     ('mul',     'ss',  ['ss',  'ss',  'xrd', 'ss',  'xos', 'ss',  'ss' ]),
     ('mul',     'tf',  ['tf',  'tf',  'xrd', 'tf',  'xos', 'tf',  'tf' ]),
@@ -84,7 +84,7 @@ def sys_dict():
     ('mul',     'ios', ['xos', 'xos', 'E',   'ios', 'ios', 'xos', 'xos']),
     ('mul',     'arr', ['ss',  'tf',  'xrd', 'xio', 'xos', 'arr', 'arr']),
     ('mul',     'flt', ['ss',  'tf',  'frd', 'xio', 'xos', 'arr', 'flt']),
-    
+
     # op        left     ss     tf    frd    lio    ios    arr    flt
     ('truediv', 'ss',  ['xs',  'tf',  'xrd', 'xio', 'xos', 'xs',  'xs' ]),
     ('truediv', 'tf',  ['tf',  'tf',  'xrd', 'tf',  'xos', 'tf',  'tf' ]),
@@ -100,7 +100,7 @@ def sys_dict():
     for rtype, expected in zip(rtype_list, expected_list):
         # Add this to the list of tests to run
         test_matrix.append([opname, ltype, rtype, expected])
-    
+
 @pytest.mark.parametrize("opname, ltype, rtype, expected", test_matrix)
 def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict):
     op = getattr(operator, opname)
@@ -110,7 +110,7 @@ def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict):
     # Get rid of warnings for InputOutputSystem objects by making a copy
     if isinstance(leftsys, ct.InputOutputSystem) and leftsys == rightsys:
         rightsys = leftsys.copy()
-            
+
     # Make sure we get the right result
     if expected == 'E' or expected[0] == 'x':
         # Exception expected
@@ -119,7 +119,7 @@ def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict):
     else:
         # Operation should work and return the given type
         result = op(leftsys, rightsys)
-                
+
         # Print out what we are testing in case something goes wrong
         assert isinstance(result, type_dict[expected])
 
@@ -138,7 +138,7 @@ def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict):
 #
 #   * For IOS/LTI, convert to IOS.  In the case of a linear I/O system (LIO),
 #     this will preserve the linear structure since the LTI system will
-#     be converted to state space.  
+#     be converted to state space.
 #
 #   * When combining state space or transfer with linear I/O systems, the
 #   * output should be of type Linear IO system, since that maintains the
@@ -149,7 +149,7 @@ def test_operator_type_conversion(opname, ltype, rtype, expected, sys_dict):
 
 type_list = ['ss', 'tf', 'tfx', 'frd', 'lio', 'ios', 'arr', 'flt']
 conversion_table = [
-    # L \ R ['ss',  'tf',  'tfx', 'frd', 'lio', 'ios', 'arr', 'flt'] 
+    # L \ R ['ss',  'tf',  'tfx', 'frd', 'lio', 'ios', 'arr', 'flt']
     ('ss',  ['ss',  'ss',  'tf'   'frd', 'lio', 'ios', 'ss',  'ss' ]),
     ('tf',  ['tf',  'tf',  'tf'   'frd', 'lio', 'ios', 'tf',  'tf' ]),
     ('tfx', ['tf',  'tf',  'tf',  'frd', 'E',   'E',   'tf',  'tf' ]),