allow saving and passing of params (for nl systems) in time responses

murrayrm · murrayrm · commit 85568a95808f · 2024-03-30T15:38:35.000-07:00
diff --git 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 ivp_rhs(t, x):
         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/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/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"):
diff --git 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 @@ def __init__(
             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 @@ def __init__(
             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 shape_matches(s_legal, s_actual):
 
 
 # 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)
