Skip to content

Commit 1dcc842

Browse files
slivingstonslivingston
committed
Merge pull request #118 "Updates to balred() and gram()"
#118 Changes are from branch `master` of https://github.com/mdclemen/python-control.git
2 parents c1f5a9e + 4815cf0 commit 1dcc842

File tree

4 files changed

+161
-48
lines changed

4 files changed

+161
-48
lines changed

‎control/modelsimp.py‎

Lines changed: 69 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -205,10 +205,17 @@ def modred(sys, ELIM, method='matchdc'):
205205
rsys = StateSpace(Ar,Br,Cr,Dr)
206206
return rsys
207207

208-
def balred(sys, orders, method='truncate'):
208+
def balred(sys, orders, method='truncate', alpha=None):
209209
"""
210210
Balanced reduced order model of sys of a given order.
211211
States are eliminated based on Hankel singular value.
212+
If sys has unstable modes, they are removed, the
213+
balanced realization is done on the stable part, then
214+
reinserted in accordance with the reference below.
215+
216+
Reference: Hsu,C.S., and Hou,D., 1991,
217+
Reducing unstable linear control systems via real Schur transformation.
218+
Electronics Letters, 27, 984-986.
212219
213220
Parameters
214221
----------
@@ -219,26 +226,46 @@ def balred(sys, orders, method='truncate'):
219226
of systems)
220227
method: string
221228
Method of removing states, either ``'truncate'`` or ``'matchdc'``.
229+
alpha: float
230+
Redefines the stability boundary for eigenvalues of the system matrix A.
231+
By default for continuous-time systems, alpha <= 0 defines the stability
232+
boundary for the real part of A's eigenvalues and for discrete-time
233+
systems, 0 <= alpha <= 1 defines the stability boundary for the modulus
234+
of A's eigenvalues. See SLICOT routines AB09MD and AB09ND for more
235+
information.
222236
223237
Returns
224238
-------
225239
rsys: StateSpace
226-
A reduced order model
240+
A reduced order model or a list of reduced order models if orders is a list
227241
228242
Raises
229243
------
230244
ValueError
231-
* if `method` is not ``'truncate'``
232-
* if eigenvalues of `sys.A` are not all in left half plane
233-
(`sys` must be stable)
245+
* if `method` is not ``'truncate'`` or ``'matchdc'``
234246
ImportError
235-
if slycot routine ab09ad is not found
247+
if slycot routine ab09md or ab09nd is not found
248+
249+
ValueError
250+
if there are more unstable modes than any value in orders
236251
237252
Examples
238253
--------
239-
>>> rsys = balred(sys, order, method='truncate')
254+
>>> rsys = balred(sys, orders, method='truncate')
240255
241256
"""
257+
if method!='truncate' and method!='matchdc':
258+
raise ValueError("supported methods are 'truncate' or 'matchdc'")
259+
elif method=='truncate':
260+
try:
261+
from slycot import ab09md, ab09ad
262+
except ImportError:
263+
raise ControlSlycot("can't find slycot subroutine ab09md or ab09ad")
264+
elif method=='matchdc':
265+
try:
266+
from slycot import ab09nd
267+
except ImportError:
268+
raise ControlSlycot("can't find slycot subroutine ab09nd")
242269

243270
#Check for ss system object, need a utility for this?
244271

@@ -250,29 +277,46 @@ def balred(sys, orders, method='truncate'):
250277
# else:
251278
dico = 'C'
252279

253-
#Check system is stable
254-
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
255-
raise ValueError("Oops, the system is unstable!")
280+
job = 'B' # balanced (B) or not (N)
281+
equil = 'N' # scale (S) or not (N)
282+
if alpha is None:
283+
if dico == 'C':
284+
alpha = 0.
285+
elif dico == 'D':
286+
alpha = 1.
256287

257-
if method=='matchdc':
258-
raise ValueError ("MatchDC not yet supported!")
259-
elif method=='truncate':
260-
try:
261-
from slycot import ab09ad
262-
except ImportError:
263-
raise ControlSlycot("can't find slycot subroutine ab09ad")
264-
job = 'B' # balanced (B) or not (N)
265-
equil = 'N' # scale (S) or not (N)
288+
rsys = [] #empty list for reduced systems
289+
290+
#check if orders is a list or a scalar
291+
try:
292+
order = iter(orders)
293+
except TypeError: #if orders is a scalar
294+
orders = [orders]
295+
296+
for i in orders:
266297
n = np.size(sys.A,0)
267298
m = np.size(sys.B,1)
268299
p = np.size(sys.C,0)
269-
Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=orders,tol=0.0)
270-
271-
rsys = StateSpace(Ar, Br, Cr, sys.D)
300+
if method == 'truncate':
301+
#check system stability
302+
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
303+
#unstable branch
304+
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)
305+
else:
306+
#stable branch
307+
Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=i,tol=0.0)
308+
rsys.append(StateSpace(Ar, Br, Cr, sys.D))
309+
310+
elif method == 'matchdc':
311+
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)
312+
rsys.append(StateSpace(Ar, Br, Cr, Dr))
313+
314+
#if orders was a scalar, just return the single reduced model, not a list
315+
if len(orders) == 1:
316+
return rsys[0]
317+
#if orders was a list/vector, return a list/vector of systems
272318
else:
273-
raise ValueError("Oops, method is not supported!")
274-
275-
return rsys
319+
return rsys
276320

277321
def minreal(sys, tol=None, verbose=True):
278322
'''

‎control/statefbk.py‎

Lines changed: 50 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -317,7 +317,8 @@ def gram(sys,type):
317317
State-space system to compute Gramian for
318318
type: String
319319
Type of desired computation.
320-
`type` is either 'c' (controllability) or 'o' (observability).
320+
`type` is either 'c' (controllability) or 'o' (observability). To compute the
321+
Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability)
321322
322323
Returns
323324
-------
@@ -328,22 +329,27 @@ def gram(sys,type):
328329
------
329330
ValueError
330331
* if system is not instance of StateSpace class
331-
* if `type` is not 'c' or 'o'
332+
* if `type` is not 'c', 'o', 'cf' or 'of'
332333
* if system is unstable (sys.A has eigenvalues not in left half plane)
333334
334335
ImportError
335-
if slycot routin sb03md cannot be found
336+
if slycot routine sb03md cannot be found
337+
if slycot routine sb03od cannot be found
336338
337339
Examples
338340
--------
339341
>>> Wc = gram(sys,'c')
340342
>>> Wo = gram(sys,'o')
343+
>>> Rc = gram(sys,'cf'), where Wc=Rc'*Rc
344+
>>> Ro = gram(sys,'of'), where Wo=Ro'*Ro
341345
342346
"""
343347

344348
#Check for ss system object
345349
if not isinstance(sys,statesp.StateSpace):
346350
raise ValueError("System must be StateSpace!")
351+
if type not in ['c', 'o', 'cf', 'of']:
352+
raise ValueError("That type is not supported!")
347353

348354
#TODO: Check for continous or discrete, only continuous supported right now
349355
# if isCont():
@@ -358,24 +364,45 @@ def gram(sys,type):
358364
if np.any(np.linalg.eigvals(sys.A).real >= 0.0):
359365
raise ValueError("Oops, the system is unstable!")
360366

361-
if type=='c':
362-
tra = 'T'
363-
C = -np.dot(sys.B,sys.B.transpose())
364-
elif type=='o':
365-
tra = 'N'
366-
C = -np.dot(sys.C.transpose(),sys.C)
367-
else:
368-
raise ValueError("Oops, neither observable, nor controllable!")
369-
370-
#Compute Gramian by the Slycot routine sb03md
367+
if type=='c' or type=='o':
368+
#Compute Gramian by the Slycot routine sb03md
371369
#make sure Slycot is installed
372-
try:
373-
from slycot import sb03md
374-
except ImportError:
375-
raise ControlSlycot("can't find slycot module 'sb03md'")
376-
n = sys.states
377-
U = np.zeros((n,n))
378-
A = np.array(sys.A) # convert to NumPy array for slycot
379-
X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
380-
gram = X
381-
return gram
370+
try:
371+
from slycot import sb03md
372+
except ImportError:
373+
raise ControlSlycot("can't find slycot module 'sb03md'")
374+
if type=='c':
375+
tra = 'T'
376+
C = -np.dot(sys.B,sys.B.transpose())
377+
elif type=='o':
378+
tra = 'N'
379+
C = -np.dot(sys.C.transpose(),sys.C)
380+
n = sys.states
381+
U = np.zeros((n,n))
382+
A = np.array(sys.A) # convert to NumPy array for slycot
383+
X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
384+
gram = X
385+
return gram
386+
387+
elif type=='cf' or type=='of':
388+
#Compute cholesky factored gramian from slycot routine sb03od
389+
try:
390+
from slycot import sb03od
391+
except ImportError:
392+
raise ControlSlycot("can't find slycot module 'sb03od'")
393+
tra='N'
394+
n = sys.states
395+
Q = np.zeros((n,n))
396+
A = np.array(sys.A) # convert to NumPy array for slycot
397+
if type=='cf':
398+
m = sys.B.shape[1]
399+
B = np.zeros_like(A)
400+
B[0:m,0:n] = sys.B.transpose()
401+
X,scale,w = sb03od(n, m, A.transpose(), Q, B, dico, fact='N', trans=tra)
402+
elif type=='of':
403+
m = sys.C.shape[0]
404+
C = np.zeros_like(A)
405+
C[0:n,0:m] = sys.C.transpose()
406+
X,scale,w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra)
407+
gram = X
408+
return gram

‎control/tests/modelsimp_test.py‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -107,6 +107,28 @@ def testBalredTruncate(self):
107107
np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4)
108108
np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=4)
109109

110+
def testBalredMatchDC(self):
111+
#controlable canonical realization computed in matlab for the transfer function:
112+
# num = [1 11 45 32], den = [1 15 60 200 60]
113+
A = np.matrix('-15., -7.5, -6.25, -1.875; \
114+
8., 0., 0., 0.; \
115+
0., 4., 0., 0.; \
116+
0., 0., 1., 0.')
117+
B = np.matrix('2.; 0.; 0.; 0.')
118+
C = np.matrix('0.5, 0.6875, 0.7031, 0.5')
119+
D = np.matrix('0.')
120+
sys = ss(A,B,C,D)
121+
orders = 2
122+
rsys = balred(sys,orders,method='matchdc')
123+
Artrue = np.matrix('-4.43094773, -4.55232904; -4.55232904, -5.36195206')
124+
Brtrue = np.matrix('1.36235673; 1.03114388')
125+
Crtrue = np.matrix('1.36235673, 1.03114388')
126+
Drtrue = np.matrix('-0.08383902')
127+
np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=2)
128+
np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=4)
129+
np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4)
130+
np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=4)
131+
110132
def suite():
111133
return unittest.TestLoader().loadTestsFromTestCase(TestModelsimp)
112134

‎control/tests/statefbk_test.py‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,16 @@ def testGramWc(self):
7171
Wc = gram(sys,'c')
7272
np.testing.assert_array_almost_equal(Wc, Wctrue)
7373

74+
def testGramRc(self):
75+
A = np.matrix("1. -2.; 3. -4.")
76+
B = np.matrix("5. 6.; 7. 8.")
77+
C = np.matrix("4. 5.; 6. 7.")
78+
D = np.matrix("13. 14.; 15. 16.")
79+
sys = ss(A, B, C, D)
80+
Rctrue = np.matrix("4.30116263 5.6961343; 0. 0.23249528")
81+
Rc = gram(sys,'cf')
82+
np.testing.assert_array_almost_equal(Rc, Rctrue)
83+
7484
@unittest.skipIf(not slycot_check(), "slycot not installed")
7585
def testGramWo(self):
7686
A = np.matrix("1. -2.; 3. -4.")
@@ -93,6 +103,16 @@ def testGramWo2(self):
93103
Wo = gram(sys,'o')
94104
np.testing.assert_array_almost_equal(Wo, Wotrue)
95105

106+
def testGramRo(self):
107+
A = np.matrix("1. -2.; 3. -4.")
108+
B = np.matrix("5. 6.; 7. 8.")
109+
C = np.matrix("4. 5.; 6. 7.")
110+
D = np.matrix("13. 14.; 15. 16.")
111+
sys = ss(A, B, C, D)
112+
Rotrue = np.matrix("16.04680654 -5.8890222; 0. 4.67112593")
113+
Ro = gram(sys,'of')
114+
np.testing.assert_array_almost_equal(Ro, Rotrue)
115+
96116
def testGramsys(self):
97117
num =[1.]
98118
den = [1., 1., 1.]

0 commit comments

Comments
 (0)