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

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

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

Log Message:
-----------
Added modred to modelsimp.  Tests pass for MatchDC and Truncate methods.

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

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

Modified: branches/control-0.4a/src/TestModelsimp.py
===================================================================
--- branches/control-0.4a/src/TestModelsimp.py	2011-02-08 22:15:21 UTC (rev 69)
+++ branches/control-0.4a/src/TestModelsimp.py	2011-02-08 22:15:26 UTC (rev 70)
@@ -35,15 +35,15 @@
         C = np.matrix('-0.9057, -0.4068, 0.3263, -0.3474')
         D = np.matrix('0.')
         sys = ss(A,B,C,D)
-        rsys = modred(sys,np.matrix('3, 4'),'matchdc')
+        rsys = modred(sys,[2, 3],'matchdc')
         Artrue = np.matrix('-4.431, -4.552; -4.552, -5.361')
         Brtrue = np.matrix('-1.362; -1.031')
         Crtrue = np.matrix('-1.362, -1.031')
         Drtrue = np.matrix('-0.08384')
-        np.testing.assert_array_almost_equal(sys.A, Artrue)
-        np.testing.assert_array_almost_equal(sys.B, Brtrue)
-        np.testing.assert_array_almost_equal(sys.C, Crtrue)
-        np.testing.assert_array_almost_equal(sys.D, Drtrue)
+        np.testing.assert_array_almost_equal(rsys.A, Artrue,decimal=3)
+        np.testing.assert_array_almost_equal(rsys.B, Brtrue,decimal=3)
+        np.testing.assert_array_almost_equal(rsys.C, Crtrue,decimal=3)
+        np.testing.assert_array_almost_equal(rsys.D, Drtrue,decimal=2)
 
     def testModredTruncate(self):
         #balanced realization computed in matlab for the transfer function:
@@ -56,15 +56,15 @@
         C = np.matrix('-0.9057, -0.4068, 0.3263, -0.3474')
         D = np.matrix('0.')
         sys = ss(A,B,C,D)
-        rsys = modred(sys,np.matrix('3, 4'),'truncate')
+        rsys = modred(sys,[2, 3],'truncate')
         Artrue = np.matrix('-1.958, -1.194; -1.194, -0.8344')
-        Brtrue = np.matrix('-0.9057, -0.4068')
+        Brtrue = np.matrix('-0.9057; -0.4068')
         Crtrue = np.matrix('-0.9057, -0.4068')
         Drtrue = np.matrix('0.')
-        np.testing.assert_array_almost_equal(sys.A, Artrue)
-        np.testing.assert_array_almost_equal(sys.B, Brtrue)
-        np.testing.assert_array_almost_equal(sys.C, Crtrue)
-        np.testing.assert_array_almost_equal(sys.D, Drtrue)
+        np.testing.assert_array_almost_equal(rsys.A, Artrue)
+        np.testing.assert_array_almost_equal(rsys.B, Brtrue)
+        np.testing.assert_array_almost_equal(rsys.C, Crtrue)
+        np.testing.assert_array_almost_equal(rsys.D, Drtrue)
 
     def testBalredMatchDC(self):
         #controlable canonical realization computed in matlab for the transfer function:

Modified: branches/control-0.4a/src/modelsimp.py
===================================================================
--- branches/control-0.4a/src/modelsimp.py	2011-02-08 22:15:21 UTC (rev 69)
+++ branches/control-0.4a/src/modelsimp.py	2011-02-08 22:15:26 UTC (rev 70)
@@ -109,27 +109,57 @@
         # else:
     dico = 'C'
 
-    #TODO: Check system is stable, perhaps a utility in ctrlutil.py
-        # or a method of the StateSpace class?
+    #Check system is stable
     D,V = np.linalg.eig(sys.A)
     for e in D:
         if e.real >= 0:
             raise ValueError, "Oops, the system is unstable!"
+    print ELIM 
+    ELIM = np.sort(ELIM)
+    print ELIM
+    NELIM = []
+    # Create list of elements not to eliminate (NELIM)
+    for i in range(0,len(sys.A)):
+        if i not in ELIM:
+            NELIM.append(i)
+    print NELIM
+    # A1 is a matrix of all columns of sys.A not to eliminate
+    A1 = sys.A[:,NELIM[0]]
+    for i in NELIM[1:]:
+        A1 = np.hstack((A1, sys.A[:,i]))
+    A11 = A1[NELIM,:]  
+    A21 = A1[ELIM,:]  
+    # A2 is a matrix of all columns of sys.A to eliminate
+    A2 = sys.A[:,ELIM[0]]
+    for i in ELIM[1:]:
+        A2 = np.hstack((A2, sys.A[:,i]))
+    A12 = A2[NELIM,:]
+    A22 = A2[ELIM,:]
 
+    C1 = sys.C[:,NELIM]
+    C2 = sys.C[:,ELIM]
+    B1 = sys.B[NELIM,:]
+    B2 = sys.B[ELIM,:]
+
+    print np.shape(A22)
+    A22I = np.linalg.inv(A22)
+
     if method=='matchdc':
-        raise ValueError, "Not supported yet!"
+        # if matchdc, residualize
+        Ar = A11 - A12*A22.I*A21
+        Br = B1 - A12*A22.I*B2
+        Cr = C1 - C2*A22.I*A21
+        Dr = sys.D - C2*A22.I*B2
     elif method=='truncate':
-        raise ValueError, "Not supported yet!"
+        # if truncate, simply discard state x2
+        Ar = A11 
+        Br = B1
+        Cr = C1
+        Dr = sys.D 
     else:
         raise ValueError, "Oops, method is not supported!"
 
-    #Compute rhs using the Slycot routine XXXXXX 
-        #make sure Slycot is installed
-    #try:
-    #    from slycot import XXXXXX
-    #except ImportError:
-    #    raise ControlSlycot("can't find slycot module 'XXXXXX'")
-    rsys = 0.
+    rsys = StateSpace(Ar,Br,Cr,Dr)
     return rsys
 
 def balred(sys,orders,method='truncate'):
@@ -162,8 +192,7 @@
         # else:
     dico = 'C'
 
-    #TODO: Check system is stable, perhaps a utility in ctrlutil.py
-        # or a method of the StateSpace class?
+    #Check system is stable
     D,V = np.linalg.eig(sys.A)
     for e in D:
         if e.real >= 0:


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