diff --git a/control/bdalg.py b/control/bdalg.py index 08300e971..f5b56fb80 100644 --- a/control/bdalg.py +++ b/control/bdalg.py @@ -334,7 +334,7 @@ def connect(sys, Q, inputv, outputv): K[inp,outp-1] = 1. elif outp < 0 and -outp >= -sys.outputs: K[inp,-outp-1] = -1. - sys = sys.feedback(sp.matrix(K), sign=1) + sys = sys.feedback(sp.array(K), sign=1) # now trim Ytrim = sp.zeros( (len(outputv), sys.outputs) ) @@ -343,4 +343,4 @@ def connect(sys, Q, inputv, outputv): Utrim[u-1,i] = 1. for i,y in enumerate(outputv): Ytrim[i,y-1] = 1. - return sp.matrix(Ytrim)*sys*sp.matrix(Utrim) + return sp.array(Ytrim) * sys * sp.array(Utrim) diff --git a/control/config.py b/control/config.py index c90b9e16b..7e9725156 100644 --- a/control/config.py +++ b/control/config.py @@ -7,6 +7,10 @@ # files. For now, you can just choose between MATLAB and FBS default # values. +import warnings +import numpy as np +from .statesp import StateSpaceMatrix + # Bode plot defaults bode_dB = False # Bode plot magnitude units bode_deg = True # Bode Plot phase units @@ -14,6 +18,9 @@ bode_number_of_samples = None # Bode plot number of samples bode_feature_periphery_decade = 1.0 # Bode plot feature periphery in decades +# State space return type (change to StateSpaceMatrix at next major revision) +ss_return_type = np.matrix + def reset_defaults(): """Reset configuration values to their default values.""" @@ -22,6 +29,7 @@ def reset_defaults(): global bode_Hz; bode_Hz = False global bode_number_of_samples; bode_number_of_samples = None global bode_feature_periphery_decade; bode_feature_periphery_decade = 1.0 + global ss_return_type; ss_return_type = np.matrix # TODO: update # Set defaults to match MATLAB @@ -52,3 +60,49 @@ def use_fbs_defaults(): global bode_deg; bode_deg = True global bode_Hz; bode_Hz = False +# +# State space function return type +# +# These functions are used to set the return type for state space functions +# that return a matrix. In the original version of python-control, these +# functions returned a numpy.matrix object. To handle the eventual +# deprecation of the numpy.matrix type, the StateSpaceMatrix object type, +# which implements matrix multiplications and related numpy.matrix operations +# was created. To avoid breaking existing code, the return type from +# functions that used to return an np.matrix object now call the +# get_ss_return_type function to obtain the return type. The return type can +# also be set using the `return_type` keyword, overriding the default state +# space matrix return type. +# + + +# Set state space return type +def set_ss_return_type(subtype, warn=True): + global ss_return_type + # If return_type is np.matrix, issue a pending deprecation warning + if (subtype is np.matrix and warn): + warnings.warn("Return type numpy.matrix is soon to be deprecated.", + stacklevel=2) + ss_return_type = subtype + + +# Get the state space return type +def get_ss_return_type(subtype=None, warn=True): + global ss_return_type + return_type = ss_return_type if subtype is None else subtype + # If return_type is np.matrix, issue a pending deprecation warning + if (return_type is np.matrix and warn): + warnings.warn("Returning numpy.matrix, soon to be deprecated; " + "make sure calling code can handle nparray.", + stacklevel=2) + return return_type + + +# Function to turn on/off use of np.matrix type +def use_numpy_matrix(flag=True, warn=True): + if flag and warn: + warnings.warn("Return type numpy.matrix is soon to be deprecated.", + stacklevel=2) + set_ss_return_type(np.matrix, warn=False) + else: + set_ss_return_type(StateSpaceMatrix) diff --git a/control/frdata.py b/control/frdata.py index 82da8420e..0ac77edc1 100644 --- a/control/frdata.py +++ b/control/frdata.py @@ -52,7 +52,7 @@ from warnings import warn import numpy as np from numpy import angle, array, empty, ones, \ - real, imag, matrix, absolute, eye, linalg, where, dot + real, imag, absolute, eye, linalg, where, dot from scipy.interpolate import splprep, splev from .lti import LTI @@ -81,6 +81,10 @@ class FRD(LTI): """ + # Allow NDarray * StateSpace to give StateSpace._rmul_() priority + # https://docs.scipy.org/doc/numpy/reference/arrays.classes.html#numpy.class.__array_priority__ + __array_priority__ = 11 # override ndarray and matrix types + epsw = 1e-8 def __init__(self, *args, **kwargs): @@ -436,12 +440,13 @@ def feedback(self, other=1, sign=-1): # TODO: vectorize this # TODO: handle omega re-mapping for k, w in enumerate(other.omega): - fresp[:, :, k] = self.fresp[:, :, k].view(type=matrix)* \ + fresp[:, :, k] = np.dot( + self.fresp[:, :, k], linalg.solve( - eye(self.inputs) + - other.fresp[:, :, k].view(type=matrix) * - self.fresp[:, :, k].view(type=matrix), - eye(self.inputs)) + eye(self.inputs) + np.dot(other.fresp[:, :, k], + self.fresp[:, :, k]), + eye(self.inputs)) + ) return FRD(fresp, other.omega, smooth=(self.ifunc is not None)) diff --git a/control/mateqn.py b/control/mateqn.py index 7e842234a..e7d53dbb8 100644 --- a/control/mateqn.py +++ b/control/mateqn.py @@ -41,7 +41,7 @@ Author: Bjorn Olofsson """ -from scipy import shape, size, asarray, asmatrix, copy, zeros, eye, dot +from scipy import shape, size, array, asarray, copy, zeros, eye, dot from scipy.linalg import eigvals, solve_discrete_are, solve from .exception import ControlSlycot, ControlArgument @@ -703,8 +703,8 @@ def dare(A, B, Q, R, S=None, E=None, stabilizing=True): if S is not None or E is not None or not stabilizing: return dare_old(A, B, Q, R, S, E, stabilizing) else: - Rmat = asmatrix(R) - Qmat = asmatrix(Q) + Rmat = array(R, ndmin=2) + Qmat = array(Q, ndmin=2) X = solve_discrete_are(A, B, Qmat, Rmat) G = solve(B.T.dot(X).dot(B) + Rmat, B.T.dot(X).dot(A)) L = eigvals(A - B.dot(G)) diff --git a/control/modelsimp.py b/control/modelsimp.py index 087b6bfaf..ef038f528 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -49,6 +49,7 @@ from .lti import isdtime, isctime from .statesp import StateSpace from .statefbk import gram +from .config import get_ss_return_type __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal'] @@ -56,7 +57,7 @@ # The following returns the Hankel singular values, which are singular values #of the matrix formed by multiplying the controllability and observability #grammians -def hsvd(sys): +def hsvd(sys, return_type=None): """Calculate the Hankel singular values. Parameters @@ -66,9 +67,12 @@ def hsvd(sys): Returns ------- - H : Matrix + H : matrix A list of Hankel singular values + return_type: nparray subtype, optional (default = numpy.matrix) + Set the ndarray subtype for the return value + See Also -------- gram @@ -96,11 +100,12 @@ def hsvd(sys): w, v = np.linalg.eig(WoWc) hsv = np.sqrt(w) - hsv = np.matrix(hsv) + hsv = np.array(hsv, ndmin=2) # was np.matrix(hsv) hsv = np.sort(hsv) hsv = np.fliplr(hsv) - # Return the Hankel singular values - return hsv + + # Return the Hankel singular values (casting type, if needed) + return hsv.view(type=get_ss_return_type(return_type)) def modred(sys, ELIM, method='matchdc'): """ @@ -409,16 +414,17 @@ def markov(Y, U, M): """ # Convert input parameters to matrices (if they aren't already) - Ymat = np.mat(Y) - Umat = np.mat(U) + Ymat = np.array(Y) + Umat = np.array(U) n = np.size(U) # Construct a matrix of control inputs to invert UU = Umat for i in range(1, M-1): - newCol = np.vstack((0, UU[0:n-1,i-2])) + #! TODO: second index on UU doesn't seem right; could be neg or pos?? + newCol = np.vstack((0, np.reshape(UU[0:n-1,i-2], (-1,1)))) UU = np.hstack((UU, newCol)) - Ulast = np.vstack((0, UU[0:n-1,M-2])) + Ulast = np.vstack((0, np.reshape(UU[0:n-1,M-2], (-1,1)))) for i in range(n-1,0,-1): Ulast[i] = np.sum(Ulast[0:i-1]) UU = np.hstack((UU, Ulast)) diff --git a/control/statefbk.py b/control/statefbk.py index d6eef0e45..4ec7f4755 100644 --- a/control/statefbk.py +++ b/control/statefbk.py @@ -43,6 +43,7 @@ import numpy as np import scipy as sp from . import statesp +from .config import get_ss_return_type from .exception import ControlSlycot, ControlArgument, ControlDimension __all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'acker'] @@ -110,6 +111,8 @@ def place(A, B, p): result = place_poles(A_mat, B_mat, placed_eigs, method='YT') K = result.gain_matrix + + # TODO: add return type conversion return K @@ -128,16 +131,17 @@ def place_varga(A, B, p, dtime=False, alpha=None): Optional Parameters --------------- - dtime: False for continuous time pole placement or True for discrete time. - The default is dtime=False. - alpha: double scalar - If DICO='C', then place_varga will leave the eigenvalues with real - real part less than alpha untouched. - If DICO='D', the place_varga will leave eigenvalues with modulus - less than alpha untouched. - - By default (alpha=None), place_varga computes alpha such that all - poles will be placed. + dtime : bool + False for continuous time pole placement or True for discrete time. + The default is dtime=False. + alpha : double scalar + If DICO='C', then place_varga will leave the eigenvalues with real + real part less than alpha untouched. + If DICO='D', the place_varga will leave eigenvalues with modulus + less than alpha untouched. + + By default (alpha=None), place_varga computes alpha such that all + poles will be placed. Returns ------- @@ -216,10 +220,11 @@ def place_varga(A, B, p, dtime=False, alpha=None): A_mat, B_mat, placed_eigs, DICO) # Return the gain matrix, with MATLAB gain convention + # TODO: add return type conversion return -F # Contributed by Roberto Bucher -def acker(A, B, poles): +def acker(A, B, poles, return_type=None): """Pole placement using Ackermann method Call: @@ -231,6 +236,9 @@ def acker(A, B, poles): State and input matrix of the system poles: 1-d list Desired eigenvalue locations + return_type : ndarray subtype, optional + Set the ndarray subtype for the return value. The default value can + be set using the `~control.use_numpy_matrix` function. Returns ------- @@ -238,12 +246,12 @@ def acker(A, B, poles): Gains such that A - B K has given eigenvalues """ - # Convert the inputs to matrices - a = np.mat(A) - b = np.mat(B) + # Convert the inputs to matrices (arrays) + a = statesp.ssmatrix(A) + b = statesp.ssmatrix(B) # Make sure the system is controllable - ct = ctrb(A, B) + ct = ctrb(A, B, return_type=np.ndarray) if np.linalg.matrix_rank(ct) != a.shape[0]: raise ValueError("System not reachable; pole placement invalid") @@ -252,13 +260,13 @@ def acker(A, B, poles): # Place the poles using Ackermann's method n = np.size(p) - pmat = p[n-1]*a**0 + pmat = p[n-1] * np.linalg.matrix_power(a, 0) for i in np.arange(1,n): - pmat = pmat + p[n-i-1]*a**i + pmat = pmat + np.dot(p[n-i-1], np.linalg.matrix_power(a, i)) K = np.linalg.solve(ct, pmat) K = K[-1][:] # Extract the last row - return K + return K.view(type=get_ss_return_type(return_type)) def lqr(*args, **keywords): """lqr(A, B, Q, R[, N]) @@ -366,19 +374,25 @@ def lqr(*args, **keywords): S = X; E = w[0:nstates]; + # TODO: add return type conversion return K, S, E -def ctrb(A, B): + +def ctrb(A, B, return_type=None): """Controllabilty matrix Parameters ---------- - A, B: array_like or string + A, B : array_like or string Dynamics and input matrix of the system + return_type : ndarray subtype, optional + Set the ndarray subtype for the return value. The default value can + be set using the `~control.use_numpy_matrix` function. + Returns ------- - C: matrix + C : matrix Controllability matrix Examples @@ -386,59 +400,72 @@ def ctrb(A, B): >>> C = ctrb(A, B) """ - # Convert input parameters to matrices (if they aren't already) - amat = np.mat(A) - bmat = np.mat(B) + amat = statesp.ssmatrix(A) + bmat = statesp.ssmatrix(B) n = np.shape(amat)[0] # Construct the controllability matrix - ctrb = np.hstack([bmat] + [amat**i*bmat for i in range(1, n)]) - return ctrb + ctrb = np.hstack([bmat] + [np.dot(np.linalg.matrix_power(amat, i), bmat) + for i in range(1, n)]) + + # Return the observability matrix in the desired type + return ctrb.view(type=get_ss_return_type(return_type)) -def obsv(A, C): + +def obsv(A, C, return_type=None): """Observability matrix Parameters ---------- - A, C: array_like or string + A, C : array_like or string Dynamics and output matrix of the system + return_type : ndarray subtype, optional + Set the ndarray subtype for the return value. The default value can + be set using the `~control.use_numpy_matrix` function. Returns ------- - O: matrix + O : matrix Observability matrix Examples -------- >>> O = obsv(A, C) - """ - + """ # Convert input parameters to matrices (if they aren't already) - amat = np.mat(A) - cmat = np.mat(C) + amat = statesp.ssmatrix(A) + cmat = statesp.ssmatrix(C) n = np.shape(amat)[0] # Construct the observability matrix - obsv = np.vstack([cmat] + [cmat*amat**i for i in range(1, n)]) - return obsv + obsv = np.vstack([cmat] + [np.dot(cmat, np.linalg.matrix_power(amat, i)) + for i in range(1, n)]) -def gram(sys,type): + # Return the observability matrix in the desired type + return obsv.view(type=get_ss_return_type(return_type)) + + +def gram(sys, type, return_type=None): """Gramian (controllability or observability) Parameters ---------- - sys: StateSpace + sys : StateSpace State-space system to compute Gramian for - type: String + type : String Type of desired computation. - `type` is either 'c' (controllability) or 'o' (observability). To compute the - Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability) + `type` is either 'c' (controllability) or 'o' (observability). To + compute the Cholesky factors of gramians use 'cf' (controllability) or + 'of' (observability) + return_type : ndarray subtype, optional + Set the ndarray subtype for the return value. The default value can + be set using the `~control.use_numpy_matrix` function. Returns ------- - gram: array + gram : array Gramian of system Raises @@ -460,7 +487,6 @@ def gram(sys,type): >>> Ro = gram(sys,'of'), where Wo=Ro'*Ro """ - #Check for ss system object if not isinstance(sys,statesp.StateSpace): raise ValueError("System must be StateSpace!") @@ -498,7 +524,7 @@ def gram(sys,type): A = np.array(sys.A) # convert to NumPy array for slycot X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra) gram = X - return gram + return gram.view(type=get_ss_return_type(return_type)) elif type=='cf' or type=='of': #Compute cholesky factored gramian from slycot routine sb03od @@ -521,4 +547,4 @@ def gram(sys,type): C[0:n,0:m] = sys.C.transpose() X,scale,w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra) gram = X - return gram + return gram.view(type=get_ss_return_type(return_type)) diff --git a/control/statesp.py b/control/statesp.py index 7c201580e..9de75b001 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -10,7 +10,7 @@ # Python 3 compatibility (needs to go here) from __future__ import print_function -from __future__ import division # for _convertToStateSpace +from __future__ import division # for _convertToStateSpace """Copyright (c) 2010 by California Institute of Technology All rights reserved. @@ -53,8 +53,8 @@ import math import numpy as np -from numpy import all, angle, any, array, asarray, concatenate, cos, delete, \ - dot, empty, exp, eye, isinf, matrix, ones, pad, shape, sin, zeros, squeeze +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 @@ -62,32 +62,177 @@ from scipy.signal import lti, cont2discrete from warnings import warn from .lti import LTI, timebase, timebaseEqual, isdtime -from .xferfcn import _convert_to_transfer_function from copy import deepcopy -__all__ = ['StateSpace', 'ss', 'rss', 'drss', 'tf2ss', 'ssdata'] +__all__ = ['StateSpaceMatrix', 'StateSpace', 'ss', 'rss', 'drss', 'tf2ss', + 'ssdata', 'ssmatrix'] -def _matrix(a): - """Wrapper around numpy.matrix that reshapes empty matrices to be 0x0 +class StateSpaceMatrix(np.ndarray): + """StateSpaceMatrix(M, [axis=0]) + + A class for representing state-space matrices + + The StateSpaceMatrix class is used to represent the state-space + matrices that define a StateSpace system (A, B, C, D). It is + mainly a wrapper for the ndarray class, but it maintains matrices + as 2-dimensional arrays. It is similar to the NDmatrix class + (which is being deprecated), but with more limited functionality. + + In addition, for empty matrices the size of the StateSpaceMatrix + instance is set to 0x0 (needed within various StateSpace methods). + + """ + # Allow ndarray * StateSpace to give StateSpace._rmul_() priority + __array_priority__ = 20 # override ndarray and matrix types + + def __new__(subtype, data=[], axis=1, dtype=float, copy=True): + """Create a StateSpaceMatrix object + + Parameters + ---------- + data: array-like object defining matrix values + axis: if data is 1D, choose its axis in the 2D representation + + Returns + ------- + M: 2D array representing the matrix. If data = [], shape = (0,0) + """ + + # Legacy code + # self = np.matrix.__new__(subtype, data, dtype=float) + # if (self.shape == (1, 0)): + # self.shape = (0, 0) + # return self + + # See if this is already a StateSpaceMatrix(from np.matrix.__new__) + if isinstance(data, StateSpaceMatrix): + dtype2 = data.dtype + if (dtype is None): + dtype = dtype2 + if (dtype2 == dtype) and (not copy): + return data + return data.astype(dtype) + + # If data is passed as a string, use (deprecated?) matrix constructor + if isinstance(data, str): + data = np.matrix(data, copy=True) + + # Convert the data into an array + arr = np.array(data, dtype=dtype, copy=copy) + 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) + + # Make sure the data representation matches the shape we used + # (from np.matrix.__new__) + order = 'C' + if (ndim == 2) and arr.flags.fortran: + order = 'F' + + # Create the actual object used to store the result + self = np.ndarray.__new__(subtype, shape, arr.dtype, + buffer=arr, order=order) + return self + + # Override multiplication operation to emulate nd.matrix (vs elementwise) + def __mul__(self, other): + """Multiply or scale state-space matrices""" + # return np.matrix.__mul__(self, other) # legacy + # Check to see if arguments are (real) scalars in disguise + # (from np.matrix.__mul__) + if isinstance(other, (np.ndarray, list, tuple)): + # return np.dot(self, np.asmatrix(other)) # legacy + # Promote 1D vectors to row matrices and return ndarray + product = np.dot(self, np.array(other, copy=False, ndmin=2)) + return product if np.isrealobj(product) else np.asarray(product) + + if np.isscalar(other) or not hasattr(other, '__rmul__'): + product = np.dot(self, other) + return product if np.isrealobj(product) else np.asarray(product) + + raise NotImplementedError("Unknown object for StateSpaceMatrix " + "multiplication.") + + def __rmul__(self, other): + """Multiply or scale state-space matrices""" + # return np.matrix.__rmul__(self, other) # legacy + product = np.dot(other, self) + return product if np.isrealobj(product) else np.asarray(product) + + def __imul__(self, other): + return self.__mul__(other) + + def __pow__(self, n): + """Raise a (square) matrix to the (integer) power `n`""" + return np.linalg.matrix_power(self, n) + + def __rpow__(self, other): + raise NotImplementedError("StateSpaceMatrix powers are not " + "implemented.") + + def __ipow__(self, n): + return self.__pow__(n) + + def __getitem__(self, index): + """Get elements of a state-space matrix and return matrix""" + # return np.matrix.__getitem__(self, index) # legacy + self._getitem = True # to get back raw item (from np.matrix.__mul__) + try: + out = np.ndarray.__getitem__(self, index) + finally: + self._getitem = False + + if out.ndim == 0: + # If we get down to a scalar, return the actual scalar + return out[()] + if out.ndim == 1: + # Got a row/column vector; figure out what to return + # (from np.matrix.__mul__) + sh = out.shape[0] + try: + n = len(index) + except Exception: + n = 0 + if n > 1 and np.isscalar(index[1]): + out.shape = (sh, 1) + else: + out.shape = (1, sh) + return out + + +def ssmatrix(object, copy=True): + """Create a state space matrix. Parameters ---------- - a: sequence passed to numpy.matrix + object : array-like object defining matrix values + copy : bool, optional + If true (default), then the object is copied. Otherwise, a + copy will be made only if numpy.array returns a copy. Returns ------- - am: result of numpy.matrix(a), except if a is empty, am will be 0x0. + out : StateSpaceMatrix object - numpy.matrix([]) has size 1x0; for empty StateSpace objects, we - need 0x0 matrices, so use this instead of numpy.matrix in this - module. """ - from numpy import matrix - am = matrix(a, dtype=float) - if (1, 0) == am.shape: - am.shape = (0, 0) - return am + return StateSpaceMatrix(object, copy=copy) class StateSpace(LTI): @@ -95,8 +240,8 @@ class StateSpace(LTI): A class for representing state-space models - The StateSpace class is used to represent state-space realizations of linear - time-invariant (LTI) systems: + 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 @@ -113,8 +258,13 @@ class StateSpace(LTI): 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 + # https://docs.scipy.org/doc/numpy/reference/arrays.classes.html#numpy.class.__array_priority__ + __array_priority__ = 11 # override ndarray and matrix types + def __init__(self, *args): """ StateSpace(A, B, C, D[, dt]) @@ -139,8 +289,9 @@ def __init__(self, *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])) + raise TypeError("The one-argument constructor can only take " + "a StateSpace object. Received %s." % + type(args[0])) A = args[0].A B = args[0].B C = args[0].C @@ -150,9 +301,10 @@ def __init__(self, *args): except NameError: dt = None else: - raise ValueError("Needs 1 or 4 arguments; received %i." % len(args)) + raise ValueError("Needs 1 or 4 arguments; received %i." % + len(args)) - A, B, C, D = [_matrix(M) for M in (A, B, C, D)] + A, B, C, D = [ssmatrix(M) for M in (A, B, C, D)] # TODO: use super here? LTI.__init__(self, inputs=D.shape[1], outputs=D.shape[0], dt=dt) @@ -166,9 +318,9 @@ def __init__(self, *args): 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) + 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]: @@ -189,8 +341,8 @@ 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. + 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. """ @@ -198,8 +350,8 @@ def _remove_useless_states(self): # as an array. 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] + ax0_A = np.where(~self.A.any(axis=0))[0] + ax0_C = np.where(~self.C.any(axis=0))[0] 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) @@ -266,11 +418,11 @@ def __add__(self, other): # 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) + 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 @@ -310,13 +462,14 @@ def __mul__(self, 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)) + 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): + if (self.dt is None and other.dt is not None): dt = other.dt # use dt from second argument - elif (other.dt == None and self.dt != None) or \ + elif (other.dt is None and self.dt is not None) or \ (timebaseEqual(self, other)): dt = self.dt # use dt from first argument else: @@ -324,11 +477,12 @@ def __mul__(self, other): # 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) + (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) + C = concatenate((self.D * other.C, self.C), axis=1) D = self.D * other.D return StateSpace(A, B, C, D, dt) @@ -353,7 +507,7 @@ def __rmul__(self, other): # try to treat this as a matrix try: - X = _matrix(other) + X = ssmatrix(other) C = X * self.C D = X * self.D return StateSpace(self.A, self.B, C, D, self.dt) @@ -548,7 +702,8 @@ def zero(self): 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) + 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 @@ -559,13 +714,14 @@ def feedback(self, other=1, sign=-1): # 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.") + 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): + 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") @@ -581,7 +737,8 @@ def feedback(self, other=1, sign=-1): F = eye(self.inputs) - sign * D2 * D1 if matrix_rank(F) != self.inputs: - raise ValueError("I - sign * D2 * D1 is singular to working precision.") + 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 @@ -595,9 +752,12 @@ def feedback(self, other=1, sign=-1): T1 = eye(self.outputs) + sign * D1 * E_D2 T2 = eye(self.inputs) + sign * E_D2 * D1 - A = concatenate((concatenate((A1 + sign * B1 * E_D2 * C1, sign * B1 * E_C2), axis=1), - concatenate((B2 * T1 * C1, A2 + sign * B2 * D1 * E_C2), axis=1)), - axis=0) + A = concatenate(( + concatenate( + (A1 + sign * B1 * E_D2 * C1, sign * B1 * E_C2), axis=1), + concatenate( + (B2 * T1 * C1, A2 + sign * B2 * D1 * E_C2), axis=1)), + axis=0) B = concatenate((B1 * T2, B2 * D1 * T2), axis=0) C = concatenate((T1 * C1, sign * D1 * E_C2), axis=1) D = D1 * T2 @@ -634,9 +794,9 @@ def lft(self, other, nu=-1, ny=-1): # TODO # Figure out the sampling time to use - if (self.dt == None and other.dt != None): + if (self.dt is None and other.dt is not None): dt = other.dt # use dt from second argument - elif (other.dt == None and self.dt != None) or \ + elif (other.dt is None and self.dt is not None) or \ timebaseEqual(self, other): dt = self.dt # use dt from first argument else: @@ -668,22 +828,26 @@ def lft(self, other, nu=-1, ny=-1): 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]] + [[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] + 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] + 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)] @@ -712,19 +876,18 @@ def minreal(self, tol=0.0): try: from slycot import tb01pd B = empty((self.states, max(self.inputs, self.outputs))) - B[:,:self.inputs] = self.B + B[:, :self.inputs] = self.B C = empty((max(self.outputs, self.inputs), self.states)) - C[:self.outputs,:] = self.C + 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) + 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. @@ -743,7 +906,7 @@ def returnScipySignalLTI(self): 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, :]), asarray(self.D[i, j])) + asarray(self.C[i, :]), self.D[i, j]) return out @@ -780,7 +943,8 @@ def __getitem__(self, indices): 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) + 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 @@ -863,10 +1027,10 @@ def dcgain(self): 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: + 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) @@ -879,8 +1043,8 @@ def _convertToStateSpace(sys, **kw): import itertools if isinstance(sys, StateSpace): if len(kw): - raise TypeError("If sys is a StateSpace, _convertToStateSpace \ -cannot take keywords.") + raise TypeError("If sys is a StateSpace object, " + "_convertToStateSpace cannot take keywords.") # Already a state space system; just return it return sys @@ -901,8 +1065,10 @@ def _convertToStateSpace(sys, **kw): 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) + 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 @@ -912,7 +1078,8 @@ def _convertToStateSpace(sys, **kw): 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)): + 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: @@ -939,18 +1106,19 @@ def _convertToStateSpace(sys, **kw): # 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))) + sys * ones((outputs, inputs))) # If this is a matrix, try to create a constant feedthrough try: - D = _matrix(sys) + D = ssmatrix(sys) return StateSpace([], [], [], D) except Exception as e: - print("Failure to assume argument is matrix-like in" \ - " _convertToStateSpace, result %s" % 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. @@ -976,13 +1144,13 @@ def _rss_generate(states, inputs, outputs, type): # Check for valid input arguments. if states < 1 or states % 1: raise ValueError("states must be a positive integer. states = %g." % - states) + states) if inputs < 1 or inputs % 1: raise ValueError("inputs must be a positive integer. inputs = %g." % - inputs) + inputs) if outputs < 1 or outputs % 1: raise ValueError("outputs must be a positive integer. outputs = %g." % - outputs) + outputs) # Make some poles for A. Preallocate a complex array. poles = zeros(states) + zeros(states) * 0.j @@ -1070,7 +1238,7 @@ def _rss_generate(states, inputs, outputs, type): # Convert a MIMO system to a SISO system # TODO: add discrete time check def _mimo2siso(sys, input, output, warn_conversion=False): - #pylint: disable=W0622 + # 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.) @@ -1110,7 +1278,7 @@ def _mimo2siso(sys, input, output, warn_conversion=False): "Selected output: {sel}, " "number of system outputs: {ext}." .format(sel=output, ext=sys.outputs)) - #Convert sys to SISO if necessary + # Convert sys to SISO if necessary if sys.inputs > 1 or sys.outputs > 1: if warn_conversion: warn("Converting MIMO system to SISO system. " diff --git a/control/tests/canonical_test.py b/control/tests/canonical_test.py index 1f1828f50..17625412d 100644 --- a/control/tests/canonical_test.py +++ b/control/tests/canonical_test.py @@ -2,6 +2,7 @@ import unittest import numpy as np +import control as ct from control import ss, tf from control.canonical import canonical_form, reachable_form, \ observable_form, modal_form @@ -17,18 +18,19 @@ def test_reachable_form(self): coeffs = [1.0, 2.0, 3.0, 4.0, 1.0] A_true = np.polynomial.polynomial.polycompanion(coeffs) A_true = np.fliplr(np.rot90(A_true)) - B_true = np.matrix("1.0 0.0 0.0 0.0").T - C_true = np.matrix("1.0 1.0 1.0 1.0") + B_true = np.array([[1.0, 0.0, 0.0, 0.0]]).T + C_true = np.array([[1.0, 1.0, 1.0, 1.0]]) D_true = 42.0 # Perform a coordinate transform with a random invertible matrix - T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], - [-0.74855725, -0.39136285, -0.18142339, -0.50356997], - [-0.40688007, 0.81416369, 0.38002113, -0.16483334], - [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) - A = np.linalg.solve(T_true, A_true)*T_true + T_true = np.array( + [[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + [-0.74855725, -0.39136285, -0.18142339, -0.50356997], + [-0.40688007, 0.81416369, 0.38002113, -0.16483334], + [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) + A = np.dot(np.linalg.solve(T_true, A_true), T_true) B = np.linalg.solve(T_true, B_true) - C = C_true*T_true + C = np.dot(C_true, T_true) D = D_true # Create a state space system and convert it to the reachable canonical form @@ -50,9 +52,9 @@ def test_unreachable_system(self): """Test reachable canonical form with an unreachable system""" # Create an unreachable system - A = np.matrix("1.0 2.0 2.0; 4.0 5.0 5.0; 7.0 8.0 8.0") - B = np.matrix("1.0 1.0 1.0").T - C = np.matrix("1.0 1.0 1.0") + A = np.array([[1.0, 2.0, 2.0], [4.0, 5.0, 5.0], [7.0, 8.0, 8.0]]) + B = np.array([[1.0, 1.0, 1.0]]).T + C = np.array([[1.0, 1.0, 1.0]]) D = 42.0 sys = ss(A, B, C, D) @@ -63,19 +65,20 @@ def test_modal_form(self): """Test the modal canonical form""" # Create a system in the modal canonical form - A_true = np.diag([4.0, 3.0, 2.0, 1.0]) # order from the largest to the smallest - B_true = np.matrix("1.1 2.2 3.3 4.4").T - C_true = np.matrix("1.3 1.4 1.5 1.6") + A_true = np.diag([4.0, 3.0, 2.0, 1.0]) # order from largest to smallest + B_true = np.array([[1.1, 2.2, 3.3, 4.4]]).T + C_true = np.array([[1.3, 1.4, 1.5, 1.6]]) D_true = 42.0 # Perform a coordinate transform with a random invertible matrix - T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], - [-0.74855725, -0.39136285, -0.18142339, -0.50356997], - [-0.40688007, 0.81416369, 0.38002113, -0.16483334], - [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) - A = np.linalg.solve(T_true, A_true)*T_true + T_true = np.array( + [[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + [-0.74855725, -0.39136285, -0.18142339, -0.50356997], + [-0.40688007, 0.81416369, 0.38002113, -0.16483334], + [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) + A = np.dot(np.linalg.solve(T_true, A_true), T_true) B = np.linalg.solve(T_true, B_true) - C = C_true*T_true + C = np.dot(C_true, T_true) D = D_true # Create a state space system and convert it to the modal canonical form @@ -98,12 +101,14 @@ def test_modal_form(self): C_true = np.array([[1, 0, 0, 1]]) D_true = np.array([[0]]) - A = np.linalg.solve(T_true, A_true) * T_true + A = np.dot(np.linalg.solve(T_true, A_true), T_true) B = np.linalg.solve(T_true, B_true) - C = C_true * T_true + C = np.dot(C_true, T_true) D = D_true # Create state space system and convert to modal canonical form + print("C = ", C, "; shape = ", C.shape) + print("D = ", D, "; shape = ", D.shape) sys_check, T_check = canonical_form(ss(A, B, C, D), 'modal') # Check A and D matrix, which are uniquely defined @@ -124,13 +129,14 @@ def test_modal_form(self): np.dot(np.dot(C, matrix_power(A, i)), B)) # Reorder rows to get complete coverage (real eigenvalue cxrtvfirst) - A_true = np.array([[-1, 0, 0, 0], + # Use StateSpaceMatrix to allow * (variant for additional testing) + A_true = ct.StateSpaceMatrix([[-1, 0, 0, 0], [ 0, -2, 1, 0], [ 0, -1, -2, 0], [ 0, 0, 0, -3]]) - B_true = np.array([[0], [0], [1], [1]]) - C_true = np.array([[0, 1, 0, 1]]) - D_true = np.array([[0]]) + B_true = ct.StateSpaceMatrix([[0], [0], [1], [1]]) + C_true = ct.StateSpaceMatrix([[0, 1, 0, 1]]) + D_true = ct.StateSpaceMatrix([[0]]) A = np.linalg.solve(T_true, A_true) * T_true B = np.linalg.solve(T_true, B_true) @@ -168,18 +174,19 @@ def test_observable_form(self): coeffs = [1.0, 2.0, 3.0, 4.0, 1.0] A_true = np.polynomial.polynomial.polycompanion(coeffs) A_true = np.fliplr(np.flipud(A_true)) - B_true = np.matrix("1.0 1.0 1.0 1.0").T - C_true = np.matrix("1.0 0.0 0.0 0.0") + B_true = np.array([[1.0, 1.0, 1.0, 1.0]]).T + C_true = np.array([[1.0, 0.0, 0.0, 0.0]]) D_true = 42.0 # Perform a coordinate transform with a random invertible matrix - T_true = np.matrix([[-0.27144004, -0.39933167, 0.75634684, 0.44135471], - [-0.74855725, -0.39136285, -0.18142339, -0.50356997], - [-0.40688007, 0.81416369, 0.38002113, -0.16483334], - [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) - A = np.linalg.solve(T_true, A_true)*T_true + T_true = np.array( + [[-0.27144004, -0.39933167, 0.75634684, 0.44135471], + [-0.74855725, -0.39136285, -0.18142339, -0.50356997], + [-0.40688007, 0.81416369, 0.38002113, -0.16483334], + [-0.44769516, 0.15654653, -0.50060858, 0.72419146]]) + A = np.dot(np.linalg.solve(T_true, A_true), T_true) B = np.linalg.solve(T_true, B_true) - C = C_true*T_true + C = np.dot(C_true, T_true) D = D_true # Create a state space system and convert it to the observable canonical form @@ -201,9 +208,9 @@ def test_unobservable_system(self): """Test observable canonical form with an unobservable system""" # Create an unobservable system - A = np.matrix("1.0 2.0 2.0; 4.0 5.0 5.0; 7.0 8.0 8.0") - B = np.matrix("1.0 1.0 1.0").T - C = np.matrix("1.0 1.0 1.0") + A = np.array([[1.0, 2.0, 2.0], [4.0, 5.0, 5.0], [7.0, 8.0, 8.0]]) + B = np.array([[1.0, 1.0, 1.0]]).T + C = np.array([[1.0, 1.0, 1.0]]) D = 42.0 sys = ss(A, B, C, D) @@ -218,8 +225,9 @@ def test_arguments(self): np.testing.assert_raises( ControlNotImplemented, canonical_form, sys, 'unknown') + def suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestFeedback) + return unittest.TestLoader().loadTestsFromTestCase(TestCanonical) if __name__ == "__main__": diff --git a/control/tests/convert_test.py b/control/tests/convert_test.py index 0340fa718..b4b81c883 100644 --- a/control/tests/convert_test.py +++ b/control/tests/convert_test.py @@ -203,7 +203,7 @@ def testTf2ssStaticMimo(self): self.assertEqual(0, gmimo.states) self.assertEqual(3, gmimo.inputs) self.assertEqual(2, gmimo.outputs) - d = np.matrix([[0.5, 30, 0.0625], [-0.5, -1.25, 101.3]]) + d = np.array([[0.5, 30, 0.0625], [-0.5, -1.25, 101.3]]) np.testing.assert_array_equal(d, gmimo.D) def testSs2tfStaticSiso(self): @@ -220,7 +220,7 @@ def testSs2tfStaticMimo(self): a = [] b = [] c = [] - d = np.matrix([[0.5, 30, 0.0625], [-0.5, -1.25, 101.3]]) + d = np.array([[0.5, 30, 0.0625], [-0.5, -1.25, 101.3]]) gtf = control.ss2tf(control.ss(a,b,c,d)) # we need a 3x2x1 array to compare with gtf.num diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py index f08a5fa5e..966150815 100644 --- a/control/tests/discrete_test.py +++ b/control/tests/discrete_test.py @@ -345,13 +345,13 @@ def test_sample_system(self): def test_sample_ss(self): # double integrators, two different ways - sys1 = StateSpace([[0.,1.],[0.,0.]], [[0.],[1.]], [[1.,0.]], 0.) - sys2 = StateSpace([[0.,0.],[1.,0.]], [[1.],[0.]], [[0.,1.]], 0.) + sys1 = StateSpace([[0., 1.], [0., 0.]], [[0.], [1.]], [[1., 0.]], 0.) + sys2 = StateSpace([[0., 0.], [1., 0.]], [[1.], [0.]], [[0., 1.]], 0.) I = np.eye(2) for sys in (sys1, sys2): for h in (0.1, 0.5, 1, 2): Ad = I + h * sys.A - Bd = h * sys.B + 0.5 * h**2 * (sys.A * sys.B) + Bd = h * sys.B + 0.5 * h**2 * np.dot(sys.A, sys.B) sysd = sample_system(sys, h, method='zoh') np.testing.assert_array_almost_equal(sysd.A, Ad) np.testing.assert_array_almost_equal(sysd.B, Bd) diff --git a/control/tests/frd_test.py b/control/tests/frd_test.py index 1a6a263f3..f3dd34282 100644 --- a/control/tests/frd_test.py +++ b/control/tests/frd_test.py @@ -169,7 +169,7 @@ def testFeedback2(self): def testAuto(self): omega = np.logspace(-1, 2, 10) f1 = _convertToFRD(1, omega) - f2 = _convertToFRD(np.matrix([[1, 0], [0.1, -1]]), omega) + f2 = _convertToFRD(np.array([[1, 0], [0.1, -1]]), omega) f2 = _convertToFRD([[1, 0], [0.1, -1]], omega) f1, f2 # reference to avoid pyflakes error @@ -216,11 +216,11 @@ def testMIMOfb(self): @unittest.skipIf(not slycot_check(), "slycot not installed") def testMIMOfb2(self): - sys = StateSpace(np.matrix('-2.0 0 0; 0 -1 1; 0 0 -3'), - np.matrix('1.0 0; 0 0; 0 1'), + sys = StateSpace(np.array([[-2.0, 0, 0], [0, -1, 1], [0, 0, -3]]), + np.array([[1.0, 0], [0, 0], [0, 1]]), np.eye(3), np.zeros((3, 2))) omega = np.logspace(-1, 2, 10) - K = np.matrix('1 0.3 0; 0.1 0 0') + K = np.array([[1, 0.3, 0], [0.1, 0, 0]]) f1 = FRD(sys, omega).feedback(K) f2 = FRD(sys.feedback(K), omega) np.testing.assert_array_almost_equal( @@ -252,7 +252,7 @@ def testMIMOSmooth(self): [[1.0, 0.0], [0.0, 1.0]], [[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]], [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]) - sys2 = np.matrix([[1, 0, 0], [0, 1, 0]]) * sys + sys2 = np.array([[1, 0, 0], [0, 1, 0]]) * sys omega = np.logspace(-1, 2, 10) f1 = FRD(sys, omega, smooth=True) f2 = FRD(sys2, omega, smooth=True) @@ -271,15 +271,15 @@ def testAgainstOctave(self): # sys = ss([-2 0 0; 0 -1 1; 0 0 -3], # [1 0; 0 0; 0 1], eye(3), zeros(3,2)) # bfr = frd(bsys, [1]) - sys = StateSpace(np.matrix('-2.0 0 0; 0 -1 1; 0 0 -3'), - np.matrix('1.0 0; 0 0; 0 1'), + sys = StateSpace(np.array([[-2.0, 0, 0], [0, -1, 1], [0, 0, -3]]), + np.array([[1.0, 0], [0, 0], [0, 1]]), np.eye(3), np.zeros((3, 2))) omega = np.logspace(-1, 2, 10) f1 = FRD(sys, omega) np.testing.assert_array_almost_equal( (f1.freqresp([1.0])[0] * np.exp(1j*f1.freqresp([1.0])[1])).reshape(3, 2), - np.matrix('0.4-0.2j 0; 0 0.1-0.2j; 0 0.3-0.1j')) + np.array([[0.4-0.2j, 0], [0, 0.1-0.2j], [0, 0.3-0.1j]])) def test_string_representation(self): sys = FRD([1, 2, 3], [4, 5, 6]) diff --git a/control/tests/freqresp_test.py b/control/tests/freqresp_test.py index 05247b46b..d2affffd4 100644 --- a/control/tests/freqresp_test.py +++ b/control/tests/freqresp_test.py @@ -18,12 +18,12 @@ class TestFreqresp(unittest.TestCase): def setUp(self): - self.A = np.matrix('1,1;0,1') - self.C = np.matrix('1,0') + self.A = np.array([[1, 1], [0, 1]]) + self.C = np.array([[1, 0]]) self.omega = np.linspace(10e-2,10e2,1000) def test_siso(self): - B = np.matrix('0;1') + B = np.array([[0], [1]]) D = 0 sys = StateSpace(self.A,B,self.C,D) @@ -83,9 +83,9 @@ def test_superimpose(self): def test_doubleint(self): # 30 May 2016, RMM: added to replicate typecast bug in freqresp.py - A = np.matrix('0, 1; 0, 0'); - B = np.matrix('0; 1'); - C = np.matrix('1, 0'); + A = np.array([[0, 1], [0, 0]]); + B = np.array([[0], [1]]); + C = np.array([[1, 0]]); D = 0; sys = ss(A, B, C, D); bode(sys); @@ -93,8 +93,8 @@ def test_doubleint(self): @unittest.skipIf(not slycot_check(), "slycot not installed") def test_mimo(self): # MIMO - B = np.matrix('1,0;0,1') - D = np.matrix('0,0') + B = np.array([[1, 0], [0, 1]]) + D = np.array([[0, 0]]) sysMIMO = ss(self.A,B,self.C,D) frqMIMO = sysMIMO.freqresp(self.omega) diff --git a/control/tests/input_element_int_test.py b/control/tests/input_element_int_test.py index c6a6f64a3..30ee3f9c7 100644 --- a/control/tests/input_element_int_test.py +++ b/control/tests/input_element_int_test.py @@ -41,12 +41,12 @@ def test_tf_input_with_int_element_works(self): self.assertAlmostEqual(1.0, ctl.dcgain(sys)) def test_ss_input_with_int_element(self): - ident = np.matrix(np.identity(2), dtype=int) - a = np.matrix([[0, 1], - [-1, -2]], dtype=int) * ident - b = np.matrix([[0], - [1]], dtype=int) - c = np.matrix([[0, 1]], dtype=int) + ident = np.array(np.identity(2), dtype=int) + a = np.array([[0, 1], + [-1, -2]], dtype=int) * ident + b = np.array([[0], + [1]], dtype=int) + c = np.array([[0, 1]], dtype=int) d = 0 sys = ctl.ss(a, b, c, d) diff --git a/control/tests/mateqn_test.py b/control/tests/mateqn_test.py index 5eadcfefa..8fd0b6e38 100644 --- a/control/tests/mateqn_test.py +++ b/control/tests/mateqn_test.py @@ -43,11 +43,11 @@ """ import unittest -from numpy import matrix from numpy.testing import assert_array_almost_equal, assert_array_less # need scipy version of eigvals for generalized eigenvalue problem from scipy.linalg import eigvals, solve from scipy import zeros,dot +from control.statesp import ssmatrix from control.mateqn import lyap,dlyap,care,dare from control.exception import slycot_check @@ -56,81 +56,81 @@ class TestMatrixEquations(unittest.TestCase): """These are tests for the matrix equation solvers in mateqn.py""" def test_lyap(self): - A = matrix([[-1, 1],[-1, 0]]) - Q = matrix([[1,0],[0,1]]) + A = ssmatrix([[-1, 1], [-1, 0]]) + Q = ssmatrix([[1, 0], [0, 1]]) X = lyap(A,Q) # print("The solution obtained is ", X) assert_array_almost_equal(A * X + X * A.T + Q, zeros((2,2))) - A = matrix([[1, 2],[-3, -4]]) - Q = matrix([[3, 1],[1, 1]]) + A = ssmatrix([[1, 2], [-3, -4]]) + Q = ssmatrix([[3, 1], [1, 1]]) X = lyap(A,Q) # print("The solution obtained is ", X) assert_array_almost_equal(A * X + X * A.T + Q, zeros((2,2))) def test_lyap_sylvester(self): A = 5 - B = matrix([[4, 3], [4, 3]]) - C = matrix([2, 1]) + B = ssmatrix([[4, 3], [4, 3]]) + C = ssmatrix([2, 1]) X = lyap(A,B,C) # print("The solution obtained is ", X) assert_array_almost_equal(A * X + X * B + C, zeros((1,2))) - A = matrix([[2,1],[1,2]]) - B = matrix([[1,2],[0.5,0.1]]) - C = matrix([[1,0],[0,1]]) + A = ssmatrix([[2, 1], [1, 2]]) + B = ssmatrix([[1, 2], [0.5, 0.1]]) + C = ssmatrix([[1, 0], [0, 1]]) X = lyap(A,B,C) # print("The solution obtained is ", X) assert_array_almost_equal(A * X + X * B + C, zeros((2,2))) def test_lyap_g(self): - A = matrix([[-1, 2],[-3, -4]]) - Q = matrix([[3, 1],[1, 1]]) - E = matrix([[1,2],[2,1]]) + A = ssmatrix([[-1, 2], [-3, -4]]) + Q = ssmatrix([[3, 1], [1, 1]]) + E = ssmatrix([[1,2], [2,1]]) X = lyap(A,Q,None,E) # print("The solution obtained is ", X) assert_array_almost_equal(A * X * E.T + E * X * A.T + Q, zeros((2,2))) def test_dlyap(self): - A = matrix([[-0.6, 0],[-0.1, -0.4]]) - Q = matrix([[1,0],[0,1]]) + A = ssmatrix([[-0.6, 0], [-0.1, -0.4]]) + Q = ssmatrix([[1,0], [0,1]]) X = dlyap(A,Q) # print("The solution obtained is ", X) assert_array_almost_equal(A * X * A.T - X + Q, zeros((2,2))) - A = matrix([[-0.6, 0],[-0.1, -0.4]]) - Q = matrix([[3, 1],[1, 1]]) + A = ssmatrix([[-0.6, 0], [-0.1, -0.4]]) + Q = ssmatrix([[3, 1], [1, 1]]) X = dlyap(A,Q) # print("The solution obtained is ", X) assert_array_almost_equal(A * X * A.T - X + Q, zeros((2,2))) def test_dlyap_g(self): - A = matrix([[-0.6, 0],[-0.1, -0.4]]) - Q = matrix([[3, 1],[1, 1]]) - E = matrix([[1, 1],[2, 1]]) + A = ssmatrix([[-0.6, 0], [-0.1, -0.4]]) + Q = ssmatrix([[3, 1], [1, 1]]) + E = ssmatrix([[1, 1], [2, 1]]) X = dlyap(A,Q,None,E) # print("The solution obtained is ", X) assert_array_almost_equal(A * X * A.T - E * X * E.T + Q, zeros((2,2))) def test_dlyap_sylvester(self): A = 5 - B = matrix([[4, 3], [4, 3]]) - C = matrix([2, 1]) + B = ssmatrix([[4, 3], [4, 3]]) + C = ssmatrix([2, 1]) X = dlyap(A,B,C) # print("The solution obtained is ", X) assert_array_almost_equal(A * X * B.T - X + C, zeros((1,2))) - A = matrix([[2,1],[1,2]]) - B = matrix([[1,2],[0.5,0.1]]) - C = matrix([[1,0],[0,1]]) + A = ssmatrix([[2,1], [1,2]]) + B = ssmatrix([[1,2], [0.5,0.1]]) + C = ssmatrix([[1,0], [0,1]]) X = dlyap(A,B,C) # print("The solution obtained is ", X) assert_array_almost_equal(A * X * B.T - X + C, zeros((2,2))) def test_care(self): - A = matrix([[-2, -1],[-1, -1]]) - Q = matrix([[0, 0],[0, 1]]) - B = matrix([[1, 0],[0, 4]]) + A = ssmatrix([[-2, -1], [-1, -1]]) + Q = ssmatrix([[0, 0], [0, 1]]) + B = ssmatrix([[1, 0], [0, 4]]) X,L,G = care(A,B,Q) # print("The solution obtained is", X) @@ -139,12 +139,12 @@ def test_care(self): assert_array_almost_equal(B.T * X, G) def test_care_g(self): - A = matrix([[-2, -1],[-1, -1]]) - Q = matrix([[0, 0],[0, 1]]) - B = matrix([[1, 0],[0, 4]]) - R = matrix([[2, 0],[0, 1]]) - S = matrix([[0, 0],[0, 0]]) - E = matrix([[2, 1],[1, 2]]) + A = ssmatrix([[-2, -1], [-1, -1]]) + Q = ssmatrix([[0, 0], [0, 1]]) + B = ssmatrix([[1, 0], [0, 4]]) + R = ssmatrix([[2, 0], [0, 1]]) + S = ssmatrix([[0, 0], [0, 0]]) + E = ssmatrix([[2, 1], [1, 2]]) X,L,G = care(A,B,Q,R,S,E) # print("The solution obtained is", X) @@ -153,12 +153,12 @@ def test_care_g(self): (E.T * X * B + S) * solve(R, B.T * X * E + S.T) + Q, zeros((2,2))) assert_array_almost_equal(solve(R, B.T * X * E + S.T), G) - A = matrix([[-2, -1],[-1, -1]]) - Q = matrix([[0, 0],[0, 1]]) - B = matrix([[1],[0]]) + A = ssmatrix([[-2, -1], [-1, -1]]) + Q = ssmatrix([[0, 0], [0, 1]]) + B = ssmatrix([[1], [0]]) R = 1 - S = matrix([[1],[0]]) - E = matrix([[2, 1],[1, 2]]) + S = ssmatrix([[1], [0]]) + E = ssmatrix([[2, 1], [1, 2]]) X,L,G = care(A,B,Q,R,S,E) # print("The solution obtained is", X) @@ -168,10 +168,10 @@ def test_care_g(self): assert_array_almost_equal(dot( 1/R , dot(B.T,dot(X,E)) + S.T) , G) def test_dare(self): - A = matrix([[-0.6, 0],[-0.1, -0.4]]) - Q = matrix([[2, 1],[1, 0]]) - B = matrix([[2, 1],[0, 1]]) - R = matrix([[1, 0],[0, 1]]) + A = ssmatrix([[-0.6, 0], [-0.1, -0.4]]) + Q = ssmatrix([[2, 1], [1, 0]]) + B = ssmatrix([[2, 1], [0, 1]]) + R = ssmatrix([[1, 0], [0, 1]]) X,L,G = dare(A,B,Q,R) # print("The solution obtained is", X) @@ -183,28 +183,29 @@ def test_dare(self): lam = eigvals(A - B * G) assert_array_less(abs(lam), 1.0) - A = matrix([[1, 0],[-1, 1]]) - Q = matrix([[0, 1],[1, 1]]) - B = matrix([[1],[0]]) + A = ssmatrix([[1, 0], [-1, 1]]) + Q = ssmatrix([[0, 1], [1, 1]]) + B = ssmatrix([[1], [0]]) R = 2 X,L,G = dare(A,B,Q,R) # print("The solution obtained is", X) assert_array_almost_equal( - A.T * X * A - X - - A.T * X * B * solve(B.T * X * B + R, B.T * X * A) + Q, zeros((2,2))) + A.T * X * A - X + - A.T * X * B * solve(B.T * X * B + R, B.T * X * A) + + Q, zeros((2,2))) assert_array_almost_equal(B.T * X * A / (B.T * X * B + R), G) # check for stable closed loop lam = eigvals(A - B * G) assert_array_less(abs(lam), 1.0) def test_dare_g(self): - A = matrix([[-0.6, 0],[-0.1, -0.4]]) - Q = matrix([[2, 1],[1, 3]]) - B = matrix([[1, 5],[2, 4]]) - R = matrix([[1, 0],[0, 1]]) - S = matrix([[1, 0],[2, 0]]) - E = matrix([[2, 1],[1, 2]]) + A = ssmatrix([[-0.6, 0], [-0.1, -0.4]]) + Q = ssmatrix([[2, 1], [1, 3]]) + B = ssmatrix([[1, 5], [2, 4]]) + R = ssmatrix([[1, 0], [0, 1]]) + S = ssmatrix([[1, 0], [2, 0]]) + E = ssmatrix([[2, 1], [1, 2]]) X,L,G = dare(A,B,Q,R,S,E) # print("The solution obtained is", X) @@ -217,12 +218,12 @@ def test_dare_g(self): lam = eigvals(A - B * G, E) assert_array_less(abs(lam), 1.0) - A = matrix([[-0.6, 0],[-0.1, -0.4]]) - Q = matrix([[2, 1],[1, 3]]) - B = matrix([[1],[2]]) + A = ssmatrix([[-0.6, 0], [-0.1, -0.4]]) + Q = ssmatrix([[2, 1], [1, 3]]) + B = ssmatrix([[1], [2]]) R = 1 - S = matrix([[1],[2]]) - E = matrix([[2, 1],[1, 2]]) + S = ssmatrix([[1], [2]]) + E = ssmatrix([[2, 1], [1, 2]]) X,L,G = dare(A,B,Q,R,S,E) # print("The solution obtained is", X) diff --git a/control/tests/matlab_test.py b/control/tests/matlab_test.py index 0e7060bea..e13a13dd0 100644 --- a/control/tests/matlab_test.py +++ b/control/tests/matlab_test.py @@ -58,10 +58,10 @@ class TestMatlab(unittest.TestCase): def setUp(self): """Set up some systems for testing out MATLAB functions""" - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5.], [7.]]) + C = np.array([[6., 8.]]) + D = np.array([[9.]]) self.siso_ss1 = ss(A,B,C,D) # Create some transfer functions @@ -75,18 +75,18 @@ def setUp(self): self.siso_tf4 = ss2tf(self.siso_ss2); #Create MIMO system, contains ``siso_ss1`` twice - A = np.matrix("1. -2. 0. 0.;" - "3. -4. 0. 0.;" - "0. 0. 1. -2.;" - "0. 0. 3. -4. ") - B = np.matrix("5. 0.;" - "7. 0.;" - "0. 5.;" - "0. 7. ") - C = np.matrix("6. 8. 0. 0.;" - "0. 0. 6. 8. ") - D = np.matrix("9. 0.;" - "0. 9. ") + A = np.array([[1., -2., 0., 0.], + [3., -4., 0., 0.], + [0., 0., 1., -2.], + [0., 0., 3., -4.]]) + B = np.array([[5., 0.], + [7., 0.], + [0., 5.], + [0., 7.]]) + C = np.array([[6., 8., 0., 0.], + [0., 0., 6., 8.]]) + D = np.array([[9., 0.], + [0., 9.]]) self.mimo_ss1 = ss(A, B, C, D) # get consistent test results @@ -157,7 +157,7 @@ def testStep(self): np.testing.assert_array_almost_equal(yout, youttrue, decimal=4) np.testing.assert_array_almost_equal(tout, t) - X0 = np.array([0, 0]); + X0 = np.array([0, 0]) yout, tout = step(sys, T=t, X0=X0) np.testing.assert_array_almost_equal(yout, youttrue, decimal=4) np.testing.assert_array_almost_equal(tout, t) @@ -188,8 +188,8 @@ def testImpulse(self): warnings.simplefilter("ignore") #Test SISO system sys = self.siso_ss1 - youttrue = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766, 31.7344, - 26.1668, 21.6292, 17.9245, 14.8945]) + youttrue = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766, + 31.7344, 26.1668, 21.6292, 17.9245, 14.8945]) yout, tout = impulse(sys, T=t) np.testing.assert_array_almost_equal(yout, youttrue, decimal=4) np.testing.assert_array_almost_equal(tout, t) @@ -220,7 +220,7 @@ def testInitial(self): #Test SISO system sys = self.siso_ss1 t = np.linspace(0, 1, 10) - x0 = np.matrix(".5; 1.") + x0 = np.array([.5, 1.]) youttrue = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092, 1.1508, 0.5833, 0.1645, -0.1391]) yout, tout = initial(sys, T=t, X0=x0) @@ -235,7 +235,7 @@ def testInitial(self): if slycot_check(): #Test MIMO system, which contains ``siso_ss1`` twice sys = self.mimo_ss1 - x0 = np.matrix(".5; 1.; .5; 1.") + x0 = np.array([[.5], [1.], [.5], [1.]]) y_00, _t = initial(sys, T=t, X0=x0, input=0, output=0) y_11, _t = initial(sys, T=t, X0=x0, input=1, output=1) np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4) @@ -257,7 +257,7 @@ def testLsim(self): #test with initial value and special algorithm for ``U=0`` u=0 - x0 = np.matrix(".5; 1.") + x0 = np.array([[.5], [1.]]) youttrue = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092, 1.1508, 0.5833, 0.1645, -0.1391]) yout, _t, _xout = lsim(self.siso_ss1, u, t, x0) @@ -268,7 +268,7 @@ def testLsim(self): #first system: initial value, second system: step response u = np.array([[0., 1.], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1]]) - x0 = np.matrix(".5; 1; 0; 0") + x0 = np.array([[.5], [1], [0], [0]]) youttrue = np.array([[11., 9.], [8.1494, 17.6457], [5.9361, 24.7072], [4.2258, 30.4855], [2.9118, 35.2234], [1.9092, 39.1165], @@ -494,14 +494,16 @@ def testSISOtfdata(self): np.testing.assert_array_almost_equal(tfdata_1[i], tfdata_2[i]) def testDamp(self): - A = np.mat('''-0.2 0.06 0 -1; - 0 0 1 0; - -17 0 -3.8 1; - 9.4 0 -0.4 -0.6''') - B = np.mat('''-0.01 0.06; - 0 0; - -32 5.4; - 2.6 -7''') + A = np.array([ + [-0.2, 0.06, 0, -1], + [ 0, 0, 1, 0], + [ -17, 0, -3.8, 1], + [ 9.4, 0, -0.4, -0.6]]) + B = np.array([ + [-0.01, 0.06], + [ 0, 0], + [ -32, 5.4], + [ 2.6, -7]]) C = np.eye(4) D = np.zeros((4,2)) sys = ss(A, B, C, D) @@ -514,20 +516,20 @@ def testDamp(self): Z, np.array([1.0, 0.07983139, 0.07983139, 1.0])) def testConnect(self): - sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.") - sys2 = ss("-1.", "1.", "1.", "0.") + sys1 = ss([[1., -2], [3., -4]], [[5.], [7]], [[6, 8]], [[9.]]) + sys2 = ss([[-1.]], [[1.]], [[1.]], [[0.]]) sys = append(sys1, sys2) - Q= np.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1 + Q= np.array([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1 sysc = connect(sys, Q, [2], [1, 2]) # print(sysc) np.testing.assert_array_almost_equal( - sysc.A, np.mat('1 -2 5; 3 -4 7; -6 -8 -10')) + sysc.A, np.array([[1, -2, 5], [3, -4, 7], [-6, -8, -10]])) np.testing.assert_array_almost_equal( - sysc.B, np.mat('0; 0; 1')) + sysc.B, np.array([[0], [0], [1]])) np.testing.assert_array_almost_equal( - sysc.C, np.mat('6 8 9; 0 0 1')) + sysc.C, np.array([[6, 8, 9], [0, 0, 1]])) np.testing.assert_array_almost_equal( - sysc.D, np.mat('0; 0')) + sysc.D, np.array([[0], [0]])) def testConnect2(self): sys = append(ss([[-5, -2.25], [4, 0]], [[2], [0]], @@ -538,18 +540,18 @@ def testConnect2(self): Q = [ [ 1, 3], [2, 1], [3, -2]] sysc = connect(sys, Q, [3], [3, 1, 2]) np.testing.assert_array_almost_equal( - sysc.A, np.mat([[-5, -2.25, 0, -6.6666], + sysc.A, np.array([[-5, -2.25, 0, -6.6666], [4, 0, 0, 0], [0, 2.25, -1.6667, 0], [0, 0, 1, 0]])) np.testing.assert_array_almost_equal( - sysc.B, np.mat([[2], [0], [0], [0]])) + sysc.B, np.array([[2], [0], [0], [0]])) np.testing.assert_array_almost_equal( - sysc.C, np.mat([[0, 0, 0, -3.3333], + sysc.C, np.array([[0, 0, 0, -3.3333], [0, 1.125, 0, 0], [0, 0, 0, 3.3333]])) np.testing.assert_array_almost_equal( - sysc.D, np.mat([[1], [0], [0]])) + sysc.D, np.array([[1], [0], [0]])) @@ -591,21 +593,22 @@ def testMinreal(self, verbose=False): def testSS2cont(self): sys = ss( - np.mat("-3 4 2; -1 -3 0; 2 5 3"), - np.mat("1 4 ; -3 -3; -2 1"), - np.mat("4 2 -3; 1 4 3"), - np.mat("-2 4; 0 1")) + np.array([[-3, 4, 2], [-1, -3, 0], [2, 5, 3]]), + np.array([[1, 4], [-3, -3], [-2, 1]]), + np.array([[4, 2, -3], [1, 4, 3]]), + np.array([[-2, 4], [0, 1]])) sysd = c2d(sys, 0.1) np.testing.assert_array_almost_equal( - np.mat( - """0.742840837331905 0.342242024293711 0.203124211149560; - -0.074130792143890 0.724553295044645 -0.009143771143630; - 0.180264783290485 0.544385612448419 1.370501013067845"""), + np.array([ + [0.742840837331905, 0.342242024293711, 0.203124211149560], + [-0.074130792143890, 0.724553295044645, -0.009143771143630], + [ 0.180264783290485, 0.544385612448419, 1.370501013067845]]), sysd.A) np.testing.assert_array_almost_equal( - np.mat(""" 0.012362066084719 0.301932197918268; - -0.260952977031384 -0.274201791021713; - -0.304617775734327 0.075182622718853"""), sysd.B) + np.array([ + [0.012362066084719, 0.301932197918268], + [-0.260952977031384, -0.274201791021713], + [-0.304617775734327, 0.075182622718853]]), sysd.B) def testCombi01(self): # test from a "real" case, combines tf, ss, connect and margin @@ -652,7 +655,7 @@ def testCombi01(self): # start with the basic satellite model sat1, and get the # payload attitude response - Hp = tf(sp.matrix([0, 0, 0, 1])*sat1) + Hp = tf(sp.array([0, 0, 0, 1])*sat1) # total open loop Hol = Hc*Hno*Hp diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py index f4bc0fd98..a30cb58b3 100644 --- a/control/tests/modelsimp_test.py +++ b/control/tests/modelsimp_test.py @@ -5,6 +5,7 @@ import unittest import numpy as np +import warnings from control.modelsimp import * from control.matlab import * from control.exception import slycot_check @@ -12,39 +13,50 @@ class TestModelsimp(unittest.TestCase): @unittest.skipIf(not slycot_check(), "slycot not installed") def testHSVD(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5.], [7.]]) + C = np.array([[6., 8.]]) + D = np.array([[9.]]) sys = ss(A,B,C,D) - hsv = hsvd(sys) - hsvtrue = np.matrix("24.42686 0.5731395") # from MATLAB + hsv = hsvd(sys, return_type=np.ndarray) + hsvtrue = np.array([[24.42686, 0.5731395]]) # from MATLAB np.testing.assert_array_almost_equal(hsv, hsvtrue) + # Make sure default type values are correct + self.assertTrue(isinstance(hsv, np.ndarray)) + + # Check that default type generates a warning + # TODO: remove this check with matrix type is deprecated + with warnings.catch_warnings(record=True) as w: + hsv = hsvd(sys) + self.assertTrue(issubclass(w[-1].category, UserWarning)) + self.assertTrue(isinstance(hsv, np.ndarray)) + def testMarkov(self): - U = np.matrix("1.; 1.; 1.; 1.; 1.") + U = np.array([[1.], [1.], [1.], [1.], [1.]]) Y = U M = 3 H = markov(Y,U,M) - Htrue = np.matrix("1.; 0.; 0.") + Htrue = np.array([[1.], [0.], [0.]]) np.testing.assert_array_almost_equal( H, Htrue ) def testModredMatchDC(self): #balanced realization computed in matlab for the transfer function: # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-1.958, -1.194, 1.824, -1.464; \ - -1.194, -0.8344, 2.563, -1.351; \ - -1.824, -2.563, -1.124, 2.704; \ - -1.464, -1.351, -2.704, -11.08') - B = np.matrix('-0.9057; -0.4068; -0.3263; -0.3474') - C = np.matrix('-0.9057, -0.4068, 0.3263, -0.3474') - D = np.matrix('0.') + A = np.array( + [[-1.958, -1.194, 1.824, -1.464], + [-1.194, -0.8344, 2.563, -1.351], + [-1.824, -2.563, -1.124, 2.704], + [-1.464, -1.351, -2.704, -11.08]]) + B = np.array([[-0.9057], [-0.4068], [-0.3263], [-0.3474]]) + C = np.array([[-0.9057, -0.4068, 0.3263, -0.3474]]) + D = np.array([[0.]]) sys = ss(A,B,C,D) rsys = modred(sys,[2, 3],'matchdc') - Artrue = np.matrix('-4.431, -4.552; -4.552, -5.361') - Brtrue = np.matrix('-1.362; -1.031') - Crtrue = np.matrix('-1.362, -1.031') - Drtrue = np.matrix('-0.08384') + Artrue = np.array([[-4.431, -4.552], [-4.552, -5.361]]) + Brtrue = np.array([[-1.362], [-1.031]]) + Crtrue = np.array([[-1.362, -1.031]]) + Drtrue = np.array([[-0.08384]]) np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=3) np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=3) np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=3) @@ -52,32 +64,34 @@ def testModredMatchDC(self): def testModredUnstable(self): # Check if an error is thrown when an unstable system is given - A = np.matrix('4.5418, 3.3999, 5.0342, 4.3808; \ - 0.3890, 0.3599, 0.4195, 0.1760; \ - -4.2117, -3.2395, -4.6760, -4.2180; \ - 0.0052, 0.0429, 0.0155, 0.2743') - B = np.matrix('1.0, 1.0; 2.0, 2.0; 3.0, 3.0; 4.0, 4.0') - C = np.matrix('1.0, 2.0, 3.0, 4.0; 1.0, 2.0, 3.0, 4.0') - D = np.matrix('0.0, 0.0; 0.0, 0.0') + A = np.array( + [[4.5418, 3.3999, 5.0342, 4.3808], + [0.3890, 0.3599, 0.4195, 0.1760], + [-4.2117, -3.2395, -4.6760, -4.2180], + [0.0052, 0.0429, 0.0155, 0.2743]]) + B = np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0], [4.0, 4.0]]) + C = np.array([[1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0]]) + D = np.array([[0.0, 0.0], [0.0, 0.0]]) sys = ss(A,B,C,D) np.testing.assert_raises(ValueError, modred, sys, [2, 3]) def testModredTruncate(self): #balanced realization computed in matlab for the transfer function: # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-1.958, -1.194, 1.824, -1.464; \ - -1.194, -0.8344, 2.563, -1.351; \ - -1.824, -2.563, -1.124, 2.704; \ - -1.464, -1.351, -2.704, -11.08') - B = np.matrix('-0.9057; -0.4068; -0.3263; -0.3474') - C = np.matrix('-0.9057, -0.4068, 0.3263, -0.3474') - D = np.matrix('0.') + A = np.array( + [[-1.958, -1.194, 1.824, -1.464], + [-1.194, -0.8344, 2.563, -1.351], + [-1.824, -2.563, -1.124, 2.704], + [-1.464, -1.351, -2.704, -11.08]]) + B = np.array([[-0.9057], [-0.4068], [-0.3263], [-0.3474]]) + C = np.array([[-0.9057, -0.4068, 0.3263, -0.3474]]) + D = np.array([[0.]]) sys = ss(A,B,C,D) rsys = modred(sys,[2, 3],'truncate') - Artrue = np.matrix('-1.958, -1.194; -1.194, -0.8344') - Brtrue = np.matrix('-0.9057; -0.4068') - Crtrue = np.matrix('-0.9057, -0.4068') - Drtrue = np.matrix('0.') + Artrue = np.array([[-1.958, -1.194], [-1.194, -0.8344]]) + Brtrue = np.array([[-0.9057], [-0.4068]]) + Crtrue = np.array([[-0.9057, -0.4068]]) + Drtrue = np.array([[0.]]) np.testing.assert_array_almost_equal(rsys.A, Artrue) np.testing.assert_array_almost_equal(rsys.B, Brtrue) np.testing.assert_array_almost_equal(rsys.C, Crtrue) @@ -88,20 +102,21 @@ def testModredTruncate(self): def testBalredTruncate(self): #controlable canonical realization computed in matlab for the transfer function: # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-15., -7.5, -6.25, -1.875; \ - 8., 0., 0., 0.; \ - 0., 4., 0., 0.; \ - 0., 0., 1., 0.') - B = np.matrix('2.; 0.; 0.; 0.') - C = np.matrix('0.5, 0.6875, 0.7031, 0.5') - D = np.matrix('0.') + A = np.array( + [[-15., -7.5, -6.25, -1.875], + [8., 0., 0., 0.], + [0., 4., 0., 0.], + [0., 0., 1., 0.]]) + B = np.array([[2.], [0.], [0.], [0.]]) + C = np.array([[0.5, 0.6875, 0.7031, 0.5]]) + D = np.array([[0.]]) sys = ss(A,B,C,D) orders = 2 rsys = balred(sys,orders,method='truncate') - Artrue = np.matrix('-1.958, -1.194; -1.194, -0.8344') - Brtrue = np.matrix('0.9057; 0.4068') - Crtrue = np.matrix('0.9057, 0.4068') - Drtrue = np.matrix('0.') + Artrue = np.array([[-1.958, -1.194], [-1.194, -0.8344]]) + Brtrue = np.array([[0.9057], [0.4068]]) + Crtrue = np.array([[0.9057, 0.4068]]) + Drtrue = np.array([[0.]]) np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=2) np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=4) np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4) @@ -111,20 +126,23 @@ def testBalredTruncate(self): def testBalredMatchDC(self): #controlable canonical realization computed in matlab for the transfer function: # num = [1 11 45 32], den = [1 15 60 200 60] - A = np.matrix('-15., -7.5, -6.25, -1.875; \ - 8., 0., 0., 0.; \ - 0., 4., 0., 0.; \ - 0., 0., 1., 0.') - B = np.matrix('2.; 0.; 0.; 0.') - C = np.matrix('0.5, 0.6875, 0.7031, 0.5') - D = np.matrix('0.') + A = np.array( + [[-15., -7.5, -6.25, -1.875], + [8., 0., 0., 0.], + [0., 4., 0., 0.], + [0., 0., 1., 0.]]) + B = np.array([[2.], [0.], [0.], [0.]]) + C = np.array([[0.5, 0.6875, 0.7031, 0.5]]) + D = np.array([[0.]]) sys = ss(A,B,C,D) orders = 2 rsys = balred(sys,orders,method='matchdc') - Artrue = np.matrix('-4.43094773, -4.55232904; -4.55232904, -5.36195206') - Brtrue = np.matrix('1.36235673; 1.03114388') - Crtrue = np.matrix('1.36235673, 1.03114388') - Drtrue = np.matrix('-0.08383902') + Artrue = np.array( + [[-4.43094773, -4.55232904], + [-4.55232904, -5.36195206]]) + Brtrue = np.array([[1.36235673], [1.03114388]]) + Crtrue = np.array([[1.36235673, 1.03114388]]) + Drtrue = np.array([[-0.08383902]]) np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=2) np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=4) np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4) diff --git a/control/tests/ssmatrix_test.py b/control/tests/ssmatrix_test.py new file mode 100644 index 000000000..02b284ce4 --- /dev/null +++ b/control/tests/ssmatrix_test.py @@ -0,0 +1,243 @@ +#!/usr/bin/env python +# +# ssmatrix_test.py - test state-space matrix class +# RMM, 13 Apr 2019 + +import unittest +import numpy as np +import control as ct + + +class TestStateSpaceMatrix(unittest.TestCase): + """Tests for the StateSpaceMatrix class.""" + + def setUp(self): + pass # Do nothing for now + + def test_constructor(self): + # Create a matrix and make sure we get back a 2D array + M = ct.StateSpaceMatrix([[1, 1], [1, 1]]) + self.assertEqual(M.shape, (2, 2)) + np.testing.assert_array_equal(M, np.array([[1, 1], [1, 1]])) + + # Passing a vector should give back a row vector (as 2D matrix) + M = ct.StateSpaceMatrix([1, 1]) + self.assertEqual(M.shape, (1, 2)) + np.testing.assert_array_equal(M, np.array([[1, 1]])) + + # Use axis to switch to a column vector + M = ct.StateSpaceMatrix([1, 1], axis=0) + self.assertEqual(M.shape, (2, 1)) + np.testing.assert_array_equal(M, np.array([[1], [1]])) + + # Scalars should get converted to 1x1 matrix + M = ct.StateSpaceMatrix(1) + self.assertEqual(M.shape, (1, 1)) + np.testing.assert_array_equal(M, np.array([[1]])) + + # Empty matrix should have shape (0, 0) + M = ct.StateSpaceMatrix([[]]) + self.assertEqual(M.shape, (0, 0)) + + # Use (deprecated?) matrix-style construction string (w/ warnings off) + import warnings + warnings.filterwarnings("ignore") + M = ct.StateSpaceMatrix("1, 1; 1, 1") + warnings.filterwarnings("default") + self.assertEqual(M.shape, (2, 2)) + self.assertTrue(isinstance(M, ct.StateSpaceMatrix)) + + def test_mul(self): + # Make sure that multiplying two matrices gives a matrix + M1 = ct.StateSpaceMatrix([[1, 1], [1, -1]]) + M2 = ct.StateSpaceMatrix([[1, 2], [2, 1]]) + Mprod = M1 * M2 + self.assertTrue(isinstance(Mprod, ct.StateSpaceMatrix)) + self.assertEqual(Mprod.shape, (2, 2)) + np.testing.assert_array_equal(Mprod, np.dot(M1, M2)) + + # Matrix times a (state-space) column vector gives a column vector + Cm = ct.StateSpaceMatrix([[1], [2]]) + MCm = M1 * Cm + self.assertTrue(isinstance(MCm, ct.StateSpaceMatrix)) + self.assertEqual(MCm.shape, (2, 1)) + np.testing.assert_array_equal(MCm, np.dot(M1, Cm)) + + # Matrix times a (ndarray) column vector gives a column vector + Ca = np.array([[1], [2]]) + MCa = M1 * Ca + self.assertTrue(isinstance(MCa, ct.StateSpaceMatrix)) + self.assertEqual(MCa.shape, (2, 1)) + np.testing.assert_array_equal(MCa, np.dot(M1, Ca)) + + # (State-space) row vector time SSMatrix gives SSMatrix row vector + Rm = ct.StateSpaceMatrix([[1, 2]]) + MRm = Rm * M1 + self.assertTrue(isinstance(MRm, ct.StateSpaceMatrix)) + self.assertEqual(MRm.shape, (1, 2)) + np.testing.assert_array_equal(MRm, np.dot(Rm, M1)) + + # (ndarray) row vector time SSMatrix gives SSMatrix row vector + Rm = np.array([[1, 2]]) + MRm = Rm * M1 + self.assertTrue(isinstance(MRm, ct.StateSpaceMatrix)) + self.assertEqual(MRm.shape, (1, 2)) + np.testing.assert_array_equal(MRm, np.dot(Rm, M1)) + + # Row vector times column vector gives 1x1 matrix + RmCm = Rm * Cm + self.assertTrue(isinstance(RmCm, ct.StateSpaceMatrix)) + self.assertEqual(RmCm.shape, (1, 1)) + np.testing.assert_array_equal(RmCm, np.dot(Rm, Cm)) + + # Column vector times row vector gives nxn matrix + CmRm = Cm * Rm + self.assertTrue(isinstance(CmRm, ct.StateSpaceMatrix)) + self.assertEqual(CmRm.shape, (2, 2)) + np.testing.assert_array_equal(CmRm, np.dot(Cm, Rm)) + + def test_multiply_other(self): + """Make sure that certain operations preserve StateSpaceMatrix type""" + M = ct.StateSpaceMatrix([[1, 1], [1, 1]]) + Mint = M * 5 + self.assertTrue(isinstance(Mint, ct.StateSpaceMatrix)) + + Mint = 5 * M + self.assertTrue(isinstance(Mint, ct.StateSpaceMatrix)) + + Mreal = M * 5.5 + self.assertTrue(isinstance(Mreal, ct.StateSpaceMatrix)) + + Mreal = 5.5 * M + self.assertTrue(isinstance(Mreal, ct.StateSpaceMatrix)) + + Mcomplex = M * 1j + self.assertFalse(isinstance(Mcomplex, ct.StateSpaceMatrix)) + + Mcomplex = 1j * M + self.assertFalse(isinstance(Mcomplex, ct.StateSpaceMatrix)) + + Mreal = M * np.array([[5.5, 0], [0, 5.5]]) + self.assertTrue(isinstance(Mreal, ct.StateSpaceMatrix)) + + Mreal = np.array([[5.5, 0], [0, 5.5]]) * M + self.assertTrue(isinstance(Mreal, ct.StateSpaceMatrix)) + + Mcomplex = M * np.array([[1j, 0], [0, 1j]]) + self.assertFalse(isinstance(Mcomplex, ct.StateSpaceMatrix)) + + Mcomplex = np.array([[1j, 0], [0, 1j]]) * M + self.assertFalse(isinstance(Mcomplex, ct.StateSpaceMatrix)) + + def test_imul(self): + # Make sure that multiplying two matrices gives a matrix + M1 = ct.StateSpaceMatrix([[1, 1], [1, -1]]) + M2 = ct.StateSpaceMatrix([[1, 2], [2, 1]]) + Mprod = M1; Mprod *= M2 + self.assertTrue(isinstance(Mprod, ct.StateSpaceMatrix)) + self.assertEqual(Mprod.shape, (2, 2)) + np.testing.assert_array_equal(Mprod, np.dot(M1, M2)) + + # Matrix times a (state-space) column vector gives a column vector + Cm = ct.StateSpaceMatrix([[1], [2]]) + MCm = M1; MCm *= Cm + self.assertTrue(isinstance(MCm, ct.StateSpaceMatrix)) + self.assertEqual(MCm.shape, (2, 1)) + np.testing.assert_array_equal(MCm, np.dot(M1, Cm)) + + # Matrix times a (ndarray) column vector gives a column vector + Ca = np.array([[1], [2]]) + MCa = M1; MCa *= Ca + self.assertTrue(isinstance(MCa, ct.StateSpaceMatrix)) + self.assertEqual(MCa.shape, (2, 1)) + np.testing.assert_array_equal(MCa, np.dot(M1, Ca)) + + def test_pow(self): + """Test matrix exponential""" + # Generate cases that should work + M = ct.ssmatrix([[1, 1], [-1, 2]]) + np.testing.assert_array_almost_equal(M**0, np.eye(2)) + np.testing.assert_array_almost_equal(M**1, M) + np.testing.assert_array_almost_equal(M**2, np.dot(M, M)) + np.testing.assert_array_almost_equal(M**3, M*M*M) + + # Make sure that we get errors if we do something wrong + self.assertRaises(TypeError, lambda: M**0.5) + self.assertRaises(np.linalg.LinAlgError, lambda: M[0, :]**2) + + def test_ipow(self): + """Test matrix exponential""" + # Generate cases that should work + M = ct.ssmatrix([[1, 1], [-1, 2]]) + M0 = M; M0 **= 0 + M1 = M; M1 **= 1 + M2 = M; M2 **= 2 + M3 = M; M3 **= 3 + np.testing.assert_array_almost_equal(M0, np.eye(2)) + np.testing.assert_array_almost_equal(M1, M) + np.testing.assert_array_almost_equal(M2, np.dot(M, M)) + np.testing.assert_array_almost_equal(M3, M*M*M) + + def test_getitem(self): + M = ct.StateSpaceMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + + # Extracting full slice gives back what we started with + S1 = M[:, :] + self.assertTrue(isinstance(S1, ct.StateSpaceMatrix)) + self.assertEqual(S1.shape, M.shape) + np.testing.assert_array_equal(S1, M) + + # Extracting multiple columns using slice + S2 = M[:, 0:2] + self.assertTrue(isinstance(S2, ct.StateSpaceMatrix)) + self.assertEqual(S2.shape, (3, 2)) + np.testing.assert_array_equal(S2, [[1, 2], [4, 5], [7, 8]]) + + # Extracting multiple columns using array + S3 = M[:, [0, 1]] + self.assertTrue(isinstance(S3, ct.StateSpaceMatrix)) + self.assertEqual(S3.shape, (3, 2)) + np.testing.assert_array_equal(S3, [[1, 2], [4, 5], [7, 8]]) + + # Extracting single column returns column matrix + S4 = M[:, 1] + self.assertTrue(isinstance(S4, ct.StateSpaceMatrix)) + self.assertEqual(S4.shape, (3, 1)) + np.testing.assert_array_equal(S4, [[2], [5], [8]]) + + # Extracting multiple rows using slice + S5 = M[0:2, :] + self.assertTrue(isinstance(S5, ct.StateSpaceMatrix)) + self.assertEqual(S5.shape, (2, 3)) + np.testing.assert_array_equal(S5, [[1, 2, 3], [4, 5, 6]]) + + # Extracting multiple rows using array + S6 = M[[0, 1], :] + self.assertTrue(isinstance(S6, ct.StateSpaceMatrix)) + self.assertEqual(S6.shape, (2, 3)) + np.testing.assert_array_equal(S6, [[1, 2, 3], [4, 5, 6]]) + + # Extracting single row returns row matrix + S6 = M[1, :] + self.assertTrue(isinstance(S6, ct.StateSpaceMatrix)) + self.assertEqual(S6.shape, (1, 3)) + np.testing.assert_array_equal(S6, [[4, 5, 6]]) + + # Extracting row and column slices returns matrix + S7 = M[0:2, 0:2] + self.assertTrue(isinstance(S7, ct.StateSpaceMatrix)) + self.assertEqual(S7.shape, (2, 2)) + np.testing.assert_array_equal(S7, [[1, 2], [4, 5]]) + + # Extracting single row and column returns scalar + S8 = M[1, 1] + self.assertTrue(np.isscalar(S8)) + self.assertEqual(S8, 5) + + +def suite(): + return unittest.TestLoader().loadTestsFromTestCase(TestStateSpaceMatrix) + + +if __name__ == "__main__": + unittest.main() diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py index 66dce2b12..515a407ef 100644 --- a/control/tests/statefbk_test.py +++ b/control/tests/statefbk_test.py @@ -5,11 +5,15 @@ from __future__ import print_function import unittest +import sys as pysys import numpy as np +import warnings from control.statefbk import ctrb, obsv, place, place_varga, lqr, gram, acker from control.matlab import * from control.exception import slycot_check, ControlDimension from control.mateqn import care, dare +from control.config import use_numpy_matrix +from control.statesp import StateSpaceMatrix class TestStatefbk(unittest.TestCase): """Test state feedback functions""" @@ -25,103 +29,164 @@ def setUp(self): np.random.seed(0) def testCtrbSISO(self): - A = np.matrix("1. 2.; 3. 4.") - B = np.matrix("5.; 7.") - Wctrue = np.matrix("5. 19.; 7. 43.") - Wc = ctrb(A,B) + A = np.array([[1., 2.], [3., 4.]]) + B = np.array([[5.], [7.]]) + Wctrue = np.array([[5., 19.], [7., 43.]]) + + use_numpy_matrix(False) + Wc = ctrb(A, B) np.testing.assert_array_almost_equal(Wc, Wctrue) + self.assertTrue(isinstance(Wc, StateSpaceMatrix)) + + # This test only works in Python 3 due to a conflict with the same + # warning type in other test modules (frd_test.py). See + # https://bugs.python.org/issue4180 for more details + @unittest.skipIf(pysys.version_info < (3, 0), "test requires Python 3+") + def test_ctrb_siso_deprecated(self): + A = np.array([[1., 2.], [3., 4.]]) + B = np.array([[5.], [7.]]) + + # Check that default using np.matrix generates a warning + # TODO: remove this check with matrix type is deprecated + use_numpy_matrix(True) + with warnings.catch_warnings(record=True) as w: + Wc = ctrb(A, B) + self.assertTrue(issubclass(w[-1].category, UserWarning)) + self.assertTrue(isinstance(Wc, np.ndarray)) def testCtrbMIMO(self): - A = np.matrix("1. 2.; 3. 4.") - B = np.matrix("5. 6.; 7. 8.") - Wctrue = np.matrix("5. 6. 19. 22.; 7. 8. 43. 50.") - Wc = ctrb(A,B) + A = np.array([[1., 2.], [3., 4.]]) + B = np.array([[5., 6.], [7., 8.]]) + Wctrue = np.array([[5., 6., 19., 22.], [7., 8., 43., 50.]]) + Wc = ctrb(A, B, return_type=np.ndarray) np.testing.assert_array_almost_equal(Wc, Wctrue) + # Make sure default type values are correct + self.assertTrue(isinstance(Wc, np.ndarray)) + def testObsvSISO(self): - A = np.matrix("1. 2.; 3. 4.") - C = np.matrix("5. 7.") - Wotrue = np.matrix("5. 7.; 26. 38.") - Wo = obsv(A,C) + A = np.array([[1., 2.], [3., 4.]]) + C = np.array([[5., 7.]]) + Wotrue = np.array([[5., 7.], [26., 38.]]) + use_numpy_matrix(True) # set default, but override in next statement + Wo = obsv(A, C, return_type=np.ndarray) np.testing.assert_array_almost_equal(Wo, Wotrue) + # Make sure default type values are correct + self.assertTrue(isinstance(Wo, np.ndarray)) + + # This test only works in Python 3 due to a conflict with the same + # warning type in other test modules (frd_test.py). See + # https://bugs.python.org/issue4180 for more details + @unittest.skipIf(pysys.version_info < (3, 0), "test requires Python 3+") + def test_obsv_siso_deprecated(self): + A = np.array([[1., 2.], [3., 4.]]) + C = np.array([[5., 7.]]) + + # Check that default type generates a warning + # TODO: remove this check with matrix type is deprecated + use_numpy_matrix(True) + with warnings.catch_warnings(record=True) as w: + Wo = obsv(A, C) + self.assertTrue(issubclass(w[-1].category, UserWarning)) + self.assertTrue(isinstance(Wo, np.ndarray)) + def testObsvMIMO(self): - A = np.matrix("1. 2.; 3. 4.") - C = np.matrix("5. 6.; 7. 8.") - Wotrue = np.matrix("5. 6.; 7. 8.; 23. 34.; 31. 46.") - Wo = obsv(A,C) + A = np.array([[1., 2.], [3., 4.]]) + C = np.array([[5., 6.], [7., 8.]]) + Wotrue = np.array([[5., 6.], [7., 8.], [23., 34.], [31., 46.]]) + Wo = obsv(A, C, return_type=np.ndarray) np.testing.assert_array_almost_equal(Wo, Wotrue) def testCtrbObsvDuality(self): - A = np.matrix("1.2 -2.3; 3.4 -4.5") - B = np.matrix("5.8 6.9; 8. 9.1") - Wc = ctrb(A,B); + A = np.array([[1.2, -2.3], [3.4, -4.5]]) + B = np.array([[5.8, 6.9], [8., 9.1]]) + Wc = ctrb(A, B, return_type=np.ndarray); A = np.transpose(A) C = np.transpose(B) - Wo = np.transpose(obsv(A,C)); + Wo = np.transpose(obsv(A, C, return_type=np.ndarray)); np.testing.assert_array_almost_equal(Wc,Wo) @unittest.skipIf(not slycot_check(), "slycot not installed") def testGramWc(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5. 6.; 7. 8.") - C = np.matrix("4. 5.; 6. 7.") - D = np.matrix("13. 14.; 15. 16.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5., 6.], [7., 8.]]) + C = np.array([[4., 5.], [6., 7.]]) + D = np.array([[13., 14.], [15., 16.]]) sys = ss(A, B, C, D) - Wctrue = np.matrix("18.5 24.5; 24.5 32.5") - Wc = gram(sys,'c') + Wctrue = np.array([[18.5, 24.5], [24.5, 32.5]]) + Wc = gram(sys, 'c', return_type=np.ndarray) np.testing.assert_array_almost_equal(Wc, Wctrue) + # This test only works in Python 3 due to a conflict with the same + # warning type in other test modules (frd_test.py). See + # https://bugs.python.org/issue4180 for more details + @unittest.skipIf(pysys.version_info < (3, 0) or not slycot_check(), + "test requires Python 3+ and slycot") + def test_gram_wc_deprecated(self): + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5., 6.], [7., 8.]]) + C = np.array([[4., 5.], [6., 7.]]) + D = np.array([[13., 14.], [15., 16.]]) + sys = ss(A, B, C, D) + # Check that default type generates a warning + # TODO: remove this check with matrix type is deprecated + use_numpy_matrix(True) + with warnings.catch_warnings(record=True) as w: + Wc = gram(sys, 'c') + self.assertTrue(issubclass(w[-1].category, UserWarning)) + self.assertTrue(isinstance(Wc, np.ndarray)) + @unittest.skipIf(not slycot_check(), "slycot not installed") def testGramRc(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5. 6.; 7. 8.") - C = np.matrix("4. 5.; 6. 7.") - D = np.matrix("13. 14.; 15. 16.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5., 6.], [7., 8.]]) + C = np.array([[4., 5.], [6., 7.]]) + D = np.array([[13., 14.], [15., 16.]]) sys = ss(A, B, C, D) - Rctrue = np.matrix("4.30116263 5.6961343; 0. 0.23249528") - Rc = gram(sys,'cf') + Rctrue = np.array([[4.30116263, 5.6961343], [0., 0.23249528]]) + Rc = gram(sys, 'cf', return_type=np.ndarray) np.testing.assert_array_almost_equal(Rc, Rctrue) @unittest.skipIf(not slycot_check(), "slycot not installed") def testGramWo(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5. 6.; 7. 8.") - C = np.matrix("4. 5.; 6. 7.") - D = np.matrix("13. 14.; 15. 16.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5., 6.], [7., 8.]]) + C = np.array([[4., 5.], [6., 7.]]) + D = np.array([[13., 14.], [15., 16.]]) sys = ss(A, B, C, D) - Wotrue = np.matrix("257.5 -94.5; -94.5 56.5") - Wo = gram(sys,'o') + Wotrue = np.array([[257.5, -94.5], [-94.5, 56.5]]) + Wo = gram(sys, 'o', return_type=np.ndarray) np.testing.assert_array_almost_equal(Wo, Wotrue) @unittest.skipIf(not slycot_check(), "slycot not installed") def testGramWo2(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5.], [7.]]) + C = np.array([[6., 8.]]) + D = np.array([[9.]]) sys = ss(A,B,C,D) - Wotrue = np.matrix("198. -72.; -72. 44.") - Wo = gram(sys,'o') + Wotrue = np.array([[198., -72.], [-72., 44.]]) + Wo = gram(sys, 'o', return_type=np.ndarray) np.testing.assert_array_almost_equal(Wo, Wotrue) @unittest.skipIf(not slycot_check(), "slycot not installed") def testGramRo(self): - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5. 6.; 7. 8.") - C = np.matrix("4. 5.; 6. 7.") - D = np.matrix("13. 14.; 15. 16.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5., 6.], [7., 8.]]) + C = np.array([[4., 5.], [6., 7.]]) + D = np.array([[13., 14.], [15., 16.]]) sys = ss(A, B, C, D) - Rotrue = np.matrix("16.04680654 -5.8890222; 0. 4.67112593") - Ro = gram(sys,'of') + Rotrue = np.array([[16.04680654, -5.8890222], [0., 4.67112593]]) + Ro = gram(sys, 'of', return_type=np.ndarray) np.testing.assert_array_almost_equal(Ro, Rotrue) def testGramsys(self): num =[1.] den = [1., 1., 1.] sys = tf(num,den) - self.assertRaises(ValueError, gram, sys, 'o') - self.assertRaises(ValueError, gram, sys, 'c') + self.assertRaises(ValueError, gram, sys, 'o', return_type=np.ndarray) + self.assertRaises(ValueError, gram, sys, 'c', return_type=np.ndarray) def testAcker(self): for states in range(1, self.maxStates): @@ -133,7 +198,7 @@ def testAcker(self): print(sys) # Make sure the system is not degenerate - Cmat = ctrb(sys.A, sys.B) + Cmat = ctrb(sys.A, sys.B, return_type=np.ndarray) if np.linalg.matrix_rank(Cmat) != states: if (self.debug): print(" skipping (not reachable or ill conditioned)") @@ -144,7 +209,7 @@ def testAcker(self): poles = pole(des) # Now place the poles using acker - K = acker(sys.A, sys.B, poles) + K = acker(sys.A, sys.B, poles, return_type=np.ndarray) new = ss(sys.A - sys.B * K, sys.B, sys.C, sys.D) placed = pole(new) @@ -331,7 +396,10 @@ def test_dare(self): def test_suite(): - return unittest.TestLoader().loadTestsFromTestCase(TestStatefbk) + + status1 = unittest.TestLoader().loadTestsFromTestCase(TestStatefbk) + status2 = unittest.TestLoader().loadTestsFromTestCase(TestStatefbk) + return status1 and status2 if __name__ == '__main__': unittest.main() diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py index 30fc2c913..8809b1eb6 100644 --- a/control/tests/statesp_test.py +++ b/control/tests/statesp_test.py @@ -300,11 +300,10 @@ def test_array_access_ss(self): np.testing.assert_array_almost_equal(sys1_11.A, sys1.A) np.testing.assert_array_almost_equal(sys1_11.B, - sys1.B[:, 1]) + sys1.B[:, [1]]) np.testing.assert_array_almost_equal(sys1_11.C, - sys1.C[0, :]) - np.testing.assert_array_almost_equal(sys1_11.D, - sys1.D[0, 1]) + sys1.C[[0], :]) + np.testing.assert_array_almost_equal(sys1_11.D, sys1.D[0,1]) assert sys1.dt == sys1_11.dt @@ -370,14 +369,14 @@ def test_scalar_static_gain(self): g4 = g1 + g2 self.assertEqual(5, g4.D[0, 0]) g5 = g1.feedback(g2) - self.assertAlmostEqual(2. / 7, g5.D[0, 0]) + np.testing.assert_array_almost_equal(2. / 7, g5.D[0, 0]) g6 = g1.append(g2) np.testing.assert_array_equal(np.diag([2, 3]), g6.D) def test_matrix_static_gain(self): """Regression: can we create matrix static gains?""" - d1 = np.matrix([[1, 2, 3], [4, 5, 6]]) - d2 = np.matrix([[7, 8], [9, 10], [11, 12]]) + d1 = np.array([[1, 2, 3], [4, 5, 6]]) + d2 = np.array([[7, 8], [9, 10], [11, 12]]) g1 = StateSpace([], [], [], d1) # _remove_useless_states was making A = [[0]] @@ -387,12 +386,12 @@ def test_matrix_static_gain(self): g3 = StateSpace([], [], [], d2.T) h1 = g1 * g2 - np.testing.assert_array_equal(d1 * d2, h1.D) + np.testing.assert_array_equal(np.dot(d1, d2), h1.D) h2 = g1 + g3 np.testing.assert_array_equal(d1 + d2.T, h2.D) h3 = g1.feedback(g2) np.testing.assert_array_almost_equal( - solve(np.eye(2) + d1 * d2, d1), h3.D) + solve(np.eye(2) + np.dot(d1, d2), d1), h3.D) h4 = g1.append(g2) np.testing.assert_array_equal(block_diag(d1, d2), h4.D) @@ -436,11 +435,11 @@ def test_empty(self): def test_matrix_to_state_space(self): """_convertToStateSpace(matrix) gives ss([],[],[],D)""" - D = np.matrix([[1, 2, 3], [4, 5, 6]]) + D = np.array([[1, 2, 3], [4, 5, 6]]) g = _convertToStateSpace(D) def empty(shape): - m = np.matrix([]) + m = np.array([]) m.shape = shape return m np.testing.assert_array_equal(empty((0, 0)), g.A) @@ -496,6 +495,17 @@ def test_lft(self): np.testing.assert_allclose(np.array(pk.C).reshape(-1), Cmatlab) np.testing.assert_allclose(np.array(pk.D).reshape(-1), Dmatlab) + def test_horner(self): + """Test horner() function""" + # Make sure we can compute the transfer function at a complex value + self.sys322.horner(1.+1.j) + + # Make sure result agrees with frequency response + mag, phase, omega = self.sys322.freqresp([1]) + np.testing.assert_array_almost_equal( + self.sys322.horner(1.j), + mag[:,:,0] * np.exp(1.j * phase[:,:,0])) + class TestRss(unittest.TestCase): """These are tests for the proper functionality of statesp.rss.""" diff --git a/control/tests/test_control_matlab.py b/control/tests/test_control_matlab.py index e45b52523..458092a66 100644 --- a/control/tests/test_control_matlab.py +++ b/control/tests/test_control_matlab.py @@ -9,7 +9,7 @@ import numpy as np import scipy.signal from numpy.testing import assert_array_almost_equal -from numpy import array, asarray, matrix, asmatrix, zeros, ones, linspace,\ +from numpy import array, asarray, zeros, ones, linspace,\ all, hstack, vstack, c_, r_ from matplotlib.pylab import show, figure, plot, legend, subplot2grid from control.matlab import ss, step, impulse, initial, lsim, dcgain, \ @@ -28,11 +28,11 @@ def plot_matrix(self): #Test: can matplotlib correctly plot matrices? #Yes, but slightly inconvenient figure() - t = matrix([[ 1.], + t = array([[ 1.], [ 2.], [ 3.], [ 4.]]) - y = matrix([[ 1., 4.], + y = array([[ 1., 4.], [ 4., 5.], [ 9., 6.], [16., 7.]]) @@ -42,11 +42,11 @@ def plot_matrix(self): def make_SISO_mats(self): """Return matrices for a SISO system""" - A = matrix([[-81.82, -45.45], + A = array([[-81.82, -45.45], [ 10., -1. ]]) - B = matrix([[9.09], + B = array([[9.09], [0. ]]) - C = matrix([[0, 0.159]]) + C = array([[0, 0.159]]) D = zeros((1, 1)) return A, B, C, D @@ -202,7 +202,7 @@ def test_initial(self): #X0=[1,1] : produces a spike subplot2grid(plot_shape, (0, 1)) - t, y = initial(sys, X0=matrix("1; 1")) + t, y = initial(sys, X0=array([1, 1])) plot(t, y) #Test MIMO system @@ -250,9 +250,10 @@ def test_check_convert_shape(self): #Convert array-like objects to arrays #Input is matrix, shape (1,3), must convert to array - arr = _check_convert_array(matrix("1. 2 3"), [(3,), (1,3)], 'Test: ') - assert isinstance(arr, np.ndarray) - assert not isinstance(arr, matrix) + #! RMM, 14 Apr 2019: removed since matrix is no longer supported + # arr = _check_convert_array(matrix("1. 2 3"), [(3,), (1,3)], 'Test: ') + # assert isinstance(arr, np.ndarray) + # assert not isinstance(arr, matrix) #Input is list, shape (1,3), must convert to array arr = _check_convert_array([[1., 2, 3]], [(3,), (1,3)], 'Test: ') @@ -320,9 +321,9 @@ def test_lsim(self): #Test with matrices subplot2grid(plot_shape, (1, 0)) - t = matrix(linspace(0, 1, 100)) - u = matrix(r_[1:1:50j, 0:0:50j]) - x0 = matrix("0.; 0") + t = array(linspace(0, 1, 100)) + u = array(r_[1:1:50j, 0:0:50j]) + x0 = array([0., 0]) y, t_out, _x = lsim(sys, u, t, x0) plot(t_out, y, label='y') plot(t_out, asarray(u/10)[0], label='u/10') @@ -332,7 +333,7 @@ def test_lsim(self): subplot2grid(plot_shape, (1, 1)) A, B, C, D = self.make_MIMO_mats() sys = ss(A, B, C, D) - t = matrix(linspace(0, 1, 100)) + t = array(linspace(0, 1, 100)) u = array([r_[1:1:50j, 0:0:50j], r_[0:1:50j, 0:0:50j]]) x0 = [0, 0, 0, 0] @@ -404,12 +405,12 @@ def test_convert_MIMO_to_SISO(self): #Test with additional systems -------------------------------------------- #They have crossed inputs and direct feedthrough #SISO system - As = matrix([[-81.82, -45.45], + As = array([[-81.82, -45.45], [ 10., -1. ]]) - Bs = matrix([[9.09], + Bs = array([[9.09], [0. ]]) - Cs = matrix([[0, 0.159]]) - Ds = matrix([[0.02]]) + Cs = array([[0, 0.159]]) + Ds = array([[0.02]]) sys_siso = ss(As, Bs, Cs, Ds) # t, y = step(sys_siso) # plot(t, y, label='sys_siso d=0.02') @@ -428,7 +429,7 @@ def test_convert_MIMO_to_SISO(self): [0 , 0 ]]) Cm = array([[0, 0, 0, 0.159], [0, 0.159, 0, 0 ]]) - Dm = matrix([[0, 0.02], + Dm = array([[0, 0.02], [0.02, 0 ]]) sys_mimo = ss(Am, Bm, Cm, Dm) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 20d715483..b1de5c246 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -19,10 +19,10 @@ class TestTimeresp(unittest.TestCase): def setUp(self): """Set up some systems for testing out MATLAB functions""" - A = np.matrix("1. -2.; 3. -4.") - B = np.matrix("5.; 7.") - C = np.matrix("6. 8.") - D = np.matrix("9.") + A = np.array([[1., -2.], [3., -4.]]) + B = np.array([[5.], [7.]]) + C = np.array([[6., 8.]]) + D = np.array([[9.]]) self.siso_ss1 = StateSpace(A, B, C, D) # Create some transfer functions @@ -30,18 +30,19 @@ def setUp(self): self.siso_tf2 = _convert_to_transfer_function(self.siso_ss1) # Create MIMO system, contains ``siso_ss1`` twice - A = np.matrix("1. -2. 0. 0.;" - "3. -4. 0. 0.;" - "0. 0. 1. -2.;" - "0. 0. 3. -4. ") - B = np.matrix("5. 0.;" - "7. 0.;" - "0. 5.;" - "0. 7. ") - C = np.matrix("6. 8. 0. 0.;" - "0. 0. 6. 8. ") - D = np.matrix("9. 0.;" - "0. 9. ") + A = np.array( + [[1., -2., 0., 0.], + [3., -4., 0., 0.], + [0., 0., 1., -2.], + [0., 0., 3., -4.]]) + B = np.array([[5., 0.], + [7., 0.], + [0., 5.], + [0., 7.]]) + C = np.array([[6., 8., 0., 0.], + [0., 0., 6., 8.,]]) + D = np.array([[9., 0.], + [0., 9.]]) self.mimo_ss1 = StateSpace(A, B, C, D) # Create discrete time systems @@ -208,7 +209,7 @@ def test_initial_response(self): # Test MIMO system, which contains ``siso_ss1`` twice sys = self.mimo_ss1 - x0 = np.matrix(".5; 1.; .5; 1.") + x0 = np.array([[.5], [1.], [.5], [1.]]) _t, y_00 = initial_response(sys, T=t, X0=x0, input=0, output=0) _t, y_11 = initial_response(sys, T=t, X0=x0, input=1, output=1) np.testing.assert_array_almost_equal(y_00, youttrue, decimal=4) @@ -217,7 +218,7 @@ def test_initial_response(self): def test_initial_response_no_trim(self): # test MIMO system without trimming t = np.linspace(0, 1, 10) - x0 = np.matrix(".5; 1.; .5; 1.") + x0 = np.array([[.5], [1.], [.5], [1.]]) youttrue = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092, 1.1508, 0.5833, 0.1645, -0.1391]) sys = self.mimo_ss1 @@ -242,7 +243,7 @@ def test_forced_response(self): # test with initial value and special algorithm for ``U=0`` u = 0 - x0 = np.matrix(".5; 1.") + x0 = np.array([[.5], [1.]]) youttrue = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092, 1.1508, 0.5833, 0.1645, -0.1391]) _t, yout, _xout = forced_response(self.siso_ss1, t, u, x0) @@ -252,7 +253,7 @@ def test_forced_response(self): # first system: initial value, second system: step response u = np.array([[0., 0, 0, 0, 0, 0, 0, 0, 0, 0], [1., 1, 1, 1, 1, 1, 1, 1, 1, 1]]) - x0 = np.matrix(".5; 1; 0; 0") + x0 = np.array([.5, 1, 0, 0]) youttrue = np.array([[11., 8.1494, 5.9361, 4.2258, 2.9118, 1.9092, 1.1508, 0.5833, 0.1645, -0.1391], [9., 17.6457, 24.7072, 30.4855, 35.2234, 39.1165, @@ -272,9 +273,9 @@ def test_forced_response(self): def test_lsim_double_integrator(self): # Note: scipy.signal.lsim fails if A is not invertible - A = np.mat("0. 1.;0. 0.") - B = np.mat("0.; 1.") - C = np.mat("1. 0.") + A = [[0., 1.], [0., 0.]] + B = [[0.], [1.]] + C = np.array([1., 0.]) D = 0. sys = StateSpace(A, B, C, D) diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py index 133cb05cf..c59aa7809 100644 --- a/control/tests/xferfcn_test.py +++ b/control/tests/xferfcn_test.py @@ -651,8 +651,8 @@ def test_matrix_multiply(self): h = (b0 + b1*s + b2*s**2)/(a0 + a1*s + a2*s**2 + a3*s**3) H = TransferFunction([[h.num[0][0]], [(h*s).num[0][0]]], [[h.den[0][0]], [h.den[0][0]]]) - H1 = (np.matrix([[1.0, 0]])*H).minreal() - H2 = (np.matrix([[0, 1.0]])*H).minreal() + H1 = (np.array([[1.0, 0]]) * H).minreal() + H2 = (np.array([[0, 1.0]]) * H).minreal() np.testing.assert_array_almost_equal(H.num[0][0], H1.num[0][0]) np.testing.assert_array_almost_equal(H.den[0][0], H1.den[0][0]) np.testing.assert_array_almost_equal(H.num[1][0], H2.num[0][0]) diff --git a/control/timeresp.py b/control/timeresp.py index 2fc794cb5..cf91615e6 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -319,7 +319,7 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, # [ u(dt) ] = exp [ 0 0 I ] [ u0 ] # [u1 - u0] [ 0 0 0 ] [u1 - u0] - M = np.bmat([[A * dt, B * dt, np.zeros((n_states, n_inputs))], + 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))]]) diff --git a/control/xferfcn.py b/control/xferfcn.py index 4598b2475..313e3f04a 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -70,7 +70,6 @@ class TransferFunction(LTI): """TransferFunction(num, den[, dt]) - A class for representing transfer functions The TransferFunction class is used to represent systems in transfer @@ -99,6 +98,10 @@ class TransferFunction(LTI): >>> G = (s + 1)/(s**2 + 2*s + 1) """ + # Allow NDarray * StateSpace to give StateSpace._rmul_() priority + # https://docs.scipy.org/doc/numpy/reference/arrays.classes.html#numpy.class.__array_priority__ + __array_priority__ = 11 # override ndarray and matrix types + def __init__(self, *args): """TransferFunction(num, den[, dt])