"""Struct of test system"""

def __init__(self, sys=None, call_kwargs=None):

self.sys = sys

self.kwargs = call_kwargs if call_kwargs else {}

def __repr__(self):

"""Show system when debugging"""

@pytest.fixture

def tsystem(self, request):

"""Define some test systems"""

"""continuous"""

A = np.array([[1., -2.], [3., -4.]])

B = np.array([[5.], [7.]])

C = np.array([[6., 8.]])

D = np.array([[9.]])

siso_ss1 = TSys(StateSpace(A, B, C, D, 0))

siso_ss1.t = np.linspace(0, 1, 10)

siso_ss1.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234,

39.1165, 42.3227, 44.9694, 47.1599,

48.9776])

siso_ss1.X0 = np.array([[.5], [1.]])

siso_ss1.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118,

1.9092, 1.1508, 0.5833, 0.1645, -0.1391])

ss1 = siso_ss1.sys

"""D=0, continuous"""

siso_ss2 = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, 0))

siso_ss2.t = siso_ss1.t

siso_ss2.ystep = siso_ss1.ystep - 9

siso_ss2.initial = siso_ss1.yinitial - 9

siso_ss2.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766,

31.7344, 26.1668, 21.6292, 17.9245,

14.8945])

"""System with unspecified timebase"""

siso_ss2_dtnone = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, None))

siso_ss2_dtnone.t = np.arange(0, 10, 1.)

siso_ss2_dtnone.ystep = np.array([0., 86., -72., 230., -360., 806.,

-1512., 3110., -6120., 12326.])

siso_tf1 = TSys(TransferFunction([1], [1, 2, 1], 0))

siso_tf2 = copy(siso_ss1)

siso_tf2.sys = ss2tf(siso_ss1.sys)

"""MIMO system, contains `siso_ss1` twice"""

mimo_ss1 = copy(siso_ss1)

A = np.zeros((4, 4))

A[:2, :2] = siso_ss1.sys.A

A[2:, 2:] = siso_ss1.sys.A

B = np.zeros((4, 2))

B[:2, :1] = siso_ss1.sys.B

B[2:, 1:] = siso_ss1.sys.B

C = np.zeros((2, 4))

C[:1, :2] = siso_ss1.sys.C

C[1:, 2:] = siso_ss1.sys.C

D = np.zeros((2, 2))

D[:1, :1] = siso_ss1.sys.D

D[1:, 1:] = siso_ss1.sys.D

mimo_ss1.sys = StateSpace(A, B, C, D)

"""MIMO system, contains `siso_ss2` twice"""

mimo_ss2 = copy(siso_ss2)

A = np.zeros((4, 4))

A[:2, :2] = siso_ss2.sys.A

A[2:, 2:] = siso_ss2.sys.A

B = np.zeros((4, 2))

B[:2, :1] = siso_ss2.sys.B

B[2:, 1:] = siso_ss2.sys.B

C = np.zeros((2, 4))

C[:1, :2] = siso_ss2.sys.C

C[1:, 2:] = siso_ss2.sys.C

D = np.zeros((2, 2))

mimo_ss2.sys = StateSpace(A, B, C, D, 0)

"""discrete"""

siso_dtf0 = TSys(TransferFunction([1.], [1., 0.], 1.))

siso_dtf0.t = np.arange(4)

siso_dtf0.yimpulse = [0., 1., 0., 0.]

siso_dtf1 = TSys(TransferFunction([1], [1, 1, 0.25], True))

siso_dtf1.t = np.arange(0, 5, 1)

siso_dtf1.ystep = np.array([0. , 0. , 1. , 0. , 0.75])

siso_dtf2 = TSys(TransferFunction([1], [1, 1, 0.25], 0.2))

siso_dtf2.t = np.arange(0, 5, 0.2)

siso_dtf2.ystep = np.array(

[0. , 0. , 1. , 0. , 0.75 , 0.25 ,

0.5625, 0.375 , 0.4844, 0.4219, 0.457 , 0.4375,

0.4482, 0.4424, 0.4456, 0.4438, 0.4448, 0.4443,

0.4445, 0.4444, 0.4445, 0.4444, 0.4445, 0.4444,

0.4444])

"""Time step which leads to rounding errors for time vector length"""

num = [-0.10966442, 0.12431949]

den = [1., -1.86789511, 0.88255018]

dt = 0.12493963338370018

siso_dtf3 = TSys(TransferFunction(num, den, dt))

siso_dtf3.t = np.linspace(0, 9*dt, 10)

siso_dtf3.ystep = np.array(

[ 0. , -0.1097, -0.1902, -0.2438, -0.2729,

-0.2799, -0.2674, -0.2377, -0.1934, -0.1368])

"""dtf1 converted statically, because Slycot and Scipy produce

siso_dss1 = copy(siso_dtf1)

siso_dss1.sys = StateSpace([[-1., -0.25],

[ 1., 0.]],

[[1.],

[0.]],

[[0., 1.]],

[[0.]],

True)

siso_dss1.X0 = [0.5, 1.]

siso_dss1.yinitial = np.array([1., 0.5, -0.75, 0.625, -0.4375])

siso_dss2 = copy(siso_dtf2)

siso_dss2.sys = tf2ss(siso_dtf2.sys)

mimo_dss1 = TSys(StateSpace(ss1.A, ss1.B, ss1.C, ss1.D, True))

mimo_dss1.t = np.arange(0, 5, 0.2)

mimo_dss2 = copy(mimo_ss1)

mimo_dss2.sys = c2d(mimo_ss1.sys, mimo_ss1.t[1]-mimo_ss1.t[0])

mimo_tf2 = copy(mimo_ss2)

tf_ = ss2tf(siso_ss2.sys)

mimo_tf2.sys = TransferFunction(

[[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]],

[[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]],

0)

mimo_dtf1 = copy(siso_dtf1)

tf_ = siso_dtf1.sys

mimo_dtf1.sys = TransferFunction(

[[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]],

[[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]],

True)

# for pole cancellation tests

pole_cancellation = TSys(TransferFunction(

[1.067e+05, 5.791e+04],

[10.67, 1.067e+05, 5.791e+04]))

no_pole_cancellation = TSys(TransferFunction(

[1.881e+06],

[188.1, 1.881e+06]))

# System Type 1 - Step response not stationary: G(s)=1/s(s+1)

siso_tf_type1 = TSys(TransferFunction(1, [1, 1, 0]))

siso_tf_type1.step_info = {

'RiseTime': np.nan,

'SettlingTime': np.nan,

'SettlingMin': np.nan,

'SettlingMax': np.nan,

'Overshoot': np.nan,

'Undershoot': np.nan,

'Peak': np.inf,

'PeakTime': np.inf,

'SteadyStateValue': np.nan}

# SISO under shoot response and positive final value

# G(s)=(-s+1)/(s²+s+1)

siso_tf_kpos = TSys(TransferFunction([-1, 1], [1, 1, 1]))

siso_tf_kpos.step_info = {

'RiseTime': 1.242,

'SettlingTime': 9.110,

'SettlingMin': 0.90,

'SettlingMax': 1.208,

'Overshoot': 20.840,

'Undershoot': 28.0,

'Peak': 1.208,

'PeakTime': 4.282,

'SteadyStateValue': 1.0}

# SISO under shoot response and negative final value

# k=-1 G(s)=-(-s+1)/(s²+s+1)

siso_tf_kneg = TSys(TransferFunction([1, -1], [1, 1, 1]))

siso_tf_kneg.step_info = {

'RiseTime': 1.242,

'SettlingTime': 9.110,

'SettlingMin': -1.208,

'SettlingMax': -0.90,

'Overshoot': 20.840,

'Undershoot': 28.0,

'Peak': 1.208,

'PeakTime': 4.282,

'SteadyStateValue': -1.0}

siso_tf_asymptotic_from_neg1 = TSys(TransferFunction([-1, 1], [1, 1]))

siso_tf_asymptotic_from_neg1.step_info = {

'RiseTime': 2.197,

'SettlingTime': 4.605,

'SettlingMin': 0.9,

'SettlingMax': 1.0,

'Overshoot': 0,

'Undershoot': 100.0,

'Peak': 1.0,

'PeakTime': 0.0,

'SteadyStateValue': 1.0}

siso_tf_asymptotic_from_neg1.kwargs = {

'step_info': {'T': np.arange(0, 5, 1e-3)}}

# example from matlab online help

# https://www.mathworks.com/help/control/ref/stepinfo.html

siso_tf_step_matlab = TSys(TransferFunction([1, 5, 5],

[1, 1.65, 5, 6.5, 2]))

siso_tf_step_matlab.step_info = {

'RiseTime': 3.8456,

'SettlingTime': 27.9762,

'SettlingMin': 2.0689,

'SettlingMax': 2.6873,

'Overshoot': 7.4915,

'Undershoot': 0,

'Peak': 2.6873,

'PeakTime': 8.0530,

'SteadyStateValue': 2.5}

A = [[0.68, -0.34],

[0.34, 0.68]]

B = [[0.18, -0.05],

[0.04, 0.11]]

C = [[0, -1.53],

[-1.12, -1.10]]

D = [[0, 0],

[0.06, -0.37]]

mimo_ss_step_matlab = TSys(StateSpace(A, B, C, D, 0.2))

mimo_ss_step_matlab.kwargs['step_info'] = {'T': 4.6}

mimo_ss_step_matlab.step_info = [[

{'RiseTime': 0.6000,

'SettlingTime': 3.0000,

'SettlingMin': -0.5999,

'SettlingMax': -0.4689,

'Overshoot': 15.5072,

'Undershoot': 0.,

'Peak': 0.5999,

'PeakTime': 1.4000,

'SteadyStateValue': -0.5193},

{'RiseTime': 0.,

'SettlingTime': 3.6000,

'SettlingMin': -0.2797,

'SettlingMax': -0.1043,

'Overshoot': 118.9918,

'Undershoot': 0,

'Peak': 0.2797,

'PeakTime': .6000,

'SteadyStateValue': -0.1277}],

[{'RiseTime': 0.4000,

'SettlingTime': 2.8000,

'SettlingMin': -0.6724,

'SettlingMax': -0.5188,

'Overshoot': 24.6476,

'Undershoot': 11.1224,

'Peak': 0.6724,

'PeakTime': 1,

'SteadyStateValue': -0.5394},

{'RiseTime': 0.0000, # (*)

'SettlingTime': 3.4000,

'SettlingMin': -0.4350, # (*)

'SettlingMax': -0.1485,

'Overshoot': 132.0170,

'Undershoot': 0.,

'Peak': 0.4350,

'PeakTime': .2,

'SteadyStateValue': -0.1875}]]

# (*): MATLAB gives 0.4 for RiseTime and -0.1034 for

# SettlingMin, but it is unclear what 10% and 90% of

# the steady state response mean, when the step for

# this channel does not start a 0.

siso_ss_step_matlab = copy(mimo_ss_step_matlab)

siso_ss_step_matlab.sys = siso_ss_step_matlab.sys[1, 0]

siso_ss_step_matlab.step_info = siso_ss_step_matlab.step_info[1][0]

Ta = [[siso_tf_kpos, siso_tf_kneg, siso_tf_step_matlab],

[siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]]

mimo_tf_step_info = TSys(TransferFunction(

[[Ti.sys.num[0][0] for Ti in Tr] for Tr in Ta],

[[Ti.sys.den[0][0] for Ti in Tr] for Tr in Ta]))

mimo_tf_step_info.step_info = [[Ti.step_info for Ti in Tr]

for Tr in Ta]

# enforce enough sample points for all channels (they have different

# characteristics)

mimo_tf_step_info.kwargs['step_info'] = {'T_num': 2000}

systems = locals()

if isinstance(request.param, str):

return systems[request.param]

else:

return [systems[sys] for sys in request.param]

@pytest.mark.parametrize(

@pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)

def test_step_response_siso(self, tsystem, kwargs):

"""Test SISO system step response"""

sys = tsystem.sys

t = tsystem.t

yref = tsystem.ystep

# SISO call

out = step_response(sys, T=t, **kwargs)

tout, yout = out[:2]

assert len(out) == 3 if ('return_x', True) in kwargs.items() else 2

np.testing.assert_array_almost_equal(tout, t)

np.testing.assert_array_almost_equal(yout, yref, decimal=4)

@pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)

def test_step_response_mimo(self, tsystem):

"""Test MIMO system, which contains `siso_ss1` twice."""

sys = tsystem.sys

t = tsystem.t

yref = tsystem.ystep

_t, y_00 = step_response(sys, T=t, input=0, output=0)

np.testing.assert_array_almost_equal(y_00, yref, decimal=4)

_t, y_11 = step_response(sys, T=t, input=1, output=1)

np.testing.assert_array_almost_equal(y_11, yref, decimal=4)

_t, y_01 = step_response(

sys, T=t, input_indices=[0], output_indices=[1])

np.testing.assert_array_almost_equal(y_01, 0 * yref, decimal=4)

# Make sure we get the same result using MIMO step response

response = step_response(sys, T=t)

np.testing.assert_allclose(response.y[0, 0, :], y_00)

np.testing.assert_allclose(response.y[1, 1, :], y_11)

np.testing.assert_allclose(response.u[0, 0, :], 1)

np.testing.assert_allclose(response.u[1, 0, :], 0)

np.testing.assert_allclose(response.u[0, 1, :], 0)

np.testing.assert_allclose(response.u[1, 1, :], 1)

# Index lists not yet implemented

with pytest.raises(NotImplementedError, match="list of .* indices"):

step_response(

sys, timepts=t, input_indices=[0, 1], output_indices=[1])

@pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)

def test_step_response_return(self, tsystem):

"""Verify continuous and discrete time use same return conventions."""

sysc = tsystem.sys

sysd = c2d(sysc, 1) # discrete-time system

Tvec = np.linspace(0, 10, 11) # make sure to use integer times 0..10

Tc, youtc = step_response(sysc, Tvec, input=0)

Td, youtd = step_response(sysd, Tvec, input=0)

np.testing.assert_array_equal(Tc.shape, Td.shape)

np.testing.assert_array_equal(youtc.shape, youtd.shape)

@pytest.mark.parametrize("dt", [0, 1], ids=["continuous", "discrete"])

def test_step_nostates(self, dt):

"""Constant system, continuous and discrete time.

sys = TransferFunction([1], [1], dt)

t, y = step_response(sys)

np.testing.assert_allclose(y, np.ones(len(t)))

def assert_step_info_match(self, sys, info, info_ref):

"""Assert reasonable step_info accuracy."""

if sys.isdtime(strict=True):

dt = sys.dt

else:

_, dt = _ideal_tfinal_and_dt(sys, is_step=True)

for k in ['RiseTime', 'SettlingTime', 'PeakTime']:

np.testing.assert_allclose(info[k], info_ref[k], atol=dt,

err_msg=f"{k} does not match")

for k in ['Overshoot', 'Undershoot', 'Peak', 'SteadyStateValue']:

np.testing.assert_allclose(info[k], info_ref[k], rtol=5e-3,

err_msg=f"{k} does not match")

# steep gradient right after RiseTime

absrefinf = np.abs(info_ref['SteadyStateValue'])

if info_ref['RiseTime'] > 0:

y_next_sample_max = 0.8*absrefinf/info_ref['RiseTime']*dt

else:

y_next_sample_max = 0

for k in ['SettlingMin', 'SettlingMax']:

if (np.abs(info_ref[k]) - 0.9 * absrefinf) > y_next_sample_max:

# local min/max peak well after signal has risen

np.testing.assert_allclose(info[k], info_ref[k], rtol=1e-3)

@pytest.mark.parametrize(

@pytest.mark.parametrize(

@pytest.mark.parametrize(

def test_step_info(self, tsystem, systype, time_2d, yfinal):

"""Test step info for SISO systems."""

step_info_kwargs = tsystem.kwargs.get('step_info', {})

if systype == "time response":

# simulate long enough for steady state value

tfinal = 3 * tsystem.step_info['SettlingTime']

if np.isnan(tfinal):

pytest.skip("test system does not settle")

t, y = step_response(tsystem.sys, T=tfinal, T_num=5000)

sysdata = y

step_info_kwargs['T'] = t[np.newaxis, :] if time_2d else t

else:

sysdata = tsystem.sys

if yfinal:

step_info_kwargs['yfinal'] = tsystem.step_info['SteadyStateValue']

info = step_info(sysdata, **step_info_kwargs)

self.assert_step_info_match(tsystem.sys, info, tsystem.step_info)

@pytest.mark.parametrize(

@pytest.mark.parametrize(

@pytest.mark.parametrize(

def test_step_info_mimo(self, tsystem, systype, yfinal):

"""Test step info for MIMO systems."""

step_info_kwargs = tsystem.kwargs.get('step_info', {})

if systype == "time response":

tfinal = 3 * max([S['SettlingTime']

for Srow in tsystem.step_info for S in Srow])

t, y = step_response(tsystem.sys, T=tfinal, T_num=5000)

sysdata = y

step_info_kwargs['T'] = t

else:

sysdata = tsystem.sys

if yfinal:

step_info_kwargs['yfinal'] = [[S['SteadyStateValue']

for S in Srow]

for Srow in tsystem.step_info]

info_dict = step_info(sysdata, **step_info_kwargs)

for i, row in enumerate(info_dict):

for j, info in enumerate(row):

self.assert_step_info_match(tsystem.sys,

info, tsystem.step_info[i][j])

def test_step_info_invalid(self):

"""Call step_info with invalid parameters."""

with pytest.raises(ValueError, match="time series data convention"):

step_info(["not numeric data"])

with pytest.raises(ValueError, match="time series data convention"):

step_info(np.ones((10, 15))) # invalid shape

with pytest.raises(ValueError, match="matching time vector"):

step_info(np.ones(15), T=np.linspace(0, 1, 20)) # time too long

with pytest.raises(ValueError, match="matching time vector"):

step_info(np.ones((2, 2, 15))) # no time vector

@pytest.mark.parametrize("tsystem",

def test_step_pole_cancellation(self, tsystem):

# confirm that pole-zero cancellation doesn't perturb results

# https://github.com/python-control/python-control/issues/440

step_info_no_cancellation = step_info(tsystem[0].sys)

step_info_cancellation = step_info(tsystem[1].sys)

self.assert_step_info_match(tsystem[0].sys,

step_info_no_cancellation,

step_info_cancellation)

@pytest.mark.parametrize(

def test_impulse_response_siso(self, tsystem, kwargs):

"""Test impulse response of SISO systems."""

sys = tsystem.sys

t = tsystem.t

yref = tsystem.yimpulse

out = impulse_response(sys, T=t, **kwargs)

tout, yout = out[:2]

assert len(out) == 3 if ('return_x', True) in kwargs.items() else 2

np.testing.assert_array_almost_equal(tout, t)

np.testing.assert_array_almost_equal(yout, yref, decimal=4)

@pytest.mark.parametrize("tsystem", ["mimo_ss2"], indirect=True)

def test_impulse_response_mimo(self, tsystem):

""""Test impulse response of MIMO systems."""

sys = tsystem.sys

t = tsystem.t

yref = tsystem.yimpulse

_t, y_00 = impulse_response(sys, T=t, input=0, output=0)

np.testing.assert_array_almost_equal(y_00, yref, decimal=4)

_t, y_11 = impulse_response(sys, T=t, input=1, output=1)

np.testing.assert_array_almost_equal(y_11, yref, decimal=4)

_t, y_01 = impulse_response(

sys, T=t, input_indices=[0], output_indices=[1])

np.testing.assert_array_almost_equal(y_01, 0 * yref, decimal=4)

yref_notrim = np.zeros((2, len(t)))

yref_notrim[:1, :] = yref

_t, yy = impulse_response(sys, T=t, input=0)

np.testing.assert_array_almost_equal(yy[:,0,:], yref_notrim, decimal=4)

# Index lists not yet implemented

with pytest.raises(NotImplementedError, match="list of .* indices"):

impulse_response(

sys, timepts=t, input_indices=[0, 1], output_indices=[1])

@pytest.mark.parametrize("tsystem", ["siso_tf1"], indirect=True)

def test_discrete_time_impulse(self, tsystem):

# discrete-time impulse sampled version should match cont time

dt = 0.1

t = np.arange(0, 3, dt)

sys = tsystem.sys

sysdt = sys.sample(dt, 'impulse')

np.testing.assert_array_almost_equal(impulse_response(sys, t)[1],

impulse_response(sysdt, t)[1])

def test_discrete_time_impulse_input(self):

# discrete-time impulse input, Only one active input for each trace

A = [[.5, 0.25],[.0, .5]]

B = [[1., 0,],[0., 1.]]

C = [[1., 0.],[0., 1.]]

D = [[0., 0.],[0., 0.]]

dt = True

sysd = ct.ss(A,B,C,D, dt=dt)

response = ct.impulse_response(sysd,T=dt*3)

Uexpected = np.zeros((2,2,4), dtype=float).astype(object)

Uexpected[0,0,0] = 1./dt

Uexpected[1,1,0] = 1./dt

np.testing.assert_array_equal(response.inputs,Uexpected)

dt = 0.5

sysd = ct.ss(A,B,C,D, dt=dt)

response = ct.impulse_response(sysd,T=dt*3)

Uexpected = np.zeros((2,2,4), dtype=float).astype(object)

Uexpected[0,0,0] = 1./dt

Uexpected[1,1,0] = 1./dt

np.testing.assert_array_equal(response.inputs,Uexpected)

@pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)

def test_impulse_response_warnD(self, tsystem):

"""Test warning about direct feedthrough"""

with pytest.warns(UserWarning, match="System has direct feedthrough"):

_ = impulse_response(tsystem.sys, tsystem.t)

@pytest.mark.parametrize(

@pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)

def test_initial_response(self, tsystem, kwargs):

"""Test initial response of SISO system"""

sys = tsystem.sys

t = tsystem.t

x0 = kwargs.get('X0', 0)

yref = tsystem.yinitial if np.any(x0) else np.zeros_like(t)

out = initial_response(sys, T=t, **kwargs)

tout, yout = out[:2]

assert len(out) == 3 if ('return_x', True) in kwargs.items() else 2

np.testing.assert_array_almost_equal(tout, t)

np.testing.assert_array_almost_equal(yout, yref, decimal=4)

@pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True)

def test_initial_response_mimo(self, tsystem):

"""Test initial response of MIMO system"""

sys = tsystem.sys

t = tsystem.t

x0 = np.array([[.5], [1.], [.5], [1.]])

yref = tsystem.yinitial

yref_notrim = np.broadcast_to(yref, (2, len(t)))

_t, y_00 = initial_response(sys, T=t, X0=x0, output=0)

np.testing.assert_array_almost_equal(y_00, yref, decimal=4)

_t, y_11 = initial_response(sys, T=t, X0=x0, output=1)

np.testing.assert_array_almost_equal(y_11, yref, decimal=4)

_t, yy = initial_response(sys, T=t, X0=x0)

np.testing.assert_array_almost_equal(yy, yref_notrim, decimal=4)

@pytest.mark.parametrize("tsystem",

def test_forced_response_step(self, tsystem):

"""Test forced response of SISO systems as step response"""

sys = tsystem.sys

t = tsystem.t

u = np.ones_like(t, dtype=float)

yref = tsystem.ystep

tout, yout = forced_response(sys, t, u)

np.testing.assert_array_almost_equal(tout, t)

np.testing.assert_array_almost_equal(yout, yref, decimal=4)

@pytest.mark.parametrize("u",

@pytest.mark.parametrize("tsystem", ["siso_ss1", "siso_tf2"],

def test_forced_response_initial(self, tsystem, u):

"""Test forced response of SISO system as intitial response."""

sys = tsystem.sys

t = tsystem.t

x0 = tsystem.X0

yref = tsystem.yinitial

if isinstance(sys, StateSpace):

tout, yout = forced_response(sys, t, u, X0=x0)

np.testing.assert_array_almost_equal(tout, t)

np.testing.assert_array_almost_equal(yout, yref, decimal=4)

else:

with pytest.warns(UserWarning, match="Non-zero initial condition "

"given for transfer function"):

tout, yout = forced_response(sys, t, u, X0=x0)

@pytest.mark.parametrize("tsystem, useT",

def test_forced_response_mimo(self, tsystem, useT):

"""Test forced response of MIMO system"""

# first system: initial value, second system: step response

sys = tsystem.sys

t = tsystem.t

u = np.array([[0., 0, 0, 0, 0, 0, 0, 0, 0, 0],

[1., 1, 1, 1, 1, 1, 1, 1, 1, 1]])

x0 = np.array([[.5], [1], [0], [0]])

yref = np.vstack([tsystem.yinitial, tsystem.ystep])

if useT:

_t, yout = forced_response(sys, t, u, x0)

else:

_t, yout = forced_response(sys, U=u, X0=x0)

np.testing.assert_array_almost_equal(yout, yref, decimal=4)

@pytest.mark.usefixtures("editsdefaults")

def test_forced_response_legacy(self):

# Define a system for testing

sys = ct.rss(2, 1, 1)

T = np.linspace(0, 10, 10)

U = np.sin(T)

"""Make sure that legacy version of forced_response works"""

with pytest.warns(

UserWarning, match="NumPy matrix class no longer"):

ct.config.use_legacy_defaults("0.8.4")

# forced_response returns x by default

t, y = ct.step_response(sys, T)

t, y, x = ct.forced_response(sys, T, U)

ct.config.use_legacy_defaults("0.9.0")

# forced_response returns input/output by default

t, y = ct.step_response(sys, T)

t, y = ct.forced_response(sys, T, U)

t, y, x = ct.forced_response(sys, T, U, return_x=True)

@pytest.mark.parametrize(

def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr):

"""Test documented forced_response behavior for parameters T and U."""

if refattr == 'yinitial':

fr_kwargs['X0'] = tsystem.X0

t, y = forced_response(tsystem.sys, **fr_kwargs)

np.testing.assert_allclose(t, tsystem.t)

np.testing.assert_allclose(y, getattr(tsystem, refattr),

rtol=1e-3, atol=1e-5)

@pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True)

def test_forced_response_invalid_c(self, tsystem):

"""Test invalid parameters."""

with pytest.raises(TypeError,

match="StateSpace.*or.*TransferFunction"):

forced_response("not a system")

with pytest.raises(ValueError, match="T.*is mandatory for continuous"):

forced_response(tsystem.sys)

with pytest.raises(ValueError, match="time values must be equally "

"spaced"):

forced_response(tsystem.sys, [0, 0.1, 0.12, 0.4])

@pytest.mark.parametrize("tsystem", ["siso_dss2"], indirect=True)

def test_forced_response_invalid_d(self, tsystem):

"""Test invalid parameters dtime with sys.dt > 0."""

with pytest.raises(ValueError, match="can't both be zero"):

forced_response(tsystem.sys)

with pytest.raises(ValueError, match="Parameter `U`: Wrong shape"):

forced_response(tsystem.sys,

T=tsystem.t, U=np.random.randn(1, 12))

with pytest.raises(ValueError, match="Parameter `U`: Wrong shape"):

forced_response(tsystem.sys,

T=tsystem.t, U=np.random.randn(12))

with pytest.raises(ValueError, match="must match sampling time"):

forced_response(tsystem.sys, T=tsystem.t*0.9)

with pytest.raises(ValueError, match="must be multiples of "

"sampling time"):

forced_response(tsystem.sys, T=tsystem.t*1.1)

# but this is ok

forced_response(tsystem.sys, T=tsystem.t*2)

@pytest.mark.parametrize("u, x0, xtrue",

def test_lsim_double_integrator(self, u, x0, xtrue):

"""Test forced response of double integrator"""

# Note: scipy.signal.lsim fails if A is not invertible

A = np.array([[0., 1.],

[0., 0.]])

B = np.array([[0.],

[1.]])

C = np.array([[1., 0.]])

D = 0.

sys = StateSpace(A, B, C, D)

t = np.linspace(0, 1, 10)

_t, yout, xout = forced_response(sys, t, u, x0, return_x=True)

np.testing.assert_array_almost_equal(xout, xtrue, decimal=6)

ytrue = np.squeeze(np.asarray(C @ xtrue))

np.testing.assert_array_almost_equal(yout, ytrue, decimal=6)

@slycotonly

def test_step_robustness(self):

"Test robustness os step_response against denomiantors: gh-240"

# Create 2 input, 2 output system

num = [[[0], [1]], [[1], [0]]]

den1 = [[[1], [1,1]], [[1, 4], [1]]]

sys1 = TransferFunction(num, den1)

den2 = [[[1], [1e-10, 1, 1]], [[1, 4], [1]]] # slight perturbation

sys2 = TransferFunction(num, den2)

t1, y1 = step_response(sys1, input=0, T=2, T_num=100)

t2, y2 = step_response(sys2, input=0, T=2, T_num=100)

np.testing.assert_array_almost_equal(y1, y2)

@pytest.mark.parametrize(

def test_auto_generated_time_vector_tfinal(self, tfsys, tfinal):

"""Confirm a TF with a pole at p simulates for tfinal seconds"""

ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(tfsys)

np.testing.assert_allclose(ideal_tfinal, tfinal, rtol=1e-4)

T = _default_time_vector(tfsys)

np.testing.assert_allclose(T[-1], tfinal, atol=0.5*ideal_dt)

@pytest.mark.parametrize("wn, zeta", [(10, 0), (100, 0), (100, .1)])

def test_auto_generated_time_vector_dt_cont1(self, wn, zeta):

"""Confirm a TF with a natural frequency of wn rad/s gets a

dtref = 0.25133 / wn

tfsys = TransferFunction(1, [1, 2*zeta*wn, wn**2])

np.testing.assert_almost_equal(_ideal_tfinal_and_dt(tfsys)[1], dtref,

decimal=5)

def test_auto_generated_time_vector_dt_cont2(self):

"""A sampled tf keeps its dt"""

wn = 100

zeta = .1

tfsys = TransferFunction(1, [1, 2*zeta*wn, wn**2]).sample(.1)

tfinal, dt = _ideal_tfinal_and_dt(tfsys)

np.testing.assert_almost_equal(dt, .1)

T, _ = initial_response(tfsys)

np.testing.assert_almost_equal(np.diff(T[:2]), [.1])

def test_default_timevector_long(self):

"""Test long time vector"""

# TF with fast oscillations simulates only 5000 time steps

# even with long tfinal

wn = 100

tfsys = TransferFunction(1, [1, 0, wn**2])

tout = _default_time_vector(tfsys, tfinal=100)

assert len(tout) == 5000

@pytest.mark.parametrize("fun", [step_response,

def test_default_timevector_functions_c(self, fun):

"""Test that functions can calculate the time vector automatically"""

sys = TransferFunction(1, [1, .5, 0])

_tfinal, _dt = _ideal_tfinal_and_dt(sys)

# test impose number of time steps

tout, _ = fun(sys, T_num=10)

assert len(tout) == 10

# test impose final time

tout, _ = fun(sys, T=100.)

np.testing.assert_allclose(tout[-1], 100., atol=0.5*_dt)

@pytest.mark.parametrize("fun", [step_response,

@pytest.mark.parametrize("dt", [0.1, 0.112])

def test_default_timevector_functions_d(self, fun, dt):

"""Test that functions can calculate the time vector automatically"""

sys = TransferFunction(1, [1, .5, 0], dt)

# test impose number of time steps is ignored with dt given

tout, _ = fun(sys, T_num=15)

assert len(tout) != 15

# test impose final time

tout, _ = fun(sys, 100)

np.testing.assert_allclose(tout[-1], 100., atol=0.5*dt)