diff --git a/control/modelsimp.py b/control/modelsimp.py index 9fd36923e..4cfcf4048 100644 --- a/control/modelsimp.py +++ b/control/modelsimp.py @@ -45,17 +45,21 @@ # External packages and modules import numpy as np -from .exception import ControlSlycot +import warnings +from .exception import ControlSlycot, ControlMIMONotImplemented, \ + ControlDimension from .lti import isdtime, isctime from .statesp import StateSpace from .statefbk import gram __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal'] + # Hankel Singular Value Decomposition -# The following returns the Hankel singular values, which are singular values -#of the matrix formed by multiplying the controllability and observability -#grammians +# +# The following returns the Hankel singular values, which are singular values +# of the matrix formed by multiplying the controllability and observability +# Gramians def hsvd(sys): """Calculate the Hankel singular values. @@ -90,8 +94,8 @@ def hsvd(sys): if (isdtime(sys, strict=True)): raise NotImplementedError("Function not implemented in discrete time") - Wc = gram(sys,'c') - Wo = gram(sys,'o') + Wc = gram(sys, 'c') + Wo = gram(sys, 'o') WoWc = np.dot(Wo, Wc) w, v = np.linalg.eig(WoWc) @@ -101,6 +105,7 @@ def hsvd(sys): # Return the Hankel singular values, high to low return hsv[::-1] + def modred(sys, ELIM, method='matchdc'): """ Model reduction of `sys` by eliminating the states in `ELIM` using a given @@ -136,21 +141,20 @@ def modred(sys, ELIM, method='matchdc'): >>> rsys = modred(sys, ELIM, method='truncate') """ - #Check for ss system object, need a utility for this? + # Check for ss system object, need a utility for this? - #TODO: Check for continous or discrete, only continuous supported right now - # if isCont(): - # dico = 'C' - # elif isDisc(): - # dico = 'D' - # else: + # TODO: Check for continous or discrete, only continuous supported for now + # if isCont(): + # dico = 'C' + # elif isDisc(): + # dico = 'D' + # else: if (isctime(sys)): dico = 'C' else: raise NotImplementedError("Function not implemented in discrete time") - - #Check system is stable + # Check system is stable if np.any(np.linalg.eigvals(sys.A).real >= 0.0): raise ValueError("Oops, the system is unstable!") @@ -160,22 +164,22 @@ def modred(sys, ELIM, method='matchdc'): # A1 is a matrix of all columns of sys.A not to eliminate A1 = sys.A[:, NELIM[0]].reshape(-1, 1) for i in NELIM[1:]: - A1 = np.hstack((A1, sys.A[:,i].reshape(-1, 1))) - A11 = A1[NELIM,:] - A21 = A1[ELIM,:] + A1 = np.hstack((A1, sys.A[:, i].reshape(-1, 1))) + A11 = A1[NELIM, :] + A21 = A1[ELIM, :] # A2 is a matrix of all columns of sys.A to eliminate A2 = sys.A[:, ELIM[0]].reshape(-1, 1) for i in ELIM[1:]: - A2 = np.hstack((A2, sys.A[:,i].reshape(-1, 1))) - A12 = A2[NELIM,:] - A22 = A2[ELIM,:] + A2 = np.hstack((A2, sys.A[:, i].reshape(-1, 1))) + A12 = A2[NELIM, :] + A22 = A2[ELIM, :] - C1 = sys.C[:,NELIM] - C2 = sys.C[:,ELIM] - B1 = sys.B[NELIM,:] - B2 = sys.B[ELIM,:] + C1 = sys.C[:, NELIM] + C2 = sys.C[:, ELIM] + B1 = sys.B[NELIM, :] + B2 = sys.B[ELIM, :] - if method=='matchdc': + if method == 'matchdc': # if matchdc, residualize # Check if the matrix A22 is invertible @@ -195,7 +199,7 @@ def modred(sys, ELIM, method='matchdc'): Br = B1 - np.dot(A12, A22I_B2) Cr = C1 - np.dot(C2, A22I_A21) Dr = sys.D - np.dot(C2, A22I_B2) - elif method=='truncate': + elif method == 'truncate': # if truncate, simply discard state x2 Ar = A11 Br = B1 @@ -204,12 +208,12 @@ def modred(sys, ELIM, method='matchdc'): else: raise ValueError("Oops, method is not supported!") - rsys = StateSpace(Ar,Br,Cr,Dr) + rsys = StateSpace(Ar, Br, Cr, Dr) return rsys + def balred(sys, orders, method='truncate', alpha=None): - """ - Balanced reduced order model of sys of a given order. + """Balanced reduced order model of sys of a given order. States are eliminated based on Hankel singular value. If sys has unstable modes, they are removed, the balanced realization is done on the stable part, then @@ -229,22 +233,23 @@ def balred(sys, orders, method='truncate', alpha=None): method: string Method of removing states, either ``'truncate'`` or ``'matchdc'``. alpha: float - Redefines the stability boundary for eigenvalues of the system matrix A. - By default for continuous-time systems, alpha <= 0 defines the stability - boundary for the real part of A's eigenvalues and for discrete-time - systems, 0 <= alpha <= 1 defines the stability boundary for the modulus - of A's eigenvalues. See SLICOT routines AB09MD and AB09ND for more - information. + Redefines the stability boundary for eigenvalues of the system + matrix A. By default for continuous-time systems, alpha <= 0 + defines the stability boundary for the real part of A's eigenvalues + and for discrete-time systems, 0 <= alpha <= 1 defines the stability + boundary for the modulus of A's eigenvalues. See SLICOT routines + AB09MD and AB09ND for more information. Returns ------- rsys: StateSpace - A reduced order model or a list of reduced order models if orders is a list + A reduced order model or a list of reduced order models if orders is + a list. Raises ------ ValueError - * if `method` is not ``'truncate'`` or ``'matchdc'`` + If `method` is not ``'truncate'`` or ``'matchdc'`` ImportError if slycot routine ab09ad, ab09md, or ab09nd is not found @@ -256,70 +261,78 @@ def balred(sys, orders, method='truncate', alpha=None): >>> rsys = balred(sys, orders, method='truncate') """ - if method!='truncate' and method!='matchdc': + if method != 'truncate' and method != 'matchdc': raise ValueError("supported methods are 'truncate' or 'matchdc'") - elif method=='truncate': + elif method == 'truncate': try: from slycot import ab09md, ab09ad except ImportError: - raise ControlSlycot("can't find slycot subroutine ab09md or ab09ad") - elif method=='matchdc': + raise ControlSlycot( + "can't find slycot subroutine ab09md or ab09ad") + elif method == 'matchdc': try: from slycot import ab09nd except ImportError: raise ControlSlycot("can't find slycot subroutine ab09nd") - #Check for ss system object, need a utility for this? + # Check for ss system object, need a utility for this? - #TODO: Check for continous or discrete, only continuous supported right now - # if isCont(): - # dico = 'C' - # elif isDisc(): - # dico = 'D' - # else: + # TODO: Check for continous or discrete, only continuous supported for now + # if isCont(): + # dico = 'C' + # elif isDisc(): + # dico = 'D' + # else: dico = 'C' - job = 'B' # balanced (B) or not (N) - equil = 'N' # scale (S) or not (N) + job = 'B' # balanced (B) or not (N) + equil = 'N' # scale (S) or not (N) if alpha is None: if dico == 'C': alpha = 0. elif dico == 'D': alpha = 1. - rsys = [] #empty list for reduced systems + rsys = [] # empty list for reduced systems - #check if orders is a list or a scalar + # check if orders is a list or a scalar try: order = iter(orders) - except TypeError: #if orders is a scalar + except TypeError: # if orders is a scalar orders = [orders] for i in orders: - n = np.size(sys.A,0) - m = np.size(sys.B,1) - p = np.size(sys.C,0) + n = np.size(sys.A, 0) + m = np.size(sys.B, 1) + p = np.size(sys.C, 0) if method == 'truncate': - #check system stability + # check system stability if np.any(np.linalg.eigvals(sys.A).real >= 0.0): - #unstable branch - Nr, Ar, Br, Cr, Ns, hsv = ab09md(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,alpha=alpha,nr=i,tol=0.0) + # unstable branch + Nr, Ar, Br, Cr, Ns, hsv = ab09md( + dico, job, equil, n, m, p, sys.A, sys.B, sys.C, + alpha=alpha, nr=i, tol=0.0) else: - #stable branch - Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=i,tol=0.0) + # stable branch + Nr, Ar, Br, Cr, hsv = ab09ad( + dico, job, equil, n, m, p, sys.A, sys.B, sys.C, + nr=i, tol=0.0) rsys.append(StateSpace(Ar, Br, Cr, sys.D)) elif method == 'matchdc': - Nr, Ar, Br, Cr, Dr, Ns, hsv = ab09nd(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,sys.D,alpha=alpha,nr=i,tol1=0.0,tol2=0.0) + Nr, Ar, Br, Cr, Dr, Ns, hsv = ab09nd( + dico, job, equil, n, m, p, sys.A, sys.B, sys.C, sys.D, + alpha=alpha, nr=i, tol1=0.0, tol2=0.0) rsys.append(StateSpace(Ar, Br, Cr, Dr)) - #if orders was a scalar, just return the single reduced model, not a list + # if orders was a scalar, just return the single reduced model, not a list if len(orders) == 1: return rsys[0] - #if orders was a list/vector, return a list/vector of systems + # if orders was a list/vector, return a list/vector of systems else: return rsys + def minreal(sys, tol=None, verbose=True): ''' Eliminates uncontrollable or unobservable states in state-space @@ -347,9 +360,10 @@ def minreal(sys, tol=None, verbose=True): nstates=len(sys.pole()) - len(sysr.pole()))) return sysr + def era(YY, m, n, nin, nout, r): - """ - Calculate an ERA model of order `r` based on the impulse-response data `YY`. + """Calculate an ERA model of order `r` based on the impulse-response data + `YY`. .. note:: This function is not implemented yet. @@ -376,54 +390,172 @@ def era(YY, m, n, nin, nout, r): Examples -------- >>> rsys = era(YY, m, n, nin, nout, r) + """ raise NotImplementedError('This function is not implemented yet.') -def markov(Y, U, m): - """ - Calculate the first `M` Markov parameters [D CB CAB ...] + +def markov(Y, U, m=None, transpose=None): + """Calculate the first `m` Markov parameters [D CB CAB ...] from input `U`, output `Y`. + This function computes the Markov parameters for a discrete time system + + .. math:: + + x[k+1] &= A x[k] + B u[k] \\\\ + y[k] &= C x[k] + D u[k] + + given data for u and y. The algorithm assumes that that C A^k B = 0 for + k > m-2 (see [1]). Note that the problem is ill-posed if the length of + the input data is less than the desired number of Markov parameters (a + warning message is generated in this case). + Parameters ---------- - Y: array_like - Output data - U: array_like - Input data - m: int - Number of Markov parameters to output + Y : array_like + Output data. If the array is 1D, the system is assumed to be single + input. If the array is 2D and transpose=False, the columns of `Y` + are taken as time points, otherwise the rows of `Y` are taken as + time points. + U : array_like + Input data, arranged in the same way as `Y`. + m : int, optional + Number of Markov parameters to output. Defaults to len(U). + transpose : bool, optional + Assume that input data is transposed relative to the standard + :ref:`time-series-convention`. The default value is true for + backward compatibility with legacy code. Returns ------- - H: ndarray - First m Markov parameters + H : ndarray + First m Markov parameters, [D CB CAB ...] + + References + ---------- + .. [1] J.-N. Juang, M. Phan, L. G. Horta, and R. W. Longman, + Identification of observer/Kalman filter Markov parameters - Theory + and experiments. Journal of Guidance Control and Dynamics, 16(2), + 320-329, 2012. http://doi.org/10.2514/3.21006 Notes ----- - Currently only works for SISO + Currently only works for SISO systems. + + This function does not currently comply with the Python Control Library + :ref:`time-series-convention` for representation of time series data. + Use `transpose=False` to make use of the standard convention (this + will be updated in a future release). Examples -------- - >>> H = markov(Y, U, m) - """ + >>> T = numpy.linspace(0, 10, 100) + >>> U = numpy.ones((1, 100)) + >>> T, Y, _ = forced_response(tf([1], [1, 0.5], True), T, U) + >>> H = markov(Y, U, 3, transpose=False) - # Convert input parameters to matrices (if they aren't already) - Ymat = np.array(Y) - Umat = np.array(U) - n = np.size(U) - - # Construct a matrix of control inputs to invert + """ + # Check on the specified format of the input + if transpose is None: + # For backwards compatibility, assume time series in rows but warn user + warnings.warn( + "Time-series data assumed to be in rows. This will change in a " + "future release. Use `transpose=True` to preserve current " + "behavior.") + transpose = True + + # Convert input parameters to 2D arrays (if they aren't already) + Umat = np.array(U, ndmin=2) + Ymat = np.array(Y, ndmin=2) + + # If data is in transposed format, switch it around + if transpose: + Umat, Ymat = np.transpose(Umat), np.transpose(Ymat) + + # Make sure the system is a SISO system + if Umat.shape[0] != 1 or Ymat.shape[0] != 1: + raise ControlMIMONotImplemented + + # Make sure the number of time points match + if Umat.shape[1] != Ymat.shape[1]: + raise ControlDimension( + "Input and output data are of differnent lengths") + n = Umat.shape[1] + + # If number of desired parameters was not given, set to size of input data + if m is None: + m = Umat.shape[1] + + # Make sure there is enough data to compute parameters + if m > n: + warn.warning("Not enough data for requested number of parameters") + + # + # Original algorithm (with mapping to standard order) + # + # RMM note, 24 Dec 2020: This algorithm sets the problem up correctly + # until the final column of the UU matrix is created, at which point it + # makes some modifications that I don't understand. This version of the + # algorithm does not seem to return the actual Markov parameters for a + # system. + # + # # Create the matrix of (shifted) inputs + # UU = np.transpose(Umat) + # for i in range(1, m-1): + # # Shift previous column down and add a zero at the top + # newCol = np.vstack((0, np.reshape(UU[0:n-1, i-1], (-1, 1)))) + # UU = np.hstack((UU, newCol)) + # + # # Shift previous column down and add a zero at the top + # Ulast = np.vstack((0, np.reshape(UU[0:n-1, m-2], (-1, 1)))) + # + # # Replace the elements of the last column new values (?) + # # Each row gets the sum of the rows above it (?) + # for i in range(n-1, 0, -1): + # Ulast[i] = np.sum(Ulast[0:i-1]) + # UU = np.hstack((UU, Ulast)) + # + # # Solve for the Markov parameters from Y = H @ UU + # # H = [[D], [CB], [CAB], ..., [C A^{m-3} B], [???]] + # H = np.linalg.lstsq(UU, np.transpose(Ymat))[0] + # + # # Markov parameters are in rows => transpose if needed + # return H if transpose else np.transpose(H) + + # + # New algorithm - Construct a matrix of control inputs to invert + # + # This algorithm sets up the following problem and solves it for + # the Markov parameters + # + # [ y(0) ] [ u(0) 0 0 ] [ D ] + # [ y(1) ] [ u(1) u(0) 0 ] [ C B ] + # [ y(2) ] = [ u(2) u(1) u(0) ] [ C A B ] + # [ : ] [ : : : : ] [ : ] + # [ y(n-1) ] [ u(n-1) u(n-2) u(n-3) ... u(n-m) ] [ C A^{m-2} B ] + # + # Note: if the number of Markov parameters (m) is less than the size of + # the input/output data (n), then this algorithm assumes C A^{j} B = 0 + # for j > m-2. See equation (3) in + # + # J.-N. Juang, M. Phan, L. G. Horta, and R. W. Longman, Identification + # of observer/Kalman filter Markov parameters - Theory and + # experiments. Journal of Guidance Control and Dynamics, 16(2), + # 320-329, 2012. http://doi.org/10.2514/3.21006 + # + + # Create matrix of (shifted) inputs UU = Umat - for i in range(1, m-1): - # 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, 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)) + for i in range(1, m): + # Shift previous column down and add a zero at the top + new_row = np.hstack((0, UU[i-1, 0:-1])) + UU = np.vstack((UU, new_row)) + UU = np.transpose(UU) # Invert and solve for Markov parameters - H = np.linalg.lstsq(UU, Y)[0] + YY = np.transpose(Ymat) + H, _, _, _ = np.linalg.lstsq(UU, YY, rcond=None) - return H + # Return the first m Markov parameters + return H if transpose else np.transpose(H) diff --git a/control/tests/modelsimp_array_test.py b/control/tests/modelsimp_array_test.py index 4a6f591e6..dbd6a5796 100644 --- a/control/tests/modelsimp_array_test.py +++ b/control/tests/modelsimp_array_test.py @@ -9,7 +9,7 @@ import control from control.modelsimp import * from control.matlab import * -from control.exception import slycot_check +from control.exception import slycot_check, ControlMIMONotImplemented class TestModelsimp(unittest.TestCase): def setUp(self): @@ -49,14 +49,91 @@ def testHSVD(self): # Go back to using the normal np.array representation control.use_numpy_matrix(False) - def testMarkov(self): - U = np.array([[1.], [1.], [1.], [1.], [1.]]) + def testMarkovSignature(self): + U = np.array([[1., 1., 1., 1., 1.]]) Y = U - M = 3 - H = markov(Y,U,M) - Htrue = np.array([[1.], [0.], [0.]]) + m = 3 + H = markov(Y, U, m, transpose=False) + Htrue = np.array([[1., 0., 0.]]) np.testing.assert_array_almost_equal( H, Htrue ) + # Make sure that transposed data also works + H = markov(np.transpose(Y), np.transpose(U), m, transpose=True) + np.testing.assert_array_almost_equal( H, np.transpose(Htrue) ) + + # Default (in v0.8.4 and below) should be transpose=True (w/ warning) + import warnings + warnings.simplefilter('always', UserWarning) # don't supress + with warnings.catch_warnings(record=True) as w: + # Set up warnings filter to only show warnings in control module + warnings.filterwarnings("ignore") + warnings.filterwarnings("always", module="control") + + # Generate Markov parameters without any arguments + H = markov(np.transpose(Y), np.transpose(U), m) + np.testing.assert_array_almost_equal( H, np.transpose(Htrue) ) + + # Make sure we got a warning + self.assertEqual(len(w), 1) + self.assertIn("assumed to be in rows", str(w[-1].message)) + self.assertIn("change in a future release", str(w[-1].message)) + + # Test example from docstring + T = np.linspace(0, 10, 100) + U = np.ones((1, 100)) + T, Y, _ = control.forced_response( + control.tf([1], [1, 0.5], True), T, U) + H = markov(Y, U, 3, transpose=False) + + # Test example from issue #395 + inp = np.array([1, 2]) + outp = np.array([2, 4]) + mrk = markov(outp, inp, 1, transpose=False) + + # Make sure MIMO generates an error + U = np.ones((2, 100)) # 2 inputs (Y unchanged, with 1 output) + np.testing.assert_raises(ControlMIMONotImplemented, markov, Y, U, m) + + # Make sure markov() returns the right answer + def testMarkovResults(self): + # + # Test over a range of parameters + # + # k = order of the system + # m = number of Markov parameters + # n = size of the data vector + # + # Values should match exactly for n = m, otherewise you get a + # close match but errors due to the assumption that C A^k B = + # 0 for k > m-2 (see modelsimp.py). + # + for k, m, n in \ + ((2, 2, 2), (2, 5, 5), (5, 2, 2), (5, 5, 5), (5, 10, 10)): + + # Generate stable continuous time system + Hc = control.rss(k, 1, 1) + + # Choose sampling time based on fastest time constant / 10 + w, _ = np.linalg.eig(Hc.A) + Ts = np.min(-np.real(w)) / 10. + + # Convert to a discrete time system via sampling + Hd = control.c2d(Hc, Ts, 'zoh') + + # Compute the Markov parameters from state space + Mtrue = np.hstack([Hd.D] + [np.dot( + Hd.C, np.dot(np.linalg.matrix_power(Hd.A, i), + Hd.B)) for i in range(m-1)]) + + # Generate input/output data + T = np.array(range(n)) * Ts + U = np.cos(T) + np.sin(T/np.pi) + _, Y, _ = control.forced_response(Hd, T, U, squeeze=True) + Mcomp = markov(Y, U, m, transpose=False) + + # Compare to results from markov() + np.testing.assert_array_almost_equal(Mtrue, Mcomp) + def testModredMatchDC(self): #balanced realization computed in matlab for the transfer function: # num = [1 11 45 32], den = [1 15 60 200 60] diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py index 2368bd92f..c0ba72a3b 100644 --- a/control/tests/modelsimp_test.py +++ b/control/tests/modelsimp_test.py @@ -25,7 +25,7 @@ def testMarkov(self): U = np.matrix("1.; 1.; 1.; 1.; 1.") Y = U M = 3 - H = markov(Y,U,M) + H = markov(Y, U, M) Htrue = np.matrix("1.; 0.; 0.") np.testing.assert_array_almost_equal( H, Htrue )