From 2915bafa80c2e4e49259e1b7c61b26bb0150b3fe Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Sat, 6 Jan 2018 08:28:00 -0800 Subject: [PATCH] Update evalfr usage for consistency with MATLAB Renamed evalfr() in the FRD, StateSpace and TransferFunction classes to _evalfr() and put a deprecation warning for use of the evalfr() method. Changed calls to evalfr() in frdata.py and margins.py to use _evalfr() Added unit tests for the the new methods + warnings These changes eliminate the inconsistency in the argument form between the (MATLAB compatible) evalfr() function, which takes a complex argument, and the evalfr() methods for LTI objects, which took a real number as argument (representing the frequency). The _evalfr() method retains the original functionality, since it is used for computing stability margins and for converting objects to FRD form. --- control/frdata.py | 33 ++++++++++++++++++++++++++++----- control/margins.py | 32 ++++++++------------------------ control/statesp.py | 28 ++++++++++++++++------------ control/tests/statesp_test.py | 13 ++++++++++++- control/tests/xferfcn_test.py | 22 ++++++++++++++++++---- control/xferfcn.py | 16 +++++++++++----- 6 files changed, 93 insertions(+), 51 deletions(-) diff --git a/control/frdata.py b/control/frdata.py index 32ba8f11b..8c4e1f2f6 100644 --- 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,11 +329,30 @@ 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 @@ -341,7 +360,11 @@ def evalfr(self, omega): 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 index 7f4f06c23..392acbc9a 100644 --- 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 index 4abc4175f..0510f87b0 100644 --- 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 index 245b396c6..0d3089c4e 100644 --- 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 index cb33a6190..7225e5323 100644 --- 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,15 +320,14 @@ 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 @@ -335,6 +335,20 @@ def testEvalFrSISO(self): 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 index 6f0b724a1..edaf19191 100644 --- 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])