Skip to content

Commit b7dc914

Browse files
slivingstonmdclemen
slivingston
authored and
mdclemen
committed
TEST: On Travis CI, fetch Miniconda using HTTPS
Prefer HTTPS, instead of HTTP, to improve security. The URLs are now https://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh and https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh for Python 2 and 3, respectively.
1 parent cdd3e73 commit b7dc914

File tree

4 files changed

+156
-53
lines changed

4 files changed

+156
-53
lines changed

‎control/modelsimp.py‎

Lines changed: 65 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@
5050
from .statesp import StateSpace
5151
from .statefbk import gram
5252

53-
__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
53+
__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal',]
5454

5555
# Hankel Singular Value Decomposition
5656
# The following returns the Hankel singular values, which are singular values
@@ -197,10 +197,17 @@ def modred(sys, ELIM, method='matchdc'):
197197
rsys = StateSpace(Ar,Br,Cr,Dr)
198198
return rsys
199199

200-
def balred(sys, orders, method='truncate'):
200+
def balred(sys, orders, method='truncate', alpha=None):
201201
"""
202202
Balanced reduced order model of sys of a given order.
203203
States are eliminated based on Hankel singular value.
204+
If sys has unstable modes, they are removed, the
205+
balanced realization is done on the stable part, then
206+
reinserted IAW reference below.
207+
208+
Reference: Hsu,C.S., and Hou,D., 1991,
209+
Reducing unstable linear control systems via real Schur transformation.
210+
Electronics Letters, 27, 984-986.
204211
205212
Parameters
206213
----------
@@ -211,26 +218,46 @@ def balred(sys, orders, method='truncate'):
211218
of systems)
212219
method: string
213220
Method of removing states, either ``'truncate'`` or ``'matchdc'``.
221+
alpha: float
222+
Specifies the alpha-stability boundary for the eigenvalues
223+
of the state dynamics matrix A. For a continuous-time
224+
system, alpha <= 0 is the boundary value for the real parts
225+
of eigenvalues, while for a discrete-time system, 0 <= alpha <= 1
226+
represents the boundary value for the moduli of eigenvalues.
227+
The alpha-stability domain does not include the boundary.
214228
215229
Returns
216230
-------
217231
rsys: StateSpace
218-
A reduced order model
232+
A reduced order model or a list of reduced order models if orders is a list
219233
220234
Raises
221235
------
222236
ValueError
223-
* if `method` is not ``'truncate'``
224-
* if eigenvalues of `sys.A` are not all in left half plane
225-
(`sys` must be stable)
237+
* if `method` is not ``'truncate'`` or ``'matchdc'``
226238
ImportError
227-
if slycot routine ab09ad is not found
239+
if slycot routine ab09md or ab09nd is not found
240+
241+
ValueError
242+
if there are more unstable modes than any value in orders
228243
229244
Examples
230245
--------
231-
>>> rsys = balred(sys, order, method='truncate')
246+
>>> rsys = balred(sys, orders, method='truncate')
232247
233248
"""
249+
if method!='truncate' and method!='matchdc':
250+
raise ValueError("supported methods are 'truncate' or 'matchdc'")
251+
elif method=='truncate':
252+
try:
253+
from slycot import ab09md
254+
except ImportError:
255+
raise ControlSlycot("can't find slycot subroutine ab09md")
256+
elif method=='matchdc':
257+
try:
258+
from slycot import ab09nd
259+
except ImportError:
260+
raise ControlSlycot("can't find slycot subroutine ab09nd")
234261

235262
#Check for ss system object, need a utility for this?
236263

@@ -241,34 +268,40 @@ def balred(sys, orders, method='truncate'):
241268
# dico = 'D'
242269
# else:
243270
dico = 'C'
244-
245-
#Check system is stable
246-
D,V = np.linalg.eig(sys.A)
247-
# print D.shape
248-
# print D
249-
for e in D:
250-
if e.real >= 0:
251-
raise ValueError("Oops, the system is unstable!")
252-
253-
if method=='matchdc':
254-
raise ValueError ("MatchDC not yet supported!")
255-
elif method=='truncate':
256-
try:
257-
from slycot import ab09ad
258-
except ImportError:
259-
raise ControlSlycot("can't find slycot subroutine ab09ad")
260-
job = 'B' # balanced (B) or not (N)
261-
equil = 'N' # scale (S) or not (N)
271+
job = 'B' # balanced (B) or not (N)
272+
equil = 'N' # scale (S) or not (N)
273+
if alpha is None:
274+
if dico == 'C':
275+
alpha = 0.
276+
elif dico == 'D':
277+
alpha = 1.
278+
279+
rsys = [] #empty list for reduced systems
280+
281+
#check if orders is a list or a scalar
282+
try:
283+
order = iter(orders)
284+
except TypeError: #if orders is a scalar
285+
orders = [orders]
286+
287+
for i in orders:
262288
n = np.size(sys.A,0)
263289
m = np.size(sys.B,1)
264290
p = np.size(sys.C,0)
265-
Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=orders,tol=0.0)
266-
267-
rsys = StateSpace(Ar, Br, Cr, sys.D)
291+
if method == 'truncate':
292+
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)
293+
rsys.append(StateSpace(Ar, Br, Cr, sys.D))
294+
295+
elif method == 'matchdc':
296+
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)
297+
rsys.append(StateSpace(Ar, Br, Cr, Dr))
298+
299+
#if orders was a scalar, just return the single reduced model, not a list
300+
if len(orders) == 1:
301+
return rsys[0]
302+
#if orders was a list/vector, return a list/vector of systems
268303
else:
269-
raise ValueError("Oops, method is not supported!")
270-
271-
return rsys
304+
return rsys
272305

273306
def minreal(sys, tol=None, verbose=True):
274307
'''

‎control/statefbk.py‎

Lines changed: 49 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -316,7 +316,8 @@ def gram(sys,type):
316316
State-space system to compute Gramian for
317317
type: String
318318
Type of desired computation.
319-
`type` is either 'c' (controllability) or 'o' (observability).
319+
`type` is either 'c' (controllability) or 'o' (observability). To compute the
320+
Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability)
320321
321322
Returns
322323
-------
@@ -327,22 +328,28 @@ def gram(sys,type):
327328
------
328329
ValueError
329330
* if system is not instance of StateSpace class
330-
* if `type` is not 'c' or 'o'
331+
* if `type` is not 'c', 'o', 'cf' or 'of'
331332
* if system is unstable (sys.A has eigenvalues not in left half plane)
332333
333334
ImportError
334-
if slycot routin sb03md cannot be found
335+
if slycot routine sb03md cannot be found
336+
if slycot routine sb03od cannot be found
335337
336338
Examples
337339
--------
338340
>>> Wc = gram(sys,'c')
339341
>>> Wo = gram(sys,'o')
342+
>>> Rc = gram(sys,'cf'), where Wc=Rc'*Rc
343+
>>> Ro = gram(sys,'of'), where Wo=Ro'*Ro
340344
341345
"""
342346

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

347354
#TODO: Check for continous or discrete, only continuous supported right now
348355
# if isCont():
@@ -358,25 +365,46 @@ def gram(sys,type):
358365
for e in D:
359366
if e.real >= 0:
360367
raise ValueError("Oops, the system is unstable!")
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!")
369368

370369
#Compute Gramian by the Slycot routine sb03md
371370
#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
371+
if type=='c' or type=='o':
372+
try:
373+
from slycot import sb03md
374+
except ImportError:
375+
raise ControlSlycot("can't find slycot module 'sb03md'")
376+
if type=='c':
377+
tra = 'T'
378+
C = -np.dot(sys.B,sys.B.transpose())
379+
elif type=='o':
380+
tra = 'N'
381+
C = -np.dot(sys.C.transpose(),sys.C)
382+
n = sys.states
383+
U = np.zeros((n,n))
384+
A = np.array(sys.A) # convert to NumPy array for slycot
385+
X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
386+
gram = X
387+
return gram
388+
#Comupte cholesky factored gramian from slycot routine sb03od
389+
elif type=='cf' or type=='of':
390+
try:
391+
from slycot import sb03od
392+
except ImportError:
393+
raise ControlSlycot("can't find slycot module 'sb03od'")
394+
tra='N'
395+
n = sys.states
396+
Q = np.zeros((n,n))
397+
A = np.array(sys.A) # convert to NumPy array for slycot
398+
if type=='cf':
399+
m = sys.B.shape[1]
400+
B = np.zeros_like(A)
401+
B[0:m,0:n] = sys.B.transpose()
402+
X,scale,w = sb03od(n, m, A.transpose(), Q, B, dico, fact='N', trans=tra)
403+
elif type=='of':
404+
m = sys.C.shape[0]
405+
C = np.zeros_like(A)
406+
C[0:n,0:m] = sys.C.transpose()
407+
X,scale,w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra)
408+
gram = X
409+
return gram
382410

‎control/tests/modelsimp_test.py‎

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

98+
def testBalredMatchDC(self):
99+
#controlable canonical realization computed in matlab for the transfer function:
100+
# num = [1 11 45 32], den = [1 15 60 200 60]
101+
A = np.matrix('-15., -7.5, -6.25, -1.875; \
102+
8., 0., 0., 0.; \
103+
0., 4., 0., 0.; \
104+
0., 0., 1., 0.')
105+
B = np.matrix('2.; 0.; 0.; 0.')
106+
C = np.matrix('0.5, 0.6875, 0.7031, 0.5')
107+
D = np.matrix('0.')
108+
sys = ss(A,B,C,D)
109+
orders = 2
110+
rsys = balred(sys,orders,method='matchdc')
111+
Artrue = np.matrix('-4.43094773, -4.55232904; -4.55232904, -5.36195206')
112+
Brtrue = np.matrix('1.36235673; 1.03114388')
113+
Crtrue = np.matrix('1.36235673, 1.03114388')
114+
Drtrue = np.matrix('-0.08383902')
115+
np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=2)
116+
np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=4)
117+
np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4)
118+
np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=4)
119+
98120
def suite():
99121
return unittest.TestLoader().loadTestsFromTestCase(TestModelsimp)
100122

‎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)