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

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

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

Log Message:
-----------
Bode for SISO systems now compatible with new statespace and transferfunction classes. Added a script for testing frequency response functions.

Lauren Padilla <lpa...@pr...>

Modified Paths:
--------------
    branches/control-0.4a/src/ctrlutil.py
    branches/control-0.4a/src/freqplot.py
    branches/control-0.4a/src/xferfcn.py

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

Added: branches/control-0.4a/src/TestFreqRsp.py
===================================================================
--- branches/control-0.4a/src/TestFreqRsp.py	                        (rev 0)
+++ branches/control-0.4a/src/TestFreqRsp.py	2011-02-08 22:16:42 UTC (rev 85)
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+
+# Script to test frequency response and frequency response plots like bode, nyquist and gang of 4.
+# Especially need to ensure that nothing SISO is broken and that MIMO at least handles exceptions and has some default to SISO in place.
+
+from statesp import StateSpace
+from matlab import ss, tf, bode
+import numpy as np
+import matplotlib.pyplot as plt
+
+# SISO
+plt.close('all')
+
+A = np.matrix('1,1;0,1')
+B = np.matrix('0;1')
+C = np.matrix('1,0')
+D = 0
+sys = StateSpace(A,B,C,D)
+#or try
+#sys = ss(A,B,C,D)
+
+# test frequency response
+omega = np.linspace(10e-2,10e2,1000)
+frq=sys.freqresp(omega)
+
+# MIMO
+B = np.matrix('1,0;0,1')
+D = np.matrix('0,0')
+sysMIMO = ss(A,B,C,D)
+frqMIMO = sysMIMO.freqresp(omega)
+
+plt.figure(1)
+bode(sys)
+
+systf = tf(sys)
+tfMIMO = tf(sysMIMO)
+
+print systf.pole()
+#print tfMIMO.pole() # - should throw not implemented exception
+#print tfMIMO.zero() # - should throw not implemented exception
+
+plt.figure(2)
+bode(systf)
+
+#bode(sysMIMO) # - should throw not implemented exception
+#bode(tfMIMO) # - should throw not implemented exception
+
+#plt.figure(3)
+#plt.semilogx(omega,20*np.log10(np.squeeze(frq[0])))
+
+#plt.figure(4)
+#bode(sysMIMO,omega)
\ No newline at end of file

Modified: branches/control-0.4a/src/ctrlutil.py
===================================================================
--- branches/control-0.4a/src/ctrlutil.py	2011-02-08 22:16:36 UTC (rev 84)
+++ branches/control-0.4a/src/ctrlutil.py	2011-02-08 22:16:42 UTC (rev 85)
@@ -68,10 +68,11 @@
     """
     wrap = 0;
     last = angle[0];
+
     for k in range(len(angle)):
         # See if we need to account for angle wrapping
         if (angle[k] - last > period/2):
-            wrap = -period
+           wrap = -period
         elif (last - angle[k] > period/2):
             wrap = period
 

Modified: branches/control-0.4a/src/freqplot.py
===================================================================
--- branches/control-0.4a/src/freqplot.py	2011-02-08 22:16:36 UTC (rev 84)
+++ branches/control-0.4a/src/freqplot.py	2011-02-08 22:16:42 UTC (rev 85)
@@ -88,50 +88,56 @@
     if (not getattr(syslist, '__iter__', False)):
         syslist = (syslist,)
 
-    # Select a default range if none is provided
-    if (omega == None):
-        omega = default_frequency_range(syslist)
-
     for sys in syslist:
-        # Get the magnitude and phase of the system
-        mag, phase, omega = sys.freqresp(omega)
-        if Hz: omega = omega/(2*sp.pi)
-        if dB: mag = 20*sp.log10(mag)
-        phase = unwrap(phase*180/sp.pi, 360)
+        if (sys.inputs > 1 or sys.outputs > 1):
+            #TODO: Add MIMO bode plots. 
+            raise NotImplementedError("Bode is currently only implemented for SISO systems.")
+        else:
+            # Select a default range if none is provided
+            if (omega == None):
+                omega = default_frequency_range(syslist)
 
-        # Get the dimensions of the current axis, which we will divide up
-        #! TODO: Not current implemented; just use subplot for now
+            # Get the magnitude and phase of the system
+            mag_tmp, phase_tmp, omega = sys.freqresp(omega)
+            mag = np.squeeze(mag_tmp)
+            phase = np.squeeze(phase_tmp)
+            if Hz: omega = omega/(2*sp.pi)
+            if dB: mag = 20*sp.log10(mag)
+            phase = unwrap(phase*180/sp.pi, 360)
 
-        # Magnitude plot
-        plt.subplot(211); 
-        if dB:
-            plt.semilogx(omega, mag)
-            plt.ylabel("Magnitude (dB)")
-        else:
-            plt.loglog(omega, mag)
-            plt.ylabel("Magnitude")
+            # Get the dimensions of the current axis, which we will divide up
+            #! TODO: Not current implemented; just use subplot for now
 
-        # Add a grid to the plot
-        plt.grid(True)
-        plt.grid(True, which='minor')
-        plt.hold(True);
+            # Magnitude plot
+            plt.subplot(211); 
+            if dB:
+                plt.semilogx(omega, mag)
+                plt.ylabel("Magnitude (dB)")
+            else:
+                plt.loglog(omega, mag)
+                plt.ylabel("Magnitude")
 
-        # Phase plot
-        plt.subplot(212);
-        plt.semilogx(omega, phase)
-        plt.hold(True)
+            # Add a grid to the plot
+            plt.grid(True)
+            plt.grid(True, which='minor')
+            plt.hold(True);
 
-        # Add a grid to the plot
-        plt.grid(True)
-        plt.grid(True, which='minor')
-        plt.ylabel("Phase (deg)")
+            # Phase plot
+            plt.subplot(212);
+            plt.semilogx(omega, phase)
+            plt.hold(True)
 
-        # Label the frequency axis
-        if Hz:
-            plt.xlabel("Frequency (Hz)")
-        else:
-            plt.xlabel("Frequency (rad/sec)")
+            # Add a grid to the plot
+            plt.grid(True)
+            plt.grid(True, which='minor')
+            plt.ylabel("Phase (deg)")
 
+            # Label the frequency axis
+            if Hz:
+                plt.xlabel("Frequency (Hz)")
+            else:
+                plt.xlabel("Frequency (rad/sec)")
+
     return (211, 212)
 
 # Nyquist plot
@@ -316,8 +322,8 @@
     features = np.array(())
     for sys in syslist:
         # Add new features to the list
-        features = np.concatenate((features, np.abs(sys.poles)))
-        features = np.concatenate((features, np.abs(sys.zeros)))
+        features = np.concatenate((features, np.abs(sys.pole())))
+        features = np.concatenate((features, np.abs(sys.zero())))
 
     # Get rid of poles and zeros at the origin
     features = features[features != 0];

Modified: branches/control-0.4a/src/xferfcn.py
===================================================================
--- branches/control-0.4a/src/xferfcn.py	2011-02-08 22:16:36 UTC (rev 84)
+++ branches/control-0.4a/src/xferfcn.py	2011-02-08 22:16:42 UTC (rev 85)
@@ -413,8 +413,12 @@
     def zero(self): 
         """Compute the zeros of a transfer function."""
         
-        raise NotImplementedError("TransferFunction.zero is not implemented \
-yet.")
+        if (self.inputs > 1 or self.outputs > 1):
+            raise NotImplementedError("TransferFunction.zero is currently \
+only implemented for SISO systems.")
+        else:
+            #for now, just give zeros of a SISO tf 
+            return roots(self.num[0][0])
 
     def feedback(self, other, sign=-1): 
         """Feedback interconnection between two LTI objects."""


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