[python-control-discuss] SF.net SVN: python-control:[40] branches/control-0.4a/src

[<kk...@us...>]

Revision: 40
          http://python-control.svn.sourceforge.net/python-control/?rev=40&view=rev
Author:   kkchen
Date:     2011-02-08 22:13:04 +0000 (Tue, 08 Feb 2011)

Log Message:
-----------
Added hsvd and markov functionality and test script TestModelsimp.py.

Steven Brunton <sbr...@pr...>

Modified Paths:
--------------
    branches/control-0.4a/src/TestStatefbk.py
    branches/control-0.4a/src/modelsimp.py

Added Paths:
-----------
    branches/control-0.4a/src/TestModelsimp.py

Added: branches/control-0.4a/src/TestModelsimp.py
===================================================================
--- branches/control-0.4a/src/TestModelsimp.py	                        (rev 0)
+++ branches/control-0.4a/src/TestModelsimp.py	2011-02-08 22:13:04 UTC (rev 40)
@@ -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()

Modified: branches/control-0.4a/src/TestStatefbk.py
===================================================================
--- branches/control-0.4a/src/TestStatefbk.py	2011-02-08 22:12:57 UTC (rev 39)
+++ branches/control-0.4a/src/TestStatefbk.py	2011-02-08 22:13:04 UTC (rev 40)
@@ -48,8 +48,7 @@
         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 @@
         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)

Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py	2011-02-08 22:12:57 UTC (rev 39)
+++ branches/control-0.4a/src/modelsimp.py	2011-02-08 22:13:04 UTC (rev 40)
@@ -53,6 +53,8 @@
     =====
     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 @@
 
     """
 
-    # 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 @@
     Usage
     =====
     H = markov(Y,U,M)
+    Currently only works for SISO
 
     Inputs
     ------
@@ -113,3 +117,23 @@
 
     """
 
+    # 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


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.