python-control · murrayrm · Jan 13, 2018 · Jan 2, 2018 · Jan 2, 2018 · Jan 3, 2018
diff --git a/control/frdata.py b/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)
 

diff --git a/control/margins.py b/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
 

diff --git a/control/statesp.py b/control/statesp.py
@@ -61,8 +61,7 @@
 from numpy.linalg.linalg import LinAlgError
 import scipy as sp
 from scipy.signal import lti, cont2discrete
-# from exceptions import Exception
-import warnings
+from warnings import warn
 from .lti import LTI, timebase, timebaseEqual, isdtime
 from .xferfcn import _convertToTransferFunction
 from copy import deepcopy
@@ -357,20 +356,26 @@ def __rdiv__(self, other):
 
         raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.")
 
-    # TODO: add discrete time check
     def evalfr(self, omega):
         """Evaluate a SS system's transfer function at a single frequency.
 
-        self.evalfr(omega) returns the value of the transfer function matrix with
-        input value s = i * omega.
+        self._evalfr(omega) returns the value of the transfer function matrix
+        with input value s = i * omega.
 
         """
+        warn("StateSpace.evalfr(omega) will be depracted in a future "
+             "release of python-control; use evalfr(sys, omega*1j) instead",
+             PendingDeprecationWarning)
+        return self._evalfr(omega)
+
+    def _evalfr(self, omega):
+        """Evaluate a SS system's transfer function at a single frequency"""
         # Figure out the point to evaluate the transfer function
         if isdtime(self, strict=True):
             dt = timebase(self)
             s = exp(1.j * omega * dt)
             if (omega * dt > math.pi):
-                warnings.warn("evalfr: frequency evaluation above Nyquist frequency")
+                warn("_evalfr: frequency evaluation above Nyquist frequency")
         else:
             s = omega * 1.j
 
@@ -900,9 +905,9 @@ def _mimo2siso(sys, input, output, warn_conversion=False):
     #Convert sys to SISO if necessary
     if sys.inputs > 1 or sys.outputs > 1:
         if warn_conversion:
-            warnings.warn("Converting MIMO system to SISO system. "
-                          "Only input {i} and output {o} are used."
-                          .format(i=input, o=output))
+            warn("Converting MIMO system to SISO system. "
+                 "Only input {i} and output {o} are used."
+                 .format(i=input, o=output))
         # $X = A*X + B*U
         #  Y = C*X + D*U
         new_B = sys.B[:, input]
@@ -950,9 +955,8 @@ def _mimo2simo(sys, input, warn_conversion=False):
     #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))
+            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]

diff --git a/control/tests/statesp_test.py b/control/tests/statesp_test.py
@@ -10,6 +10,7 @@
 from control import matlab
 from control.statesp import StateSpace, _convertToStateSpace
 from control.xferfcn import TransferFunction
+from control.lti import evalfr
 from control.exception import slycot_check
 
 class TestStateSpace(unittest.TestCase):
@@ -113,7 +114,17 @@ def testEvalFr(self):
                 [-0.331544857768052 + 0.0576105032822757j,
                  0.128919037199125 - 0.143824945295405j]]
 
-        np.testing.assert_almost_equal(sys.evalfr(1.), resp)
+        # Correct versions of the call
+        np.testing.assert_almost_equal(evalfr(sys, 1j), resp)
+        np.testing.assert_almost_equal(sys._evalfr(1.), resp)
+
+        # Deprecated version of the call (should generate warning)
+        import warnings
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter("always")
+            sys.evalfr(1.)
+            assert len(w) == 1
+            assert issubclass(w[-1].category, PendingDeprecationWarning)
 
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testFreqResp(self):

diff --git a/control/tests/xferfcn_test.py b/control/tests/xferfcn_test.py
@@ -7,6 +7,7 @@
 import numpy as np
 from control.statesp import StateSpace, _convertToStateSpace
 from control.xferfcn import TransferFunction, _convertToTransferFunction
+from control.lti import evalfr
 from control.exception import slycot_check
 # from control.lti import isdtime
 
@@ -319,22 +320,35 @@ def testDivSISO(self):
         np.testing.assert_array_equal(sys4.den, sys3.num)
 
     # Tests for TransferFunction.evalfr.
-
     def testEvalFrSISO(self):
         """Evaluate the frequency response of a SISO system at one frequency."""
 
         sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
 
-        np.testing.assert_array_almost_equal(sys.evalfr(1.),
+        np.testing.assert_array_almost_equal(evalfr(sys, 1j),
             np.array([[-0.5 - 0.5j]]))
-        np.testing.assert_array_almost_equal(sys.evalfr(32.),
+        np.testing.assert_array_almost_equal(evalfr(sys, 32j),
             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)
 
+        # Test internal version (with real argument)
+        np.testing.assert_array_almost_equal(sys._evalfr(1.),
+            np.array([[-0.5 - 0.5j]]))
+        np.testing.assert_array_almost_equal(sys._evalfr(32.),
+            np.array([[0.00281959302585077 - 0.030628473607392j]]))
+
+        # Deprecated version of the call (should generate warning)
+        import warnings
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter("always")
+            sys.evalfr(1.)
+            assert len(w) == 1
+            assert issubclass(w[-1].category, PendingDeprecationWarning)
+
     @unittest.skipIf(not slycot_check(), "slycot not installed")
     def testEvalFrMIMO(self):
         """Evaluate the frequency response of a MIMO system at one frequency."""
@@ -348,7 +362,7 @@ def testEvalFrMIMO(self):
                 [-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
                  -1. - 8.j]]
 
-        np.testing.assert_array_almost_equal(sys.evalfr(2.), resp)
+        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)

diff --git a/control/xferfcn.py b/control/xferfcn.py
@@ -59,6 +59,7 @@
 import scipy as sp
 from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete
 from copy import deepcopy
+import warnings
 from warnings import warn
 from .lti import LTI, timebaseEqual, timebase, isdtime
 
@@ -491,19 +492,24 @@ def __pow__(self, other):
     def evalfr(self, omega):
         """Evaluate a transfer function at a single angular frequency.
 
-        self.evalfr(omega) returns the value of the
-        transfer function matrix with
-        input value s = i * omega.
+        self._evalfr(omega) returns the value of the transfer function
+        matrix with input value s = i * omega.
 
         """
+        warn("TransferFunction.evalfr(omega) will be deprecated in a "
+             "future release of python-control; use evalfr(sys, omega*1j) "
+             "instead", PendingDeprecationWarning)
+        return self._evalfr(omega)
 
+    def _evalfr(self, omega):
+        """Evaluate a transfer function at a single angular frequency."""
         # TODO: implement for discrete time systems
         if isdtime(self, strict=True):
             # Convert the frequency to discrete time
             dt = timebase(self)
             s = exp(1.j * omega * dt)
             if (omega * dt > pi):
-                warn("evalfr: frequency evaluation above Nyquist frequency")
+                warn("_evalfr: frequency evaluation above Nyquist frequency")
         else:
             s = 1.j * omega
 
@@ -552,7 +558,7 @@ def freqresp(self, omega):
             dt = timebase(self)
             slist = np.array([exp(1.j * w * dt) for w in omega])
             if (max(omega) * dt > pi):
-                warn("evalfr: frequency evaluation above Nyquist frequency")
+                warn("freqresp: frequency evaluation above Nyquist frequency")
         else:
             slist = np.array([1j * w for w in omega])