transpose=False):

"""

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

"""Simulate the output of a linear system.

if not isinstance(sys, LTI):

raise TypeError('Parameter ``sys``: must be a ``LTI`` object. '

'(For example ``StateSpace`` or ``TransferFunction``)')

sys = _convertToStateSpace(sys)

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

np.asarray(sys.D)

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 = np.asarray(T)

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

if isdtime(sys, strict=True):

if T is None:

if U is None:

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]

T = np.array(range(n_steps)) * (1 if sys.dt is True else sys.dt)

else:

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

# TODO: allow interpolation of the input vector

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

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

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

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

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

# T must be array-like and values must be increasing.

# The length of T determines the length of the input vector.

if T is None:

raise ValueError('Parameter ``T``: must be array-like, and contain '

'(strictly monotonic) increasing numbers.')

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

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

transpose=transpose)

dt = T[1] - T[0]

if not np.allclose(T[1:] - T[:-1], dt):

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

"equally spaced.")

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

# 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):

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

dot = np.dot # Faster and shorter code

# Faster algorithm if U is zero

if U is None or (isinstance(U, (int, float)) and 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

if sys.dt is not True:

# 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

elif 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

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), T, X0)

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)

# Get rid of unneeded dimensions

if squeeze:

yout = np.squeeze(yout)

xout = np.squeeze(xout)

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

if transpose:

tout = np.transpose(tout)

yout = np.transpose(yout)

xout = np.transpose(xout)

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

sys_ss = _convertToStateSpace(sys)

if sys_ss.issiso():

return sys_ss

warn = False

if input is None:

# issue warning if input is not given

warn = True

input = 0

if output is None:

return _mimo2simo(sys_ss, input, warn_conversion=warn)

else:

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

# pylint: disable=W0622

"""Step response of a linear system

sys = _get_ss_simo(sys, input, output)

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)

T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose,

squeeze=squeeze)

if return_x:

return T, yout, xout

RiseTimeLimits=(0.1, 0.9)):

'''

sys = _get_ss_simo(sys)

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

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

T, yout = step_response(sys, T)

# Steady state value

InfValue = yout[-1]

# RiseTime

tr_lower_index = (np.where(yout >= RiseTimeLimits[0] * InfValue)[0])[0]

tr_upper_index = (np.where(yout >= RiseTimeLimits[1] * InfValue)[0])[0]

RiseTime = T[tr_upper_index] - T[tr_lower_index]

# SettlingTime

sup_margin = (1. + SettlingTimeThreshold) * InfValue

inf_margin = (1. - SettlingTimeThreshold) * InfValue

# find Steady State looking for the first point out of specified limits

for i in reversed(range(T.size)):

if((yout[i] <= inf_margin) | (yout[i] >= sup_margin)):

SettlingTime = T[i + 1]

break

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

return {

'RiseTime': RiseTime,

'SettlingTime': SettlingTime,

'SettlingMin': yout[tr_upper_index:].min(),

'SettlingMax': yout.max(),

'Overshoot': 100. * (yout.max() - InfValue) / (InfValue - yout[0]),

'Undershoot': yout.min(), # not very confident about this

'Peak': yout[PeakIndex],

'PeakTime': T[PeakIndex],

'SteadyStateValue': InfValue

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

# pylint: disable=W0622

"""Initial condition response of a linear system

sys = _get_ss_simo(sys, input, output)

# Create time and input vectors; checking is done in forced_response(...)

# The initial vector X0 is created in forced_response(...) if necessary

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

T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)

U = np.zeros_like(T)

T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,

squeeze=squeeze)

if return_x:

return T, yout, _xout

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

# pylint: disable=W0622

"""Impulse response of a linear system

sys = _get_ss_simo(sys, input, output)

# if system has direct feedthrough, can't simulate impulse response

# numerically

if np.any(sys.D != 0) and isctime(sys):

warnings.warn("System has direct feedthrough: ``D != 0``. The "

"infinite impulse at ``t=0`` does not appear in the "

"output.\n"

"Results may be meaningless!")

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

# Must be done here because it is used for computations below.

n_states = sys.A.shape[0]

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

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

# Compute T and U, no checks necessary, will be checked in forced_response

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

T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False)

U = np.zeros_like(T)

# Compute new X0 that contains the impulse

# We can't put the impulse into U because there is no numerical

# representation for it (infinitesimally short, infinitely high).

# See also: http://www.mathworks.com/support/tech-notes/1900/1901.html

if isctime(sys):

B = np.asarray(sys.B).squeeze()

new_X0 = B + X0

else:

new_X0 = X0

U[0] = 1.

T, yout, _xout = forced_response(sys, T, U, new_X0, transpose=transpose,

squeeze=squeeze)

if return_x:

return T, yout, _xout

"""helper function to compute ideal simulation duration tfinal and dt, the

sqrt_eps = np.sqrt(np.spacing(1.))

default_tfinal = 5 # Default simulation horizon

default_dt = 0.1

total_cycles = 5 # number of cycles for oscillating modes

pts_per_cycle = 25 # Number of points divide a period of oscillation

log_decay_percent = np.log(100) # Factor of reduction for real pole decays

if sys.is_static_gain():

tfinal = default_tfinal

dt = sys.dt if isdtime(sys, strict=True) else default_dt

elif isdtime(sys, strict=True):

dt = sys.dt

A = _convertToStateSpace(sys).A

tfinal = default_tfinal

p = eigvals(A)

# Array Masks

# unstable

m_u = (np.abs(p) >= 1 + sqrt_eps)

p_u, p = p[m_u], p[~m_u]

if p_u.size > 0:

m_u = (p_u.real < 0) & (np.abs(p_u.imag) < sqrt_eps)

t_emp = np.max(log_decay_percent / np.abs(np.log(p_u[~m_u])/dt))

tfinal = max(tfinal, t_emp)

# zero - negligible effect on tfinal

m_z = np.abs(p) < sqrt_eps

p = p[~m_z]

# Negative reals- treated as oscillary mode

m_nr = (p.real < 0) & (np.abs(p.imag) < sqrt_eps)

p_nr, p = p[m_nr], p[~m_nr]

if p_nr.size > 0:

t_emp = np.max(log_decay_percent / np.abs((np.log(p_nr)/dt).real))

tfinal = max(tfinal, t_emp)

# discrete integrators

m_int = (p.real - 1 < sqrt_eps) & (np.abs(p.imag) < sqrt_eps)

p_int, p = p[m_int], p[~m_int]

# pure oscillatory modes

m_w = (np.abs(np.abs(p) - 1) < sqrt_eps)

p_w, p = p[m_w], p[~m_w]

if p_w.size > 0:

t_emp = total_cycles * 2 * np.pi / np.abs(np.log(p_w)/dt).min()

tfinal = max(tfinal, t_emp)

if p.size > 0:

t_emp = log_decay_percent / np.abs((np.log(p)/dt).real).min()

tfinal = max(tfinal, t_emp)

if p_int.size > 0:

tfinal = tfinal * 5

else: # cont time

sys_ss = _convertToStateSpace(sys)

# Improve conditioning via balancing and zeroing tiny entries

# See <w,v> for [[1,2,0], [9,1,0.01], [1,2,10*np.pi]] before/after balance

b, (sca, perm) = matrix_balance(sys_ss.A, separate=True)

p, l, r = eig(b, left=True, right=True)

# Reciprocal of inner product <w,v> for each eigval, (bound the ~infs by 1e12)

# G = Transfer([1], [1,0,1]) gives zero sensitivity (bound by 1e-12)

eig_sens = np.reciprocal(maximum(1e-12, einsum('ij,ij->j', l, r).real))

eig_sens = minimum(1e12, eig_sens)

# Tolerances

p[np.abs(p) < np.spacing(eig_sens * norm(b, 1))] = 0.

# Incorporate balancing to outer factors

l[perm, :] *= np.reciprocal(sca)[:, None]

r[perm, :] *= sca[:, None]

w, v = sys_ss.C.dot(r), l.T.conj().dot(sys_ss.B)

origin = False

# Computing the "size" of the response of each simple mode

wn = np.abs(p)

if np.any(wn == 0.):

origin = True

dc = np.zeros_like(p, dtype=float)

# well-conditioned nonzero poles, np.abs just in case

ok = np.abs(eig_sens) <= 1/sqrt_eps

# the averaged t->inf response of each simple eigval on each i/o channel

# See, A = [[-1, k], [0, -2]], response sizes are k-dependent (that is

# R/L eigenvector dependent)

dc[ok] = norm(v[ok, :], axis=1)*norm(w[:, ok], axis=0)*eig_sens[ok]

dc[wn != 0.] /= wn[wn != 0] if is_step else 1.

dc[wn == 0.] = 0.

# double the oscillating mode magnitude for the conjugate

dc[p.imag != 0.] *= 2

# Now get rid of noncontributing integrators and simple modes if any

relevance = (dc > 0.1*dc.max()) | ~ok

psub = p[relevance]

wnsub = wn[relevance]

tfinal, dt = [], []

ints = wnsub == 0.

iw = (psub.imag != 0.) & (np.abs(psub.real) <= sqrt_eps)

# Pure imaginary?

if np.any(iw):

tfinal += (total_cycles * 2 * np.pi / wnsub[iw]).tolist()

dt += (2 * np.pi / pts_per_cycle / wnsub[iw]).tolist()

# The rest ~ts = log(%ss value) / exp(Re(eigval)t)

texp_mode = log_decay_percent / np.abs(psub[~iw & ~ints].real)

tfinal += texp_mode.tolist()

dt += minimum(texp_mode / 50,

(2 * np.pi / pts_per_cycle / wnsub[~iw & ~ints])).tolist()

# All integrators?

if len(tfinal) == 0:

return default_tfinal*5, default_dt*5

tfinal = np.max(tfinal)*(5 if origin else 1)

dt = np.min(dt)

"""Returns a time vector that has a reasonable number of points.

N_max = 5000

N_min_ct = 100 # min points for cont time systems

N_min_dt = 20 # more common to see just a few samples in discrete-time