Added hsvd and markov functionality and test script TestModelsimp.py.

Kevin Chen · Kevin Chen · commit d1b6b809393e · 2011-02-08T22:13:04.000Z
Steven Brunton &lt;sbrunton@princeton.edu&gt;
diff --git a/src/TestModelsimp.py b/src/TestModelsimp.py
@@ -0,0 +1,28 @@
+#!/usr/bin/env python
+
+from modelsimp import *
+from matlab import *
+import numpy as np
+import unittest
+
+class TestModelsimp(unittest.TestCase):
+    def testHSVD(self):
+        A = np.matrix("1. -2.; 3. -4.")
+        B = np.matrix("5.; 7.")
+        C = np.matrix("6. 8.")
+        D = np.matrix("9.")
+        sys = ss(A,B,C,D)
+        hsv = hsvd(sys)
+        hsvtrue = np.matrix("24.42686 0.5731395")
+        np.testing.assert_array_almost_equal(hsv, hsvtrue)
+
+    def testMarkov(self):
+        U = np.matrix("1.; 1.; 1.; 1.; 1.")
+        Y = U
+        M = 3
+        H = markov(Y,U,M)
+        Htrue = np.matrix("1.; 0.; 0.")
+        np.testing.assert_array_almost_equal( H, Htrue )
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/TestStatefbk.py b/src/TestStatefbk.py
@@ -48,8 +48,7 @@ def testGramWc(self):
         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)
-        sys = 1.
+        sys = ss(A, B, C, D)
         Wctrue = np.matrix("18.5 24.5; 24.5 32.5")
         Wc = gram(sys,'c')
         np.testing.assert_array_almost_equal(Wc, Wctrue)
@@ -60,7 +59,6 @@ def testGramWo(self):
         C = np.matrix("4. 5.; 6. 7.")
         D = np.matrix("13. 14.; 15. 16.")
         sys = ss(A, B, C, D)
-        sys = 1.
         Wotrue = np.matrix("257.5 -94.5; -94.5 56.5")
         Wo = gram(sys,'o')
         np.testing.assert_array_almost_equal(Wo, Wotrue)
diff --git a/src/modelsimp.py b/src/modelsimp.py
@@ -53,6 +53,8 @@ def hsvd(sys):
     =====
     H = hsvd(sys)
 
+    The Hankel singular values are the singular values of the Hankel operator.  In practice, we compute the square root of the eigenvalues of the matrix formed by taking the product of the observability and controllability gramians.  There are other (more efficient) methods based on solving the Lyapunov equation in a particular way (more details soon).  
+
     Inputs
     ------
     sys : a state space system 
@@ -63,15 +65,16 @@ def hsvd(sys):
 
     """
 
-    # Make sure that SLICOT is installed
-    try:
-        from slycot import sb01bd
-    except ImportError:
-        raise ControlSlycot("can't find slycot module 'sb01bd'")
- 
-    H = 1.
+    Wc = gram(sys,'c')
+    Wo = gram(sys,'o')
+
+    WoWc = np.dot(Wo, Wc)
+    w, v = LA.eig(WoWc)
+
+    hsv = np.sqrt(w)
+    hsv = np.matrix(hsv)
     # Return the Hankel singular values
-    return H
+    return hsv
 
 def era(YY,m,n,nin,nout,r):
     """Calculate an ERA model of order r based on the impulse-response data YY
@@ -100,6 +103,7 @@ def markov(Y,U,M):
     Usage
     =====
     H = markov(Y,U,M)
+    Currently only works for SISO
 
     Inputs
     ------
@@ -113,3 +117,23 @@ def markov(Y,U,M):
 
     """
 
+    # Convert input parameters to matrices (if they aren't already)
+    Ymat = np.mat(Y)
+    Umat = np.mat(U)
+    n = np.size(U)
+
+    # Construct a matrix of control inputs to invert
+    UU = Umat
+    for i in range(1, M-1):
+        newCol = np.vstack((0, UU[0:n-1,i-2]))
+        UU = np.hstack((UU, newCol))
+    Ulast = np.vstack((0, UU[0:n-1,M-2]))
+    for i in range(n-1,0,-1):
+        Ulast[i] = np.sum(Ulast[0:i-1])
+    UU = np.hstack((UU, Ulast))
+
+    # Invert and solve for Markov parameters
+    H = UU.I
+    H = np.dot(H, Y)
+
+    return H
