balred() now has 'matchdc' option and may handle unstable systems and will accept a list of 'orders'. gram() now does Cholesky factored gramians if desired. Requires slycot routines AB09MD, AB09ND, SB03OD.

mdclemen · mdclemen · commit f2bdde4cfcdd · 2016-10-20T17:04:00.000-07:00
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -50,7 +50,7 @@
 from .statesp import StateSpace
 from .statefbk import gram
 
-__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
+__all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal',]
 
 # Hankel Singular Value Decomposition
 #   The following returns the Hankel singular values, which are singular values
@@ -197,38 +197,7 @@ def modred(sys, ELIM, method='matchdc'):
     rsys = StateSpace(Ar,Br,Cr,Dr)
     return rsys
 
-def stabsep(T_schur, Z_schur, sys, ldim, no_in, no_out):
-    """
-    Performs stable/unstabe decomposition of sys after Schur forms have been computed for system matrix.
-
-    Reference: Hsu,C.S., and Hou,D., 1991,
-    Reducing unstable linear control systems via real Schur transformation.
-    Electronics Letters, 27, 984-986.
-
-    """
-    #Author: M. Clement (mdclemen@eng.ucsd.edu) 2016
-    As = np.asmatrix(T_schur)
-    Bs = Z_schur.T*sys.B
-    Cs = sys.C*Z_schur
-    #from ref 1 eq(1) As = [A_ Ac], Bs = [B_], and Cs = [C_ C+]; _ denotes stable subsystem
-    #                      [0  A+]       [B+]
-    A_ = As[0:ldim,0:ldim]
-    Ac = As[0:ldim,ldim::]
-    Ap = As[ldim::,ldim::]
-    
-    B_ = Bs[0:ldim,:]
-    Bp = Bs[ldim::,:]
-    
-    C_ = Cs[:,0:ldim]
-    Cp = Cs[:,ldim::]
-    #do some more tricky math IAW ref 1 eq(3)
-    B_tilde = np.bmat([[B_, Ac]])
-    D_tilde = np.bmat([[np.zeros((no_out, no_in)), Cp]])
-
-    return A_, B_tilde, C_, D_tilde, Ap, Bp, Cp
-
-
-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.
@@ -249,6 +218,13 @@ def balred(sys, orders, method='truncate'):
         of systems)
     method: string
         Method of removing states, either ``'truncate'`` or ``'matchdc'``.
+    alpha: float
+        Specifies the alpha-stability boundary for the eigenvalues
+        of the state dynamics matrix A. For a continuous-time
+        system, alpha <= 0 is the boundary value for the real parts
+        of eigenvalues, while for a discrete-time system, 0 <= alpha <= 1
+        represents the boundary value for the moduli of eigenvalues.
+        The alpha-stability domain does not include the boundary.
 
     Returns
     -------
@@ -260,7 +236,7 @@ def balred(sys, orders, method='truncate'):
     ValueError
         * if `method` is not ``'truncate'`` or ``'matchdc'``
     ImportError
-        if slycot routine ab09ad or ab09bd is not found
+        if slycot routine ab09md or ab09nd is not found
 
     ValueError
         if there are more unstable modes than any value in orders
@@ -274,18 +250,15 @@ def balred(sys, orders, method='truncate'):
         raise ValueError("supported methods are 'truncate' or 'matchdc'")
     elif method=='truncate':
         try:
-            from slycot import ab09ad
+            from slycot import ab09md
         except ImportError:
-            raise ControlSlycot("can't find slycot subroutine ab09ad") 
+            raise ControlSlycot("can't find slycot subroutine ab09md") 
     elif method=='matchdc':
         try:
-            from slycot import ab09bd
+            from slycot import ab09nd
         except ImportError:
-            raise ControlSlycot("can't find slycot subroutine ab09bd") 
-
+            raise ControlSlycot("can't find slycot subroutine ab09nd") 
 
-    from scipy.linalg import schur#, cholesky, svd
-    from numpy.linalg import cholesky, svd
     #Check for ss system object, need a utility for this?
 
     #TODO: Check for continous or discrete, only continuous supported right now
@@ -297,6 +270,11 @@ def balred(sys, orders, method='truncate'):
     dico = 'C'
     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
 
@@ -305,64 +283,18 @@ def balred(sys, orders, method='truncate'):
         order = iter(orders)
     except TypeError: #if orders is a scalar
         orders = [orders]
-    
-    #first get original system order
-    nn = sys.A.shape[0] #no. of states
-    mm = sys.B.shape[1] #no. of inputs
-    rr = sys.C.shape[0] #no. of outputs
-    #first do the schur decomposition
-    T, V, l = schur(sys.A, sort = 'lhp') #l will contain the number of eigenvalues in the open left half plane, i.e. no. of stable eigenvalues
 
     for i in orders:
-        rorder = i - (nn - l)
-        if rorder <= 0:
-            raise ValueError("System has %i unstable states which is more than ORDER(%i)" % (nn-l, i))
-
-    for i in orders:
-        if (nn - l) > 0: #handles the stable/unstable decomposition if unstable eigenvalues are found, nn - l is the number of ustable eigenvalues
-            #Author: M. Clement (mdclemen@eng.ucsd.edu) 2016
-            print("Unstable eigenvalues found, performing stable/unstable decomposition")
-
-            rorder = i - (nn - l)
-            A_, B_tilde, C_, D_tilde, Ap, Bp, Cp = stabsep(T, V, sys, l, mm, rr)
-
-            subSys = StateSpace(A_, B_tilde, C_, D_tilde)
-            n = np.size(subSys.A,0)
-            m = np.size(subSys.B,1)
-            p = np.size(subSys.C,0)
-
-            if method == 'truncate':
-                Nr, Ar, Br, Cr, hsv = ab09ad(dico,job,equil,n,m,p,subSys.A,subSys.B,subSys.C,nr=rorder,tol=0.0)
-                rsubSys = StateSpace(Ar, Br, Cr, np.zeros((p,m)))
-
-            elif method == 'matchdc':
-                Nr, Ar, Br, Cr, Dr, hsv = ab09bd(dico,job,equil,n,m,p,subSys.A,subSys.B,subSys.C,subSys.D,nr=rorder,tol1=0.0,tol2=0.0)
-                rsubSys = StateSpace(Ar, Br, Cr, Dr)
-                
-            A_r = rsubSys.A
-            #IAW ref 1 eq(4) B^{tilde}_r = [B_r, Acr]
-            B_r = rsubSys.B[:,0:mm]
-            Acr = rsubSys.B[:,mm:mm+(nn-l)]
-            C_r = rsubSys.C
-
-            #now put the unstable subsystem back in
-            Ar = np.bmat([[A_r, Acr], [np.zeros((nn-l,rorder)), Ap]])
-            Br = np.bmat([[B_r], [Bp]])
-            Cr = np.bmat([[C_r, Cp]])
-
+        n = np.size(sys.A,0)
+        m = np.size(sys.B,1)
+        p = np.size(sys.C,0)
+        if method == 'truncate':
+            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)
             rsys.append(StateSpace(Ar, Br, Cr, sys.D))
 
-        else: #stable system branch
-            n = np.size(sys.A,0)
-            m = np.size(sys.B,1)
-            p = np.size(sys.C,0)
-            if method == 'truncate':
-                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, hsv = ab09bd(dico,job,equil,n,m,p,sys.A,sys.B,sys.C,sys.D,nr=rorder,tol1=0.0,tol2=0.0)
-                rsys.append(StateSpace(Ar, Br, Cr, Dr))
+        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:
diff --git a/control/statefbk.py b/control/statefbk.py
@@ -347,6 +347,9 @@ def gram(sys,type):
     #Check for ss system object
     if not isinstance(sys,statesp.StateSpace):
         raise ValueError("System must be StateSpace!")
+    gramtypes = ['c', 'o', 'cf', 'of']
+    if type not in gramtypes:
+        raise ValueError("That type is not supported!")
 
     #TODO: Check for continous or discrete, only continuous supported right now
         # if isCont():
@@ -362,16 +365,6 @@ def gram(sys,type):
     for e in D:
         if e.real >= 0:
             raise ValueError("Oops, the system is unstable!")
-    if type=='c' or type=='cf':
-        tra = 'T'
-        if type=='c':
-            C = -np.dot(sys.B,sys.B.transpose())
-    elif type=='o' or type=='of':
-        tra = 'N'
-        if type=='o':
-            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
@@ -380,24 +373,38 @@ def gram(sys,type):
             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
-        B = np.zeros_like(A)
-        B[0:sys.B.shape[0],0:sys.B.shape[1]] = sys.B
-        m = sys.B.shape[0]
-        X,scale,w = sb03od(n, m, A, Q, B, dico, fact='N', trans=tra)
+        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
 
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
@@ -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)
 
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
@@ -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.")
@@ -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.]
