"""statesp.py State space representation and functions. This file contains the StateSpace class, which is used to represent linear systems in state space. This is the primary representation for the python-control library. """ # Python 3 compatibility (needs to go here) from __future__ import print_function from __future__ import division # for _convertToStateSpace """Copyright (c) 2010 by California Institute of Technology All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the California Institute of Technology nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CALTECH OR THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. Author: Richard M. Murray Date: 24 May 09 Revised: Kevin K. Chen, Dec 10 $Id$ """ import math import numpy as np from numpy import any, array, asarray, concatenate, cos, delete, \ dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze from numpy.random import rand, randn from numpy.linalg import solve, eigvals, matrix_rank from numpy.linalg.linalg import LinAlgError import scipy as sp from scipy.signal import lti, cont2discrete from warnings import warn from .lti import LTI, timebase, timebaseEqual, isdtime from . import config from copy import deepcopy __all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] # Define module default parameter values _statesp_defaults = { 'statesp.use_numpy_matrix': True, } def _ssmatrix(data, axis=1): """Convert argument to a (possibly empty) state space matrix. Parameters ---------- data : array, list, or string Input data defining the contents of the 2D array axis : 0 or 1 If input data is 1D, which axis to use for return object. The default is 1, corresponding to a row matrix. Returns ------- arr : 2D array, with shape (0, 0) if a is empty """ # Convert the data into an array or matrix, as configured # If data is passed as a string, use (deprecated?) matrix constructor if config.defaults['statesp.use_numpy_matrix'] or isinstance(data, str): arr = np.matrix(data, dtype=float) else: arr = np.array(data, dtype=float) ndim = arr.ndim shape = arr.shape # Change the shape of the array into a 2D array if (ndim > 2): raise ValueError("state-space matrix must be 2-dimensional") elif (ndim == 2 and shape == (1, 0)) or \ (ndim == 1 and shape == (0, )): # Passed an empty matrix or empty vector; change shape to (0, 0) shape = (0, 0) elif ndim == 1: # Passed a row or column vector shape = (1, shape[0]) if axis == 1 else (shape[0], 1) elif ndim == 0: # Passed a constant; turn into a matrix shape = (1, 1) # Create the actual object used to store the result return arr.reshape(shape) class StateSpace(LTI): """StateSpace(A, B, C, D[, dt]) A class for representing state-space models The StateSpace class is used to represent state-space realizations of linear time-invariant (LTI) systems: dx/dt = A x + B u y = C x + D u where u is the input, y is the output, and x is the state. The main data members are the A, B, C, and D matrices. The class also keeps track of the number of states (i.e., the size of A). The data format used to store state space matrices is set using the value of `config.defaults['use_numpy_matrix']`. If True (default), the state space elements are stored as `numpy.matrix` objects; otherwise they are `numpy.ndarray` objects. The :func:`~control.use_numpy_matrix` function can be used to set the storage type. Discrete-time state space system are implemented by using the 'dt' instance variable and setting it to the sampling period. If 'dt' is not None, then it must match whenever two state space systems are combined. Setting dt = 0 specifies a continuous system, while leaving dt = None means the system timebase is not specified. If 'dt' is set to True, the system will be treated as a discrete time system with unspecified sampling time. """ # Allow ndarray * StateSpace to give StateSpace._rmul_() priority __array_priority__ = 11 # override ndarray and matrix types def __init__(self, *args, **kw): """ StateSpace(A, B, C, D[, dt]) Construct a state space object. The default constructor is StateSpace(A, B, C, D), where A, B, C, D are matrices or equivalent objects. To create a discrete time system, use StateSpace(A, B, C, D, dt) where 'dt' is the sampling time (or True for unspecified sampling time). To call the copy constructor, call StateSpace(sys), where sys is a StateSpace object. """ if len(args) == 4: # The user provided A, B, C, and D matrices. (A, B, C, D) = args dt = None elif len(args) == 5: # Discrete time system (A, B, C, D, dt) = args 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 try: dt = args[0].dt except NameError: dt = None else: raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) # Process keyword arguments remove_useless = kw.get('remove_useless', True) # Convert all matrices to standard form A = _ssmatrix(A) B = _ssmatrix(B, axis=0) C = _ssmatrix(C, axis=1) 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) # TODO: use super here? LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0], dt=dt) self.A = A self.B = B self.C = C self.D = D self.states = A.shape[1] if 0 == self.states: # static gain # matrix's default "empty" shape is 1x0 A.shape = (0,0) B.shape = (0,self.inputs) C.shape = (self.outputs,0) # Check that the matrix sizes are consistent. if self.states != A.shape[0]: raise ValueError("A must be square.") if self.states != B.shape[0]: raise ValueError("A and B must have the same number of rows.") if self.states != C.shape[1]: raise ValueError("A and C must have the same number of columns.") if self.inputs != B.shape[1]: raise ValueError("B and D must have the same number of columns.") if self.outputs != C.shape[0]: raise ValueError("C and D must have the same number of rows.") # Check for states that don't do anything, and remove them. if remove_useless: self._remove_useless_states() def _remove_useless_states(self): """Check for states that don't do anything, and remove them. Scan the A, B, and C matrices for rows or columns of zeros. If the zeros are such that a particular state has no effect on the input-output dynamics, then remove that state from the A, B, and C matrices. """ # Search for useless states and get indices of these states. # # Note: shape from np.where depends on whether we are storing state # space objects as np.matrix or np.array. Code below will work # correctly in either case. 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) self.states = self.A.shape[0] self.inputs = self.B.shape[1] self.outputs = self.C.shape[0] def __str__(self): """String representation of the state space.""" str = "A = " + self.A.__str__() + "\n\n" str += "B = " + self.B.__str__() + "\n\n" str += "C = " + self.C.__str__() + "\n\n" str += "D = " + self.D.__str__() + "\n" # TODO: replace with standard calls to lti functions if (type(self.dt) == bool and self.dt is True): str += "\ndt unspecified\n" elif (not (self.dt is None) and type(self.dt) != bool and self.dt > 0): str += "\ndt = " + self.dt.__str__() + "\n" return str # represent as string, makes display work for IPython __repr__ = __str__ # 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).""" # 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 else: other = _convertToStateSpace(other) # Check to make sure the dimensions are OK if ((self.inputs != other.inputs) or (self.outputs != other.outputs)): raise ValueError("Systems have different shapes.") # Figure out the sampling time to use if self.dt is None and other.dt is not None: dt = other.dt # use dt from second argument elif (other.dt is None and self.dt is not None) or \ (timebaseEqual(self, other)): dt = self.dt # use dt from first argument else: raise ValueError("Systems have different sampling times") # 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).""" # Check for a couple of special cases if isinstance(other, (int, float, complex, np.number)): # Just multiplying by a scalar; change the output A, B = self.A, self.B C = self.C * other D = self.D * other dt = self.dt else: other = _convertToStateSpace(other) # Check to make sure the dimensions are OK if self.inputs != other.outputs: raise ValueError("C = A * B: A has %i column(s) (input(s)), \ but B has %i row(s)\n(output(s))." % (self.inputs, other.outputs)) # Figure out the sampling time to use if (self.dt == None and other.dt != None): dt = other.dt # use dt from second argument elif (other.dt == None and self.dt != None) or \ (timebaseEqual(self, other)): dt = self.dt # use dt from first argument else: raise ValueError("Systems have different sampling times") # Concatenate the various arrays A = concatenate( (concatenate((other.A, zeros((other.A.shape[0], self.A.shape[1]))), axis=1), concatenate((np.dot(self.B, other.C), self.A), axis=1)), axis=0) B = concatenate((other.B, np.dot(self.B, other.D)), axis=0) C = concatenate((np.dot(self.D, other.C), self.C),axis=1) D = np.dot(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 # TODO: __rmul__ only works for special cases (??) def __rmul__(self, other): """Right multiply two LTI objects (serial connection).""" # Check for a couple of special cases if isinstance(other, (int, float, complex, np.number)): # Just multiplying by a scalar; change the input A, C = self.A, self.C B = self.B * other D = self.D * other return StateSpace(A, B, C, D, self.dt) # is lti, and convertible? if isinstance(other, LTI): return _convertToStateSpace(other) * self # try to treat this as a matrix try: X = _ssmatrix(other) C = np.dot(X, self.C) D = np.dot(X, self.D) return StateSpace(self.A, self.B, C, D, self.dt) except Exception as e: print(e) pass raise TypeError("can't interconnect systems") # TODO: __div__ and __rdiv__ are not written yet. def __div__(self, other): """Divide two LTI systems.""" raise NotImplementedError("StateSpace.__div__ is not implemented yet.") def __rdiv__(self, other): """Right divide two LTI systems.""" raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.") def evalfr(self, omega): """Evaluate a SS system's transfer function at a single frequency. self._evalfr(omega) returns the value of the transfer function matrix with input value s = i * omega. """ warn("StateSpace.evalfr(omega) will be deprecated in a future " "release of python-control; use evalfr(sys, omega*1j) instead", PendingDeprecationWarning) return self._evalfr(omega) def _evalfr(self, omega): """Evaluate a SS system's transfer function at a single frequency""" # Figure out the point to evaluate the transfer function if isdtime(self, strict=True): dt = timebase(self) s = exp(1.j * omega * dt) if omega * dt > math.pi: warn("_evalfr: frequency evaluation above Nyquist frequency") else: s = omega * 1.j return self.horner(s) def horner(self, s): """Evaluate the systems's transfer function for a complex variable Returns a matrix of values evaluated at complex variable s. """ resp = np.dot(self.C, solve(s * eye(self.states) - self.A, self.B)) + self.D return array(resp) # Method for generating the frequency response of the system def freqresp(self, omega): """Evaluate the system's transfer func. at a list of freqs, omega. mag, phase, omega = self.freqresp(omega) Reports the frequency response of the system, G(j*omega) = mag*exp(j*phase) for continuous time. For discrete time systems, the response is evaluated around the unit circle such that G(exp(j*omega*dt)) = mag*exp(j*phase). Parameters ---------- omega : array A list of frequencies in radians/sec at which the system should be evaluated. The list can be either a python list or a numpy array and will be sorted before evaluation. Returns ------- mag : float The magnitude (absolute value, not dB or log10) of the system frequency response. phase : float The wrapped phase in radians of the system frequency response. omega : array The list of sorted frequencies at which the response was evaluated. """ # In case omega is passed in as a list, rather than a proper array. omega = np.asarray(omega) numFreqs = len(omega) Gfrf = np.empty((self.outputs, self.inputs, numFreqs), dtype=np.complex128) # Sort frequency and calculate complex frequencies on either imaginary # axis (continuous time) or unit circle (discrete time). omega.sort() if isdtime(self, strict=True): dt = timebase(self) cmplx_freqs = exp(1.j * omega * dt) if max(np.abs(omega)) * dt > math.pi: warn("freqresp: frequency evaluation above Nyquist frequency") else: cmplx_freqs = omega * 1.j # Do the frequency response evaluation. Use TB05AD from Slycot # if it's available, otherwise use the built-in horners function. try: from slycot import tb05ad n = np.shape(self.A)[0] m = self.inputs p = self.outputs # The first call both evaluates C(sI-A)^-1 B and also returns # Hessenberg transformed matrices at, bt, ct. result = tb05ad(n, m, p, cmplx_freqs[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. Gfrf[:, :, 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, cmplx_freqs_kk in enumerate(cmplx_freqs[1:numFreqs]): result = tb05ad(n, m, p, cmplx_freqs_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. Gfrf[:, :, kk+1] = result[0] + self.D except ImportError: # Slycot unavailable. Fall back to horner. for kk, cmplx_freqs_kk in enumerate(cmplx_freqs): Gfrf[:, :, kk] = self.horner(cmplx_freqs_kk) # mag phase omega return np.abs(Gfrf), np.angle(Gfrf), omega # Compute poles and zeros def pole(self): """Compute the poles of a state space system.""" return eigvals(self.A) if self.states else np.array([]) def zero(self): """Compute the zeros of a state space system.""" if not self.states: 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: return sp.linalg.eigvals(out[8][0:nu, 0:nu], out[9][0:nu, 0:nu]) 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)]) # Feedback around a state space system def feedback(self, other=1, sign=-1): """Feedback interconnection between two LTI systems.""" other = _convertToStateSpace(other) # Check to make sure the dimensions are OK if (self.inputs != other.outputs) or (self.outputs != other.inputs): raise ValueError("State space systems don't have compatible inputs/outputs for " "feedback.") # Figure out the sampling time to use if self.dt is None and other.dt is not None: dt = other.dt # use dt from second argument elif other.dt is None and self.dt is not None or timebaseEqual(self, other): dt = self.dt # use dt from first argument else: raise ValueError("Systems have different sampling times") A1 = self.A B1 = self.B C1 = self.C D1 = self.D A2 = other.A B2 = other.B C2 = other.C D2 = other.D F = eye(self.inputs) - sign * np.dot(D2, D1) if matrix_rank(F) != self.inputs: raise ValueError("I - sign * D2 * D1 is singular to working precision.") # Precompute F\D2 and F\C2 (E = inv(F)) # We can solve two linear systems in one pass, since the # coefficients matrix F is the same. Thus, we perform the LU # decomposition (cubic runtime complexity) of F only once! # The remaining back substitutions are only quadratic in runtime. E_D2_C2 = solve(F, concatenate((D2, C2), axis=1)) E_D2 = E_D2_C2[:, :other.inputs] E_C2 = E_D2_C2[:, other.inputs:] T1 = eye(self.outputs) + sign * np.dot(D1, E_D2) T2 = eye(self.inputs) + sign * np.dot(E_D2, D1) A = concatenate( (concatenate( (A1 + sign * np.dot(np.dot(B1, E_D2), C1), sign * np.dot(B1, E_C2)), axis=1), concatenate( (np.dot(B2, np.dot(T1, C1)), A2 + sign * np.dot(np.dot(B2, D1), E_C2)), axis=1)), axis=0) B = concatenate((np.dot(B1, T2), np.dot(np.dot(B2, D1), T2)), axis=0) C = concatenate((np.dot(T1, C1), sign * np.dot(D1, E_C2)), axis=1) D = np.dot(D1, T2) return StateSpace(A, B, C, D, dt) def lft(self, other, nu=-1, ny=-1): """Return the Linear Fractional Transformation. A definition of the LFT operator can be found in Appendix A.7, page 512 in the 2nd Edition, Multivariable Feedback Control by Sigurd Skogestad. An alternative definition can be found here: https://www.mathworks.com/help/control/ref/lft.html Parameters ---------- other : LTI The lower LTI system ny : int, optional Dimension of (plant) measurement output. nu : int, optional Dimension of (plant) control input. """ other = _convertToStateSpace(other) # maximal values for nu, ny if ny == -1: ny = min(other.inputs, self.outputs) if nu == -1: nu = min(other.outputs, self.inputs) # dimension check # TODO # Figure out the sampling time to use if (self.dt == None and other.dt != None): dt = other.dt # use dt from second argument elif (other.dt == None and self.dt != None) or \ timebaseEqual(self, other): dt = self.dt # use dt from first argument else: raise ValueError("Systems have different time bases") # submatrices A = self.A B1 = self.B[:, :self.inputs - nu] B2 = self.B[:, self.inputs - nu:] C1 = self.C[:self.outputs - ny, :] C2 = self.C[self.outputs - ny:, :] D11 = self.D[:self.outputs - ny, :self.inputs - nu] D12 = self.D[:self.outputs - ny, self.inputs - nu:] D21 = self.D[self.outputs - ny:, :self.inputs - nu] D22 = self.D[self.outputs - ny:, self.inputs - nu:] # submatrices Abar = other.A Bbar1 = other.B[:, :ny] Bbar2 = other.B[:, ny:] Cbar1 = other.C[:nu, :] Cbar2 = other.C[nu:, :] Dbar11 = other.D[:nu, :ny] Dbar12 = other.D[:nu, ny:] Dbar21 = other.D[nu:, :ny] Dbar22 = other.D[nu:, ny:] # well-posed check F = np.block([[np.eye(ny), -D22], [-Dbar11, np.eye(nu)]]) if matrix_rank(F) != ny + nu: raise ValueError("lft not well-posed to working precision.") # solve for the resulting ss by solving for [y, u] using [x, # xbar] and [w1, w2]. TH = np.linalg.solve(F, np.block( [[C2, np.zeros((ny, other.states)), D21, np.zeros((ny, other.inputs - ny))], [np.zeros((nu, self.states)), Cbar1, np.zeros((nu, self.inputs - nu)), Dbar12]] )) T11 = TH[:ny, :self.states] T12 = TH[:ny, self.states: self.states + other.states] T21 = TH[ny:, :self.states] T22 = TH[ny:, self.states: self.states + other.states] H11 = TH[:ny, self.states + other.states: self.states + other.states + self.inputs - nu] H12 = TH[:ny, self.states + other.states + self.inputs - nu:] H21 = TH[ny:, self.states + other.states: self.states + other.states + self.inputs - nu] H22 = TH[ny:, self.states + other.states + self.inputs - nu:] Ares = np.block([ [A + B2.dot(T21), B2.dot(T22)], [Bbar1.dot(T11), Abar + Bbar1.dot(T12)] ]) Bres = np.block([ [B1 + B2.dot(H21), B2.dot(H22)], [Bbar1.dot(H11), Bbar2 + Bbar1.dot(H12)] ]) Cres = np.block([ [C1 + D12.dot(T21), D12.dot(T22)], [Dbar21.dot(T11), Cbar2 + Dbar21.dot(T12)] ]) Dres = np.block([ [D11 + D12.dot(H21), D12.dot(H22)], [Dbar21.dot(H11), Dbar22 + Dbar21.dot(H12)] ]) return StateSpace(Ares, Bres, Cres, Dres, dt) def minreal(self, tol=0.0): """Calculate a minimal realization, removes unobservable and uncontrollable states""" if self.states: try: from slycot import tb01pd B = empty((self.states, max(self.inputs, self.outputs))) B[:,:self.inputs] = self.B C = empty((max(self.outputs, self.inputs), self.states)) C[:self.outputs,:] = self.C A, B, C, nr = tb01pd(self.states, self.inputs, self.outputs, self.A, B, C, tol=tol) return StateSpace(A[:nr,:nr], B[:nr,:self.inputs], C[:self.outputs,:nr], self.D) except ImportError: raise TypeError("minreal requires slycot tb01pd") else: return StateSpace(self) # TODO: add discrete time check def returnScipySignalLTI(self): """Return a list of a list of scipy.signal.lti objects. For instance, >>> out = ssobject.returnScipySignalLTI() >>> out[3][5] is a signal.scipy.lti object corresponding to the transfer function from the 6th input to the 4th output.""" # Preallocate the output. out = [[[] for _ in range(self.inputs)] for _ in range(self.outputs)] for i in range(self.outputs): for j in range(self.inputs): out[i][j] = lti(asarray(self.A), asarray(self.B[:, j]), asarray(self.C[i, :]), self.D[i, j]) return out def append(self, other): """Append a second model to the present model. The second model is converted to state-space if necessary, inputs and outputs are appended and their order is preserved""" if not isinstance(other, StateSpace): other = _convertToStateSpace(other) if self.dt != other.dt: raise ValueError("Systems must have the same time step") n = self.states + other.states m = self.inputs + other.inputs p = self.outputs + other.outputs A = zeros((n, n)) B = zeros((n, m)) C = zeros((p, n)) D = zeros((p, m)) A[:self.states, :self.states] = self.A A[self.states:, self.states:] = other.A B[:self.states, :self.inputs] = self.B B[self.states:, self.inputs:] = other.B C[:self.outputs, :self.states] = self.C C[self.outputs:, self.states:] = other.C D[:self.outputs, :self.inputs] = self.D D[self.outputs:, self.inputs:] = other.D return StateSpace(A, B, C, D, self.dt) def __getitem__(self, indices): """Array style access""" if len(indices) != 2: raise IOError('must provide indices of length 2 for state space') i = indices[0] j = indices[1] return StateSpace(self.A, self.B[:, j], self.C[i, :], self.D[i, j], self.dt) def sample(self, Ts, method='zoh', alpha=None): """Convert a continuous time system to discrete time Creates a discrete-time system from a continuous-time system by sampling. Multiple methods of conversion are supported. Parameters ---------- Ts : float Sampling period method : {"gbt", "bilinear", "euler", "backward_diff", "zoh"} Which method to use: * gbt: generalized bilinear transformation * bilinear: Tustin's approximation ("gbt" with alpha=0.5) * euler: Euler (or forward differencing) method ("gbt" with alpha=0) * backward_diff: Backwards differencing ("gbt" with alpha=1.0) * zoh: zero-order hold (default) alpha : float within [0, 1] The generalized bilinear transformation weighting parameter, which should only be specified with method="gbt", and is ignored otherwise Returns ------- sysd : StateSpace Discrete time system, with sampling rate Ts Notes ----- Uses the command 'cont2discrete' from scipy.signal Examples -------- >>> sys = StateSpace(0, 1, 1, 0) >>> sysd = sys.sample(0.5, method='bilinear') """ if not self.isctime(): raise ValueError("System must be continuous time system") sys = (self.A, self.B, self.C, self.D) Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha) return StateSpace(Ad, Bd, C, D, dt) def dcgain(self): """Return the zero-frequency gain The zero-frequency gain of a continuous-time state-space system is given by: .. math: G(0) = - C A^{-1} B + D and of a discrete-time state-space system by: .. math: G(1) = C (I - A)^{-1} B + D Returns ------- gain : ndarray An array of shape (outputs,inputs); the array will either be the zero-frequency (or DC) gain, or, if the frequency response is singular, the array will be filled with np.nan. """ try: if self.isctime(): gain = np.asarray(self.D-self.C.dot(np.linalg.solve(self.A, self.B))) else: gain = self.horner(1) except LinAlgError: # eigenvalue at DC gain = np.tile(np.nan, (self.outputs, self.inputs)) return np.squeeze(gain) # TODO: add discrete time check def _convertToStateSpace(sys, **kw): """Convert a system to state space form (if needed). If sys is already a state space, then it is returned. If sys is a transfer function object, then it is converted to a state space and returned. If sys is a scalar, then the number of inputs and outputs can be specified manually, as in: >>> sys = _convertToStateSpace(3.) # Assumes inputs = outputs = 1 >>> sys = _convertToStateSpace(1., inputs=3, outputs=2) In the latter example, A = B = C = 0 and D = [[1., 1., 1.] [1., 1., 1.]]. """ from .xferfcn import TransferFunction import itertools if isinstance(sys, StateSpace): if len(kw): raise TypeError("If sys is a StateSpace, _convertToStateSpace \ cannot take keywords.") # Already a state space system; just return it return sys elif isinstance(sys, TransferFunction): try: from slycot import td04ad if len(kw): raise TypeError("If sys is a TransferFunction, " "_convertToStateSpace cannot take keywords.") # Change the numerator and denominator arrays so that the transfer # function matrix has a common denominator. # matrices are also sized/padded to fit td04ad num, den, denorder = sys.minreal()._common_den() # transfer function to state space conversion now should work! ssout = td04ad('C', sys.inputs, sys.outputs, denorder, den, num, tol=0) states = ssout[0] return StateSpace(ssout[1][:states, :states], ssout[2][:states, :sys.inputs], ssout[3][:sys.outputs, :states], ssout[4], sys.dt) except ImportError: # No Slycot. Scipy tf->ss can't handle MIMO, but static # MIMO is an easy special case we can check for here maxn = max(max(len(n) for n in nrow) for nrow in sys.num) maxd = max(max(len(d) for d in drow) for drow in sys.den) if 1 == maxn and 1 == maxd: D = empty((sys.outputs, sys.inputs), dtype=float) for i, j in itertools.product(range(sys.outputs), range(sys.inputs)): D[i, j] = sys.num[i][j][0] / sys.den[i][j][0] return StateSpace([], [], [], D, sys.dt) else: if sys.inputs != 1 or sys.outputs != 1: raise TypeError("No support for MIMO without slycot") # TODO: do we want to squeeze first and check dimenations? # I think this will fail if num and den aren't 1-D after # the squeeze A, B, C, D = sp.signal.tf2ss(squeeze(sys.num), squeeze(sys.den)) return StateSpace(A, B, C, D, sys.dt) elif isinstance(sys, (int, float, complex, np.number)): if "inputs" in kw: inputs = kw["inputs"] else: inputs = 1 if "outputs" in kw: outputs = kw["outputs"] else: outputs = 1 # Generate a simple state space system of the desired dimension # The following Doesn't work due to inconsistencies in ltisys: # return StateSpace([[]], [[]], [[]], eye(outputs, inputs)) return StateSpace(0., zeros((1, inputs)), zeros((outputs, 1)), sys * ones((outputs, inputs))) # If this is a matrix, try to create a constant feedthrough try: D = _ssmatrix(sys) return StateSpace([], [], [], D) except Exception as e: print("Failure to assume argument is matrix-like in" \ " _convertToStateSpace, result %s" % e) raise TypeError("Can't convert given type to StateSpace system.") # TODO: add discrete time option def _rss_generate(states, inputs, outputs, type): """Generate a random state space. This does the actual random state space generation expected from rss and drss. type is 'c' for continuous systems and 'd' for discrete systems. """ # Probability of repeating a previous root. pRepeat = 0.05 # Probability of choosing a real root. Note that when choosing a complex # root, the conjugate gets chosen as well. So the expected proportion of # real roots is pReal / (pReal + 2 * (1 - pReal)). pReal = 0.6 # Probability that an element in B or C will not be masked out. pBCmask = 0.8 # Probability that an element in D will not be masked out. pDmask = 0.3 # Probability that D = 0. pDzero = 0.5 # Check for valid input arguments. if states < 1 or states % 1: raise ValueError("states must be a positive integer. states = %g." % states) if inputs < 1 or inputs % 1: raise ValueError("inputs must be a positive integer. inputs = %g." % inputs) if outputs < 1 or outputs % 1: raise ValueError("outputs must be a positive integer. outputs = %g." % outputs) # Make some poles for A. Preallocate a complex array. poles = zeros(states) + zeros(states) * 0.j i = 0 while i < states: if rand() < pRepeat and i != 0 and i != states - 1: # Small chance of copying poles, if we're not at the first or last # element. if poles[i-1].imag == 0: # Copy previous real pole. poles[i] = poles[i-1] i += 1 else: # Copy previous complex conjugate pair of poles. poles[i:i+2] = poles[i-2:i] i += 2 elif rand() < pReal or i == states - 1: # No-oscillation pole. if type == 'c': poles[i] = -exp(randn()) + 0.j elif type == 'd': poles[i] = 2. * rand() - 1. i += 1 else: # Complex conjugate pair of oscillating poles. if type == 'c': poles[i] = complex(-exp(randn()), 3. * exp(randn())) elif type == 'd': mag = rand() phase = 2. * math.pi * rand() poles[i] = complex(mag * cos(phase), mag * sin(phase)) poles[i+1] = complex(poles[i].real, -poles[i].imag) i += 2 # Now put the poles in A as real blocks on the diagonal. A = zeros((states, states)) i = 0 while i < states: if poles[i].imag == 0: A[i, i] = poles[i].real i += 1 else: A[i, i] = A[i+1, i+1] = poles[i].real A[i, i+1] = poles[i].imag A[i+1, i] = -poles[i].imag i += 2 # Finally, apply a transformation so that A is not block-diagonal. while True: T = randn(states, states) try: A = dot(solve(T, A), T) # A = T \ A * T break except LinAlgError: # In the unlikely event that T is rank-deficient, iterate again. pass # Make the remaining matrices. B = randn(states, inputs) C = randn(outputs, states) D = randn(outputs, inputs) # Make masks to zero out some of the elements. while True: Bmask = rand(states, inputs) < pBCmask if any(Bmask): # Retry if we get all zeros. break while True: Cmask = rand(outputs, states) < pBCmask if any(Cmask): # Retry if we get all zeros. break if rand() < pDzero: Dmask = zeros((outputs, inputs)) else: Dmask = rand(outputs, inputs) < pDmask # Apply masks. B = B * Bmask C = C * Cmask D = D * Dmask return StateSpace(A, B, C, D) # Convert a MIMO system to a SISO system # TODO: add discrete time check def _mimo2siso(sys, input, output, warn_conversion=False): #pylint: disable=W0622 """ Convert a MIMO system to a SISO system. (Convert a system with multiple inputs and/or outputs, to a system with a single input and output.) The input and output that are used in the SISO system can be selected with the parameters ``input`` and ``output``. All other inputs are set to 0, all other outputs are ignored. If ``sys`` is already a SISO system, it will be returned unaltered. Parameters ---------- sys : StateSpace Linear (MIMO) system that should be converted. input : int Index of the input that will become the SISO system's only input. output : int Index of the output that will become the SISO system's only output. warn_conversion : bool, optional If `True`, print a message when sys is a MIMO system, warning that a conversion will take place. Default is False. Returns sys : StateSpace The converted (SISO) system. """ if not (isinstance(input, int) and isinstance(output, int)): raise TypeError("Parameters ``input`` and ``output`` must both " "be integer numbers.") if not (0 <= input < sys.inputs): raise ValueError("Selected input does not exist. " "Selected input: {sel}, " "number of system inputs: {ext}." .format(sel=input, ext=sys.inputs)) if not (0 <= output < sys.outputs): raise ValueError("Selected output does not exist. " "Selected output: {sel}, " "number of system outputs: {ext}." .format(sel=output, ext=sys.outputs)) #Convert sys to SISO if necessary if sys.inputs > 1 or sys.outputs > 1: if warn_conversion: warn("Converting MIMO system to SISO system. " "Only input {i} and output {o} are used." .format(i=input, o=output)) # $X = A*X + B*U # Y = C*X + D*U new_B = sys.B[:, input] new_C = sys.C[output, :] new_D = sys.D[output, input] sys = StateSpace(sys.A, new_B, new_C, new_D, sys.dt) return sys def _mimo2simo(sys, input, warn_conversion=False): # pylint: disable=W0622 """ Convert a MIMO system to a SIMO system. (Convert a system with multiple inputs and/or outputs, to a system with a single input but possibly multiple outputs.) The input that is used in the SIMO system can be selected with the parameter ``input``. All other inputs are set to 0, all other outputs are ignored. If ``sys`` is already a SIMO system, it will be returned unaltered. Parameters ---------- sys: StateSpace Linear (MIMO) system that should be converted. input: int Index of the input that will become the SIMO system's only input. warn_conversion: bool If True: print a warning message when sys is a MIMO system. Warn that a conversion will take place. Returns ------- sys: StateSpace The converted (SIMO) system. """ if not (isinstance(input, int)): raise TypeError("Parameter ``input`` be an integer number.") if not (0 <= input < sys.inputs): raise ValueError("Selected input does not exist. " "Selected input: {sel}, " "number of system inputs: {ext}." .format(sel=input, ext=sys.inputs)) # Convert sys to SISO if necessary if sys.inputs > 1: if warn_conversion: warn("Converting MIMO system to SIMO system. " "Only input {i} is used." .format(i=input)) # $X = A*X + B*U # Y = C*X + D*U new_B = sys.B[:, input] new_D = sys.D[:, input] sys = StateSpace(sys.A, new_B, sys.C, new_D, sys.dt) return sys def ss(*args): """ss(A, B, C, D[, dt]) Create a state space system. The function accepts either 1, 4 or 5 parameters: ``ss(sys)`` Convert a linear system into space system form. Always creates a new system, even if sys is already a StateSpace object. ``ss(A, B, C, D)`` Create a state space system from the matrices of its state and output equations: .. math:: \\dot x = A \\cdot x + B \\cdot u y = C \\cdot x + D \\cdot u ``ss(A, B, C, D, dt)`` Create a discrete-time state space system from the matrices of its state and output equations: .. math:: x[k+1] = A \\cdot x[k] + B \\cdot u[k] y[k] = C \\cdot x[k] + D \\cdot u[ki] The matrices can be given as *array like* data types or strings. Everything that the constructor of :class:`numpy.matrix` accepts is permissible here too. Parameters ---------- sys: StateSpace or TransferFunction A linear system A: array_like or string System matrix B: array_like or string Control matrix C: array_like or string Output matrix D: array_like or string Feed forward matrix dt: If present, specifies the sampling period and a discrete time system is created Returns ------- out: :class:`StateSpace` The new linear system Raises ------ ValueError if matrix sizes are not self-consistent See Also -------- StateSpace tf ss2tf tf2ss Examples -------- >>> # Create a StateSpace object from four "matrices". >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.") >>> # Convert a TransferFunction to a StateSpace object. >>> sys_tf = tf([2.], [1., 3]) >>> sys2 = ss(sys_tf) """ if len(args) == 4 or len(args) == 5: return StateSpace(*args) elif len(args) == 1: from .xferfcn import TransferFunction sys = args[0] if isinstance(sys, StateSpace): return deepcopy(sys) elif isinstance(sys, TransferFunction): return tf2ss(sys) else: raise TypeError("ss(sys): sys must be a StateSpace or \ TransferFunction object. It is %s." % type(sys)) else: raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) def tf2ss(*args): """tf2ss(sys) Transform a transfer function to a state space system. The function accepts either 1 or 2 parameters: ``tf2ss(sys)`` Convert a linear system into transfer function form. Always creates a new system, even if sys is already a TransferFunction object. ``tf2ss(num, den)`` Create a transfer function system from its numerator and denominator polynomial coefficients. For details see: :func:`tf` Parameters ---------- sys: LTI (StateSpace or TransferFunction) A linear system num: array_like, or list of list of array_like Polynomial coefficients of the numerator den: array_like, or list of list of array_like Polynomial coefficients of the denominator Returns ------- out: StateSpace New linear system in state space form Raises ------ ValueError if `num` and `den` have invalid or unequal dimensions, or if an invalid number of arguments is passed in TypeError if `num` or `den` are of incorrect type, or if sys is not a TransferFunction object See Also -------- ss tf ss2tf Examples -------- >>> num = [[[1., 2.], [3., 4.]], [[5., 6.], [7., 8.]]] >>> den = [[[9., 8., 7.], [6., 5., 4.]], [[3., 2., 1.], [-1., -2., -3.]]] >>> sys1 = tf2ss(num, den) >>> sys_tf = tf(num, den) >>> sys2 = tf2ss(sys_tf) """ from .xferfcn import TransferFunction if len(args) == 2 or len(args) == 3: # Assume we were given the num, den return _convertToStateSpace(TransferFunction(*args)) elif len(args) == 1: sys = args[0] if not isinstance(sys, TransferFunction): raise TypeError("tf2ss(sys): sys must be a TransferFunction \ object.") return _convertToStateSpace(sys) else: raise ValueError("Needs 1 or 2 arguments; received %i." % len(args)) def rss(states=1, outputs=1, inputs=1): """ Create a stable *continuous* random state space object. Parameters ---------- states : integer Number of state variables inputs : integer Number of system inputs outputs : integer Number of system outputs Returns ------- sys : StateSpace The randomly created linear system Raises ------ ValueError if any input is not a positive integer See Also -------- drss Notes ----- If the number of states, inputs, or outputs is not specified, then the missing numbers are assumed to be 1. The poles of the returned system will always have a negative real part. """ return _rss_generate(states, inputs, outputs, 'c') def drss(states=1, outputs=1, inputs=1): """ Create a stable *discrete* random state space object. Parameters ---------- states : integer Number of state variables inputs : integer Number of system inputs outputs : integer Number of system outputs Returns ------- sys : StateSpace The randomly created linear system Raises ------ ValueError if any input is not a positive integer See Also -------- rss Notes ----- If the number of states, inputs, or outputs is not specified, then the missing numbers are assumed to be 1. The poles of the returned system will always have a magnitude less than 1. """ return _rss_generate(states, inputs, outputs, 'd') def ssdata(sys): """ Return state space data objects for a system Parameters ---------- sys : LTI (StateSpace, or TransferFunction) LTI system whose data will be returned Returns ------- (A, B, C, D): list of matrices State space data for the system """ ss = _convertToStateSpace(sys) return ss.A, ss.B, ss.C, ss.D