'statesp.remove_useless_states': False,

'statesp.latex_num_format': '.3g',

'statesp.latex_repr_type': 'partitioned',

'statesp.latex_maxsize': 10,

r"""StateSpace(A, B, C, D[, dt])

def __init__(self, *args, **kwargs):

"""StateSpace(A, B, C, D[, dt])

#

# Process positional arguments

#

if len(args) == 4:

# The user provided A, B, C, and D matrices.

A, B, C, D = args

elif len(args) == 5:

# Discrete time system

A, B, C, D, dt = args

if 'dt' in kwargs:

warn("received multiple dt arguments, "

"using positional arg dt = %s" % dt)

kwargs['dt'] = dt

args = args[:-1]

elif len(args) == 1:

# Use the copy constructor

if not isinstance(args[0], StateSpace):

raise TypeError(

"the one-argument constructor can only take in a "

"StateSpace object; received %s" % type(args[0]))

A = args[0].A

B = args[0].B

C = args[0].C

D = args[0].D

if 'dt' not in kwargs:

kwargs['dt'] = args[0].dt

else:

raise TypeError(

"Expected 1, 4, or 5 arguments; received %i." % len(args))

# Convert all matrices to standard form (sizes checked later)

A = _ssmatrix(A, square=True, name="A")

B = _ssmatrix(

B, axis=0 if np.asarray(B).ndim == 1 and len(B) == A.shape[0]

else 1, name="B")

C = _ssmatrix(

C, axis=1 if np.asarray(C).ndim == 1 and len(C) == A.shape[0]

else 0, name="C")

if np.isscalar(D) and D == 0 and B.shape[1] > 0 and C.shape[0] > 0:

# If D is a scalar zero, broadcast it to the proper size

D = np.zeros((C.shape[0], B.shape[1]))

D = _ssmatrix(D, name="D")

# If only direct term is present, adjust sizes of C and D if needed

if D.size > 0 and B.size == 0:

B = np.zeros((0, D.shape[1]))

if D.size > 0 and C.size == 0:

C = np.zeros((D.shape[0], 0))

# Matrices defining the linear system

self.A = A

self.B = B

self.C = C

self.D = D

# Determine if the system is static (memoryless)

static = (A.size == 0)

#

# Process keyword arguments

#

remove_useless_states = kwargs.pop(

'remove_useless_states',

config.defaults['statesp.remove_useless_states'])

# Process iosys keywords

defaults = args[0] if len(args) == 1 else \

{'inputs': B.shape[1], 'outputs': C.shape[0],

'states': A.shape[0]}

name, inputs, outputs, states, dt = _process_iosys_keywords(

kwargs, defaults, static=static)

# Create updfcn and outfcn

updfcn = lambda t, x, u, params: \

self.A @ np.atleast_1d(x) + self.B @ np.atleast_1d(u)

outfcn = lambda t, x, u, params: \

self.C @ np.atleast_1d(x) + self.D @ np.atleast_1d(u)

# Initialize NonlinearIOSystem object

super().__init__(

updfcn, outfcn,

name=name, inputs=inputs, outputs=outputs,

states=states, dt=dt, **kwargs)

# Reset shapes if the system is static

if static:

A.shape = (0, 0)

B.shape = (0, self.ninputs)

C.shape = (self.noutputs, 0)

# Check to make sure everything is consistent

_check_shape(A, self.nstates, self.nstates, name="A")

_check_shape(B, self.nstates, self.ninputs, name="B")

_check_shape(C, self.noutputs, self.nstates, name="C")

_check_shape(D, self.noutputs, self.ninputs, name="D")

#

# Final processing

#

# Check for states that don't do anything, and remove them

if remove_useless_states:

self._remove_useless_states()

#

# Class attributes

#

# These attributes are defined as class attributes so that they are

# documented properly. They are "overwritten" in __init__.

#

#: Number of system inputs.

#:

#: :meta hide-value:

ninputs = 0

#: Number of system outputs.

#:

#: :meta hide-value:

noutputs = 0

#: Number of system states.

#:

#: :meta hide-value:

nstates = 0

#: Dynamics matrix.

#:

#: :meta hide-value:

A = []

#: Input matrix.

#:

#: :meta hide-value:

B = []

#: Output matrix.

#:

#: :meta hide-value:

C = []

#: Direct term.

#:

#: :meta hide-value:

D = []

#

# Getter and setter functions for legacy state attributes

#

# For this iteration, generate a deprecation warning whenever the

# getter/setter is called. For a future iteration, turn it into a

# future warning, so that users will see it.

#

def _get_states(self):

warn("The StateSpace `states` attribute will be deprecated in a "

"future release. Use `nstates` instead.",

FutureWarning, stacklevel=2)

return self.nstates

def _set_states(self, value):

warn("The StateSpace `states` attribute will be deprecated in a "

"future release. Use `nstates` instead.",

FutureWarning, stacklevel=2)

self.nstates = value

#: Deprecated attribute; use `nstates` instead.

#:

#: The `state` attribute was used to store the number of states for : a

#: state space system. It is no longer used. If you need to access the

#: number of states, use `nstates`.

states = property(_get_states, _set_states)

def _remove_useless_states(self):

"""Check for states that don't do anything, and remove them.

# Search for useless states and get indices of these states.

ax1_A = np.where(~self.A.any(axis=1))[0]

ax1_B = np.where(~self.B.any(axis=1))[0]

ax0_A = np.where(~self.A.any(axis=0))[-1]

ax0_C = np.where(~self.C.any(axis=0))[-1]

useless_1 = np.intersect1d(ax1_A, ax1_B, assume_unique=True)

useless_2 = np.intersect1d(ax0_A, ax0_C, assume_unique=True)

useless = np.union1d(useless_1, useless_2)

# Remove the useless states.

self.A = delete(self.A, useless, 0)

self.A = delete(self.A, useless, 1)

self.B = delete(self.B, useless, 0)

self.C = delete(self.C, useless, 1)

# Remove any state names that we don't need

self.set_states(

[self.state_labels[i] for i in range(self.nstates)

if i not in useless])

def __str__(self):

"""Return string representation of the state space system."""

string = f"{InputOutputSystem.__str__(self)}\n\n"

string += "\n\n".join([

"{} = {}".format(Mvar,

"\n ".join(str(M).splitlines()))

for Mvar, M in zip(["A", "B", "C", "D"],

[self.A, self.B, self.C, self.D])])

return string

def _repr_eval_(self):

# Loadable format

out = "StateSpace(\n{A},\n{B},\n{C},\n{D}".format(

A=self.A.__repr__(), B=self.B.__repr__(),

C=self.C.__repr__(), D=self.D.__repr__())

out += super()._dt_repr(separator=",\n", space="")

if len(labels := super()._label_repr()) > 0:

out += ",\n" + labels

out += ")"

return out

def _repr_html_(self):

"""HTML representation of state-space model.

syssize = self.nstates + max(self.noutputs, self.ninputs)

if syssize > config.defaults['statesp.latex_maxsize']:

return None

elif config.defaults['statesp.latex_repr_type'] == 'partitioned':

return super()._repr_info_(html=True) + \

"\n" + self._latex_partitioned()

elif config.defaults['statesp.latex_repr_type'] == 'separate':

return super()._repr_info_(html=True) + \

"\n" + self._latex_separate()

else:

raise ValueError(

"Unknown statesp.latex_repr_type '{cfg}'".format(

cfg=config.defaults['statesp.latex_repr_type']))

def _latex_partitioned_stateless(self):

"""`Partitioned` matrix LaTeX representation for stateless systems

# Apply NumPy formatting

with np.printoptions(threshold=sys.maxsize):

D = eval(repr(self.D))

lines = [

r'$$',

(r'\left['

+ r'\begin{array}'

+ r'{' + 'rll' * self.ninputs + '}')

]

for Di in asarray(D):

lines.append('&'.join(_f2s(Dij) for Dij in Di)

+ '\\\\')

lines.extend([

r'\end{array}'

r'\right]',

r'$$'])

return '\n'.join(lines)

def _latex_partitioned(self):

"""Partitioned matrix LaTeX representation of state-space model

if self.nstates == 0:

return self._latex_partitioned_stateless()

# Apply NumPy formatting

with np.printoptions(threshold=sys.maxsize):

A, B, C, D = (

eval(repr(getattr(self, M))) for M in ['A', 'B', 'C', 'D'])

lines = [

r'$$',

(r'\left['

+ r'\begin{array}'

+ r'{' + 'rll' * self.nstates + '|' + 'rll' * self.ninputs + '}')

]

for Ai, Bi in zip(asarray(A), asarray(B)):

lines.append('&'.join([_f2s(Aij) for Aij in Ai]

+ [_f2s(Bij) for Bij in Bi])

+ '\\\\')

lines.append(r'\hline')

for Ci, Di in zip(asarray(C), asarray(D)):

lines.append('&'.join([_f2s(Cij) for Cij in Ci]

+ [_f2s(Dij) for Dij in Di])

+ '\\\\')

lines.extend([

r'\end{array}'

+ r'\right]',

r'$$'])

return '\n'.join(lines)

def _latex_separate(self):

"""Separate matrices LaTeX representation of state-space model

lines = [

r'$$',

r'\begin{array}{ll}',

]

def fmt_matrix(matrix, name):

matlines = [name

+ r' = \left[\begin{array}{'

+ 'rll' * matrix.shape[1]

+ '}']

for row in asarray(matrix):

matlines.append('&'.join(_f2s(entry) for entry in row)

+ '\\\\')

matlines.extend([

r'\end{array}'

r'\right]'])

return matlines

if self.nstates > 0:

lines.extend(fmt_matrix(self.A, 'A'))

lines.append('&')

lines.extend(fmt_matrix(self.B, 'B'))

lines.append('\\\\')

lines.extend(fmt_matrix(self.C, 'C'))

lines.append('&')

lines.extend(fmt_matrix(self.D, 'D'))

lines.extend([

r'\end{array}',

r'$$'])

return '\n'.join(lines)

# Negation of a system

def __neg__(self):

"""Negate a state space system."""

return StateSpace(self.A, self.B, -self.C, -self.D, self.dt)

# Addition of two state space systems (parallel interconnection)

def __add__(self, other):

"""Add two LTI systems (parallel connection)."""

from .xferfcn import TransferFunction

# Convert transfer functions to state space

if isinstance(other, TransferFunction):

# Convert the other argument to state space

other = _convert_to_statespace(other)

# Check for a couple of special cases

if isinstance(other, (int, float, complex, np.number)):

# Just adding a scalar; put it in the D matrix

A, B, C = self.A, self.B, self.C

D = self.D + other

dt = self.dt

elif isinstance(other, np.ndarray):

other = np.atleast_2d(other)

# Special case for SISO

if self.issiso():

self = np.ones_like(other) * self

if self.ninputs != other.shape[0]:

raise ValueError("array has incompatible shape")

A, B, C = self.A, self.B, self.C

D = self.D + other

dt = self.dt

elif not isinstance(other, StateSpace):

return NotImplemented # let other.__rmul__ handle it

else:

# Promote SISO object to compatible dimension

if self.issiso() and not other.issiso():

self = np.ones((other.noutputs, other.ninputs)) * self

elif not self.issiso() and other.issiso():

other = np.ones((self.noutputs, self.ninputs)) * other

# Check to make sure the dimensions are OK

if ((self.ninputs != other.ninputs) or

(self.noutputs != other.noutputs)):

raise ValueError(

"can't add systems with incompatible inputs and outputs")

dt = common_timebase(self.dt, other.dt)

# Concatenate the various arrays

A = concatenate((

concatenate((self.A, zeros((self.A.shape[0],

other.A.shape[-1]))), axis=1),

concatenate((zeros((other.A.shape[0], self.A.shape[-1])),

other.A), axis=1)), axis=0)

B = concatenate((self.B, other.B), axis=0)

C = concatenate((self.C, other.C), axis=1)

D = self.D + other.D

return StateSpace(A, B, C, D, dt)

# Right addition - just switch the arguments

def __radd__(self, other):

"""Right add two LTI systems (parallel connection)."""

return self + other

# Subtraction of two state space systems (parallel interconnection)

def __sub__(self, other):

"""Subtract two LTI systems."""

return self + (-other)

def __rsub__(self, other):

"""Right subtract two LTI systems."""

return other + (-self)

# Multiplication of two state space systems (series interconnection)

def __mul__(self, other):

"""Multiply two LTI objects (serial connection)."""

from .xferfcn import TransferFunction

# Convert transfer functions to state space

if isinstance(other, TransferFunction):

# Convert the other argument to state space

other = _convert_to_statespace(other)

# Check for a couple of special cases

if isinstance(other, (int, float, complex, np.number)):

# Just multiplying by a scalar; change the output

A, C = self.A, self.C

B = self.B * other

D = self.D * other

dt = self.dt

elif isinstance(other, np.ndarray):

other = np.atleast_2d(other)

# Special case for SISO

if self.issiso():

self = bdalg.append(*([self] * other.shape[0]))

# Dimension check after broadcasting

if self.ninputs != other.shape[0]:

raise ValueError("array has incompatible shape")

A, C = self.A, self.C

B = self.B @ other

D = self.D @ other

dt = self.dt

elif not isinstance(other, StateSpace):

return NotImplemented # let other.__rmul__ handle it

else:

# Promote SISO object to compatible dimension

if self.issiso() and not other.issiso():

self = bdalg.append(*([self] * other.noutputs))

elif not self.issiso() and other.issiso():

other = bdalg.append(*([other] * self.ninputs))

# Check to make sure the dimensions are OK

if self.ninputs != other.noutputs:

raise ValueError(

"can't multiply systems with incompatible"

" inputs and outputs")

dt = common_timebase(self.dt, other.dt)

# Concatenate the various arrays

A = concatenate(

(concatenate((other.A,

zeros((other.A.shape[0], self.A.shape[1]))),

axis=1),

concatenate((self.B @ other.C, self.A), axis=1)),

axis=0)

B = concatenate((other.B, self.B @ other.D), axis=0)

C = concatenate((self.D @ other.C, self.C), axis=1)

D = self.D @ other.D

return StateSpace(A, B, C, D, dt)

# Right multiplication of two state space systems (series interconnection)

# Just need to convert LH argument to a state space object

def __rmul__(self, other):

"""Right multiply two LTI objects (serial connection)."""

from .xferfcn import TransferFunction

# Convert transfer functions to state space

if isinstance(other, TransferFunction):

# Convert the other argument to state space

other = _convert_to_statespace(other)

# Check for a couple of special cases

if isinstance(other, (int, float, complex, np.number)):

# Just multiplying by a scalar; change the input

B = other * self.B

D = other * self.D

return StateSpace(self.A, B, self.C, D, self.dt)

elif isinstance(other, np.ndarray):

other = np.atleast_2d(other)

# Special case for SISO transfer function

if self.issiso():

self = bdalg.append(*([self] * other.shape[1]))

# Dimension check after broadcasting

if self.noutputs != other.shape[1]:

raise ValueError("array has incompatible shape")

C = other @ self.C

D = other @ self.D

return StateSpace(self.A, self.B, C, D, self.dt)

if not isinstance(other, StateSpace):

return NotImplemented

# Promote SISO object to compatible dimension

if self.issiso() and not other.issiso():

self = bdalg.append(*([self] * other.ninputs))

elif not self.issiso() and other.issiso():

other = bdalg.append(*([other] * self.noutputs))

return other * self

# TODO: general __truediv__ requires descriptor system support

def __truediv__(self, other):

"""Division of state space systems by TFs, FRDs, scalars, and arrays"""

# Let ``other.__rtruediv__`` handle it

try:

return self * (1 / other)

except ValueError:

return NotImplemented

def __rtruediv__(self, other):

"""Division by state space system"""

return other * self**-1

def __pow__(self, other):

"""Power of a state space system"""

if not type(other) == int:

raise ValueError("Exponent must be an integer")

if self.ninputs != self.noutputs:

# System must have same number of inputs and outputs

return NotImplemented

if other < -1:

return (self**-1)**(-other)

elif other == -1:

try:

Di = scipy.linalg.inv(self.D)

except scipy.linalg.LinAlgError:

# D matrix must be nonsingular

return NotImplemented

Ai = self.A - self.B @ Di @ self.C

Bi = self.B @ Di

Ci = -Di @ self.C

return StateSpace(Ai, Bi, Ci, Di, self.dt)

elif other == 0:

return StateSpace([], [], [], np.eye(self.ninputs), self.dt)

elif other == 1:

return self

elif other > 1:

return self * (self**(other - 1))

def __call__(self, x, squeeze=None, warn_infinite=True):

"""Evaluate system transfer function at point in complex plane.

# Use Slycot if available

out = self.horner(x, warn_infinite=warn_infinite)

return _process_frequency_response(self, x, out, squeeze=squeeze)

def slycot_laub(self, x):

"""Laub's method to evaluate response at complex frequency.

from slycot import tb05ad

# Make sure the argument is a 1D array of complex numbers

x_arr = np.atleast_1d(x).astype(complex, copy=False)

# Make sure that we are operating on a simple list

if len(x_arr.shape) > 1:

raise ValueError("input list must be 1D")

# preallocate

n = self.nstates

m = self.ninputs

p = self.noutputs

out = np.empty((p, m, len(x_arr)), dtype=complex)

# The first call both evaluates C(sI-A)^-1 B and also returns

# Hessenberg transformed matrices at, bt, ct.

result = tb05ad(n, m, p, x_arr[0], self.A, self.B, self.C, job='NG')

# When job='NG', result = (at, bt, ct, g_i, hinvb, info)

at = result[0]

bt = result[1]

ct = result[2]

# TB05AD frequency evaluation does not include direct feedthrough.

out[:, :, 0] = result[3] + self.D

# Now, iterate through the remaining frequencies using the

# transformed state matrices, at, bt, ct.

# Start at the second frequency, already have the first.

for kk, x_kk in enumerate(x_arr[1:]):

result = tb05ad(n, m, p, x_kk, at, bt, ct, job='NH')

# When job='NH', result = (g_i, hinvb, info)

# kk+1 because enumerate starts at kk = 0.

# but zero-th spot is already filled.

out[:, :, kk+1] = result[0] + self.D

return out

def horner(self, x, warn_infinite=True):

"""Evaluate value of transfer function using Horner's method.

# Make sure the argument is a 1D array of complex numbers

x_arr = np.atleast_1d(x).astype(complex, copy=False)

# return fast on systems with 0 or 1 state

if self.nstates == 0:

return self.D[:, :, np.newaxis] \

* np.ones_like(x_arr, dtype=complex)

elif self.nstates == 1:

with np.errstate(divide='ignore', invalid='ignore'):

out = self.C[:, :, np.newaxis] \

/ (x_arr - self.A[0, 0]) \

* self.B[:, :, np.newaxis] \

+ self.D[:, :, np.newaxis]

out[np.isnan(out)] = complex(np.inf, np.nan)

return out

try:

out = self.slycot_laub(x_arr)

except (ImportError, Exception):

# Fall back because either Slycot unavailable or cannot handle

# certain cases.

# Make sure that we are operating on a simple list

if len(x_arr.shape) > 1:

raise ValueError("input list must be 1D")

# Preallocate

out = empty((self.noutputs, self.ninputs, len(x_arr)),

dtype=complex)

# TODO: can this be vectorized?

for idx, x_idx in enumerate(x_arr):

try:

xr = solve(x_idx * eye(self.nstates) - self.A, self.B)

out[:, :, idx] = self.C @ xr + self.D

except LinAlgError:

# Issue a warning message, for consistency with xferfcn

if warn_infinite:

warn("singular matrix in frequency response",

RuntimeWarning)

# Evaluating at a pole. Return value depends if there

# is a zero at the same point or not.

if x_idx in self.zeros():

out[:, :, idx] = complex(np.nan, np.nan)

else:

out[:, :, idx] = complex(np.inf, np.nan)

return out

def freqresp(self, omega):

"""(deprecated) Evaluate transfer function at complex frequencies.

warn("StateSpace.freqresp(omega) will be removed in a "

"future release of python-control; use "

"sys.frequency_response(omega), or freqresp(sys, omega) in the "

"MATLAB compatibility module instead", FutureWarning)

return self.frequency_response(omega)

# Compute poles and zeros

def poles(self):

"""Compute the poles of a state space system."""

return eigvals(self.A).astype(complex) if self.nstates \

else np.array([])

def zeros(self):

"""Compute the zeros of a state space system."""

if not self.nstates:

return np.array([])

# Use AB08ND from Slycot if it's available, otherwise use

# scipy.lingalg.eigvals().

try:

from slycot import ab08nd

out = ab08nd(self.A.shape[0], self.B.shape[1], self.C.shape[0],

self.A, self.B, self.C, self.D)

nu = out[0]

if nu == 0:

return np.array([])

else:

# Use SciPy generalized eigenvalue function

return sp.linalg.eigvals(out[8][0:nu, 0:nu],

out[9][0:nu, 0:nu]).astype(complex)

except ImportError: # Slycot unavailable. Fall back to SciPy.

if self.C.shape[0] != self.D.shape[1]:

raise NotImplementedError(

"StateSpace.zero only supports systems with the same "

"number of inputs as outputs.")

# This implements the QZ algorithm for finding transmission zeros

# from

# https://dspace.mit.edu/bitstream/handle/1721.1/841/P-0802-06587335.pdf.

# The QZ algorithm solves the generalized eigenvalue problem: given

# `L = [A, B; C, D]` and `M = [I_nxn 0]`, find all finite lambda

# for which there exist nontrivial solutions of the equation

# `Lz - lamba Mz`.

#

# The generalized eigenvalue problem is only solvable if its

# arguments are square matrices.

L = concatenate((concatenate((self.A, self.B), axis=1),

concatenate((self.C, self.D), axis=1)), axis=0)

M = pad(eye(self.A.shape[0]), ((0, self.C.shape[0]),

(0, self.B.shape[1])), "constant")

return np.array([x for x in sp.linalg.eigvals(L, M,

overwrite_a=True)

if not isinf(x)], dtype=complex)

# Feedback around a state space system

def feedback(self, other=1, sign=-1):

"""Feedback interconnection between two LTI objects.