python-control-discuss — Discussion list for Python Control Systems Library

Re: [python-control-discuss] __repr__

[Rene v. P. <ren...@gm...>]

Numpy does this:

Imported NumPy 1.6.2, SciPy 0.10.1, Matplotlib 1.0.1
Type "scientific" for more details.
>>> import numpy as np
>>> x = np.arange(0, 9)
>>> x
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> y = array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> y
array([0, 1, 2, 3, 4, 5, 6, 7, 8])
>>> print x
[0 1 2 3 4 5 6 7 8]

But for big arrays

>>> y = np.arange(0, 100, 0.01)


>>> y

array([ 0.00000000e+00, 1.00000000e-02, 2.00000000e-02, ...,

         9.99700000e+01, 9.99800000e+01, 9.99900000e+01])


I guess for python-control the proper thing is something like:


>>> h

TransferFunction([1], [1, 2, 1])


I wouldn't have known the distinction between __repr__ and __str__ either,
but stumbled upon it since I am currently teaching python to 1st year
undergrads. I also teach control theory to 2nd year, of which a quarter is
using Python (from last year's Python programming trial) and the rest is
using Matlab.


René


Associate Professor

Aerospace Engineering

Delft University of Technology


On 12 June 2013 04:26, Ryan Krauss <rk...@si...> wrote:

> > __repr__ is ideally for printing a representation that can directly be
> fed back into Python to > re-create the object.
>
> I guess I didn't realize this convention.  Like you said, I don't think it
> is typically followed.  When I want to look at a variable at the IPython
> prompt, I am not used to having to type print in front of it.
>
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
>
>
> On Tue, Jun 11, 2013 at 2:37 PM, Rene van Paassen <
> ren...@gm...> wrote:
>
>> __str__ is for pretty printing
>> __repr__ is ideally for printing a representation that can directly be
>> fed back into Python to re-create the object. Many objects don't implement
>> that, though, it is basically limited to the standard Python objects.
>>
>>
>> On 10 June 2013 17:53, Ryan Krauss <rk...@si...> wrote:
>>
>>> I probably won't get to this one for a few weeks (I have a conference
>>> and another paper coming up), but I would like to contribute a __repr__
>>> method for TransferFunctions if it would be welcome.  I personally would
>>> prefer that when I type the name of a TransferFunction instance and hit
>>> enter that it pretty print an ascii numerator and denominator and a
>>> horizontal line.  It might be tricky to do it really nicely with exponents,
>>> but this is what happend now:
>>>
>>> In [2]: G = control.TransferFunction(1,[1,2])
>>>
>>> In [3]: G
>>> Out[3]: <control.xferfcn.TransferFunction instance at 0x10e94cf80>
>>>
>>> I then have to ask ro G.num and G.den separately, which I personally
>>> don't like.
>>>
>>> I see that there is a __str__ method which pretty much does what I would
>>> want, so that "print G" gives me an ascii representation.
>>>
>>> Is there any reason not to have __repr__ simply call __str__?  Is there
>>> some historical reason to have one and not the other?
>>>
>>> --
>>> Ryan Krauss, Ph.D.
>>> Associate Professor
>>> Mechanical Engineering
>>> Southern Illinois University Edwardsville
>>>
>>>
>>> ------------------------------------------------------------------------------
>>> How ServiceNow helps IT people transform IT departments:
>>> 1. A cloud service to automate IT design, transition and operations
>>> 2. Dashboards that offer high-level views of enterprise services
>>> 3. A single system of record for all IT processes
>>> http://p.sf.net/sfu/servicenow-d2d-j
>>> _______________________________________________
>>> python-control-discuss mailing list
>>> pyt...@li...
>>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>>>
>>>
>>
>>
>> --
>> René van Paassen        | ______o____/_|     Ren...@gm...
>>                        <[___\_\_-----<               t: +31 15 2628685
>>                         |  o'                   mobile: +31 6 39846891
>>
>>
>>
>> ------------------------------------------------------------------------------
>> This SF.net email is sponsored by Windows:
>>
>> Build for Windows Store.
>>
>> http://p.sf.net/sfu/windows-dev2dev
>>
>> _______________________________________________
>> python-control-discuss mailing list
>> pyt...@li...
>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>>
>>
>
>
> ------------------------------------------------------------------------------
> This SF.net email is sponsored by Windows:
>
> Build for Windows Store.
>
> http://p.sf.net/sfu/windows-dev2dev
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>


-- 
René van Paassen        | ______o____/_|     Ren...@gm...
                       <[___\_\_-----<               t: +31 15 2628685
                        |  o'                   mobile: +31 6 39846891

Re: [python-control-discuss] __repr__

[Ryan K. <rk...@si...>]

I guess we are already following the scipy.signal.ltisys convention:

In [4]: test = signal.lti([1,2,3],[1,4,5,6,7])

In [5]: test
Out[5]: <scipy.signal.ltisys.lti at 0x1129c5510>

In [6]: test.num
Out[6]: array([ 1.,  2.,  3.])

In [7]: test.den
Out[7]: array([ 1.,  4.,  5.,  6.,  7.])

We are even a step ahead of them in user friendliness:

In [8]: print test
<scipy.signal.ltisys.lti object at 0x1129c5510>


That just doesn't seem user friendly to me.

--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville


On Tue, Jun 11, 2013 at 9:26 PM, Ryan Krauss <rk...@si...> wrote:

> > __repr__ is ideally for printing a representation that can directly be
> fed back into Python to > re-create the object.
>
> I guess I didn't realize this convention.  Like you said, I don't think it
> is typically followed.  When I want to look at a variable at the IPython
> prompt, I am not used to having to type print in front of it.
>
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
>
>
> On Tue, Jun 11, 2013 at 2:37 PM, Rene van Paassen <
> ren...@gm...> wrote:
>
>> __str__ is for pretty printing
>> __repr__ is ideally for printing a representation that can directly be
>> fed back into Python to re-create the object. Many objects don't implement
>> that, though, it is basically limited to the standard Python objects.
>>
>>
>> On 10 June 2013 17:53, Ryan Krauss <rk...@si...> wrote:
>>
>>> I probably won't get to this one for a few weeks (I have a conference
>>> and another paper coming up), but I would like to contribute a __repr__
>>> method for TransferFunctions if it would be welcome.  I personally would
>>> prefer that when I type the name of a TransferFunction instance and hit
>>> enter that it pretty print an ascii numerator and denominator and a
>>> horizontal line.  It might be tricky to do it really nicely with exponents,
>>> but this is what happend now:
>>>
>>> In [2]: G = control.TransferFunction(1,[1,2])
>>>
>>> In [3]: G
>>> Out[3]: <control.xferfcn.TransferFunction instance at 0x10e94cf80>
>>>
>>> I then have to ask ro G.num and G.den separately, which I personally
>>> don't like.
>>>
>>> I see that there is a __str__ method which pretty much does what I would
>>> want, so that "print G" gives me an ascii representation.
>>>
>>> Is there any reason not to have __repr__ simply call __str__?  Is there
>>> some historical reason to have one and not the other?
>>>
>>> --
>>> Ryan Krauss, Ph.D.
>>> Associate Professor
>>> Mechanical Engineering
>>> Southern Illinois University Edwardsville
>>>
>>>
>>> ------------------------------------------------------------------------------
>>> How ServiceNow helps IT people transform IT departments:
>>> 1. A cloud service to automate IT design, transition and operations
>>> 2. Dashboards that offer high-level views of enterprise services
>>> 3. A single system of record for all IT processes
>>> http://p.sf.net/sfu/servicenow-d2d-j
>>> _______________________________________________
>>> python-control-discuss mailing list
>>> pyt...@li...
>>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>>>
>>>
>>
>>
>> --
>> René van Paassen        | ______o____/_|     Ren...@gm...
>>                        <[___\_\_-----<               t: +31 15 2628685
>>                         |  o'                   mobile: +31 6 39846891
>>
>>
>>
>> ------------------------------------------------------------------------------
>> This SF.net email is sponsored by Windows:
>>
>> Build for Windows Store.
>>
>> http://p.sf.net/sfu/windows-dev2dev
>>
>> _______________________________________________
>> python-control-discuss mailing list
>> pyt...@li...
>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>>
>>
>

Re: [python-control-discuss] __repr__

[Ryan K. <rk...@si...>]

> __repr__ is ideally for printing a representation that can directly be
fed back into Python to > re-create the object.

I guess I didn't realize this convention.  Like you said, I don't think it
is typically followed.  When I want to look at a variable at the IPython
prompt, I am not used to having to type print in front of it.


--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville


On Tue, Jun 11, 2013 at 2:37 PM, Rene van Paassen <ren...@gm...
> wrote:

> __str__ is for pretty printing
> __repr__ is ideally for printing a representation that can directly be fed
> back into Python to re-create the object. Many objects don't implement
> that, though, it is basically limited to the standard Python objects.
>
>
> On 10 June 2013 17:53, Ryan Krauss <rk...@si...> wrote:
>
>> I probably won't get to this one for a few weeks (I have a conference and
>> another paper coming up), but I would like to contribute a __repr__ method
>> for TransferFunctions if it would be welcome.  I personally would prefer
>> that when I type the name of a TransferFunction instance and hit enter that
>> it pretty print an ascii numerator and denominator and a horizontal line.
>>  It might be tricky to do it really nicely with exponents, but this is what
>> happend now:
>>
>> In [2]: G = control.TransferFunction(1,[1,2])
>>
>> In [3]: G
>> Out[3]: <control.xferfcn.TransferFunction instance at 0x10e94cf80>
>>
>> I then have to ask ro G.num and G.den separately, which I personally
>> don't like.
>>
>> I see that there is a __str__ method which pretty much does what I would
>> want, so that "print G" gives me an ascii representation.
>>
>> Is there any reason not to have __repr__ simply call __str__?  Is there
>> some historical reason to have one and not the other?
>>
>> --
>> Ryan Krauss, Ph.D.
>> Associate Professor
>> Mechanical Engineering
>> Southern Illinois University Edwardsville
>>
>>
>> ------------------------------------------------------------------------------
>> How ServiceNow helps IT people transform IT departments:
>> 1. A cloud service to automate IT design, transition and operations
>> 2. Dashboards that offer high-level views of enterprise services
>> 3. A single system of record for all IT processes
>> http://p.sf.net/sfu/servicenow-d2d-j
>> _______________________________________________
>> python-control-discuss mailing list
>> pyt...@li...
>> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>>
>>
>
>
> --
> René van Paassen        | ______o____/_|     Ren...@gm...
>                        <[___\_\_-----<               t: +31 15 2628685
>                         |  o'                   mobile: +31 6 39846891
>
>
>
> ------------------------------------------------------------------------------
> This SF.net email is sponsored by Windows:
>
> Build for Windows Store.
>
> http://p.sf.net/sfu/windows-dev2dev
>
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>

Re: [python-control-discuss] __repr__

[Rene v. P. <ren...@gm...>]

__str__ is for pretty printing
__repr__ is ideally for printing a representation that can directly be fed
back into Python to re-create the object. Many objects don't implement
that, though, it is basically limited to the standard Python objects.


On 10 June 2013 17:53, Ryan Krauss <rk...@si...> wrote:

> I probably won't get to this one for a few weeks (I have a conference and
> another paper coming up), but I would like to contribute a __repr__ method
> for TransferFunctions if it would be welcome.  I personally would prefer
> that when I type the name of a TransferFunction instance and hit enter that
> it pretty print an ascii numerator and denominator and a horizontal line.
>  It might be tricky to do it really nicely with exponents, but this is what
> happend now:
>
> In [2]: G = control.TransferFunction(1,[1,2])
>
> In [3]: G
> Out[3]: <control.xferfcn.TransferFunction instance at 0x10e94cf80>
>
> I then have to ask ro G.num and G.den separately, which I personally don't
> like.
>
> I see that there is a __str__ method which pretty much does what I would
> want, so that "print G" gives me an ascii representation.
>
> Is there any reason not to have __repr__ simply call __str__?  Is there
> some historical reason to have one and not the other?
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
>
>
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>


-- 
René van Paassen        | ______o____/_|     Ren...@gm...
                       <[___\_\_-----<               t: +31 15 2628685
                        |  o'                   mobile: +31 6 39846891

Re: [python-control-discuss] __repr__

[Richard M. <mu...@cd...>]

We should definitely do something nicer than "<control.xferfcn.TransferFunction instance at 0x10e94cf80>".  What's standard for scipy objects (or scipy.signal objects)?  I can imagine this eventually becoming part of scipy, so perhaps we use whatever standards exist there as a guide?

-richard

On 10 Jun 2013, at 8:53 , Ryan Krauss <rk...@si...> wrote:

> I probably won't get to this one for a few weeks (I have a conference and another paper coming up), but I would like to contribute a __repr__ method for TransferFunctions if it would be welcome.  I personally would prefer that when I type the name of a TransferFunction instance and hit enter that it pretty print an ascii numerator and denominator and a horizontal line.  It might be tricky to do it really nicely with exponents, but this is what happend now:
> 
> In [2]: G = control.TransferFunction(1,[1,2])
> 
> In [3]: G
> Out[3]: <control.xferfcn.TransferFunction instance at 0x10e94cf80>
> 
> I then have to ask ro G.num and G.den separately, which I personally don't like.
> 
> I see that there is a __str__ method which pretty much does what I would want, so that "print G" gives me an ascii representation.
> 
> Is there any reason not to have __repr__ simply call __str__?  Is there some historical reason to have one and not the other?
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

[python-control-discuss] __repr__

[Ryan K. <rk...@si...>]

I probably won't get to this one for a few weeks (I have a conference and
another paper coming up), but I would like to contribute a __repr__ method
for TransferFunctions if it would be welcome.  I personally would prefer
that when I type the name of a TransferFunction instance and hit enter that
it pretty print an ascii numerator and denominator and a horizontal line.
 It might be tricky to do it really nicely with exponents, but this is what
happend now:

In [2]: G = control.TransferFunction(1,[1,2])

In [3]: G
Out[3]: <control.xferfcn.TransferFunction instance at 0x10e94cf80>

I then have to ask ro G.num and G.den separately, which I personally don't
like.

I see that there is a __str__ method which pretty much does what I would
want, so that "print G" gives me an ascii representation.

Is there any reason not to have __repr__ simply call __str__?  Is there
some historical reason to have one and not the other?

--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville

[python-control-discuss] SF.net python-control commit:[281] trunk

[<mur...@us...>]

Revision: 281
          http://sourceforge.net/p/python-control/code/281
Author:   murrayrm
Date:     2013-06-10 15:38:47 +0000 (Mon, 10 Jun 2013)
Log Message:
-----------
docstring updates

Modified Paths:
--------------
    trunk/ChangeLog
    trunk/src/freqplot.py
    trunk/src/xferfcn.py

Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog	2013-06-10 05:05:46 UTC (rev 280)
+++ trunk/ChangeLog	2013-06-10 15:38:47 UTC (rev 281)
@@ -1,3 +1,9 @@
+2013-06-10  Richard Murray  <murray@altura-2.local>
+
+	* src/xferfcn.py (TransferFunction.horner): small fix to docstring
+
+	* src/freqplot.py (nyquist_plot): small fix to docstring
+
 2013-06-09  Richard Murray  <murray@altura-2.local>
 
 	* tests/bdalg_test.py (TestFeedback.testScalarSS)

Modified: trunk/src/freqplot.py
===================================================================
--- trunk/src/freqplot.py	2013-06-10 05:05:46 UTC (rev 280)
+++ trunk/src/freqplot.py	2013-06-10 15:38:47 UTC (rev 281)
@@ -204,7 +204,7 @@
     Examples
     --------
     >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
-    >>> real, imag, freq = nyquist(sys)
+    >>> real, imag, freq = nyquist_plot(sys)
     """
     # If argument was a singleton, turn it into a list
     if (not getattr(syslist, '__iter__', False)):

Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py	2013-06-10 05:05:46 UTC (rev 280)
+++ trunk/src/xferfcn.py	2013-06-10 15:38:47 UTC (rev 281)
@@ -545,10 +545,10 @@
         return self.horner(s)
 
     def horner(self, s):
-        '''Evaluate the systems's transfer function for a complex variable
+        """Evaluate the systems's transfer function for a complex variable
         
         Returns a matrix of values evaluated at complex variable s.
-        '''
+        """
 
         # Preallocate the output.
         if getattr(s, '__iter__', False):

[python-control-discuss] SF.net python-control commit:[275] trunk

[<mur...@us...>]

Revision: 275
          http://sourceforge.net/p/python-control/code/275
Author:   murrayrm
Date:     2013-06-09 15:50:50 +0000 (Sun, 09 Jun 2013)
Log Message:
-----------
Added in patches from Ryan Krauss <rk...@si...>: 

* now possible to call a transfer function with a complex argument and get
  back the value of the transfer function (returns a matrix in MIMO case)

* default argument for second transfer function in feedback() call is now 1.
  Note that this will only work for SISO systems (wasn't sure extending to
  MIMO was that useful).

Modified Paths:
--------------
    trunk/ChangeLog
    trunk/src/bdalg.py
    trunk/src/frdata.py
    trunk/src/statesp.py
    trunk/src/xferfcn.py
    trunk/tests/bdalg_test.py
    trunk/tests/frd_test.py
    trunk/tests/xferfcn_test.py

Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/ChangeLog	2013-06-09 15:50:50 UTC (rev 275)
@@ -1,3 +1,28 @@
+2013-06-09  Richard Murray  <murray@altura-2.local>
+
+	* tests/bdalg_test.py (TestFeedback.testScalarSS)
+	(TestFeedback.testScalarTF): added simple test of default argument
+	to feedback() for state space and transfer function systems
+
+	* tests/frd_test.py (TestFRD.testFeedback): added test of default
+	argument to feedback() for FRD type
+
+	* src/bdalg.py (feedback): 
+	* src/frdata.py (FRD.feedback):
+	* src/statesp.py (StateSpace.feedback):
+	patched in default value of 1 for other system in feedback().
+	Contributed by Ryan Krauss <rk...@si...>
+	
+	* tests/xferfcn_test.py (TestXferFcn.testEvalFrMIMO)
+	(TestXferFcn.testEvalFrSISO): added tests of __call__ version of
+	evalfr. 
+
+	* src/xferfcn.py (TransferFunction.__call__): patched in new
+	__call__ functionality from Ryan Krauss <rk...@si...>, but
+	modified to use evalfr (works with MIMO).  SISO transfer functions
+	return a single complex number, MIMO transfer functions return a
+	matrix of complex numbers.
+
 2013-05-23  Rene van Paassen <ren...@gm...>
 
 	* src/margin.py: re-vamped stability_margins function. Now

Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/bdalg.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -172,7 +172,8 @@
     
     return -sys;
 
-def feedback(sys1, sys2, sign=-1):
+#! TODO: expand to allow sys2 default to work in MIMO case?
+def feedback(sys1, sys2=1, sign=-1):
     """
     Feedback interconnection between two I/O systems.
 

Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/frdata.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -393,7 +393,7 @@
 
         return mag, phase, omega
 
-    def feedback(self, other, sign=-1): 
+    def feedback(self, other=1, sign=-1): 
         """Feedback interconnection between two FRD objects."""
         
         other = _convertToFRD(other, omega=self.omega)

Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/statesp.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -447,7 +447,7 @@
         return roots(num)
 
     # Feedback around a state space system
-    def feedback(self, other, sign=-1):
+    def feedback(self, other=1, sign=-1):
         """Feedback interconnection between two LTI systems."""
  
         other = _convertToStateSpace(other)

Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/src/xferfcn.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -81,7 +81,7 @@
 # External function declarations
 from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
     polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, exp, pi, \
-    where, delete, real, poly
+    where, delete, real, poly, poly1d
 from scipy.signal import lti
 from copy import deepcopy
 from warnings import warn
@@ -231,7 +231,21 @@
         self.den = den
         
         self._truncatecoeff()
-        
+
+    def __call__(self, s):
+        """Evaluate the system's transfer function for a complex vairable
+
+        For a SISO transfer function, returns the value of the
+        transfer function.  For a MIMO transfer fuction, returns a
+        matrix of values evaluated at complex variable s."""
+
+        if (self.inputs > 1 or self.outputs > 1):
+            # MIMO transfer function, return a matrix
+            return self.horner(s)
+        else:
+            # SISO transfer function, return a scalar
+            return self.horner(s)[0][0]
+
     def _truncatecoeff(self):
         """Remove extraneous zero coefficients from num and den.
 
@@ -604,7 +618,7 @@
             #for now, just give zeros of a SISO tf 
             return roots(self.num[0][0])
 
-    def feedback(self, other, sign=-1): 
+    def feedback(self, other=1, sign=-1): 
         """Feedback interconnection between two LTI objects."""
         
         other = _convertToTransferFunction(other)

Modified: trunk/tests/bdalg_test.py
===================================================================
--- trunk/tests/bdalg_test.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/tests/bdalg_test.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -50,6 +50,14 @@
         np.testing.assert_array_almost_equal(ans2.C, [[2.5, 0.]])
         np.testing.assert_array_almost_equal(ans2.D, [[2.5]])
 
+        # Make sure default arugments work as well
+        ans3 = feedback(self.sys2, 1)
+        ans4 = feedback(self.sys2)
+        np.testing.assert_array_almost_equal(ans3.A, ans4.A)
+        np.testing.assert_array_almost_equal(ans3.B, ans4.B)
+        np.testing.assert_array_almost_equal(ans3.C, ans4.C)
+        np.testing.assert_array_almost_equal(ans3.D, ans4.D)
+
     def testScalarTF(self):
         """Scalar system with transfer function feedback block."""
 
@@ -61,6 +69,12 @@
         np.testing.assert_array_almost_equal(ans2.num, [[[2.5, 5., 7.5]]])
         np.testing.assert_array_almost_equal(ans2.den, [[[1., -0.5, -2.]]])
 
+        # Make sure default arugments work as well
+        ans3 = feedback(self.sys1, 1)
+        ans4 = feedback(self.sys1)
+        np.testing.assert_array_almost_equal(ans3.num, ans4.num)
+        np.testing.assert_array_almost_equal(ans3.den, ans4.den)
+
     def testSSScalar(self):
         """State space system with scalar feedback block."""
         

Modified: trunk/tests/frd_test.py
===================================================================
--- trunk/tests/frd_test.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/tests/frd_test.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -84,6 +84,11 @@
         np.testing.assert_array_almost_equal(
             f1.feedback(1).freqresp([0.1, 1.0, 10])[0],
             h1.feedback(1).freqresp([0.1, 1.0, 10])[0])
+
+        # Make sure default argument also works
+        np.testing.assert_array_almost_equal(
+            f1.feedback().freqresp([0.1, 1.0, 10])[0],
+            h1.feedback().freqresp([0.1, 1.0, 10])[0])
         
     def testFeedback2(self):
         h2 = StateSpace([[-1.0, 0], [0, -2.0]], [[0.4], [0.1]],

Modified: trunk/tests/xferfcn_test.py
===================================================================
--- trunk/tests/xferfcn_test.py	2013-06-09 14:57:59 UTC (rev 274)
+++ trunk/tests/xferfcn_test.py	2013-06-09 15:50:50 UTC (rev 275)
@@ -324,6 +324,11 @@
         np.testing.assert_array_almost_equal(sys.evalfr(32.),
             np.array([[0.00281959302585077 - 0.030628473607392j]]))
 
+        # Test call version as well
+        np.testing.assert_almost_equal(sys(1.j), -0.5 - 0.5j)
+        np.testing.assert_almost_equal(sys(32.j), 
+            0.00281959302585077 - 0.030628473607392j)
+
     def testEvalFrMIMO(self):
         """Evaluate the frequency response of a MIMO system at one frequency."""
 
@@ -338,6 +343,9 @@
         
         np.testing.assert_array_almost_equal(sys.evalfr(2.), resp)
 
+        # Test call version as well
+        np.testing.assert_array_almost_equal(sys(2.j), resp)
+
     # Tests for TransferFunction.freqresp.
 
     def testFreqRespSISO(self):

[python-control-discuss] SF.net python-control commit:[274] trunk

[<mur...@us...>]

Revision: 274
          http://sourceforge.net/p/python-control/code/274
Author:   murrayrm
Date:     2013-06-09 14:57:59 +0000 (Sun, 09 Jun 2013)
Log Message:
-----------
merging changes from Rane van Paassen - see 2013-05-23 ChangeLog entry for details

Modified Paths:
--------------
    trunk/ChangeLog
    trunk/src/bdalg.py
    trunk/src/frdata.py
    trunk/src/lti.py
    trunk/src/margins.py
    trunk/src/matlab.py
    trunk/src/modelsimp.py
    trunk/src/rlocus.py
    trunk/src/statesp.py
    trunk/src/timeresp.py
    trunk/src/xferfcn.py
    trunk/tests/convert_test.py
    trunk/tests/matlab_test.py
    trunk/tests/slycot_convert_test.py
    trunk/tests/statesp_test.py
    trunk/tests/timeresp_test.py
    trunk/tests/xferfcn_test.py

Removed Paths:
-------------
    branches/add-frd/
    branches/mix-matrices-and-sys/

Property Changed:
----------------
    trunk/
    trunk/src/ctrlutil.py

Index: trunk
===================================================================
--- trunk	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk	2013-06-09 14:57:59 UTC (rev 274)

Property changes on: trunk
___________________________________________________________________
Modified: svn:mergeinfo
## -1 +1,2 ##
 /branches/add-frd:188-219
+/branches/mix-matrices-and-sys:249-273
\ No newline at end of property
Modified: trunk/ChangeLog
===================================================================
--- trunk/ChangeLog	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/ChangeLog	2013-06-09 14:57:59 UTC (rev 274)
@@ -1,3 +1,24 @@
+2013-05-23  Rene van Paassen <ren...@gm...>
+
+	* src/margin.py: re-vamped stability_margins function. Now
+          uses analytical method for TF and state-space, works on FRD
+          objects too, and can optionally return all margins found for
+          the analytical method.
+
+	* src/xferfcn.py and src/frdata.py: Now return array result
+          for fresp when called with an iterable frequency vector
+
+	* src/xferfcn.py and src/statesp.py: added minreal functions
+
+	* src/bdalg.py: implemented append and connect functions (old
+          version, with indexes, not with names)
+	  
+	* src/matlab.py: added connect, append, minreal
+
+	* src/xferfcn.py and src/statesp.py: added/corrected
+          __rmul__. In most cases a gain matrix (np.matrix) may now be
+          used in combination with ss or tf systems.
+
 2012-11-10  Richard Murray  <murray@altura.local>
 
 	* src/canonical.py: new module implementing conversions to selected

Modified: trunk/src/bdalg.py
===================================================================
--- trunk/src/bdalg.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/bdalg.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -4,10 +4,12 @@
 
 Routines in this module:
 
+append
 series
 parallel
 negate
 feedback
+connect
 
 """
 
@@ -236,3 +238,97 @@
             sys2 = tf._convertToTransferFunction(sys2)
 
     return sys1.feedback(sys2, sign)
+
+def append(*sys):
+    '''
+    Group models by appending their inputs and outputs
+
+    Forms an augmented system model, and appends the inputs and
+    outputs together. The system type will be the type of the first
+    system given; if you mix state-space systems and gain matrices,
+    make sure the gain matrices are not first.
+
+    Parameters.
+    -----------
+    sys1, sys2, ... sysn: StateSpace or Transferfunction
+        LTI systems to combine
+
+        
+    Returns
+    -------
+    sys: LTI system
+        Combined LTI system, with input/output vectors consisting of all 
+        input/output vectors appended
+        
+    Examples
+    --------
+    >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
+    >>> sys2 = ss("-1.", "1.", "1.", "0.")
+    >>> sys = append(sys1, sys2)
+    
+    .. todo::
+        also implement for transfer function, zpk, etc.
+    '''
+    s1 = sys[0]
+    for s in sys[1:]:
+        s1 = s1.append(s)
+    return s1
+
+def connect(sys, Q, inputv, outputv):
+    '''
+    Index-base interconnection of system
+
+    The system sys is a system typically constructed with append, with
+    multiple inputs and outputs. The inputs and outputs are connected
+    according to the interconnection matrix Q, and then the final
+    inputs and outputs are trimmed according to the inputs and outputs
+    listed in inputv and outputv.
+
+    Note: to have this work, inputs and outputs start counting at 1!!!!
+
+    Parameters.
+    -----------
+    sys: StateSpace Transferfunction
+        System to be connected
+    Q: 2d array
+        Interconnection matrix. First column gives the input to be connected
+        second column gives the output to be fed into this input. Negative
+        values for the second column mean the feedback is negative, 0 means
+        no connection is made
+    inputv: 1d array
+        list of final external inputs
+    outputv: 1d array
+        list of final external outputs
+
+    Returns
+    -------
+    sys: LTI system
+        Connected and trimmed LTI system
+
+    Examples
+    --------
+    >>> sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
+    >>> sys2 = ss("-1.", "1.", "1.", "0.")
+    >>> sys = append(sys1, sys2)
+    >>> Q = sp.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1
+    >>> sysc = connect(sys, Q, [2], [1, 2])
+    '''
+    # first connect
+    K = sp.zeros( (sys.inputs, sys.outputs) )
+    for r in sp.array(Q).astype(int):
+        inp = r[0]-1
+        for outp in r[1:]:
+            if outp > 0 and outp <= sys.outputs:
+                K[inp,outp-1] = 1.
+            elif outp < 0 and -outp >= -sys.outputs:
+                K[inp,-outp-1] = -1.
+    sys = sys.feedback(sp.matrix(K), sign=1)
+    
+    # now trim
+    Ytrim = sp.zeros( (len(outputv), sys.outputs) )
+    Utrim = sp.zeros( (sys.inputs, len(inputv)) )
+    for i,u in enumerate(inputv):
+        Utrim[u-1,i] = 1.
+    for i,y in enumerate(outputv):
+        Ytrim[i,y-1] = 1.
+    return sp.matrix(Ytrim)*sys*sp.matrix(Utrim)  

Index: trunk/src/ctrlutil.py
===================================================================
--- trunk/src/ctrlutil.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/ctrlutil.py	2013-06-09 14:57:59 UTC (rev 274)

Property changes on: trunk/src/ctrlutil.py
___________________________________________________________________
Modified: svn:mergeinfo
## -1,2 +1,3 ##
 /branches/add-frd/src/ctrlutil.py:188-219
+/branches/mix-matrices-and-sys/src/ctrlutil.py:249-273
 /trunk/src/ctrlutil.py:140-142
\ No newline at end of property
Modified: trunk/src/frdata.py
===================================================================
--- trunk/src/frdata.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/frdata.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -107,7 +107,7 @@
     def __init__(self, *args, **kwargs):
         """Construct a transfer function.
         
-        The default constructor is FRD(w, d), where w is an iterable of 
+        The default constructor is FRD(d, w), where w is an iterable of 
         frequency points, and d is the matching frequency data. 
         
         If d is a single list, 1d array, or tuple, a SISO system description 
@@ -142,8 +142,8 @@
                 # The user provided a response and a freq vector
                 self.fresp = array(args[0], dtype=complex)
                 if len(self.fresp.shape) == 1:
-                    self.fresp.reshape(1, 1, len(args[0]))
-                self.omega = array(args[1])
+                    self.fresp = self.fresp.reshape(1, 1, len(args[0]))
+                self.omega = array(args[1], dtype=float)
                 if len(self.fresp.shape) != 3 or \
                         self.fresp.shape[-1] != self.omega.shape[-1] or \
                         len(self.omega.shape) != 1:
@@ -189,9 +189,9 @@
                 if mimo:
                     outstr.append("Input %i to output %i:" % (i + 1, j + 1))
                 outstr.append('Freq [rad/s]  Response   ')
-                outstr.append('------------  ------------------------')
+                outstr.append('------------  ---------------------')
                 outstr.extend(
-                    [ '%12.3f  %10.4g + %10.4g' % (w, m, p)
+                    [ '%12.3f  %10.4g%+10.4gj' % (w, m, p)
                       for m, p, w in zip(real(self.fresp[j,i,:]), imag(self.fresp[j,i,:]), wt) ])
 
 
@@ -340,7 +340,10 @@
         """
 
         # Preallocate the output.
-        out = empty((self.outputs, self.inputs), dtype=complex)
+        if getattr(omega, '__iter__', False):
+            out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
+        else:
+            out = empty((self.outputs, self.inputs), dtype=complex)
 
         if self.ifunc is None:
             try:
@@ -350,11 +353,18 @@
                     "Frequency %f not in frequency list, try an interpolating"
                     " FRD if you want additional points")
         else:
-            for i in range(self.outputs):
-                for j in range(self.inputs):
-                    frraw = splev(omega, self.ifunc[i,j], der=0)
-                    out[i,j] = frraw[0] + 1.0j*frraw[1]
-
+            if getattr(omega, '__iter__', False):
+                for i in range(self.outputs):
+                    for j in range(self.inputs):
+                        for k,w in enumerate(omega): 
+                            frraw = splev(w, self.ifunc[i,j], der=0)
+                            out[i,j,k] = frraw[0] + 1.0j*frraw[1]
+            else:
+                for i in range(self.outputs):
+                    for j in range(self.inputs):
+                        frraw = splev(omega, self.ifunc[i,j], der=0)
+                        out[i,j] = frraw[0] + 1.0j*frraw[1]
+                
         return out
 
     # Method for generating the frequency response of the system

Modified: trunk/src/lti.py
===================================================================
--- trunk/src/lti.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/lti.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -12,6 +12,8 @@
 timebaseEqual()
 """
 
+from numpy import absolute, real
+
 class Lti:
     """Lti is a parent class to linear time invariant control (LTI) objects.
     
@@ -64,6 +66,12 @@
         self.outputs = outputs
         self.dt = dt
 
+    def damp(self):
+        poles = self.pole()
+        wn = absolute(poles)
+        Z = -real(poles)/wn
+        return wn, Z, poles
+
 # Test to see if a system is SISO
 def issiso(sys, strict=False):
     if isinstance(sys, (int, float, complex)) and not strict:

Modified: trunk/src/margins.py
===================================================================
--- trunk/src/margins.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/margins.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -53,15 +53,37 @@
 import numpy as np
 import control.xferfcn as xferfcn
 from control.freqplot import bode
-from control.lti import isdtime
+from control.lti import isdtime, issiso
+import control.frdata as frdata
+import scipy as sp
 
-# gain and phase margins
-# contributed by Sawyer B. Fuller <mi...@ca...>
-#! TODO - need to add unit test functions
-def stability_margins(sysdata, deg=True):
+# helper functions for stability_margins
+def _polyimsplit(pol):
+    """split a polynomial with (iw) applied into a real and an
+       imaginary part with w applied"""
+    rpencil = np.zeros_like(pol)
+    ipencil = np.zeros_like(pol)
+    rpencil[-1::-4] = 1.
+    rpencil[-3::-4] = -1.
+    ipencil[-2::-4] = 1.
+    ipencil[-4::-4] = -1.
+    return pol * rpencil, pol*ipencil
+
+def _polysqr(pol):
+    """return a polynomial squared"""
+    return np.polymul(pol, pol)
+
+# Took the framework for the old function by 
+# Sawyer B. Fuller <mi...@ca...>, removed a lot of the innards
+# and replaced with analytical polynomial functions for Lti systems.
+#
+# idea for the frequency data solution copied/adapted from
+# https://github.com/alchemyst/Skogestad-Python/blob/master/BODE.py
+# Rene van Paassen <ren...@gm...>
+def stability_margins(sysdata, deg=True, returnall=False):
     """Calculate gain, phase and stability margins and associated
     crossover frequencies.
-
+    
     Usage
     -----
     gm, pm, sm, wg, wp, ws = stability_margins(sysdata, deg=True)
@@ -76,86 +98,117 @@
             bode frequency response data 
     deg=True: boolean  
         If true, all input and output phases in degrees, else in radians
-        
+    returnall=False: boolean
+        If true, return all margins found. Note that for frequency data or
+        FRD systems, only one margin is found and returned. 
+       
     Returns
     -------
-    gm, pm, sm, wg, wp, ws: float
+    gm, pm, sm, wg, wp, ws: float or array_like
         Gain margin gm, phase margin pm, stability margin sm, and 
         associated crossover
         frequencies wg, wp, and ws of SISO open-loop. If more than
         one crossover frequency is detected, returns the lowest corresponding
         margin. 
-    """
-    #TODO do this precisely without the effects of discretization of frequencies?
-    #TODO assumes SISO
-    #TODO unit tests, margin plot
+        When requesting all margins, the return values are array_like, 
+        and all margins are returns for linear systems not equal to FRD
+        """
+    
+    try:
+        if isinstance(sysdata, frdata.FRD):
+            sys = frdata.FRD(sysdata, smooth=True) 
+        elif isinstance(sysdata, xferfcn.TransferFunction):
+            sys = sysdata
+        elif getattr(sysdata, '__iter__', False) and len(sysdata) == 3:
+            mag, phase, omega = sysdata
+            sys = frdata.FRD(mag*np.exp((1j/360.)*phase), omega, smooth=True)
+        else:
+            sys = xferfcn._convertToTransferFunction(sysdata)
+    except Exception, e:
+        print (e)
+        raise ValueError("Margin sysdata must be either a linear system or "
+                         "a 3-sequence of mag, phase, omega.")
 
-    if (not getattr(sysdata, '__iter__', False)):
-        sys = sysdata
+    # calculate gain of system
+    if isinstance(sys, xferfcn.TransferFunction):
+        
+        # check for siso
+        if not issiso(sys):
+            raise ValueError("Can only do margins for SISO system")
+        
+        # real and imaginary part polynomials in omega:
+        rnum, inum = _polyimsplit(sys.num[0][0])
+        rden, iden = _polyimsplit(sys.den[0][0])
 
-        # TODO: implement for discrete time systems
-        if (isdtime(sys, strict=True)):
-            raise NotImplementedError("Function not implemented in discrete time")
+        # test imaginary part of tf == 0, for phase crossover/gain margins
+        test_w_180 = np.polyadd(np.polymul(inum, rden), np.polymul(rnum, -iden))
+        w_180 = np.roots(test_w_180)
+        w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 > 0.)])
+        w_180.sort()
 
-        mag, phase, omega = bode(sys, deg=deg, Plot=False)
-    elif len(sysdata) == 3:
-        # TODO: replace with FRD object type?
-        mag, phase, omega = sysdata
-    else: 
-        raise ValueError("Margin sysdata must be either a linear system or a 3-sequence of mag, phase, omega.")
-        
-    if deg:
-        cycle = 360. 
-        crossover = 180. 
+        # test magnitude is 1 for gain crossover/phase margins
+        test_wc = np.polysub(np.polyadd(_polysqr(rnum), _polysqr(inum)), 
+                             np.polyadd(_polysqr(rden), _polysqr(iden)))
+        wc = np.roots(test_wc)
+        wc = np.real(wc[(np.imag(wc) == 0) * (wc > 0.)])
+        wc.sort()
+
+        # stability margin was a bitch to elaborate, relies on magnitude to
+        # point -1, then take the derivative. Second derivative needs to be >0
+        # to have a minimum
+        test_wstabn = np.polyadd(_polysqr(rnum), _polysqr(inum))
+        test_wstabd = np.polyadd(_polysqr(np.polyadd(rnum,rden)), 
+                                 _polysqr(np.polyadd(inum,iden)))
+        test_wstab = np.polysub(
+            np.polymul(np.polyder(test_wstabn),test_wstabd), 
+            np.polymul(np.polyder(test_wstabd),test_wstabn))
+
+        # find the solutions
+        wstab = np.roots(test_wstab)
+
+        # and find the value of the 2nd derivative there, needs to be positive
+        wstabplus = np.polyval(np.polyder(test_wstab), wstab)
+        wstab = np.real(wstab[(np.imag(wstab) == 0) * (wstab > 0.) *
+                              (np.abs(wstabplus) > 0.)])
+        wstab.sort()
+
     else:
-        cycle = 2 * np.pi
-        crossover = np.pi
-        
-    wrapped_phase = -np.mod(phase, cycle)
-    
-    # phase margin from minimum phase among all gain crossovers
-    neg_mag_crossings_i = np.nonzero(np.diff(mag < 1) > 0)[0]
-    mag_crossings_p = wrapped_phase[neg_mag_crossings_i]
-    if len(neg_mag_crossings_i) == 0:
-        if mag[0] < 1: # gain always less than one
-            wp = np.nan
-            pm = np.inf
-        else: # gain always greater than one
-            print("margin: no magnitude crossings found")
-            wp = np.nan
-            pm = np.nan
-    else:
-        min_mag_crossing_i = neg_mag_crossings_i[np.argmin(mag_crossings_p)]
-        wp = omega[min_mag_crossing_i]
-        pm = crossover + phase[min_mag_crossing_i] 
-        if pm < 0:
-            print("warning: system unstable: negative phase margin")
-    
-    # gain margin from minimum gain margin among all phase crossovers
-    neg_phase_crossings_i = np.nonzero(np.diff(wrapped_phase < -crossover) > 0)[0]
-    neg_phase_crossings_g = mag[neg_phase_crossings_i]
-    if len(neg_phase_crossings_i) == 0:
-        wg = np.nan
-        gm = np.inf
-    else:
-        min_phase_crossing_i = neg_phase_crossings_i[
-            np.argmax(neg_phase_crossings_g)]
-        wg = omega[min_phase_crossing_i]
-        gm = abs(1/mag[min_phase_crossing_i])
-        if gm < 1: 
-            print("warning: system unstable: gain margin < 1")
+        # a bit coarse, have the interpolated frd evaluated again
+        def mod(w):
+            """to give the function to calculate |G(jw)| = 1"""
+            return [np.abs(sys.evalfr(w[0])[0][0]) - 1]
 
-    # stability margin from minimum abs distance from -1 point
-    if deg:
-        phase_rad = phase * np.pi / 180.
+        def arg(w):
+            """function to calculate the phase angle at -180 deg"""
+            return [np.angle(sys.evalfr(w[0])[0][0]) + np.pi]
+
+        def dstab(w):
+            """function to calculate the distance from -1 point"""
+            return np.abs(sys.evalfr(w[0])[0][0] + 1.)
+
+        # how to calculate the frequency at which |G(jw)| = 1
+        wc = np.array([sp.optimize.fsolve(mod, sys.omega[0])])[0]
+        w_180 = np.array([sp.optimize.fsolve(arg, sys.omega[0])])[0]
+        wstab = np.real(
+            np.array([sp.optimize.fmin(dstab, sys.omega[0], disp=0)])[0])
+
+    # margins, as iterables, converted frdata and xferfcn calculations to
+    # vector for this
+    PM = np.angle(sys.evalfr(wc)[0][0], deg=True) + 180
+    GM = 1/(np.abs(sys.evalfr(w_180)[0][0]))
+    SM = np.abs(sys.evalfr(wstab)[0][0]+1)
+
+    if returnall:
+        return GM, PM, SM, wc, w_180, wstab
     else:
-        phase_rad = phase
-    L = mag * np.exp(1j * phase_rad) # complex loop response to -1 pt
-    min_Lplus1_i = np.argmin(np.abs(L + 1))
-    sm = np.abs(L[min_Lplus1_i] + 1)
-    ws = phase[min_Lplus1_i]
+        return (
+            (GM.shape[0] or None) and GM[0], 
+            (PM.shape[0] or None) and PM[0], 
+            (SM.shape[0] or None) and SM[0], 
+            (wc.shape[0] or None) and wc[0],
+            (w_180.shape[0] or None) and w_180[0],
+            (wstab.shape[0] or None) and wstab[0])
 
-    return gm, pm, sm, wg, wp, ws 
 
 # Contributed by Steffen Waldherr <wal...@is...>
 #! TODO - need to add test functions

Modified: trunk/src/matlab.py
===================================================================
--- trunk/src/matlab.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/matlab.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -89,6 +89,7 @@
 from control.statesp import StateSpace, _rss_generate, _convertToStateSpace
 from control.xferfcn import TransferFunction, _convertToTransferFunction
 from control.lti import Lti #base class of StateSpace, TransferFunction
+from control.frdata import FRD
 from control.dtime import sample_system
 from control.exception import ControlArgument
 
@@ -96,11 +97,11 @@
 from control.ctrlutil import unwrap
 from control.freqplot import nyquist, gangof4
 from control.nichols import nichols
-from control.bdalg import series, parallel, negate, feedback
+from control.bdalg import series, parallel, negate, feedback, append, connect
 from control.pzmap import pzmap
 from control.statefbk import ctrb, obsv, gram, place, lqr
 from control.delay import pade
-from control.modelsimp import hsvd, balred, modred
+from control.modelsimp import hsvd, balred, modred, minreal
 from control.mateqn import lyap, dlyap, dare, care
 
 __doc__ += r"""
@@ -125,7 +126,7 @@
 \*  :func:`ss`                  create state-space (SS) models
 \   dss                         create descriptor state-space models
 \   delayss                     create state-space models with delayed terms
-\   frd                         create frequency response data (FRD) models
+\*  :func:`frd`                 create frequency response data (FRD) models
 \   lti/exp                     create pure continuous-time delays (TF and
                                 ZPK only)
 \   filt                        specify digital filters
@@ -156,7 +157,7 @@
 \*  :func:`tf`                  conversion to transfer function
 \   zpk                         conversion to zero/pole/gain
 \*  :func:`ss`                  conversion to state space
-\   frd                         conversion to frequency data
+\*  :func:`frd`                 conversion to frequency data
 \   c2d                         continuous to discrete conversion
 \   d2c                         discrete to continuous conversion
 \   d2d                         resample discrete-time model
@@ -174,7 +175,7 @@
 ----------------------------------------------------------------------------
 
 ==  ==========================  ============================================
-\   append                      group LTI models by appending inputs/outputs
+\*  :func:`~bdalg.append`       group LTI models by appending inputs/outputs
 \*  :func:`~bdalg.parallel`     connect LTI models in parallel 
                                 (see also overloaded ``+``)
 \*  :func:`~bdalg.series`       connect LTI models in series 
@@ -200,7 +201,7 @@
 \   lti/order                   model order (number of states)
 \*  :func:`~pzmap.pzmap`        pole-zero map (TF only)
 \   lti/iopzmap                 input/output pole-zero map
-\   damp                        natural frequency, damping of system poles
+\*  :func:`damp`                natural frequency, damping of system poles
 \   esort                       sort continuous poles by real part
 \   dsort                       sort discrete poles by magnitude
 \   lti/stabsep                 stable/unstable decomposition
@@ -243,7 +244,7 @@
 ----------------------------------------------------------------------------
 
 ==  ==========================  ============================================
-\   minreal                     minimal realization; pole/zero cancellation
+\*  :func:`~modelsimp.minreal`  minimal realization; pole/zero cancellation
 \   ss/sminreal                 structurally minimal realization
 \*  :func:`~modelsimp.hsvd`     hankel singular values (state contributions)
 \*  :func:`~modelsimp.balred`   reduced-order approximations of LTI models
@@ -569,7 +570,43 @@
     else:
         raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
+def frd(*args):
+    '''
+    Construct a Frequency Response Data model, or convert a system
 
+    frd models store the (measured) frequency response of a system.
+
+    This function can be called in different ways:
+
+    ``frd(response, freqs)``
+        Create an frd model with the given response data, in the form of
+        complex response vector, at matching frequency freqs [in rad/s]
+
+    ``frd(sys, freqs)``
+        Convert an Lti system into an frd model with data at frequencies
+        freqs. 
+
+    Parameters
+    ----------
+    response: array_like, or list
+        complex vector with the system response
+    freq: array_lik or lis
+        vector with frequencies
+    sys: Lti (StateSpace or TransferFunction)
+        A linear system
+
+    Returns
+    -------
+    sys: FRD
+        New frequency response system
+
+    See Also
+    --------
+    ss, tf
+    '''
+    return FRD(*args)
+
+
 def ss2tf(*args):
     """
     Transform a state space system to a transfer function.
@@ -714,7 +751,7 @@
     else:
         raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))
 
-def rss(states=1, inputs=1, outputs=1):
+def rss(states=1, outputs=1, inputs=1):
     """
     Create a stable **continuous** random state space object.
     
@@ -751,7 +788,7 @@
     
     return _rss_generate(states, inputs, outputs, 'c')
     
-def drss(states=1, inputs=1, outputs=1):
+def drss(states=1, outputs=1, inputs=1):
     """
     Create a stable **discrete** random state space object.
     
@@ -852,16 +889,20 @@
 
     return sys.zero()
 
-def evalfr(sys, omega):
+def evalfr(sys, x):
     """
-    Evaluate the transfer function of an LTI system at an angular frequency.
+    Evaluate the transfer function of an LTI system for a single complex 
+    number x.
+    
+    To evaluate at a frequency, enter x = omega*j, where omega is the 
+    frequency in radians
 
     Parameters
     ----------
     sys: StateSpace or TransferFunction
         Linear system
-    omega: scalar
-        Frequency 
+    x: scalar
+        Complex number 
 
     Returns
     -------
@@ -880,15 +921,17 @@
     Examples
     --------
     >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
-    >>> evalfr(sys, 1.)
+    >>> evalfr(sys, 1j)
     array([[ 44.8-21.4j]])
     >>> # This is the transfer function matrix evaluated at s = i.
 
     .. todo:: Add example with MIMO system
     """
+    if issiso(sys):
+        return sys.horner(x)[0][0]
+    return sys.horner(x)
+        
 
-    return sys.evalfr(omega)
-
 def freqresp(sys, omega): 
     """
     Frequency response of an LTI system at multiple angular frequencies.
@@ -1057,10 +1100,11 @@
 
     Returns
     -------
-    gm, pm, wg, wp : float
-        Gain margin gm, phase margin pm (in deg), and associated crossover
-        frequencies wg and wp (in rad/sec) of SISO open-loop. If more than
-        one crossover frequency is detected, returns the lowest
+    gm, pm, Wcg, Wcp : float
+        Gain margin gm, phase margin pm (in deg), gain crossover frequency 
+        (corresponding to phase margin) and phase crossover frequency
+        (corresponding to gain margin), in rad/sec of SISO open-loop.
+        If more than one crossover frequency is detected, returns the lowest
         corresponding margin.
 
     Examples
@@ -1083,7 +1127,7 @@
         raise ValueError("Margin needs 1 or 3 arguments; received %i." 
             % len(args))
             
-    return margin[0], margin[1], margin[3], margin[4]
+    return margin[0], margin[1], margin[4], margin[3]
 
 def dcgain(*args):
     '''
@@ -1139,10 +1183,48 @@
     gain = sys.D - sys.C * sys.A.I * sys.B
     return gain
 
+def damp(sys, doprint=True):
+    '''
+    Compute natural frequency, damping and poles of a system
+    
+    The function takes 1 or 2 parameters
+
+    Parameters
+    ----------
+    sys: Lti (StateSpace or TransferFunction)
+        A linear system object
+    doprint: 
+        if true, print table with values
+
+    Returns
+    -------
+    wn: array
+        Natural frequencies of the poles
+    damping: array
+        Damping values
+    poles: array
+        Pole locations
+
+    See Also
+    --------
+    pole        
+    '''
+    wn, damping, poles = sys.damp()
+    if doprint:
+        print('_____Eigenvalue______ Damping___ Frequency_')
+        for p, d, w in zip(poles, damping, wn) :
+            if abs(p.imag) < 1e-12:
+                print("%10.4g            %10.4g %10.4g" % 
+                      (p.real, 1.0, -p.real)) 
+            else:
+                print("%10.4g%+10.4gj %10.4g %10.4g" % 
+                      (p.real, p.imag, d, w)) 
+    return wn, damping, poles
+
 # Simulation routines 
 # Call corresponding functions in timeresp, with arguments transposed
 
-def step(sys, T=None, X0=0., input=0, output=0, **keywords):
+def step(sys, T=None, X0=0., input=0, output=None, **keywords):
     '''
     Step response of a linear system
     
@@ -1192,7 +1274,7 @@
 
     Examples
     --------
-    >>> T, yout = step(sys, T, X0)
+    >>> yout, T = step(sys, T, X0)
     '''
     T, yout = timeresp.step_response(sys, T, X0, input, output, 
                                    transpose = True, **keywords)
@@ -1405,3 +1487,4 @@
     # TODO: add docstring
     #  Call the sample_system() function to do the work
     return sample_system(sysc, Ts, method)
+

Modified: trunk/src/modelsimp.py
===================================================================
--- trunk/src/modelsimp.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/modelsimp.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -269,6 +269,33 @@
 
     return rsys
 
+def minreal(sys, tol=None, verbose=True):
+    '''
+    Eliminates uncontrollable or unobservable states in state-space
+    models or cancelling pole-zero pairs in transfer functions. The
+    output sysr has minimal order and the same response
+    characteristics as the original model sys.
+
+    Parameters
+    ----------
+    sys: StateSpace or TransferFunction
+        Original system
+    tol: real
+        Tolerance
+    verbose: bool
+        Print results if True
+
+    Returns
+    -------
+    rsys: StateSpace or TransferFunction
+        Cleaned model
+    '''
+    sysr = sys.minreal(tol)
+    if verbose:
+        print("{nstates} states have been removed from the model".format(
+                nstates=len(sys.pole()) - len(sysr.pole())))
+    return sysr
+
 def era(YY, m, n, nin, nout, r):
     """
     Calculate an ERA model of order `r` based on the impulse-response data `YY`.

Modified: trunk/src/rlocus.py
===================================================================
--- trunk/src/rlocus.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/rlocus.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -50,9 +50,11 @@
 import scipy.signal             # signal processing toolbox
 import pylab                    # plotting routines
 import control.xferfcn as xferfcn
+from functools import partial
 
 # Main function: compute a root locus diagram
-def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True):
+def root_locus(sys, kvect, xlim=None, ylim=None, plotstr='-', Plot=True, 
+               PrintGain=True):
     """Calculate the root locus by finding the roots of 1+k*TF(s)
     where TF is self.num(s)/self.den(s) and each k is an element
     of kvect.
@@ -65,7 +67,9 @@
         List of gains to use in computing diagram
     Plot : boolean (default = True)
         If True, plot magnitude and phase
-
+    PrintGain: boolean (default = True)
+        If True, report mouse clicks when close to the root-locus branches,
+        calculate gain, damping and print
     Return values
     -------------
     rlist : list of computed root locations
@@ -80,6 +84,10 @@
 
     # Create the plot
     if (Plot):
+        f = pylab.figure()
+        if PrintGain:
+            cid = f.canvas.mpl_connect(
+                'button_release_event', partial(_RLFeedbackClicks, sys=sys))
         ax = pylab.axes();
 
         # plot open loop poles
@@ -165,3 +173,12 @@
                 sorted[n,ind] = elem
         prevrow = sorted[n,:]
     return sorted
+
+def _RLFeedbackClicks(event, sys):
+    """Print root-locus gain feedback for clicks on the root-locus plot
+    """
+    s = complex(event.xdata, event.ydata)
+    K = -1./sys.horner(s)
+    if abs(K.real) > 1e-8 and abs(K.imag/K.real) < 0.04:
+        print("Clicked at %10.4g%+10.4gj gain %10.4g damp %10.4g" %
+              (s.real, s.imag, K.real, -1*s.real/abs(s)))

Modified: trunk/src/statesp.py
===================================================================
--- trunk/src/statesp.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/statesp.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -27,6 +27,7 @@
 StateSpace.zero
 StateSpace.feedback
 StateSpace.returnScipySignalLti
+StateSpace.append
 _convertToStateSpace
 _rss_generate
 
@@ -338,10 +339,23 @@
             B = self.B * other;
             D = self.D * other;
             return StateSpace(A, B, C, D, self.dt)
+        
+        # is lti, and convertible?
+        if isinstance(other, Lti):
+            return _convertToStateSpace(other) * self
 
-        else:
-            raise TypeError("can't interconnect systems")
+        # try to treat this as a matrix
+        try:
+            X = matrix(other)
+            C = X * self.C
+            D = X * self.D
+            return StateSpace(self.A, self.B, C, D, self.dt)
 
+        except Exception, e:
+            print(e)
+            pass
+        raise TypeError("can't interconnect systems")
+
     # TODO: __div__ and __rdiv__ are not written yet.
     def __div__(self, other):
         """Divide two LTI systems."""
@@ -370,10 +384,16 @@
         else:
             s = omega * 1.j
 
-        fresp = self.C * solve(s * eye(self.states) - self.A,
-            self.B) + self.D
+        return self.horner(s)
 
-        return array(fresp)
+    def horner(self, s):
+        '''Evaluate the systems's transfer function for a complex variable
+        
+        Returns a matrix of values evaluated at complex variable s.
+        '''
+        resp = self.C * solve(s * eye(self.states) - self.A,
+                              self.B) + self.D
+        return array(resp)
 
     # Method for generating the frequency response of the system
     # TODO: add discrete time check
@@ -475,6 +495,17 @@
 
         return StateSpace(A, B, C, D, dt)
 
+    def minreal(self, tol=None):
+        """Calculate a minimal realization, removes unobservable and
+        uncontrollable states"""
+        try:
+            from slycot import tb01pd
+            A, B, C, nr = tb01pd(self.states, self.inputs, self.outputs, 
+                                    self.A, self.B, self.C, tol=tol)
+            return StateSpace(A[:nr,:nr], B[:nr,:], C[:,:nr], self.D)
+        except ImportError:
+            raise TypeError("minreal requires slycot tb01pd")
+       
     # TODO: add discrete time check
     def returnScipySignalLti(self):
         """Return a list of a list of scipy.signal.lti objects.
@@ -497,6 +528,33 @@
 
         return out
 
+    def append(self, other):
+        """Append a second model to the present model. The second
+        model is converted to state-space if necessary, inputs and
+        outputs are appended and their order is preserved"""
+        if not isinstance(other, StateSpace):
+            other = _convertToStateSpace(other)
+        
+        if self.dt != other.dt:
+            raise ValueError("Systems must have the same time step")
+
+        n = self.states + other.states
+        m = self.inputs + other.inputs
+        p = self.outputs + other.outputs
+        A = zeros( (n, n) )
+        B = zeros( (n, m) )
+        C = zeros( (p, n) )
+        D = zeros( (p, m) )
+        A[:self.states,:self.states] = self.A
+        A[self.states:,self.states:] = other.A
+        B[:self.states,:self.inputs] = self.B
+        B[self.states:,self.inputs:] = other.B
+        C[:self.outputs,:self.states] = self.C
+        C[self.outputs:,self.states:] = other.C
+        D[:self.outputs,:self.inputs] = self.D
+        D[self.outputs:,self.inputs:] = other.D
+        return StateSpace(A, B, C, D, self.dt)
+
 # TODO: add discrete time check
 def _convertToStateSpace(sys, **kw):
     """Convert a system to state space form (if needed).
@@ -757,3 +815,51 @@
         
     return sys
 
+def _mimo2simo(sys, input, warn_conversion=False):
+    #pylint: disable=W0622
+    """
+    Convert a MIMO system to a SIMO system. (Convert a system with multiple
+    inputs and/or outputs, to a system with a single input but possibly
+    multiple outputs.)
+    
+    The input that is used in the SIMO system can be selected with the
+    parameter ``input``. All other inputs are set to 0, all other
+    outputs are ignored.
+    
+    If ``sys`` is already a SIMO system, it will be returned unaltered. 
+    
+    Parameters
+    ----------
+    sys: StateSpace
+        Linear (MIMO) system that should be converted.
+    input: int
+        Index of the input that will become the SIMO system's only input.
+    warn_conversion: bool
+        If True: print a warning message when sys is a MIMO system. 
+        Warn that a conversion will take place.
+        
+    Returns:
+    --------
+    sys: StateSpace
+        The converted (SIMO) system.
+    """
+    if not (isinstance(input, int)):
+        raise TypeError("Parameter ``input`` be an integer number.")
+    if not (0 <= input < sys.inputs):
+        raise ValueError("Selected input does not exist. "
+                         "Selected input: {sel}, "
+                         "number of system inputs: {ext}."
+                         .format(sel=input, ext=sys.inputs))
+    #Convert sys to SISO if necessary
+    if sys.inputs > 1:
+        if warn_conversion:
+            warnings.warn("Converting MIMO system to SIMO system. "
+                          "Only input {i} is used."
+                          .format(i=input))
+        # $X = A*X + B*U
+        #  Y = C*X + D*U
+        new_B = sys.B[:, input]
+        new_D = sys.D[:, input]
+        sys = StateSpace(sys.A, new_B, sys.C, new_D, sys.dt)
+        
+    return sys

Modified: trunk/src/timeresp.py
===================================================================
--- trunk/src/timeresp.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/timeresp.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -121,7 +121,7 @@
 from copy import deepcopy
 import warnings
 from control.lti import Lti     # base class of StateSpace, TransferFunction
-from control. statesp import StateSpace, _rss_generate, _convertToStateSpace, _mimo2siso
+from control. statesp import StateSpace, _rss_generate, _convertToStateSpace, _mimo2simo, _mimo2siso
 from control.lti import isdtime, isctime
 
 # Helper function for checking array-like parameters
@@ -386,15 +386,15 @@
 
     return T, yout, xout
 
-def step_response(sys, T=None, X0=0., input=0, output=0, \
-                     transpose = False, **keywords):
+def step_response(sys, T=None, X0=0., input=0, output=None,
+                  transpose = False, **keywords):
     #pylint: disable=W0622
     """Step response of a linear system
     
-    If the system has multiple inputs or outputs (MIMO), one input and one 
-    output have to be selected for the simulation. The parameters `input` 
-    and `output` do this. All other inputs are set to 0, all other outputs 
-    are ignored.
+    If the system has multiple inputs or outputs (MIMO), one input has
+    to be selected for the simulation. Optionally, one output may be
+    selected. The parameters `input` and `output` do this. All other
+    inputs are set to 0, all other outputs are ignored.
     
     For information on the **shape** of parameters `T`, `X0` and 
     return values `T`, `yout` see: :ref:`time-series-convention`
@@ -416,7 +416,8 @@
         Index of the input that will be used in this simulation.
 
     output: int
-        Index of the output that will be used in this simulation.
+        Index of the output that will be used in this simulation. Set to None
+        to not trim outputs
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -447,7 +448,10 @@
     >>> T, yout = step_response(sys, T, X0)
     """
     sys = _convertToStateSpace(sys)
-    sys = _mimo2siso(sys, input, output, warn_conversion=True)
+    if output == None:
+        sys = _mimo2simo(sys, input, warn_conversion=True)
+    else:
+        sys = _mimo2siso(sys, input, output, warn_conversion=True)
     if T is None:
         if isctime(sys):
             T = _default_response_times(sys.A, 100)
@@ -459,13 +463,13 @@
     U = np.ones_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, 
-                                    transpose=transpose, **keywords)
+                                     transpose=transpose, **keywords)
 
     return T, yout
 
 
-def initial_response(sys, T=None, X0=0., input=0, output=0, transpose=False,
-                    **keywords):
+def initial_response(sys, T=None, X0=0., input=0, output=None, transpose=False,
+                     **keywords):
     #pylint: disable=W0622
     """Initial condition response of a linear system
     
@@ -494,7 +498,8 @@
         Index of the input that will be used in this simulation.
 
     output: int
-        Index of the output that will be used in this simulation.
+        Index of the output that will be used in this simulation. Set to None
+        to not trim outputs
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -525,7 +530,10 @@
     >>> T, yout = initial_response(sys, T, X0)
     """
     sys = _convertToStateSpace(sys) 
-    sys = _mimo2siso(sys, input, output, warn_conversion=True)
+    if output == None:
+        sys = _mimo2simo(sys, input, warn_conversion=False)
+    else:
+        sys = _mimo2siso(sys, input, output, warn_conversion=False)
 
     # Create time and input vectors; checking is done in forced_response(...)
     # The initial vector X0 is created in forced_response(...) if necessary
@@ -534,11 +542,11 @@
     U = np.zeros_like(T)
 
     T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose,
-                                    **keywords)
+                                     **keywords)
     return T, yout
 
 
-def impulse_response(sys, T=None, X0=0., input=0, output=0,
+def impulse_response(sys, T=None, X0=0., input=0, output=None,
                     transpose=False, **keywords):
     #pylint: disable=W0622
     """Impulse response of a linear system
@@ -568,7 +576,8 @@
         Index of the input that will be used in this simulation.
 
     output: int
-        Index of the output that will be used in this simulation.
+        Index of the output that will be used in this simulation. Set to None
+        to not trim outputs
 
     transpose: bool
         If True, transpose all input and output arrays (for backward
@@ -599,7 +608,10 @@
     >>> T, yout = impulse_response(sys, T, X0) 
     """
     sys = _convertToStateSpace(sys) 
-    sys = _mimo2siso(sys, input, output, warn_conversion=True)
+    if output == None:
+        sys = _mimo2simo(sys, input, warn_conversion=True)
+    else:
+        sys = _mimo2siso(sys, input, output, warn_conversion=True)
     
     # System has direct feedthrough, can't simulate impulse response numerically
     if np.any(sys.D != 0) and isctime(sys):

Modified: trunk/src/xferfcn.py
===================================================================
--- trunk/src/xferfcn.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/src/xferfcn.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -404,7 +404,8 @@
                 for k in range(self.inputs): # Multiply & add.
                     num_summand[k] = polymul(self.num[i][k], other.num[k][j])
                     den_summand[k] = polymul(self.den[i][k], other.den[k][j])
-                    num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
+                    num[i][j], den[i][j] = _addSISO(
+                        num[i][j], den[i][j],
                         num_summand[k], den_summand[k])
         
         return TransferFunction(num, den, dt)
@@ -412,8 +413,48 @@
     def __rmul__(self, other): 
         """Right multiply two LTI objects (serial connection)."""
         
-        return self * other
+        # Convert the second argument to a transfer function.
+        if isinstance(other, (int, float, complex)):
+            other = _convertToTransferFunction(other, inputs=self.inputs, 
+                outputs=self.inputs)
+        else:
+            other = _convertToTransferFunction(other)
+            
+        # Check that the input-output sizes are consistent.
+        if other.inputs != self.outputs:
+            raise ValueError("C = A * B: A has %i column(s) (input(s)), but B \
+has %i row(s)\n(output(s))." % (other.inputs, self.outputs))
 
+        inputs = self.inputs
+        outputs = other.outputs
+        
+        # Figure out the sampling time to use
+        if (self.dt == None and other.dt != None):
+            dt = other.dt       # use dt from second argument
+        elif (other.dt == None and self.dt != None) or (self.dt == other.dt):
+            dt = self.dt        # use dt from first argument
+        else:
+            raise ValueError("Systems have different sampling times")
+
+        # Preallocate the numerator and denominator of the sum.
+        num = [[[0] for j in range(inputs)] for i in range(outputs)]
+        den = [[[1] for j in range(inputs)] for i in range(outputs)]
+        
+        # Temporary storage for the summands needed to find the (i, j)th element
+        # of the product.
+        num_summand = [[] for k in range(other.inputs)]
+        den_summand = [[] for k in range(other.inputs)]
+        
+        for i in range(outputs): # Iterate through rows of product.
+            for j in range(inputs): # Iterate through columns of product.
+                for k in range(other.inputs): # Multiply & add.
+                    num_summand[k] = polymul(other.num[i][k], self.num[k][j])
+                    den_summand[k] = polymul(other.den[i][k], self.den[k][j])
+                    num[i][j], den[i][j] = _addSISO(num[i][j], den[i][j],
+                        num_summand[k], den_summand[k])
+       
+        return TransferFunction(num, den, dt)
+
     # TODO: Division of MIMO transfer function objects is not written yet.
     def __div__(self, other):
         """Divide two LTI objects."""
@@ -486,9 +527,20 @@
                 warn("evalfr: frequency evaluation above Nyquist frequency")
         else:
             s = 1.j * omega
+            
+        return self.horner(s)
 
+    def horner(self, s):
+        '''Evaluate the systems's transfer function for a complex variable
+        
+        Returns a matrix of values evaluated at complex variable s.
+        '''
+
         # Preallocate the output.
-        out = empty((self.outputs, self.inputs), dtype=complex)
+        if getattr(s, '__iter__', False):
+            out = empty((self.outputs, self.inputs, len(s)), dtype=complex)
+        else:
+            out = empty((self.outputs, self.inputs), dtype=complex)
 
         for i in range(self.outputs):
             for j in range(self.inputs):
@@ -620,8 +672,12 @@
                         newzeros.append(z)
                         
                 # keep result
-                num[i][j] = gain * real(poly(newzeros))
+                if len(newzeros):
+                    num[i][j] = gain * real(poly(newzeros))
+                else:
+                    num[i][j] = array([gain])
                 den[i][j] = real(poly(poles))
+                    
 
         # end result
         return TransferFunction(num, den)
@@ -906,6 +962,15 @@
 
     In the latter example, sys's matrix transfer function is [[1., 1., 1.]
                                                               [1., 1., 1.]].
+
+    If sys is an array-like type, then it is converted to a constant-gain
+    transfer function.
+
+    >>> sys = _convertToTransferFunction([[1. 0.], [2. 3.]])
+
+    In this example, the numerator matrix will be 
+       [[[1.0], [0.0]], [[2.0], [3.0]]]
+    and the denominator matrix [[[1.0], [1.0]], [[1.0], [1.0]]]
     
     """
     from control.statesp import StateSpace
@@ -966,5 +1031,16 @@
         den = [[[1] for j in range(inputs)] for i in range(outputs)]
         
         return TransferFunction(num, den)
-    else:
-        raise TypeError("Can't convert given type to TransferFunction system.")
+
+    # If this is array-like, try to create a constant feedthrough
+    try:
+        D = array(sys)
+        outputs, inputs = D.shape
+        num = [[[D[i,j]] for j in range(inputs)] for i in range(outputs)]
+        den = [[[1] for j in range(inputs)] for i in range(outputs)]
+        return TransferFunction(num, den)
+    except Exception, e: 
+        print("Failure to assume argument is matrix-like in" 
+              " _convertToTransferFunction, result %s" % e)
+        
+    raise TypeError("Can't convert given type to TransferFunction system.")

Modified: trunk/tests/convert_test.py
===================================================================
--- trunk/tests/convert_test.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/convert_test.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -57,7 +57,7 @@
                 for outputs in range(1, self.maxIO):
                     # start with a random SS system and transform to TF then
                     # back to SS, check that the matrices are the same.
-                    ssOriginal = matlab.rss(states, inputs, outputs)
+                    ssOriginal = matlab.rss(states, outputs, inputs)
                     if (verbose):
                         self.printSys(ssOriginal, 1)
 

Modified: trunk/tests/matlab_test.py
===================================================================
--- trunk/tests/matlab_test.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/matlab_test.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -11,9 +11,48 @@
 from __future__ import print_function
 import unittest
 import numpy as np
+from scipy.linalg import eigvals
 import scipy as sp
 from control.matlab import *
+from control.frdata import FRD
 
+# for running these through Matlab or Octave
+'''
+siso_ss1 = ss([1. -2.; 3. -4.], [5.; 7.], [6. 8.], [0])
+
+siso_tf1 = tf([1], [1, 2, 1])
+siso_tf2 = tf([1, 1], [1, 2, 3, 1])
+
+siso_tf3 = tf(siso_ss1)
+siso_ss2 = ss(siso_tf2)
+siso_ss3 = ss(siso_tf3)
+siso_tf4 = tf(siso_ss2)
+
+A =[ 1. -2. 0.  0.;
+     3. -4. 0.  0.;
+     0.  0. 1. -2.;
+     0.  0. 3. -4. ]
+B = [ 5. 0.;
+      7. 0.;
+      0. 5.;
+      0. 7. ]
+C = [ 6. 8. 0. 0.;
+      0. 0. 6. 8. ]
+D = [ 9. 0.;
+      0. 9. ]
+mimo_ss1 = ss(A, B, C, D)
+
+% all boring, since no cross-over
+margin(siso_tf1)
+margin(siso_tf2)
+margin(siso_ss1)
+margin(siso_ss2)
+
+% make a bit better
+[gm, pm, gmc, pmc] = margin(siso_ss2*siso_ss2*2)
+
+'''
+
 class TestMatlab(unittest.TestCase):
     def setUp(self):
         """Set up some systems for testing out MATLAB functions"""
@@ -194,6 +233,9 @@
         gm, pm, wg, wp = margin(self.siso_tf2);
         gm, pm, wg, wp = margin(self.siso_ss1);
         gm, pm, wg, wp = margin(self.siso_ss2);
+        gm, pm, wg, wp = margin(self.siso_ss2*self.siso_ss2*2);
+        np.testing.assert_array_almost_equal(
+            [gm, pm, wg, wp], [1.5451, 75.9933, 1.2720, 0.6559], decimal=3)
         
     def testDcgain(self):
         #Create different forms of a SISO system
@@ -274,13 +316,16 @@
         freqresp(self.siso_tf3, w)
 
     def testEvalfr(self):
-        w = 1
-        evalfr(self.siso_ss1, w)
+        w = 1j
+        self.assertEqual(evalfr(self.siso_ss1, w), 44.8-21.4j)
         evalfr(self.siso_ss2, w)
         evalfr(self.siso_ss3, w)
         evalfr(self.siso_tf1, w)
         evalfr(self.siso_tf2, w)
         evalfr(self.siso_tf3, w)
+        np.testing.assert_array_almost_equal(
+            evalfr(self.mimo_ss1, w), 
+            np.array( [[44.8-21.4j, 0.], [0., 44.8-21.4j]]))
 
     def testHsvd(self):
         hsvd(self.siso_ss1)
@@ -309,12 +354,12 @@
     def testRss(self):
         rss(1)
         rss(2)
-        rss(2, 3, 1)
+        rss(2, 1, 3)
 
     def testDrss(self):
         drss(1)
         drss(2)
-        drss(2, 3, 1)
+        drss(2, 1, 3)
 
     def testCtrb(self):
         ctrb(self.siso_ss1.A, self.siso_ss1.B)
@@ -372,6 +417,102 @@
         for i in range(len(tfdata_1)):
             np.testing.assert_array_almost_equal(tfdata_1[i], tfdata_2[i])
 
+    def testDamp(self):
+        A = np.mat('''-0.2  0.06 0    -1;
+               0    0    1     0;
+             -17    0   -3.8   1;
+               9.4  0   -0.4  -0.6''')
+        B = np.mat('''-0.01  0.06;
+               0     0;
+             -32     5.4;
+               2.6  -7''')
+        C = np.eye(4)
+        D = np.zeros((4,2))
+        sys = ss(A, B, C, D)
+        wn, Z, p = damp(sys, False)
+        print (wn)
+        np.testing.assert_array_almost_equal(
+            wn, np.array([4.07381994,   3.28874827,   3.28874827,
+                          1.08937685e-03]))
+        np.testing.assert_array_almost_equal(
+            Z, np.array([1.0, 0.07983139,  0.07983139, 1.0]))
+        
+    def testConnect(self):
+        sys1 = ss("1. -2; 3. -4", "5.; 7", "6, 8", "9.")
+        sys2 = ss("-1.", "1.", "1.", "0.")
+        sys = append(sys1, sys2)
+        Q= np.mat([ [ 1, 2], [2, -1] ]) # basically feedback, output 2 in 1
+        sysc = connect(sys, Q, [2], [1, 2])
+        print(sysc)
+        np.testing.assert_array_almost_equal(
+            sysc.A, np.mat('1 -2 5; 3 -4 7; -6 -8 -10'))
+        np.testing.assert_array_almost_equal(
+            sysc.B, np.mat('0; 0; 1'))
+        np.testing.assert_array_almost_equal(
+            sysc.C, np.mat('6 8 9; 0 0 1'))
+        np.testing.assert_array_almost_equal(
+            sysc.D, np.mat('0; 0'))
+
+    def testConnect2(self):
+        sys = append(ss([[-5, -2.25], [4, 0]], [[2], [0]], 
+                          [[0, 1.125]], [[0]]), 
+                       ss([[-1.6667, 0], [1, 0]], [[2], [0]],
+                          [[0, 3.3333]], [[0]]),
+                       1)
+        Q = [ [ 1, 3], [2, 1], [3, -2]]
+        sysc = connect(sys, Q, [3], [3, 1, 2])
+        np.testing.assert_array_almost_equal(
+            sysc.A, np.mat([[-5, -2.25, 0, -6.6666], 
+                            [4, 0, 0, 0],
+                            [0, 2.25, -1.6667, 0], 
+                            [0, 0, 1, 0]]))
+        np.testing.assert_array_almost_equal(
+            sysc.B, np.mat([[2], [0], [0], [0]]))
+        np.testing.assert_array_almost_equal(
+            sysc.C, np.mat([[0, 0, 0, -3.3333], 
+                            [0, 1.125, 0, 0],
+                            [0, 0, 0, 3.3333]]))
+        np.testing.assert_array_almost_equal(
+            sysc.D, np.mat([[1], [0], [0]]))
+            
+                          
+                                 
+    def testFRD(self):
+        h = tf([1], [1, 2, 2])
+        omega = np.logspace(-1, 2, 10)
+        frd1 = frd(h, omega)
+        assert isinstance(frd1, FRD)
+        frd2 = frd(frd1.fresp[0,0,:], omega)
+        assert isinstance(frd2, FRD)
+
+    def testMinreal(self, verbose=False):
+        """Test a minreal model reduction"""
+        #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
+        A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+        #B = [0.3, -1.3; 0.1, 0; 1, 0]
+        B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+        #C = [0, 0.1, 0; -0.3, -0.2, 0]
+        C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+        #D = [0 -0.8; -0.3 0]
+        D = [[0., -0.8], [-0.3, 0.]]
+        # sys = ss(A, B, C, D)
+        
+        sys = ss(A, B, C, D)
+        sysr = minreal(sys)
+        self.assertEqual(sysr.states, 2)
+        self.assertEqual(sysr.inputs, sys.inputs)
+        self.assertEqual(sysr.outputs, sys.outputs)
+        np.testing.assert_array_almost_equal(
+            eigvals(sysr.A), [-2.136154, -0.1638459])
+
+        s = tf([1, 0], [1])
+        h = (s+1)*(s+2.00000000001)/(s+2)/(s**2+s+1)
+        hm = minreal(h)
+        hr = (s+1)/(s**2+s+1)
+        np.testing.assert_array_almost_equal(hm.num[0][0], hr.num[0][0])
+        np.testing.assert_array_almost_equal(hm.den[0][0], hr.den[0][0])
+
+
 #! TODO: not yet implemented
 #    def testMIMOtfdata(self):
 #        sisotf = ss2tf(self.siso_ss1)

Modified: trunk/tests/slycot_convert_test.py
===================================================================
--- trunk/tests/slycot_convert_test.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/slycot_convert_test.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -37,7 +37,7 @@
             for inputs in range(1, self.maxI+1):
                 for outputs in range(1, self.maxO+1):
                     for testNum in range(self.numTests):
-                        ssOriginal = matlab.rss(states, inputs, outputs)
+                        ssOriginal = matlab.rss(states, outputs, inputs)
                         if (verbose):
                             print('====== Original SS ==========')
                             print(ssOriginal)
@@ -82,7 +82,7 @@
             for testNum in range(self.numTests):                       
                 for inputs in range(1,1):
                     for outputs in range(1,1):       
-                        ssOriginal = matlab.rss(states, inputs, outputs)
+                        ssOriginal = matlab.rss(states, outputs, inputs)
                         
                         tfOriginal_Actrb, tfOriginal_Bctrb, tfOriginal_Cctrb, tfOrigingal_nctrb, tfOriginal_index,\
                             tfOriginal_dcoeff, tfOriginal_ucoeff = tb04ad(states,inputs,outputs,\

Modified: trunk/tests/statesp_test.py
===================================================================
--- trunk/tests/statesp_test.py	2013-05-31 11:53:29 UTC (rev 273)
+++ trunk/tests/statesp_test.py	2013-06-09 14:57:59 UTC (rev 274)
@@ -5,8 +5,10 @@
 
 import unittest
 import numpy as np
+from scipy.linalg import eigvals
 import control.matlab as matlab
-from control.statesp import StateSpace
+from control.statesp import StateSpace, _convertToStateSpace
+from control.xferfcn import TransferFunction
 
 class TestStateSpace(unittest.TestCase):
     """Tests for the StateSpace class."""
@@ -134,7 +136,74 @@
         np.testing.assert_almost_equal(mag, truemag)
         np.testing.assert_almost_equal(phase, truephase)
         np.testing.assert_equal(omega, trueomega)
+                
+    def testMinreal(self):
+        """Test a minreal model reduction"""
+        #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
+        A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+        #B = [0.3, -1.3; 0.1, 0; 1, 0]
+        B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+        #C = [0, 0.1, 0; -0.3, -0.2, 0]
+        C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+        #D = [0 -0.8; -0.3 0]
+        D = [[0., -0.8], [-0.3, 0.]]
+        # sys = ss(A, B, C, D)
+        
+        sys = StateSpace(A, B, C, D)
+        sysr = sys.minreal()
+        self.assertEqual(sysr.states, 2)
+        self.assertEqual(sysr.inputs, sys.inputs)
+        self.assertEqual(sysr.outputs, sys.outputs)
+        np.testing.assert_array_almost_equal(
+            eigvals(sysr.A), [-2.136154, -0.1638459])
+                                    
+    def testAppendSS(self):
+        """Test appending two state-space systems"""
+        A1 = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+        B1 = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
+        C1 = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
+        D1 = [[0., -0.8], [-0.3, 0.]]
+        A2 = [[-1.]]
+        B2 = [[1.2]]
+        C2 = [[0.5]]
+        D2 = [[0.4]]
+        A3 = [[-2, 0.5, 0, 0], [0.5, -0.3, 0, 0], [0, 0, -0.1, 0], 
+              [0, 0, 0., -1.]]
+        B3 = [[0.3, -1.3, 0], [0.1, 0., 0], [1.0, 0.0, 0], [0., 0, 1.2]]
+        C3 = [[0., 0.1, 0.0, 0.0], [-0.3, -0.2, 0.0, 0.0], [0., 0., 0., 0.5]]
+        D3 = [[0., -0.8, 0.], [-0.3, 0., 0.], [0., 0., 0.4]]
+        sys1 = StateSpace(A1, B1, C1, D1)
+        sys2 = StateSpace(A2, B2, C2, D2)
+        sys3 = StateSpace(A3, B3, C3, D3)
+        sys3c = sys1.append(sys2)
+        np.testing.assert_array_almost_equal(sys3.A, sys3c.A)
+        np.testing.assert_array_almost_equal(sys3.B, sys3c.B)
+        np.testing.assert_array_almost_equal(sys3.C, sys3c.C)
+        np.testing.assert_array_almost_equal(sys3.D, sys3c.D)
+  
+    def testAppendTF(self):
+        """Test appending a state-space system with a tf"""
+        A1 = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
+        B1 = [[0.3...
 
[truncated message content]

Re: [python-control-discuss] control.feedback(sys1, sys2=1, sign=-1)

[Richard M. <mu...@cd...>]

No problem e-mailing me a patch and I'll put it in.

I should get Rene's changes integrated this weekend (just need to stop getting on and off airplanes!).

-richard

On 7 Jun 2013, at 7:42 , Ryan Krauss <rk...@si...> wrote:

> We probably need to make the change in a few places because I see that TransferFuction also has a feedback method.
> 
> My SVN skills are pretty minimal.  Do I just email you a patch?
> 
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> 
> 
> On Thu, Jun 6, 2013 at 4:15 PM, Richard Murray <mu...@cd...> wrote:
> I don't see any reason not to do this.  I agree that unity feedback is a good default.
> 
> I've got some changes pending from Rene van Paassen that I need merge into trunk.  Will try to do that in the next few days.
> 
> -richard
> 
> On 6 Jun 2013, at 12:41 , Ryan Krauss <rk...@si...> wrote:
> 
> > It would save me a lot of hassle in my conversion if the default for sys2 in control.feedback was 1 (i.e. unity feedback is assumed).
> >
> > Are there any objections to this?
> >
> > Thanks,
> >
> > Ryan
> >
> > --
> > Ryan Krauss, Ph.D.
> > Associate Professor
> > Mechanical Engineering
> > Southern Illinois University Edwardsville
> > ------------------------------------------------------------------------------
> > How ServiceNow helps IT people transform IT departments:
> > 1. A cloud service to automate IT design, transition and operations
> > 2. Dashboards that offer high-level views of enterprise services
> > 3. A single system of record for all IT processes
> > http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> > python-control-discuss mailing list
> > pyt...@li...
> > https://lists.sourceforge.net/lists/listinfo/python-control-discuss
> 
> 
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
> 
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

Re: [python-control-discuss] control.feedback(sys1, sys2=1, sign=-1)

[Ryan K. <rk...@si...>]

We probably need to make the change in a few places because I see that
TransferFuction also has a feedback method.

My SVN skills are pretty minimal.  Do I just email you a patch?


--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville


On Thu, Jun 6, 2013 at 4:15 PM, Richard Murray <mu...@cd...>wrote:

> I don't see any reason not to do this.  I agree that unity feedback is a
> good default.
>
> I've got some changes pending from Rene van Paassen that I need merge into
> trunk.  Will try to do that in the next few days.
>
> -richard
>
> On 6 Jun 2013, at 12:41 , Ryan Krauss <rk...@si...> wrote:
>
> > It would save me a lot of hassle in my conversion if the default for
> sys2 in control.feedback was 1 (i.e. unity feedback is assumed).
> >
> > Are there any objections to this?
> >
> > Thanks,
> >
> > Ryan
> >
> > --
> > Ryan Krauss, Ph.D.
> > Associate Professor
> > Mechanical Engineering
> > Southern Illinois University Edwardsville
> >
> ------------------------------------------------------------------------------
> > How ServiceNow helps IT people transform IT departments:
> > 1. A cloud service to automate IT design, transition and operations
> > 2. Dashboards that offer high-level views of enterprise services
> > 3. A single system of record for all IT processes
> >
> http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> > python-control-discuss mailing list
> > pyt...@li...
> > https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>
>
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>

Re: [python-control-discuss] call method

[Rene v. P. <ren...@gm...>]

On my experimental branch, I implemented something similar, (but maybe
yours is better), with the method horner (after Scilab's horner). But I
agree that __call__ is a perfectly fine name for this, and I think that we
should simply implement this for both StateSpace and TransferFunction.

I want to vote for keeping it mimo, though

René


On 6 June 2013 23:16, Richard Murray <mu...@cd...> wrote:

> Good call on __call__.  This seems like a fine thing to integrate in.
>
> -richard
>
> On 6 Jun 2013, at 12:36 , Ryan Krauss <rk...@si...> wrote:
>
> > poly1d makes it naturally vectorize so that the __call__ method could
> also be used in Bode-type calculations:
> >
> > In [16]: w = arange(1,5,0.5)
> >
> > In [17]: control.freqresp(test,w)
> > Out[17]:
> > (array([[[ 0.31622777,  0.2981424 ,  0.2773501 ,  0.25607376,
>  0.23570226,
> >           0.21693046,  0.2       ,  0.18490007]]]),
> >  array([[[-0.32175055, -0.46364761, -0.5880026 , -0.69473828,
> -0.78539816,
> >          -0.86217005, -0.92729522, -0.98279372]]]),
> >  array([ 1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5]))
> >
> > In [18]: s_vect = 1.0j*w
> >
> > In [18.5]: test(s_vect)
> > Out[18.5]:
> > array([ 0.30000000-0.1j       ,  0.26666667-0.13333333j,
> >         0.23076923-0.15384615j,  0.19672131-0.16393443j,
> >         0.16666667-0.16666667j,  0.14117647-0.16470588j,
> >         0.12000000-0.16j      ,  0.10256410-0.15384615j])
> >
> >
> > In [19]: freq_test = test(s_vect)
> >
> > In [20]: abs(freq_test)
> > Out[20]:
> > array([ 0.31622777,  0.2981424 ,  0.2773501 ,  0.25607376,  0.23570226,
> >         0.21693046,  0.2       ,  0.18490007])
> >
> > In [21]: angle(freq_test)
> > Out[21]:
> > array([-0.32175055, -0.46364761, -0.5880026 , -0.69473828, -0.78539816,
> >        -0.86217005, -0.92729522, -0.98279372])
> >
> >
> >
> > --
> > Ryan Krauss, Ph.D.
> > Associate Professor
> > Mechanical Engineering
> > Southern Illinois University Edwardsville
> >
> >
> > On Thu, Jun 6, 2013 at 2:31 PM, Ryan Krauss <rk...@si...> wrote:
> > FYI, I made this change in my local checkout code and ran a simple test
> and it is working:
> >
> > In [8]: test = control.TransferFunction(1,[1,3])
> >
> > In [9]: s1 = -1+2.0j
> >
> > In [10]: test(s1)
> > Out[10]: (0.25-0.25j)
> >
> > In [11]: 1.0/(s1+3)
> > Out[11]: (0.25-0.25j)
> >
> > Nothing too fancy, but solves my issue.
> >
> >
> > --
> > Ryan Krauss, Ph.D.
> > Associate Professor
> > Mechanical Engineering
> > Southern Illinois University Edwardsville
> >
> >
> > On Thu, Jun 6, 2013 at 1:47 PM, Ryan Krauss <rk...@si...> wrote:
> > I would like to submit some minor patches over the course of the summer
> to help me completely get rid of my old controls module and switch over to
> python-control for everything.  My most immediate need is for a __call__
> method.  This method might not make any sense for a state-space system, but
> for certain aspects of my modeling work, I need to evaluate transfer
> functions at certain numerical values of s.
> >
> > I would propose something like this:
> >
> > def __call__(self, s):
> >     """Evaluate a transfer function at s."""
> >     N_poly = poly1d(squeeze(self.num))
> >     D_poly = poly1d(squeeze(self.den))
> >     return N_poly(s)/D_poly(s)
> >
> >
> > It could also include a check that self is SISO and possibly throw an
> exception if it is not (pretty sure that poly1d will throw one if we pass
> in a vector of coefficients that isn't 1d).
> >
> > Would this be welcome?  Is there anything I should tweak?
> >
> > Thanks,
> >
> > Ryan
> >
> > --
> > Ryan Krauss, Ph.D.
> > Associate Professor
> > Mechanical Engineering
> > Southern Illinois University Edwardsville
> >
> >
> >
> ------------------------------------------------------------------------------
> > How ServiceNow helps IT people transform IT departments:
> > 1. A cloud service to automate IT design, transition and operations
> > 2. Dashboards that offer high-level views of enterprise services
> > 3. A single system of record for all IT processes
> >
> http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> > python-control-discuss mailing list
> > pyt...@li...
> > https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>
>
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>



-- 
René van Paassen        | ______o____/_|     Ren...@gm...
                       <[___\_\_-----<               t: +31 15 2628685
                        |  o'                   mobile: +31 6 39846891

Re: [python-control-discuss] call method

[Richard M. <mu...@cd...>]

Good call on __call__.  This seems like a fine thing to integrate in.

-richard

On 6 Jun 2013, at 12:36 , Ryan Krauss <rk...@si...> wrote:

> poly1d makes it naturally vectorize so that the __call__ method could also be used in Bode-type calculations:
> 
> In [16]: w = arange(1,5,0.5)
> 
> In [17]: control.freqresp(test,w)
> Out[17]: 
> (array([[[ 0.31622777,  0.2981424 ,  0.2773501 ,  0.25607376,  0.23570226,
>           0.21693046,  0.2       ,  0.18490007]]]),
>  array([[[-0.32175055, -0.46364761, -0.5880026 , -0.69473828, -0.78539816,
>          -0.86217005, -0.92729522, -0.98279372]]]),
>  array([ 1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5]))
> 
> In [18]: s_vect = 1.0j*w
> 
> In [18.5]: test(s_vect)
> Out[18.5]: 
> array([ 0.30000000-0.1j       ,  0.26666667-0.13333333j,
>         0.23076923-0.15384615j,  0.19672131-0.16393443j,
>         0.16666667-0.16666667j,  0.14117647-0.16470588j,
>         0.12000000-0.16j      ,  0.10256410-0.15384615j])
> 
> 
> In [19]: freq_test = test(s_vect)
> 
> In [20]: abs(freq_test)
> Out[20]: 
> array([ 0.31622777,  0.2981424 ,  0.2773501 ,  0.25607376,  0.23570226,
>         0.21693046,  0.2       ,  0.18490007])
> 
> In [21]: angle(freq_test)
> Out[21]: 
> array([-0.32175055, -0.46364761, -0.5880026 , -0.69473828, -0.78539816,
>        -0.86217005, -0.92729522, -0.98279372])
> 
> 
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> 
> 
> On Thu, Jun 6, 2013 at 2:31 PM, Ryan Krauss <rk...@si...> wrote:
> FYI, I made this change in my local checkout code and ran a simple test and it is working:
> 
> In [8]: test = control.TransferFunction(1,[1,3])
> 
> In [9]: s1 = -1+2.0j
> 
> In [10]: test(s1)
> Out[10]: (0.25-0.25j)
> 
> In [11]: 1.0/(s1+3)
> Out[11]: (0.25-0.25j)
> 
> Nothing too fancy, but solves my issue.
> 
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> 
> 
> On Thu, Jun 6, 2013 at 1:47 PM, Ryan Krauss <rk...@si...> wrote:
> I would like to submit some minor patches over the course of the summer to help me completely get rid of my old controls module and switch over to python-control for everything.  My most immediate need is for a __call__ method.  This method might not make any sense for a state-space system, but for certain aspects of my modeling work, I need to evaluate transfer functions at certain numerical values of s.  
> 
> I would propose something like this:
> 
> def __call__(self, s):
>     """Evaluate a transfer function at s."""
>     N_poly = poly1d(squeeze(self.num))
>     D_poly = poly1d(squeeze(self.den))
>     return N_poly(s)/D_poly(s)
> 
> 
> It could also include a check that self is SISO and possibly throw an exception if it is not (pretty sure that poly1d will throw one if we pass in a vector of coefficients that isn't 1d).
> 
> Would this be welcome?  Is there anything I should tweak?
> 
> Thanks,
> 
> Ryan
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> 
> 
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

Re: [python-control-discuss] control.feedback(sys1, sys2=1, sign=-1)

[Richard M. <mu...@cd...>]

I don't see any reason not to do this.  I agree that unity feedback is a good default.

I've got some changes pending from Rene van Paassen that I need merge into trunk.  Will try to do that in the next few days.  

-richard

On 6 Jun 2013, at 12:41 , Ryan Krauss <rk...@si...> wrote:

> It would save me a lot of hassle in my conversion if the default for sys2 in control.feedback was 1 (i.e. unity feedback is assumed).
> 
> Are there any objections to this?
> 
> Thanks,
> 
> Ryan
> 
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
> ------------------------------------------------------------------------------
> How ServiceNow helps IT people transform IT departments:
> 1. A cloud service to automate IT design, transition and operations
> 2. Dashboards that offer high-level views of enterprise services
> 3. A single system of record for all IT processes
> http://p.sf.net/sfu/servicenow-d2d-j_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

[python-control-discuss] control.feedback(sys1, sys2=1, sign=-1)

[Ryan K. <rk...@si...>]

It would save me a lot of hassle in my conversion if the default for sys2
in control.feedback was 1 (i.e. unity feedback is assumed).

Are there any objections to this?

Thanks,

Ryan

--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville

Re: [python-control-discuss] call method

[Ryan K. <rk...@si...>]

poly1d makes it naturally vectorize so that the __call__ method could also
be used in Bode-type calculations:

In [16]: w = arange(1,5,0.5)

In [17]: control.freqresp(test,w)
Out[17]:
(array([[[ 0.31622777,  0.2981424 ,  0.2773501 ,  0.25607376,  0.23570226,
          0.21693046,  0.2       ,  0.18490007]]]),
 array([[[-0.32175055, -0.46364761, -0.5880026 , -0.69473828, -0.78539816,
         -0.86217005, -0.92729522, -0.98279372]]]),
 array([ 1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5]))

In [18]: s_vect = 1.0j*w

In [18.5]: test(s_vect)
Out[18.5]:
array([ 0.30000000-0.1j       ,  0.26666667-0.13333333j,
        0.23076923-0.15384615j,  0.19672131-0.16393443j,
        0.16666667-0.16666667j,  0.14117647-0.16470588j,
        0.12000000-0.16j      ,  0.10256410-0.15384615j])


In [19]: freq_test = test(s_vect)

In [20]: abs(freq_test)
Out[20]:
array([ 0.31622777,  0.2981424 ,  0.2773501 ,  0.25607376,  0.23570226,
        0.21693046,  0.2       ,  0.18490007])

In [21]: angle(freq_test)
Out[21]:
array([-0.32175055, -0.46364761, -0.5880026 , -0.69473828, -0.78539816,
       -0.86217005, -0.92729522, -0.98279372])



--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville


On Thu, Jun 6, 2013 at 2:31 PM, Ryan Krauss <rk...@si...> wrote:

> FYI, I made this change in my local checkout code and ran a simple test
> and it is working:
>
> In [8]: test = control.TransferFunction(1,[1,3])
>
> In [9]: s1 = -1+2.0j
>
> In [10]: test(s1)
> Out[10]: (0.25-0.25j)
>
> In [11]: 1.0/(s1+3)
> Out[11]: (0.25-0.25j)
>
> Nothing too fancy, but solves my issue.
>
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
>
>
> On Thu, Jun 6, 2013 at 1:47 PM, Ryan Krauss <rk...@si...> wrote:
>
>> I would like to submit some minor patches over the course of the summer
>> to help me completely get rid of my old controls module and switch over to
>> python-control for everything.  My most immediate need is for a __call__
>> method.  This method might not make any sense for a state-space system, but
>> for certain aspects of my modeling work, I need to evaluate transfer
>> functions at certain numerical values of s.
>>
>> I would propose something like this:
>>
>> def __call__(self, s):
>>     """Evaluate a transfer function at s."""
>>     N_poly = poly1d(squeeze(self.num))
>>     D_poly = poly1d(squeeze(self.den))
>>     return N_poly(s)/D_poly(s)
>>
>>
>> It could also include a check that self is SISO and possibly throw an
>> exception if it is not (pretty sure that poly1d will throw one if we pass
>> in a vector of coefficients that isn't 1d).
>>
>> Would this be welcome?  Is there anything I should tweak?
>>
>> Thanks,
>>
>> Ryan
>>
>> --
>> Ryan Krauss, Ph.D.
>> Associate Professor
>> Mechanical Engineering
>> Southern Illinois University Edwardsville
>>
>
>

Re: [python-control-discuss] call method

[Ryan K. <rk...@si...>]

FYI, I made this change in my local checkout code and ran a simple test and
it is working:

In [8]: test = control.TransferFunction(1,[1,3])

In [9]: s1 = -1+2.0j

In [10]: test(s1)
Out[10]: (0.25-0.25j)

In [11]: 1.0/(s1+3)
Out[11]: (0.25-0.25j)

Nothing too fancy, but solves my issue.


--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville


On Thu, Jun 6, 2013 at 1:47 PM, Ryan Krauss <rk...@si...> wrote:

> I would like to submit some minor patches over the course of the summer to
> help me completely get rid of my old controls module and switch over to
> python-control for everything.  My most immediate need is for a __call__
> method.  This method might not make any sense for a state-space system, but
> for certain aspects of my modeling work, I need to evaluate transfer
> functions at certain numerical values of s.
>
> I would propose something like this:
>
> def __call__(self, s):
>     """Evaluate a transfer function at s."""
>     N_poly = poly1d(squeeze(self.num))
>     D_poly = poly1d(squeeze(self.den))
>     return N_poly(s)/D_poly(s)
>
>
> It could also include a check that self is SISO and possibly throw an
> exception if it is not (pretty sure that poly1d will throw one if we pass
> in a vector of coefficients that isn't 1d).
>
> Would this be welcome?  Is there anything I should tweak?
>
> Thanks,
>
> Ryan
>
> --
> Ryan Krauss, Ph.D.
> Associate Professor
> Mechanical Engineering
> Southern Illinois University Edwardsville
>

[python-control-discuss] call method

[Ryan K. <rk...@si...>]

I would like to submit some minor patches over the course of the summer to
help me completely get rid of my old controls module and switch over to
python-control for everything.  My most immediate need is for a __call__
method.  This method might not make any sense for a state-space system, but
for certain aspects of my modeling work, I need to evaluate transfer
functions at certain numerical values of s.

I would propose something like this:

def __call__(self, s):
    """Evaluate a transfer function at s."""
    N_poly = poly1d(squeeze(self.num))
    D_poly = poly1d(squeeze(self.den))
    return N_poly(s)/D_poly(s)


It could also include a check that self is SISO and possibly throw an
exception if it is not (pretty sure that poly1d will throw one if we pass
in a vector of coefficients that isn't 1d).

Would this be welcome?  Is there anything I should tweak?

Thanks,

Ryan

--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville

Re: [python-control-discuss] ODE integrator problems

[Ryan K. <rk...@si...>]

I have run into this same problem with pulse inputs.  The issue is that the
underlying Fortran numeric integration algorithm tries to determine the
largest possible step-size that can be used and still get accurate results.
 If the input starts and ends at zero but is non-zero in the middle, it is
possible for the adaptive step-size algorithm to essentially step over the
input.  There is no way for the algorithm to correctly choose the correct
step-size in that situation without exhaustively trying every possible
point on the input vector.  So, setting hmax is a perfectly valid solution
in this case.  If you input also starts at zero, you could try truncating
the input vector to start at the first non-zero element.


--
Ryan Krauss, Ph.D.
Associate Professor
Mechanical Engineering
Southern Illinois University Edwardsville


On Tue, Mar 26, 2013 at 10:55 AM, Sawyer Fuller <mi...@se...>wrote:

> Is anybody experiencing terrible performance of the integrator, used in
> functions like "step_response"? I see glitches and other random behavior
> when running my own code, and even for example the pvtol-nested.py example
> code. The integrator it uses, scipy.integrate.odeint, is derived from some
> fortran code somewhere and  appears to be doing optimization that is giving
> unwanted behavior.
>
> One can force a short simulation interval dt for a fixed-step integration
> by passing "hmax=dt" to "step-respose" and get a response that looks like
> Matlab's, but I'd rather just have a working integrator that optimizes
> correctly.
>
> Anybody know a good workaround, eg settings for the integrator that make
> it behave normally? Or why it is giving odd behavior?
>
> S
>
>
> ------------------------
> Sawyer B. Fuller, Ph.D.
> Harvard School of Engineering and Applied Science
> Wyss Institute for Biologically Inspired Engineering
> http://people.seas.harvard.edu/~minster/
>
>
>
> ------------------------------------------------------------------------------
> Own the Future-Intel&reg; Level Up Game Demo Contest 2013
> Rise to greatness in Intel's independent game demo contest.
> Compete for recognition, cash, and the chance to get your game
> on Steam. $5K grand prize plus 10 genre and skill prizes.
> Submit your demo by 6/6/13. http://p.sf.net/sfu/intel_levelupd2d
>
> _______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss
>
>

[python-control-discuss] ODE integrator problems

[Sawyer F. <mi...@se...>]

Is anybody experiencing terrible performance of the integrator, used in
functions like "step_response"? I see glitches and other random behavior
when running my own code, and even for example the pvtol-nested.py example
code. The integrator it uses, scipy.integrate.odeint, is derived from some
fortran code somewhere and  appears to be doing optimization that is giving
unwanted behavior.

One can force a short simulation interval dt for a fixed-step integration
by passing "hmax=dt" to "step-respose" and get a response that looks like
Matlab's, but I'd rather just have a working integrator that optimizes
correctly.

Anybody know a good workaround, eg settings for the integrator that make it
behave normally? Or why it is giving odd behavior?

S


------------------------
Sawyer B. Fuller, Ph.D.
Harvard School of Engineering and Applied Science
Wyss Institute for Biologically Inspired Engineering
http://people.seas.harvard.edu/~minster/

Re: [python-control-discuss] exceptions import error

[Richard M. <mu...@cd...>]

Fixed in trunk.  Thanks.

-richard

On 25 Feb 2013, at 11:33 , Dale Lukas Peterson <haz...@gm...> wrote:

> Line 83 of src/__init__.py doesn't follow the same convention for
> importing and causes an import error for me.  Attached is a patch that
> fixes the problem.
> 
> Dale Peterson
> 
> -- 
> “People call me a perfectionist, but I'm not. I'm a rightist. I do
> something until it's right, and then I move on to the next thing.”
> ― James Cameron
> <0001-Fixed-exception-import.patch>------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_d2d_feb_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

[python-control-discuss] SF.net python-control commit:[248] trunk/src/__init__.py

[<mur...@us...>]

Revision: 248
          http://sourceforge.net/p/python-control/code/248
Author:   murrayrm
Date:     2013-03-03 01:12:34 +0000 (Sun, 03 Mar 2013)
Log Message:
-----------
import error fix from Dale Peterson

Modified Paths:
--------------
    trunk/src/__init__.py

Modified: trunk/src/__init__.py
===================================================================
--- trunk/src/__init__.py	2012-11-19 05:28:05 UTC (rev 247)
+++ trunk/src/__init__.py	2013-03-03 01:12:34 UTC (rev 248)
@@ -80,7 +80,7 @@
 from control.canonical import canonical_form, reachable_form
 
 # Exceptions
-from exception import *
+from control.exception import *
 
 # Import some of the more common (and benign) MATLAB shortcuts
 # By default, don't import conflicting commands here

[python-control-discuss] exceptions import error

[Dale L. P. <haz...@gm...>]

Line 83 of src/__init__.py doesn't follow the same convention for
importing and causes an import error for me.  Attached is a patch that
fixes the problem.

Dale Peterson

-- 
“People call me a perfectionist, but I'm not. I'm a rightist. I do
something until it's right, and then I move on to the next thing.”
― James Cameron

Re: [python-control-discuss] Controlling an inverted n-link pendulum

[Richard M. <mu...@cd...>]

Thanks for sending this, Jason.

-richard

On 27 Jan 2013, at 17:47 , Jason Moore <moo...@gm...> wrote:

> You all might like this demonstration of a classic control problem all done it Python (including python-control):
> 
> http://www.moorepants.info/blog/npendulum.html
> 
> Jason
> moorepants.info
> +01 530-601-9791
> ------------------------------------------------------------------------------
> Master Visual Studio, SharePoint, SQL, ASP.NET, C# 2012, HTML5, CSS,
> MVC, Windows 8 Apps, JavaScript and much more. Keep your skills current
> with LearnDevNow - 3,200 step-by-step video tutorials by Microsoft
> MVPs and experts. ON SALE this month only -- learn more at:
> http://p.sf.net/sfu/learnnow-d2d_______________________________________________
> python-control-discuss mailing list
> pyt...@li...
> https://lists.sourceforge.net/lists/listinfo/python-control-discuss

[python-control-discuss] Controlling an inverted n-link pendulum

[Jason M. <moo...@gm...>]

You all might like this demonstration of a classic control problem all done
it Python (including python-control):

http://www.moorepants.info/blog/npendulum.html

Jason
moorepants.info
+01 530-601-9791