diff --git a/control/nlsys.py b/control/nlsys.py index fd6e207fc..c154c0818 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 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 index 3f8114c3b..1c2976c56 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/trdata_test.py b/control/tests/trdata_test.py index 15ce75121..7d0c20e7a 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"): diff --git a/control/timeresp.py b/control/timeresp.py index b45dbc0ea..58207e88e 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 @@ 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)