Skip to content

Commit 7ed68ec

Browse files
murrayrmmurrayrm
committed
updated markov() to add tranpose keyword + default warning; tests, PEP8
1 parent 120a926 commit 7ed68ec

File tree

3 files changed

+164
-90
lines changed

3 files changed

+164
-90
lines changed

‎control/modelsimp.py

Lines changed: 133 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -45,17 +45,20 @@
4545

4646
# External packages and modules
4747
import numpy as np
48-
from .exception import ControlSlycot
48+
import warnings
49+
from .exception import ControlSlycot, ControlMIMONotImplemented
4950
from .lti import isdtime, isctime
5051
from .statesp import StateSpace
5152
from .statefbk import gram
5253

5354
__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
5455

56+
5557
# Hankel Singular Value Decomposition
56-
# The following returns the Hankel singular values, which are singular values
57-
#of the matrix formed by multiplying the controllability and observability
58-
#grammians
58+
#
59+
# The following returns the Hankel singular values, which are singular values
60+
# of the matrix formed by multiplying the controllability and observability
61+
# Gramians
5962
def hsvd(sys):
6063
"""Calculate the Hankel singular values.
6164
@@ -90,8 +93,8 @@ def hsvd(sys):
9093
if (isdtime(sys, strict=True)):
9194
raise NotImplementedError("Function not implemented in discrete time")
9295

93-
Wc = gram(sys,'c')
94-
Wo = gram(sys,'o')
96+
Wc = gram(sys, 'c')
97+
Wo = gram(sys, 'o')
9598
WoWc = np.dot(Wo, Wc)
9699
w, v = np.linalg.eig(WoWc)
97100

@@ -101,6 +104,7 @@ def hsvd(sys):
101104
# Return the Hankel singular values, high to low
102105
return hsv[::-1]
103106

107+
104108
def modred(sys, ELIM, method='matchdc'):
105109
"""
106110
Model reduction of `sys` by eliminating the states in `ELIM` using a given
@@ -136,21 +140,20 @@ def modred(sys, ELIM, method='matchdc'):
136140
>>> rsys = modred(sys, ELIM, method='truncate')
137141
"""
138142

139-
#Check for ss system object, need a utility for this?
143+
# Check for ss system object, need a utility for this?
140144

141-
#TODO: Check for continous or discrete, only continuous supported right now
142-
# if isCont():
143-
# dico = 'C'
144-
# elif isDisc():
145-
# dico = 'D'
146-
# else:
145+
# TODO: Check for continous or discrete, only continuous supported for now
146+
# if isCont():
147+
# dico = 'C'
148+
# elif isDisc():
149+
# dico = 'D'
150+
# else:
147151
if (isctime(sys)):
148152
dico = 'C'
149153
else:
150154
raise NotImplementedError("Function not implemented in discrete time")
151155

152-
153-
#Check system is stable
156+
# Check system is stable
154157
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
155158
raise ValueError("Oops, the system is unstable!")
156159

@@ -160,22 +163,22 @@ def modred(sys, ELIM, method='matchdc'):
160163
# A1 is a matrix of all columns of sys.A not to eliminate
161164
A1 = sys.A[:, NELIM[0]].reshape(-1, 1)
162165
for i in NELIM[1:]:
163-
A1 = np.hstack((A1, sys.A[:,i].reshape(-1, 1)))
164-
A11 = A1[NELIM,:]
165-
A21 = A1[ELIM,:]
166+
A1 = np.hstack((A1, sys.A[:, i].reshape(-1, 1)))
167+
A11 = A1[NELIM, :]
168+
A21 = A1[ELIM, :]
166169
# A2 is a matrix of all columns of sys.A to eliminate
167170
A2 = sys.A[:, ELIM[0]].reshape(-1, 1)
168171
for i in ELIM[1:]:
169-
A2 = np.hstack((A2, sys.A[:,i].reshape(-1, 1)))
170-
A12 = A2[NELIM,:]
171-
A22 = A2[ELIM,:]
172+
A2 = np.hstack((A2, sys.A[:, i].reshape(-1, 1)))
173+
A12 = A2[NELIM, :]
174+
A22 = A2[ELIM, :]
172175

173-
C1 = sys.C[:,NELIM]
174-
C2 = sys.C[:,ELIM]
175-
B1 = sys.B[NELIM,:]
176-
B2 = sys.B[ELIM,:]
176+
C1 = sys.C[:, NELIM]
177+
C2 = sys.C[:, ELIM]
178+
B1 = sys.B[NELIM, :]
179+
B2 = sys.B[ELIM, :]
177180

178-
if method=='matchdc':
181+
if method == 'matchdc':
179182
# if matchdc, residualize
180183

181184
# Check if the matrix A22 is invertible
@@ -195,7 +198,7 @@ def modred(sys, ELIM, method='matchdc'):
195198
Br = B1 - np.dot(A12, A22I_B2)
196199
Cr = C1 - np.dot(C2, A22I_A21)
197200
Dr = sys.D - np.dot(C2, A22I_B2)
198-
elif method=='truncate':
201+
elif method == 'truncate':
199202
# if truncate, simply discard state x2
200203
Ar = A11
201204
Br = B1
@@ -204,12 +207,12 @@ def modred(sys, ELIM, method='matchdc'):
204207
else:
205208
raise ValueError("Oops, method is not supported!")
206209

207-
rsys = StateSpace(Ar,Br,Cr,Dr)
210+
rsys = StateSpace(Ar, Br, Cr, Dr)
208211
return rsys
209212

213+
210214
def balred(sys, orders, method='truncate', alpha=None):
211-
"""
212-
Balanced reduced order model of sys of a given order.
215+
"""Balanced reduced order model of sys of a given order.
213216
States are eliminated based on Hankel singular value.
214217
If sys has unstable modes, they are removed, the
215218
balanced realization is done on the stable part, then
@@ -229,22 +232,23 @@ def balred(sys, orders, method='truncate', alpha=None):
229232
method: string
230233
Method of removing states, either ``'truncate'`` or ``'matchdc'``.
231234
alpha: float
232-
Redefines the stability boundary for eigenvalues of the system matrix A.
233-
By default for continuous-time systems, alpha <= 0 defines the stability
234-
boundary for the real part of A's eigenvalues and for discrete-time
235-
systems, 0 <= alpha <= 1 defines the stability boundary for the modulus
236-
of A's eigenvalues. See SLICOT routines AB09MD and AB09ND for more
237-
information.
235+
Redefines the stability boundary for eigenvalues of the system
236+
matrix A. By default for continuous-time systems, alpha <= 0
237+
defines the stability boundary for the real part of A's eigenvalues
238+
and for discrete-time systems, 0 <= alpha <= 1 defines the stability
239+
boundary for the modulus of A's eigenvalues. See SLICOT routines
240+
AB09MD and AB09ND for more information.
238241
239242
Returns
240243
-------
241244
rsys: StateSpace
242-
A reduced order model or a list of reduced order models if orders is a list
245+
A reduced order model or a list of reduced order models if orders is
246+
a list.
243247
244248
Raises
245249
------
246250
ValueError
247-
* if `method` is not ``'truncate'`` or ``'matchdc'``
251+
If `method` is not ``'truncate'`` or ``'matchdc'``
248252
ImportError
249253
if slycot routine ab09ad, ab09md, or ab09nd is not found
250254
@@ -256,70 +260,78 @@ def balred(sys, orders, method='truncate', alpha=None):
256260
>>> rsys = balred(sys, orders, method='truncate')
257261
258262
"""
259-
if method!='truncate' and method!='matchdc':
263+
if method != 'truncate' and method != 'matchdc':
260264
raise ValueError("supported methods are 'truncate' or 'matchdc'")
261-
elif method=='truncate':
265+
elif method == 'truncate':
262266
try:
263267
from slycot import ab09md, ab09ad
264268
except ImportError:
265-
raise ControlSlycot("can't find slycot subroutine ab09md or ab09ad")
266-
elif method=='matchdc':
269+
raise ControlSlycot(
270+
"can't find slycot subroutine ab09md or ab09ad")
271+
elif method == 'matchdc':
267272
try:
268273
from slycot import ab09nd
269274
except ImportError:
270275
raise ControlSlycot("can't find slycot subroutine ab09nd")
271276

272-
#Check for ss system object, need a utility for this?
277+
# Check for ss system object, need a utility for this?
273278

274-
#TODO: Check for continous or discrete, only continuous supported right now
275-
# if isCont():
276-
# dico = 'C'
277-
# elif isDisc():
278-
# dico = 'D'
279-
# else:
279+
# TODO: Check for continous or discrete, only continuous supported for now
280+
# if isCont():
281+
# dico = 'C'
282+
# elif isDisc():
283+
# dico = 'D'
284+
# else:
280285
dico = 'C'
281286

282-
job = 'B' # balanced (B) or not (N)
283-
equil = 'N' # scale (S) or not (N)
287+
job = 'B' # balanced (B) or not (N)
288+
equil = 'N' # scale (S) or not (N)
284289
if alpha is None:
285290
if dico == 'C':
286291
alpha = 0.
287292
elif dico == 'D':
288293
alpha = 1.
289294

290-
rsys = [] #empty list for reduced systems
295+
rsys = [] # empty list for reduced systems
291296

292-
#check if orders is a list or a scalar
297+
# check if orders is a list or a scalar
293298
try:
294299
order = iter(orders)
295-
except TypeError: #if orders is a scalar
300+
except TypeError: # if orders is a scalar
296301
orders = [orders]
297302

298303
for i in orders:
299-
n = np.size(sys.A,0)
300-
m = np.size(sys.B,1)
301-
p = np.size(sys.C,0)
304+
n = np.size(sys.A, 0)
305+
m = np.size(sys.B, 1)
306+
p = np.size(sys.C, 0)
302307
if method == 'truncate':
303-
#check system stability
308+
# check system stability
304309
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
305-
#unstable branch
306-
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)
310+
# unstable branch
311+
Nr, Ar, Br, Cr, Ns, hsv = ab09md(
312+
dico, job, equil, n, m, p, sys.A, sys.B, sys.C,
313+
alpha=alpha, nr=i, tol=0.0)
307314
else:
308-
#stable branch
309-
Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=i,tol=0.0)
315+
# stable branch
316+
Nr, Ar, Br, Cr, hsv = ab09ad(
317+
dico, job, equil, n, m, p, sys.A, sys.B, sys.C,
318+
nr=i, tol=0.0)
310319
rsys.append(StateSpace(Ar, Br, Cr, sys.D))
311320

312321
elif method == 'matchdc':
313-
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)
322+
Nr, Ar, Br, Cr, Dr, Ns, hsv = ab09nd(
323+
dico, job, equil, n, m, p, sys.A, sys.B, sys.C, sys.D,
324+
alpha=alpha, nr=i, tol1=0.0, tol2=0.0)
314325
rsys.append(StateSpace(Ar, Br, Cr, Dr))
315326

316-
#if orders was a scalar, just return the single reduced model, not a list
327+
# if orders was a scalar, just return the single reduced model, not a list
317328
if len(orders) == 1:
318329
return rsys[0]
319-
#if orders was a list/vector, return a list/vector of systems
330+
# if orders was a list/vector, return a list/vector of systems
320331
else:
321332
return rsys
322333

334+
323335
def minreal(sys, tol=None, verbose=True):
324336
'''
325337
Eliminates uncontrollable or unobservable states in state-space
@@ -347,9 +359,10 @@ def minreal(sys, tol=None, verbose=True):
347359
nstates=len(sys.pole()) - len(sysr.pole())))
348360
return sysr
349361

362+
350363
def era(YY, m, n, nin, nout, r):
351-
"""
352-
Calculate an ERA model of order `r` based on the impulse-response data `YY`.
364+
"""Calculate an ERA model of order `r` based on the impulse-response data
365+
`YY`.
353366
354367
.. note:: This function is not implemented yet.
355368
@@ -376,41 +389,76 @@ def era(YY, m, n, nin, nout, r):
376389
Examples
377390
--------
378391
>>> rsys = era(YY, m, n, nin, nout, r)
392+
379393
"""
380394
raise NotImplementedError('This function is not implemented yet.')
381395

382-
def markov(Y, U, m):
383-
"""
384-
Calculate the first `M` Markov parameters [D CB CAB ...]
396+
397+
def markov(Y, U, m, transpose=None):
398+
"""Calculate the first `M` Markov parameters [D CB CAB ...]
385399
from input `U`, output `Y`.
386400
387401
Parameters
388402
----------
389-
Y: array_like
390-
Output data
391-
U: array_like
392-
Input data
393-
m: int
394-
Number of Markov parameters to output
403+
Y : array_like
404+
Output data. If the array is 1D, the system is assumed to be single
405+
input. If the array is 2D and transpose=False, the columns of `Y`
406+
are taken as time points, otherwise the rows of `Y` are taken as
407+
time points.
408+
U : array_like
409+
Input data, arranged in the same way as `Y`.
410+
m : int
411+
Number of Markov parameters to output.
412+
transpose : bool, optional
413+
Assume that input data is transposed relative to the standard
414+
:ref:`time-series-convention`. The default value is true for
415+
backward compatibility with legacy code.
395416
396417
Returns
397418
-------
398-
H: ndarray
419+
H : ndarray
399420
First m Markov parameters
400421
401422
Notes
402423
-----
403-
Currently only works for SISO
424+
Currently only works for SISO systems.
425+
426+
This function does not currently comply with the Python Control Library
427+
:ref:`time-series-convention` for representation of time series data.
428+
Use `transpose=False` to make use of the standard convention (this
429+
will be updated in a future release).
404430
405431
Examples
406432
--------
407-
>>> H = markov(Y, U, m)
408-
"""
433+
>>> T = numpy.linspace(0, 10, 100)
434+
>>> U = numpy.ones((1, 100))
435+
>>> Y = forced_response(tf([1], [1, 1]), T, U)
436+
>>> H = markov(Y, U, m, transpose=False)
409437
410-
# Convert input parameters to matrices (if they aren't already)
411-
Ymat = np.array(Y)
412-
Umat = np.array(U)
413-
n = np.size(U)
438+
"""
439+
# Check on the specified format of the input
440+
if transpose is None:
441+
# For backwards compatibility, assume time series in rows but warn user
442+
warnings.warn(
443+
"Time-series data assumed to be in rows. This will change in a "
444+
"future release. Use `transpose=True` to preserve current "
445+
"behavior.")
446+
transpose = True
447+
448+
# Convert input parameters to 2D arrays (if they aren't already)
449+
Umat = np.array(U, ndmin=2)
450+
Ymat = np.array(Y, ndmin=2)
451+
452+
# Transpose the data from our normal convention to match algorithm
453+
# TODO: rewrite the algorithm to use standard conventions
454+
if not transpose:
455+
Umat = np.transpose(Umat)
456+
Ymat = np.transpose(Ymat)
457+
458+
# Make sure the system is a SISO system
459+
if Umat.shape[1] != 1 or Ymat.shape[1] != 1:
460+
raise ControlMIMONotImplemented
461+
n = Umat.shape[0]
414462

415463
# Construct a matrix of control inputs to invert
416464
UU = Umat
@@ -424,6 +472,6 @@ def markov(Y, U, m):
424472
UU = np.hstack((UU, Ulast))
425473

426474
# Invert and solve for Markov parameters
427-
H = np.linalg.lstsq(UU, Y)[0]
475+
H = np.linalg.lstsq(UU, Ymat)[0]
428476

429477
return H

0 commit comments

Comments
 (0)