Skip to content

Updates to balred() and gram() #118

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Jan 20, 2017
Merged

Updates to balred() and gram() #118

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9d840b8
balred() now has 'matchdc' option and may handle unstable systems and…
mdclemen Jan 6, 2017
4815cf0
Updated docstring for balred to remove abbreviations and better descr…
mdclemen Jan 19, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 70 additions & 31 deletions control/modelsimp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -197,10 +197,17 @@ def modred(sys, ELIM, method='matchdc'):
rsys = StateSpace(Ar,Br,Cr,Dr)
return rsys

def balred(sys, orders, method='truncate'):
def balred(sys, orders, method='truncate', alpha=None):
"""
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
reinserted in accordance with the reference below.

Reference: Hsu,C.S., and Hou,D., 1991,
Reducing unstable linear control systems via real Schur transformation.
Electronics Letters, 27, 984-986.

Parameters
----------
Expand All @@ -211,26 +218,46 @@ def balred(sys, orders, method='truncate'):
of systems)
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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While resolving the merge conflict, I added a space here after the .. So, now it reads "...eigenvalues. See SLICOT...".

information.

Returns
-------
rsys: StateSpace
A reduced order model
A reduced order model or a list of reduced order models if orders is a list

Raises
------
ValueError
* if `method` is not ``'truncate'``
* if eigenvalues of `sys.A` are not all in left half plane
(`sys` must be stable)
* if `method` is not ``'truncate'`` or ``'matchdc'``
ImportError
if slycot routine ab09ad is not found
if slycot routine ab09md or ab09nd is not found

ValueError
if there are more unstable modes than any value in orders

Examples
--------
>>> rsys = balred(sys, order, method='truncate')
>>> rsys = balred(sys, orders, method='truncate')

"""
if method!='truncate' and method!='matchdc':
raise ValueError("supported methods are 'truncate' or 'matchdc'")
elif method=='truncate':
try:
from slycot import ab09md, ab09ad
except ImportError:
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?

Expand All @@ -241,34 +268,46 @@ def balred(sys, orders, method='truncate'):
# dico = 'D'
# else:
dico = 'C'

#Check system is stable
D,V = np.linalg.eig(sys.A)
# print D.shape
# print D
for e in D:
if e.real >= 0:
raise ValueError("Oops, the system is unstable!")

if method=='matchdc':
raise ValueError ("MatchDC not yet supported!")
elif method=='truncate':
try:
from slycot import ab09ad
except ImportError:
raise ControlSlycot("can't find slycot subroutine ab09ad")
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

#check if orders is a list or a scalar
try:
order = iter(orders)
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)
Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,nr=orders,tol=0.0)

rsys = StateSpace(Ar, Br, Cr, sys.D)
if method == 'truncate':
#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)
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)
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)
rsys.append(StateSpace(Ar, Br, Cr, Dr))

#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
else:
raise ValueError("Oops, method is not supported!")

return rsys
return rsys

def minreal(sys, tol=None, verbose=True):
'''
Expand Down
69 changes: 48 additions & 21 deletions control/statefbk.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -316,7 +316,8 @@ def gram(sys,type):
State-space system to compute Gramian for
type: String
Type of desired computation.
`type` is either 'c' (controllability) or 'o' (observability).
`type` is either 'c' (controllability) or 'o' (observability). To compute the
Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability)

Returns
-------
Expand All @@ -327,22 +328,27 @@ def gram(sys,type):
------
ValueError
* if system is not instance of StateSpace class
* if `type` is not 'c' or 'o'
* if `type` is not 'c', 'o', 'cf' or 'of'
* if system is unstable (sys.A has eigenvalues not in left half plane)

ImportError
if slycot routin sb03md cannot be found
if slycot routine sb03md cannot be found
if slycot routine sb03od cannot be found

Examples
--------
>>> Wc = gram(sys,'c')
>>> Wo = gram(sys,'o')
>>> Rc = gram(sys,'cf'), where Wc=Rc'*Rc
>>> 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!")
if type not in ['c', 'o', 'cf', 'of']:
raise ValueError("That type is not supported!")

#TODO: Check for continous or discrete, only continuous supported right now
# if isCont():
Expand All @@ -358,25 +364,46 @@ def gram(sys,type):
for e in D:
if e.real >= 0:
raise ValueError("Oops, the system is unstable!")
if type=='c':
tra = 'T'
C = -np.dot(sys.B,sys.B.transpose())
elif type=='o':
tra = 'N'
C = -np.dot(sys.C.transpose(),sys.C)
else:
raise ValueError("Oops, neither observable, nor controllable!")

#Compute Gramian by the Slycot routine sb03md
#make sure Slycot is installed
try:
from slycot import sb03md
except ImportError:
raise ControlSlycot("can't find slycot module 'sb03md'")
n = sys.states
U = np.zeros((n,n))
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
if type=='c' or type=='o':
try:
from slycot import sb03md
except ImportError:
raise ControlSlycot("can't find slycot module 'sb03md'")
if type=='c':
tra = 'T'
C = -np.dot(sys.B,sys.B.transpose())
elif type=='o':
tra = 'N'
C = -np.dot(sys.C.transpose(),sys.C)
n = sys.states
U = np.zeros((n,n))
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
#Comupte cholesky factored gramian from slycot routine sb03od
elif type=='cf' or type=='of':
try:
from slycot import sb03od
except ImportError:
raise ControlSlycot("can't find slycot module 'sb03od'")
tra='N'
n = sys.states
Q = np.zeros((n,n))
A = np.array(sys.A) # convert to NumPy array for slycot
if type=='cf':
m = sys.B.shape[1]
B = np.zeros_like(A)
B[0:m,0:n] = sys.B.transpose()
X,scale,w = sb03od(n, m, A.transpose(), Q, B, dico, fact='N', trans=tra)
elif type=='of':
m = sys.C.shape[0]
C = np.zeros_like(A)
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

22 changes: 22 additions & 0 deletions control/tests/modelsimp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,28 @@ def testBalredTruncate(self):
np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=4)
np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=4)

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.')
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')
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)
np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=4)

def suite():
return unittest.TestLoader().loadTestsFromTestCase(TestModelsimp)

Expand Down
20 changes: 20 additions & 0 deletions control/tests/statefbk_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,16 @@ def testGramWc(self):
Wc = gram(sys,'c')
np.testing.assert_array_almost_equal(Wc, Wctrue)

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.")
sys = ss(A, B, C, D)
Rctrue = np.matrix("4.30116263 5.6961343; 0. 0.23249528")
Rc = gram(sys,'cf')
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.")
Expand All @@ -93,6 +103,16 @@ def testGramWo2(self):
Wo = gram(sys,'o')
np.testing.assert_array_almost_equal(Wo, Wotrue)

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.")
sys = ss(A, B, C, D)
Rotrue = np.matrix("16.04680654 -5.8890222; 0. 4.67112593")
Ro = gram(sys,'of')
np.testing.assert_array_almost_equal(Ro, Rotrue)

def testGramsys(self):
num =[1.]
den = [1., 1., 1.]
Expand Down