Repository: python-control/python-control
Pull request #982 (2024-03-30 22:42:39)

Allow passing and saving of params in time responses
This PR allows parameter values to be saved in time responses as well as passing parameter values (via the `params` keyword) for the `initial_response`, ,`step_response`, and `final_response` functions (which call through to `input_output_response` for nonlinear I/O systems).

With this functionality, you can do things like
```
    nlsys = ct.nlsys(
        lambda t, x, u, params: params['a'] * x[0] + u[0],
        states = 1, inputs = 1, outputs = 1, params={'a': 0})
    timevec = np.linspace(0, 1)
    resp = step_response(nlsys, timevec, params={'a': -0.5})
```
The simulation will run with params['a'] = -0.5 and this information will be stored in `resp.params`.



Patch:
diff --git a/control/nlsys.py b/control/nlsys.py
index fd6e207f..c154c081 100644
--- a/control/nlsys.py
+++ b/control/nlsys.py
@@ -1367,6 +1367,8 @@ def input_output_response(
 
         * inputs (array): Input(s) to the system, indexed by input and time.
 
+        * params (dict): Parameters values used for the simulation.
+
         The return value of the system can also be accessed by assigning the
         function to a tuple of length 2 (time, output) or of length 3 (time,
         output, state) if ``return_x`` is ``True``.  If the input/output
@@ -1643,7 +1645,7 @@ def input_output_response(
         raise TypeError("Can't determine system type")
 
     return TimeResponseData(
-        soln.t, y, soln.y, u, issiso=sys.issiso(),
+        soln.t, y, soln.y, u, params=params, issiso=sys.issiso(),
         output_labels=sys.output_labels, input_labels=sys.input_labels,
         state_labels=sys.state_labels, sysname=sys.name,
         title="Input/output response for " + sys.name,
diff --git a/control/statefbk.py b/control/statefbk.py
index 15bba545..3d2ef6d9 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -1119,21 +1119,18 @@ def gram(sys, type):
     if type not in ['c', 'o', 'cf', 'of']:
         raise ValueError("That type is not supported!")
 
-    # Check if system is continuous or discrete
-    if sys.isctime():
-        dico = 'C'
-
-        # TODO: Check system is stable, perhaps a utility in ctrlutil.py
-        # or a method of the StateSpace class?
-        if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
-            raise ValueError("Oops, the system is unstable!")
-
-    else:
-        assert sys.isdtime()
-        dico = 'D'
-
-        if np.any(np.abs(sys.poles()) >= 1.):
-            raise ValueError("Oops, the system is unstable!")
+    # TODO: Check for continuous or discrete, only continuous supported for now
+        # if isCont():
+        #    dico = 'C'
+        # elif isDisc():
+        #    dico = 'D'
+        # else:
+    dico = 'C'
+
+    # TODO: Check system is stable, perhaps a utility in ctrlutil.py
+    # or a method of the StateSpace class?
+    if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
+        raise ValueError("Oops, the system is unstable!")
 
     if type == 'c' or type == 'o':
         # Compute Gramian by the Slycot routine sb03md
diff --git a/control/sysnorm.py b/control/sysnorm.py
index f5e583dc..49873dba 100644
--- a/control/sysnorm.py
+++ b/control/sysnorm.py
@@ -94,13 +94,13 @@ def norm(system, p=2, tol=1e-6, print_warning=True, method=None):
         Type of norm to be computed. ``p=2`` gives the H2 norm, and ``p='inf'`` gives the L-infinity norm.
     tol : float
         Relative tolerance for accuracy of L-infinity norm computation. Ignored
-        unless ``p='inf'``.
+        unless p='inf'.
     print_warning : bool
         Print warning message in case norm value may be uncertain.
     method : str, optional
         Set the method used for computing the result.  Current methods are
-        ``'slycot'`` and ``'scipy'``. If set to ``None`` (default), try ``'slycot'`` first
-        and then ``'scipy'``.
+        'slycot' and 'scipy'. If set to None (default), try 'slycot' first
+        and then 'scipy'.
     
     Returns
     -------
@@ -114,10 +114,10 @@ def norm(system, p=2, tol=1e-6, print_warning=True, method=None):
     Examples
     --------
     >>> Gc = ct.tf([1], [1, 2, 1])
-    >>> round(ct.norm(Gc, 2), 3)
-    0.5
-    >>> round(ct.norm(Gc, 'inf', tol=1e-5, method='scipy'), 3)
-    1.0
+    >>> ct.norm(Gc, 2)
+    0.5000000000000001
+    >>> ct.norm(Gc, 'inf', tol=1e-11, method='scipy')
+    1.000000000007276
     """
            
     if not isinstance(system, (ct.StateSpace, ct.TransferFunction)):
diff --git a/control/timeresp.py b/control/timeresp.py
index b45dbc0e..58207e88 100644
--- a/control/timeresp.py
+++ b/control/timeresp.py
@@ -173,6 +173,9 @@ class TimeResponseData:
     sysname : str, optional
         Name of the system that created the data.
 
+    params : dict, optional
+        If system is a nonlinear I/O system, set parameter values.
+
     plot_inputs : bool, optional
         Whether or not to plot the inputs by default (can be overridden in
         the plot() method)
@@ -227,7 +230,7 @@ class TimeResponseData:
             output_labels=None, state_labels=None, input_labels=None,
             title=None, transpose=False, return_x=False, squeeze=None,
             multi_trace=False, trace_labels=None, trace_types=None,
-            plot_inputs=True, sysname=None
+            plot_inputs=True, sysname=None, params=None
     ):
         """Create an input/output time response object.
 
@@ -315,6 +318,7 @@ class TimeResponseData:
             raise ValueError("Time vector must be 1D array")
         self.title = title
         self.sysname = sysname
+        self.params = params
 
         #
         # Output vector (and number of traces)
@@ -856,7 +860,7 @@ def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False,
 
 
 # Forced response of a linear system
-def forced_response(sys, T=None, U=0., X0=0., transpose=False,
+def forced_response(sys, T=None, U=0., X0=0., transpose=False, params=None,
                     interpolate=False, return_x=None, squeeze=None):
     """Compute the output of a linear system given the input.
 
@@ -889,6 +893,9 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
     X0 : array_like or float, default=0.
         Initial condition.
 
+    params : dict, optional
+        If system is a nonlinear I/O system, set parameter values.
+
     transpose : bool, default=False
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`).
@@ -981,7 +988,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
                 warnings.warn(
                     "interpolation not supported for nonlinear I/O systems")
             return input_output_response(
-                sys, T, U, X0, transpose=transpose,
+                sys, T, U, X0, params=params, transpose=transpose,
                 return_x=return_x, squeeze=squeeze)
         else:
             raise TypeError('Parameter ``sys``: must be a ``StateSpace`` or'
@@ -1169,7 +1176,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False,
         yout = np.transpose(yout)
 
     return TimeResponseData(
-        tout, yout, xout, U, issiso=sys.issiso(),
+        tout, yout, xout, U, params=params, issiso=sys.issiso(),
         output_labels=sys.output_labels, input_labels=sys.input_labels,
         state_labels=sys.state_labels, sysname=sys.name, plot_inputs=True,
         title="Forced response for " + sys.name, trace_types=['forced'],
@@ -1256,7 +1263,7 @@ def _process_time_response(
 
 
 def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
-                  transpose=False, return_x=False, squeeze=None):
+                  transpose=False, return_x=False, squeeze=None, params=None):
     # pylint: disable=W0622
     """Compute the step response for a linear system.
 
@@ -1298,6 +1305,9 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
         Only report the step response for the listed output.  If not
         specified, all outputs are reported.
 
+    params : dict, optional
+        If system is a nonlinear I/O system, set parameter values.
+
     T_num : int, optional
         Number of time steps to use in simulation if T is not provided as an
         array (autocomputed if not given); ignored if sys is discrete-time.
@@ -1403,7 +1413,7 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
         U = np.zeros((sys.ninputs, T.size))
         U[i, :] = np.ones_like(T)
 
-        response = forced_response(sys, T, U, X0, squeeze=True)
+        response = forced_response(sys, T, U, X0, squeeze=True, params=params)
         inpidx = i if input is None else 0
         yout[:, inpidx, :] = response.y if output is None \
             else response.y[output]
@@ -1424,11 +1434,11 @@ def step_response(sys, T=None, X0=0, input=None, output=None, T_num=None,
         output_labels=output_labels, input_labels=input_labels,
         state_labels=sys.state_labels, title="Step response for " + sys.name,
         transpose=transpose, return_x=return_x, squeeze=squeeze,
-        sysname=sys.name, trace_labels=trace_labels,
+        sysname=sys.name, params=params, trace_labels=trace_labels,
         trace_types=trace_types, plot_inputs=False)
 
 
-def step_info(sysdata, T=None, T_num=None, yfinal=None,
+def step_info(sysdata, T=None, T_num=None, yfinal=None, params=None,
               SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)):
     """
     Step response characteristics (Rise time, Settling Time, Peak and others).
@@ -1451,6 +1461,8 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None,
         systems to simulate and the last value of the the response data is
         used for a given time series of response data. Scalar for SISO,
         (noutputs, ninputs) array_like for MIMO systems.
+    params : dict, optional
+        If system is a nonlinear I/O system, set parameter values.
     SettlingTimeThreshold : float, optional
         Defines the error to compute settling time (default = 0.02)
     RiseTimeLimits : tuple (lower_threshold, upper_theshold)
@@ -1536,7 +1548,7 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None,
     from .nlsys import NonlinearIOSystem
 
     if isinstance(sysdata, (StateSpace, TransferFunction, NonlinearIOSystem)):
-        T, Yout = step_response(sysdata, T, squeeze=False)
+        T, Yout = step_response(sysdata, T, squeeze=False, params=params)
         if yfinal:
             InfValues = np.atleast_2d(yfinal)
         else:
@@ -1651,7 +1663,7 @@ def step_info(sysdata, T=None, T_num=None, yfinal=None,
     return ret[0][0] if retsiso else ret
 
 
-def initial_response(sys, T=None, X0=0, output=None, T_num=None,
+def initial_response(sys, T=None, X0=0, output=None, T_num=None, params=None,
                      transpose=False, return_x=False, squeeze=None):
     # pylint: disable=W0622
     """Compute the initial condition response for a linear system.
@@ -1684,6 +1696,9 @@ def initial_response(sys, T=None, X0=0, output=None, T_num=None,
         Number of time steps to use in simulation if T is not provided as an
         array (autocomputed if not given); ignored if sys is discrete-time.
 
+    params : dict, optional
+        If system is a nonlinear I/O system, set parameter values.
+
     transpose : bool, optional
         If True, transpose all input and output arrays (for backward
         compatibility with MATLAB and :func:`scipy.signal.lsim`).  Default
@@ -1748,7 +1763,7 @@ def initial_response(sys, T=None, X0=0, output=None, T_num=None,
         raise ValueError("invalid value of T for this type of system")
 
     # Compute the forced response
-    response = forced_response(sys, T, 0, X0)
+    response = forced_response(sys, T, 0, X0, params=params)
 
     # Figure out if the system is SISO or not
     issiso = sys.issiso() or output is not None
@@ -1760,7 +1775,7 @@ def initial_response(sys, T=None, X0=0, output=None, T_num=None,
 
     # Store the response without an input
     return TimeResponseData(
-        response.t, yout, response.x, None, issiso=issiso,
+        response.t, yout, response.x, None, params=params, issiso=issiso,
         output_labels=output_labels, input_labels=None,
         state_labels=sys.state_labels, sysname=sys.name,
         title="Initial response for " + sys.name, trace_types=['initial'],
@@ -1863,6 +1878,10 @@ def impulse_response(sys, T=None, input=None, output=None, T_num=None,
     from .statesp import _convert_to_statespace
     from .lti import LTI
 
+    # Make sure we have an LTI system
+    if not isinstance(sys, LTI):
+        raise ValueError("system must be LTI system for impulse response")
+
     # Create the time and input vectors
     if T is None or np.asarray(T).size == 1:
         T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)
diff --git a/doc/control.rst b/doc/control.rst
index 1b1b7406..96714bf7 100644
--- a/doc/control.rst
+++ b/doc/control.rst
@@ -179,7 +179,6 @@ Utility functions and conversions
     issys
     mag2db
     modal_form
-    norm
     observable_form
     pade
     reachable_form


Test patch:
diff --git a/control/tests/nlsys_test.py b/control/tests/nlsys_test.py
index 3f8114c3..1c2976c5 100644
--- a/control/tests/nlsys_test.py
+++ b/control/tests/nlsys_test.py
@@ -75,3 +75,20 @@ def test_lti_nlsys_response(nin, nout, input, output):
     resp_nl = ct.forced_response(sys_nl, timepts, U, X0=X0)
     np.testing.assert_equal(resp_ss.time, resp_nl.time)
     np.testing.assert_allclose(resp_ss.states, resp_nl.states, atol=0.05)
+
+
+# Test to make sure that impulse responses are not allowed
+def test_nlsys_impulse():
+    sys_ss = ct.rss(4, 1, 1, strictly_proper=True)
+    sys_nl = ct.nlsys(
+        lambda t, x, u, params: sys_ss.A @ x + sys_ss.B @ u,
+        lambda t, x, u, params: sys_ss.C @ x + sys_ss.D @ u,
+        inputs=1, outputs=1, states=4)
+
+    # Figure out the time to use from the linear impulse response
+    resp_ss = ct.impulse_response(sys_ss)
+    timepts = np.linspace(0, resp_ss.time[-1]/10, 100)
+
+    # Impulse_response (not implemented)
+    with pytest.raises(ValueError, match="system must be LTI"):
+        resp_nl = ct.impulse_response(sys_nl, timepts)
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
index 4a0472de..cb677b40 100644
--- a/control/tests/statefbk_test.py
+++ b/control/tests/statefbk_test.py
@@ -98,28 +98,6 @@ class TestStatefbk:
         Wctrue = np.array([[18.5, 24.5], [24.5, 32.5]])
         Wc = gram(sys, 'c')
         np.testing.assert_array_almost_equal(Wc, Wctrue)
-        sysd = ct.c2d(sys, 0.2)
-        Wctrue = np.array([[3.666767, 4.853625],
-                           [4.853625, 6.435233]])
-        Wc = gram(sysd, 'c')
-        np.testing.assert_array_almost_equal(Wc, Wctrue)
-
-    @slycotonly
-    def testGramWc2(self):
-        A = np.array([[1., -2.], [3., -4.]])
-        B = np.array([[5.], [7.]])
-        C = np.array([[6., 8.]])
-        D = np.array([[9.]])
-        sys = ss(A,B,C,D)
-        Wctrue = np.array([[ 7.166667,  9.833333],
-                           [ 9.833333,  13.5]])
-        Wc = gram(sys, 'c')
-        np.testing.assert_array_almost_equal(Wc, Wctrue)
-        sysd = ct.c2d(sys, 0.2)
-        Wctrue = np.array([[1.418978, 1.946180],
-                           [1.946180, 2.670758]])
-        Wc = gram(sysd, 'c')
-        np.testing.assert_array_almost_equal(Wc, Wctrue)
 
     @slycotonly
     def testGramRc(self):
@@ -128,15 +106,9 @@ class TestStatefbk:
         C = np.array([[4., 5.], [6., 7.]])
         D = np.array([[13., 14.], [15., 16.]])
         sys = ss(A, B, C, D)
-        Rctrue = np.array([[4.30116263, 5.6961343], 
-                           [0., 0.23249528]])
+        Rctrue = np.array([[4.30116263, 5.6961343], [0., 0.23249528]])
         Rc = gram(sys, 'cf')
         np.testing.assert_array_almost_equal(Rc, Rctrue)
-        sysd = ct.c2d(sys, 0.2)
-        Rctrue = np.array([[1.91488054, 2.53468814],
-                           [0.        , 0.10290372]])
-        Rc = gram(sysd, 'cf')
-        np.testing.assert_array_almost_equal(Rc, Rctrue)        
 
     @slycotonly
     def testGramWo(self):
@@ -148,11 +120,6 @@ class TestStatefbk:
         Wotrue = np.array([[257.5, -94.5], [-94.5, 56.5]])
         Wo = gram(sys, 'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
-        sysd = ct.c2d(sys, 0.2)
-        Wotrue = np.array([[ 1305.369179, -440.046414], 
-                           [ -440.046414,  333.034844]])
-        Wo = gram(sysd, 'o')
-        np.testing.assert_array_almost_equal(Wo, Wotrue)
 
     @slycotonly
     def testGramWo2(self):
@@ -164,11 +131,6 @@ class TestStatefbk:
         Wotrue = np.array([[198., -72.], [-72., 44.]])
         Wo = gram(sys, 'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
-        sysd = ct.c2d(sys, 0.2)
-        Wotrue = np.array([[ 1001.835511, -335.337663],
-                           [ -335.337663,  263.355793]])
-        Wo = gram(sysd, 'o')
-        np.testing.assert_array_almost_equal(Wo, Wotrue)
 
     @slycotonly
     def testGramRo(self):
@@ -180,34 +142,15 @@ class TestStatefbk:
         Rotrue = np.array([[16.04680654, -5.8890222], [0., 4.67112593]])
         Ro = gram(sys, 'of')
         np.testing.assert_array_almost_equal(Ro, Rotrue)
-        sysd = ct.c2d(sys, 0.2)
-        Rotrue = np.array([[ 36.12989315, -12.17956588],
-                           [  0.        ,  13.59018097]])
-        Ro = gram(sysd, 'of')
-        np.testing.assert_array_almost_equal(Ro, Rotrue) 
 
     def testGramsys(self):
-        sys = tf([1.], [1., 1., 1.])
-        with pytest.raises(ValueError) as excinfo:
-            gram(sys, 'o')
-        assert "must be StateSpace" in str(excinfo.value)
-        with pytest.raises(ValueError) as excinfo:
-            gram(sys, 'c')
-        assert "must be StateSpace" in str(excinfo.value)
-        sys = tf([1], [1, -1], 0.5)
-        with pytest.raises(ValueError) as excinfo:
-            gram(sys, 'o')
-        assert "must be StateSpace" in str(excinfo.value)
-        with pytest.raises(ValueError) as excinfo:
-            gram(sys, 'c')
-        assert "must be StateSpace" in str(excinfo.value)
-        sys = ct.ss(sys)  # this system is unstable
-        with pytest.raises(ValueError) as excinfo:
+        num =[1.]
+        den = [1., 1., 1.]
+        sys = tf(num,den)
+        with pytest.raises(ValueError):
             gram(sys, 'o')
-        assert "is unstable" in str(excinfo.value)
-        with pytest.raises(ValueError) as excinfo:
+        with pytest.raises(ValueError):
             gram(sys, 'c')
-        assert "is unstable" in str(excinfo.value)
 
     def testAcker(self, fixedseed):
         for states in range(1, self.maxStates):
diff --git a/control/tests/trdata_test.py b/control/tests/trdata_test.py
index 15ce7512..7d0c20e7 100644
--- a/control/tests/trdata_test.py
+++ b/control/tests/trdata_test.py
@@ -236,7 +236,8 @@ def test_trdata_labels():
     np.testing.assert_equal(
         response.input_labels, ["u[%d]" % i for i in range(sys.ninputs)])
 
-    # Make sure the selected input and output are both correctly transferred to the response
+    # Make sure the selected input and output are both correctly
+    # transferred to the response
     for nu in range(sys.ninputs):
         for ny in range(sys.noutputs):
             step_response = ct.step_response(sys, T, input=nu, output=ny)
@@ -339,6 +340,37 @@ def test_trdata_multitrace():
             np.zeros(5), np.ones(5), np.zeros((1, 5)), np.zeros(6))
 
 
+@pytest.mark.parametrize("func, args", [
+    (ct.step_response, ()),
+    (ct.initial_response, (1, )),
+    (ct.forced_response, (0, 1)),
+    (ct.input_output_response, (0, 1)),
+])
+@pytest.mark.parametrize("dt", [0, 1])
+def test_trdata_params(func, args, dt):
+    # Create a nonlinear system with parameters, neutrally stable
+    nlsys = ct.nlsys(
+        lambda t, x, u, params: params['a'] * x[0] + u[0],
+        states = 1, inputs = 1, outputs = 1, params={'a': 0}, dt=dt)
+    lnsys = ct.ss([[-0.5]], [[1]], [[1]], 0, dt=dt)
+
+    # Compute the response, setting parameters to make things stable
+    timevec = np.linspace(0, 1) if dt == 0 else np.arange(0, 10, 1)
+    nlresp = func(nlsys, timevec, *args, params={'a': -0.5})
+    lnresp = func(lnsys, timevec, *args)
+
+    # Make sure the modified system was stable
+    np.testing.assert_allclose(
+        nlresp.states, lnresp.states, rtol=1e-3, atol=1e-5)
+    assert lnresp.params == None
+    assert nlresp.params['a'] == -0.5
+
+    # Make sure the match was not accidental
+    bdresp = func(nlsys, timevec, *args)
+    assert not np.allclose(
+        bdresp.states, nlresp.states, rtol=1e-3, atol=1e-5)
+
+
 def test_trdata_exceptions():
     # Incorrect dimension for time vector
     with pytest.raises(ValueError, match="Time vector must be 1D"):