Merge branch 'master' of https://github.com/python-control/python-control

Author: murrayrm

control/bdalg.py

@@ -61,13 +61,13 @@
 
 __all__ = ['series', 'parallel', 'negate', 'feedback', 'append', 'connect']
 
-def series(sys1, sys2):
-    """Return the series connection sys2 * sys1 for --> sys1 --> sys2 -->.
+def series(sys1, *sysn):
+    """Return the series connection (... * sys3 *) sys2 * sys1
 
     Parameters
     ----------
     sys1: scalar, StateSpace, TransferFunction, or FRD
-    sys2: scalar, StateSpace, TransferFunction, or FRD
+    *sysn: other scalers, StateSpaces, TransferFunctions, or FRDs
 
     Returns
     -------
@@ -97,20 +97,22 @@ def series(sys1, sys2):
 
     Examples
     --------
-    >>> sys3 = series(sys1, sys2) # Same as sys3 = sys2 * sys1.
+    >>> sys3 = series(sys1, sys2) # Same as sys3 = sys2 * sys1
 
-    """
+    >>> sys5 = series(sys1, sys2, sys3, sys4) # More systems
 
-    return sys2 * sys1
+    """
+    from functools import reduce
+    return reduce(lambda x, y:x*y, sysn, sys1)
 
-def parallel(sys1, sys2):
+def parallel(sys1, *sysn):
     """
-    Return the parallel connection sys1 + sys2.
+    Return the parallel connection sys1 + sys2 (+ sys3 + ...)
 
     Parameters
     ----------
     sys1: scalar, StateSpace, TransferFunction, or FRD
-    sys2: scalar, StateSpace, TransferFunction, or FRD
+    *sysn: other scalers, StateSpaces, TransferFunctions, or FRDs
 
     Returns
     -------
@@ -140,11 +142,13 @@ def parallel(sys1, sys2):
 
     Examples
     --------
-    >>> sys3 = parallel(sys1, sys2) # Same as sys3 = sys1 + sys2.
+    >>> sys3 = parallel(sys1, sys2) # Same as sys3 = sys1 + sys2
 
-    """
+    >>> sys5 = parallel(sys1, sys2, sys3, sys4) # More systems
 
-    return sys1 + sys2
+    """
+    from functools import reduce
+    return reduce(lambda x, y:x+y, sysn, sys1)
 
 def negate(sys):
     """
@@ -247,7 +251,7 @@ def feedback(sys1, sys2=1, sign=-1):
     return sys1.feedback(sys2, sign)
 
 def append(*sys):
-    '''append(sys1, sys2, ... sysn)
+    '''append(sys1, sys2, ..., sysn)
 
     Group models by appending their inputs and outputs
 

control/frdata.py

@@ -116,7 +116,7 @@ def __init__(self, *args, **kwargs):
                     (otherlti.outputs, otherlti.inputs, numfreq),
                     dtype=complex)
                 for k, w in enumerate(self.omega):
-                    self.fresp[:, :, k] = otherlti.evalfr(w)
+                    self.fresp[:, :, k] = otherlti._evalfr(w)
 
             else:
                 # The user provided a response and a freq vector
@@ -329,19 +329,42 @@ def __pow__(self,other):
             return (FRD(ones(self.fresp.shape), self.omega) / self) * \
                 (self**(other+1))
 
-
     def evalfr(self, omega):
         """Evaluate a transfer function at a single angular frequency.
 
-        self.evalfr(omega) returns the value of the frequency response
+        self._evalfr(omega) returns the value of the frequency response
+        at frequency omega.
+
+        Note that a "normal" FRD only returns values for which there is an
+        entry in the omega vector. An interpolating FRD can return
+        intermediate values.
+
+        """
+        warn("FRD.evalfr(omega) will be deprecated in a future release of python-control; use sys.eval(omega) instead", 
+             PendingDeprecationWarning)
+        return self._evalfr(omega)
+
+    # Define the `eval` function to evaluate an FRD at a given (real)
+    # frequency.  Note that we choose to use `eval` instead of `evalfr` to
+    # avoid confusion with :func:`evalfr`, which takes a complex number as its
+    # argument.  Similarly, we don't use `__call__` to avoid confusion between
+    # G(s) for a transfer function and G(omega) for an FRD object.
+    def eval(self, omega):
+        """Evaluate a transfer function at a single angular frequency.
+
+        self._evalfr(omega) returns the value of the frequency response
         at frequency omega.
 
         Note that a "normal" FRD only returns values for which there is an
         entry in the omega vector. An interpolating FRD can return
         intermediate values.
 
         """
+        return self._evalfr(omega)
 
+    # Internal function to evaluate the frequency responses
+    def _evalfr(self, omega):
+        """Evaluate a transfer function at a single angular frequency."""
         # Preallocate the output.
         if getattr(omega, '__iter__', False):
             out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
@@ -390,7 +413,7 @@ def freqresp(self, omega):
         omega.sort()
 
         for k, w in enumerate(omega):
-            fresp = self.evalfr(w)
+            fresp = self._evalfr(w)
             mag[:, :, k] = abs(fresp)
             phase[:, :, k] = angle(fresp)
 
@@ -450,7 +473,7 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
         omega.sort()
         fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
         for k, w in enumerate(omega):
-            fresp[:, :, k] = sys.evalfr(w)
+            fresp[:, :, k] = sys._evalfr(w)
 
         return FRD(fresp, omega, smooth=True)
 

control/margins.py

@@ -164,12 +164,10 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
         # 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)
-        #print ('1:w_180', w_180)
 
         # first remove imaginary and negative frequencies, epsw removes the
         # "0" frequency for type-2 systems
         w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 >= epsw)])
-        #print ('2:w_180', w_180)
 
         # evaluate response at remaining frequencies, to test for phase 180 vs 0
         with np.errstate(all='ignore'):
@@ -182,7 +180,6 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
 
         # and sort
         w_180.sort()
-        #print ('3:w_180', w_180)
 
         # test magnitude is 1 for gain crossover/phase margins
         test_wc = np.polysub(np.polyadd(_polysqr(rnum), _polysqr(inum)),
@@ -203,32 +200,28 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
 
         # find the solutions, for positive omega, and only real ones
         wstab = np.roots(test_wstab)
-        #print('wstabr', wstab)
         wstab = np.real(wstab[(np.imag(wstab) == 0) *
                         (np.real(wstab) >= 0)])
-        #print('wstab', wstab)
 
         # and find the value of the 2nd derivative there, needs to be positive
         wstabplus = np.polyval(np.polyder(test_wstab), wstab)
-        #print('wstabplus', wstabplus)
         wstab = np.real(wstab[(np.imag(wstab) == 0) * (wstab > epsw) *
                               (wstabplus > 0.)])
-        #print('wstab', wstab)
         wstab.sort()
 
     else:
         # 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]) - 1
+            return np.abs(sys._evalfr(w)[0][0]) - 1
 
         def arg(w):
             """function to calculate the phase angle at -180 deg"""
-            return np.angle(-sys.evalfr(w)[0][0])
+            return np.angle(-sys._evalfr(w)[0][0])
 
         def dstab(w):
             """function to calculate the distance from -1 point"""
-            return np.abs(sys.evalfr(w)[0][0] + 1.)
+            return np.abs(sys._evalfr(w)[0][0] + 1.)
 
         # Find all crossings, note that this depends on omega having
         # a correct range
@@ -239,37 +232,28 @@ def dstab(w):
 
         # find the phase crossings ang(H(jw) == -180
         widx = np.where(np.diff(np.sign(arg(sys.omega))))[0]
-        #print('widx (180)', widx, sys.omega[widx])
-        #print('x', sys.evalfr(sys.omega[widx])[0][0])
-        widx = widx[np.real(sys.evalfr(sys.omega[widx])[0][0]) <= 0]
-        #print('widx (180,2)', widx)
+        widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0]
         w_180 = np.array(
             [ sp.optimize.brentq(arg, sys.omega[i], sys.omega[i+1])
               for i in widx if i+1 < len(sys.omega) ])
-        #print('x', sys.evalfr(w_180)[0][0])
-        #print('w_180', w_180)
 
         # find all stab margins?
         widx = np.where(np.diff(np.sign(np.diff(dstab(sys.omega)))))[0]
-        #print('widx', widx)
-        #print('wstabx', sys.omega[widx])
         wstab = np.array([ sp.optimize.minimize_scalar(
                   dstab, bracket=(sys.omega[i], sys.omega[i+1])).x
               for i in widx if i+1 < len(sys.omega) and
               np.diff(np.diff(dstab(sys.omega[i-1:i+2])))[0] > 0 ])
-        #print('wstabf0', wstab)
         wstab = wstab[(wstab >= sys.omega[0]) *
                       (wstab <= sys.omega[-1])]
-        #print ('wstabf', wstab)
 
 
     # margins, as iterables, converted frdata and xferfcn calculations to
     # vector for this
     with np.errstate(all='ignore'):
-        gain_w_180 = np.abs(sys.evalfr(w_180)[0][0])
+        gain_w_180 = np.abs(sys._evalfr(w_180)[0][0])
         GM = 1.0/gain_w_180
-    SM = np.abs(sys.evalfr(wstab)[0][0]+1)
-    PM = np.remainder(np.angle(sys.evalfr(wc)[0][0], deg=True), 360.0) - 180.0
+    SM = np.abs(sys._evalfr(wstab)[0][0]+1)
+    PM = np.remainder(np.angle(sys._evalfr(wc)[0][0], deg=True), 360.0) - 180.0
 
     if returnall:
         return GM, PM, SM, w_180, wc, wstab
@@ -331,7 +315,7 @@ def phase_crossover_frequencies(sys):
 
     # using real() to avoid rounding errors and results like 1+0j
     # it would be nice to have a vectorized version of self.evalfr here
-    gain = np.real(np.asarray([tf.evalfr(f)[0][0] for f in realposfreq]))
+    gain = np.real(np.asarray([tf._evalfr(f)[0][0] for f in realposfreq]))
 
     return realposfreq, gain