transpose=False):

"""Helper function for checking array_like parameters.

# convert nearly everything to an array.

out_array = np.asarray(in_obj)

if (transpose):

out_array = np.transpose(out_array)

# Test element data type, elements must be numbers

legal_kinds = set(("i", "f", "c")) # integer, float, complex

if out_array.dtype.kind not in legal_kinds:

err_msg = "Wrong element data type: '{d}'. Array elements " \

"must be numbers.".format(d=str(out_array.dtype))

raise TypeError(err_msg_start + err_msg)

# If array is zero dimensional (in_obj is scalar):

# create array with legal shape filled with the original value.

if out_array.ndim == 0:

for s_legal in legal_shapes:

# search for shape that does not contain the special symbol any.

if "any" in s_legal:

continue

the_val = out_array[()]

out_array = np.empty(s_legal, 'd')

out_array.fill(the_val)

break

# Test shape

def shape_matches(s_legal, s_actual):

"""Test if two shape tuples match"""

# Array must have required number of dimensions

if len(s_legal) != len(s_actual):

return False

# All dimensions must contain required number of elements. Joker: "all"

for n_legal, n_actual in zip(s_legal, s_actual):

if n_legal == "any":

continue

if n_legal != n_actual:

return False

return True

# Iterate over legal shapes, and see if any matches out_array's shape.

for s_legal in legal_shapes:

if shape_matches(s_legal, out_array.shape):

break

else:

legal_shape_str = " or ".join([str(s) for s in legal_shapes])

err_msg = "Wrong shape (rows, columns): {a}. Expected: {e}." \

.format(e=legal_shape_str, a=str(out_array.shape))

raise ValueError(err_msg_start + err_msg)

# Convert shape

if squeeze:

out_array = np.squeeze(out_array)

# We don't want zero dimensional arrays

if out_array.shape == tuple():

out_array = out_array.reshape((1,))

interpolate=False, return_x=None, squeeze=None):

"""Simulate the output of a linear system.

if not isinstance(sys, (StateSpace, TransferFunction)):

raise TypeError('Parameter ``sys``: must be a ``StateSpace`` or'

' ``TransferFunction``)')

# If return_x was not specified, figure out the default

if return_x is None:

return_x = config.defaults['forced_response.return_x']

# If return_x is used for TransferFunction, issue a warning

if return_x and isinstance(sys, TransferFunction):

warnings.warn(

"return_x specified for a transfer function system. Internal "

"conversion to state space used; results may meaningless.")

# If we are passed a transfer function and X0 is non-zero, warn the user

if isinstance(sys, TransferFunction) and np.any(X0 != 0):

warnings.warn(

"Non-zero initial condition given for transfer function system. "

"Internal conversion to state space used; may not be consistent "

"with given X0.")

sys = _convert_to_statespace(sys)

A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \

np.asarray(sys.D)

# d_type = A.dtype

n_states = A.shape[0]

n_inputs = B.shape[1]

n_outputs = C.shape[0]

# Convert inputs to numpy arrays for easier shape checking

if U is not None:

U = np.asarray(U)

if T is not None:

# T must be array-like

T = np.asarray(T)

# Set and/or check time vector in discrete time case

if isdtime(sys):

if T is None:

if U is None or (U.ndim == 0 and U == 0.):

raise ValueError('Parameters ``T`` and ``U`` can\'t both be '

'zero for discrete-time simulation')

# Set T to equally spaced samples with same length as U

if U.ndim == 1:

n_steps = U.shape[0]

else:

n_steps = U.shape[1]

dt = 1. if sys.dt in [True, None] else sys.dt

T = np.array(range(n_steps)) * dt

else:

# Make sure the input vector and time vector have same length

if (U.ndim == 1 and U.shape[0] != T.shape[0]) or \

(U.ndim > 1 and U.shape[1] != T.shape[0]):

raise ValueError('Pamameter ``T`` must have same elements as'

' the number of columns in input array ``U``')

if U.ndim == 0:

U = np.full((n_inputs, T.shape[0]), U)

else:

if T is None:

raise ValueError('Parameter ``T`` is mandatory for continuous '

'time systems.')

# Test if T has shape (n,) or (1, n);

T = _check_convert_array(T, [('any',), (1, 'any')],

'Parameter ``T``: ', squeeze=True,

transpose=transpose)

n_steps = T.shape[0] # number of simulation steps

# equally spaced also implies strictly monotonic increase,

dt = (T[-1] - T[0]) / (n_steps - 1)

if not np.allclose(np.diff(T), dt):

raise ValueError("Parameter ``T``: time values must be equally "

"spaced.")

# create X0 if not given, test if X0 has correct shape

X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)],

'Parameter ``X0``: ', squeeze=True)

xout = np.zeros((n_states, n_steps))

xout[:, 0] = X0

yout = np.zeros((n_outputs, n_steps))

# Separate out the discrete and continuous time cases

if isctime(sys, strict=True):

# Solve the differential equation, copied from scipy.signal.ltisys.

dot = np.dot # Faster and shorter code

# Faster algorithm if U is zero

# (if not None, it was converted to array above)

if U is None or np.all(U == 0):

# Solve using matrix exponential

expAdt = sp.linalg.expm(A * dt)

for i in range(1, n_steps):

xout[:, i] = dot(expAdt, xout[:, i-1])

yout = dot(C, xout)

# General algorithm that interpolates U in between output points

else:

# Test if U has correct shape and type

legal_shapes = [(n_steps,), (1, n_steps)] if n_inputs == 1 else \

[(n_inputs, n_steps)]

U = _check_convert_array(U, legal_shapes,

'Parameter ``U``: ', squeeze=False,

transpose=transpose)

# convert 1D array to 2D array with only one row

if len(U.shape) == 1:

U = U.reshape(1, -1) # pylint: disable=E1103

# Algorithm: to integrate from time 0 to time dt, with linear

# interpolation between inputs u(0) = u0 and u(dt) = u1, we solve

# xdot = A x + B u, x(0) = x0

# udot = (u1 - u0) / dt, u(0) = u0.

#

# Solution is

# [ x(dt) ] [ A*dt B*dt 0 ] [ x0 ]

# [ u(dt) ] = exp [ 0 0 I ] [ u0 ]

# [u1 - u0] [ 0 0 0 ] [u1 - u0]

M = np.block([[A * dt, B * dt, np.zeros((n_states, n_inputs))],

[np.zeros((n_inputs, n_states + n_inputs)),

np.identity(n_inputs)],

[np.zeros((n_inputs, n_states + 2 * n_inputs))]])

expM = sp.linalg.expm(M)

Ad = expM[:n_states, :n_states]

Bd1 = expM[:n_states, n_states+n_inputs:]

Bd0 = expM[:n_states, n_states:n_states + n_inputs] - Bd1

for i in range(1, n_steps):

xout[:, i] = (dot(Ad, xout[:, i-1]) + dot(Bd0, U[:, i-1]) +

dot(Bd1, U[:, i]))

yout = dot(C, xout) + dot(D, U)

tout = T

else:

# Discrete type system => use SciPy signal processing toolbox

# sp.signal.dlsim assumes T[0] == 0

spT = T - T[0]

if sys.dt is not True and sys.dt is not None:

# Make sure that the time increment is a multiple of sampling time

# First make sure that time increment is bigger than sampling time

# (with allowance for small precision errors)

if dt < sys.dt and not np.isclose(dt, sys.dt):

raise ValueError("Time steps ``T`` must match sampling time")

# Now check to make sure it is a multiple (with check against

# sys.dt because floating point mod can have small errors

if not (np.isclose(dt % sys.dt, 0) or

np.isclose(dt % sys.dt, sys.dt)):

raise ValueError("Time steps ``T`` must be multiples of "

"sampling time")

sys_dt = sys.dt

# sp.signal.dlsim returns not enough samples if

# T[-1] - T[0] < sys_dt * decimation * (n_steps - 1)

# due to rounding errors.

# https://github.com/scipyscipy/blob/v1.6.1/scipy/signal/ltisys.py#L3462

scipy_out_samples = int(np.floor(spT[-1] / sys_dt)) + 1

if scipy_out_samples < n_steps:

# parantheses: order of evaluation is important

spT[-1] = spT[-1] * (n_steps / (spT[-1] / sys_dt + 1))

else:

sys_dt = dt # For unspecified sampling time, use time incr

# Discrete time simulation using signal processing toolbox

dsys = (A, B, C, D, sys_dt)

# Use signal processing toolbox for the discrete time simulation

# Transpose the input to match toolbox convention

tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), spT, X0)

tout = tout + T[0]

if not interpolate:

# If dt is different from sys.dt, resample the output

inc = int(round(dt / sys_dt))

tout = T # Return exact list of time steps

yout = yout[::inc, :]

xout = xout[::inc, :]

# Transpose the output and state vectors to match local convention

xout = np.transpose(xout)

yout = np.transpose(yout)

return _process_time_response(sys, tout, yout, xout, transpose=transpose,

sys, tout, yout, xout, transpose=None, return_x=False,

squeeze=None, input=None, output=None):

"""Process time response signals.

# If squeeze was not specified, figure out the default (might remain None)

if squeeze is None:

squeeze = config.defaults['control.squeeze_time_response']

# Determine if the system is SISO

issiso = sys.issiso() or (input is not None and output is not None)

# Figure out whether and how to squeeze output data

if squeeze is True: # squeeze all dimensions

yout = np.squeeze(yout)

elif squeeze is False: # squeeze no dimensions

pass

elif squeeze is None: # squeeze signals if SISO

if issiso:

if len(yout.shape) == 3:

yout = yout[0][0] # remove input and output

else:

yout = yout[0] # remove input

else:

raise ValueError("unknown squeeze value")

# Figure out whether and how to squeeze the state data

if issiso and xout is not None and len(xout.shape) > 2:

xout = xout[:, 0, :] # remove input

# See if we need to transpose the data back into MATLAB form

if transpose:

# Transpose time vector in case we are using np.matrix

tout = np.transpose(tout)

# For signals, put the last index (time) into the first slot

yout = np.transpose(yout, np.roll(range(yout.ndim), 1))

if xout is not None:

xout = np.transpose(xout, np.roll(range(xout.ndim), 1))

# Return time, output, and (optionally) state

"""Return a SISO or SIMO state-space version of sys.

# If squeeze was not specified, figure out the default

if squeeze is None:

squeeze = config.defaults['control.squeeze_time_response']

sys_ss = _convert_to_statespace(sys)

if sys_ss.issiso():

return squeeze, sys_ss

elif squeeze == None and (input is None or output is None):

# Don't squeeze outputs if resulting system turns out to be siso

# Note: if we expand input to allow a tuple, need to update this check

squeeze = False

warn = False

if input is None:

# issue warning if input is not given

warn = True

input = 0

if output is None:

return squeeze, _mimo2simo(sys_ss, input, warn_conversion=warn)

else:

transpose=False, return_x=False, squeeze=None):

# pylint: disable=W0622

"""Compute the step response for a linear system.

# 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=True)

U = np.ones_like(T)

# If we are passed a transfer function and X0 is non-zero, warn the user

if isinstance(sys, TransferFunction) and np.any(X0 != 0):

warnings.warn(

"Non-zero initial condition given for transfer function system. "

"Internal conversion to state space used; may not be consistent "

"with given X0.")

# Convert to state space so that we can simulate

sys = _convert_to_statespace(sys)

# Set up arrays to handle the output

ninputs = sys.ninputs if input is None else 1

noutputs = sys.noutputs if output is None else 1

yout = np.empty((noutputs, ninputs, np.asarray(T).size))

xout = np.empty((sys.nstates, ninputs, np.asarray(T).size))

# Simulate the response for each input

for i in range(sys.ninputs):

# If input keyword was specified, only simulate for that input

if isinstance(input, int) and i != input:

continue

# Create a set of single inputs system for simulation

squeeze, simo = _get_ss_simo(sys, i, output, squeeze=squeeze)

out = forced_response(simo, T, U, X0, transpose=False,

return_x=return_x, squeeze=True)

inpidx = i if input is None else 0

yout[:, inpidx, :] = out[1]

if return_x:

xout[:, i, :] = out[2]

return _process_time_response(

sys, out[0], yout, xout, transpose=transpose, return_x=return_x,

SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1, 0.9)):

"""

if isinstance(sysdata, (StateSpace, TransferFunction)):

if T is None or np.asarray(T).size == 1:

T = _default_time_vector(sysdata, N=T_num, tfinal=T, is_step=True)

T, Yout = step_response(sysdata, T, squeeze=False)

if yfinal:

InfValues = np.atleast_2d(yfinal)

else:

InfValues = np.atleast_2d(sysdata.dcgain())

retsiso = sysdata.issiso()

noutputs = sysdata.noutputs

ninputs = sysdata.ninputs

else:

# Time series of response data

errmsg = ("`sys` must be a LTI system, or time response data"

" with a shape following the python-control"

" time series data convention.")

try:

Yout = np.array(sysdata, dtype=float)

except ValueError:

raise ValueError(errmsg)

if Yout.ndim == 1 or (Yout.ndim == 2 and Yout.shape[0] == 1):

Yout = Yout[np.newaxis, np.newaxis, :]

retsiso = True

elif Yout.ndim == 3:

retsiso = False

else:

raise ValueError(errmsg)

if T is None or Yout.shape[2] != len(np.squeeze(T)):

raise ValueError("For time response data, a matching time vector"

" must be given")

T = np.squeeze(T)

noutputs = Yout.shape[0]

ninputs = Yout.shape[1]

InfValues = np.atleast_2d(yfinal) if yfinal else Yout[:, :, -1]

ret = []

for i in range(noutputs):

retrow = []

for j in range(ninputs):

yout = Yout[i, j, :]

# Steady state value

InfValue = InfValues[i, j]

sgnInf = np.sign(InfValue.real)

rise_time: float = np.NaN

settling_time: float = np.NaN

settling_min: float = np.NaN

settling_max: float = np.NaN

peak_value: float = np.Inf

peak_time: float = np.Inf

undershoot: float = np.NaN

overshoot: float = np.NaN

steady_state_value: complex = np.NaN

if not np.isnan(InfValue) and not np.isinf(InfValue):

# RiseTime

tr_lower_index = np.where(

sgnInf * (yout - RiseTimeLimits[0] * InfValue) >= 0

)[0][0]

tr_upper_index = np.where(

sgnInf * (yout - RiseTimeLimits[1] * InfValue) >= 0

)[0][0]

rise_time = T[tr_upper_index] - T[tr_lower_index]

# SettlingTime

settled = np.where(

np.abs(yout/InfValue-1) >= SettlingTimeThreshold)[0][-1]+1

# MIMO systems can have unsettled channels without infinite

# InfValue

if settled < len(T):

settling_time = T[settled]

settling_min = min((yout[tr_upper_index:]).min(), InfValue)

settling_max = max((yout[tr_upper_index:]).max(), InfValue)

# Overshoot

y_os = (sgnInf * yout).max()

dy_os = np.abs(y_os) - np.abs(InfValue)

if dy_os > 0:

overshoot = np.abs(100. * dy_os / InfValue)

else:

overshoot = 0

# Undershoot : InfValue and undershoot must have opposite sign

y_us_index = (sgnInf * yout).argmin()

y_us = yout[y_us_index]

if (sgnInf * y_us) < 0:

undershoot = (-100. * y_us / InfValue)

else:

undershoot = 0

# Peak

peak_index = np.abs(yout).argmax()

peak_value = np.abs(yout[peak_index])

peak_time = T[peak_index]

# SteadyStateValue

steady_state_value = InfValue

retij = {

'RiseTime': rise_time,

'SettlingTime': settling_time,

'SettlingMin': settling_min,

'SettlingMax': settling_max,

'Overshoot': overshoot,

'Undershoot': undershoot,

'Peak': peak_value,

'PeakTime': peak_time,

'SteadyStateValue': steady_state_value

}

retrow.append(retij)

ret.append(retrow)

transpose=False, return_x=False, squeeze=None):

# pylint: disable=W0622

"""Initial condition response of a linear system