diff --git a/control/modelsimp.py b/control/modelsimp.py
index 1bdeb7e99..70b086a5e 100644
--- a/control/modelsimp.py
+++ b/control/modelsimp.py
@@ -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
     ----------
@@ -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
+        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?
 
@@ -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):
     '''
diff --git a/control/statefbk.py b/control/statefbk.py
index 26778960e..c7681f883 100644
--- a/control/statefbk.py
+++ b/control/statefbk.py
@@ -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
     -------
@@ -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():
@@ -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
 
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
index 8db556045..6bf028ce2 100644
--- 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
index e94ba0d0e..ba38ba391 100644
--- 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.]