Skip to content

Commit bb65f31

Browse files
murrayrmmurrayrm
authored
Merge pull request #179 from murrayrm/update_evalfr
Update evalfr() to be consistent with MATLAB usage
2 parents 7795a13 + 2915baf commit bb65f31

File tree

6 files changed

+93
-51
lines changed

6 files changed

+93
-51
lines changed

‎control/frdata.py‎

Lines changed: 28 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ def __init__(self, *args, **kwargs):
116116
(otherlti.outputs, otherlti.inputs, numfreq),
117117
dtype=complex)
118118
for k, w in enumerate(self.omega):
119-
self.fresp[:, :, k] = otherlti.evalfr(w)
119+
self.fresp[:, :, k] = otherlti._evalfr(w)
120120

121121
else:
122122
# The user provided a response and a freq vector
@@ -329,19 +329,42 @@ def __pow__(self,other):
329329
return (FRD(ones(self.fresp.shape), self.omega) / self) * \
330330
(self**(other+1))
331331

332-
333332
def evalfr(self, omega):
334333
"""Evaluate a transfer function at a single angular frequency.
335334
336-
self.evalfr(omega) returns the value of the frequency response
335+
self._evalfr(omega) returns the value of the frequency response
336+
at frequency omega.
337+
338+
Note that a "normal" FRD only returns values for which there is an
339+
entry in the omega vector. An interpolating FRD can return
340+
intermediate values.
341+
342+
"""
343+
warn("FRD.evalfr(omega) will be deprecated in a future release of python-control; use sys.eval(omega) instead",
344+
PendingDeprecationWarning)
345+
return self._evalfr(omega)
346+
347+
# Define the `eval` function to evaluate an FRD at a given (real)
348+
# frequency. Note that we choose to use `eval` instead of `evalfr` to
349+
# avoid confusion with :func:`evalfr`, which takes a complex number as its
350+
# argument. Similarly, we don't use `__call__` to avoid confusion between
351+
# G(s) for a transfer function and G(omega) for an FRD object.
352+
def eval(self, omega):
353+
"""Evaluate a transfer function at a single angular frequency.
354+
355+
self._evalfr(omega) returns the value of the frequency response
337356
at frequency omega.
338357
339358
Note that a "normal" FRD only returns values for which there is an
340359
entry in the omega vector. An interpolating FRD can return
341360
intermediate values.
342361
343362
"""
363+
return self._evalfr(omega)
344364

365+
# Internal function to evaluate the frequency responses
366+
def _evalfr(self, omega):
367+
"""Evaluate a transfer function at a single angular frequency."""
345368
# Preallocate the output.
346369
if getattr(omega, '__iter__', False):
347370
out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
@@ -390,7 +413,7 @@ def freqresp(self, omega):
390413
omega.sort()
391414

392415
for k, w in enumerate(omega):
393-
fresp = self.evalfr(w)
416+
fresp = self._evalfr(w)
394417
mag[:, :, k] = abs(fresp)
395418
phase[:, :, k] = angle(fresp)
396419

@@ -450,7 +473,7 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
450473
omega.sort()
451474
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
452475
for k, w in enumerate(omega):
453-
fresp[:, :, k] = sys.evalfr(w)
476+
fresp[:, :, k] = sys._evalfr(w)
454477

455478
return FRD(fresp, omega, smooth=True)
456479

‎control/margins.py‎

Lines changed: 8 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -164,12 +164,10 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
164164
# test (imaginary part of tf) == 0, for phase crossover/gain margins
165165
test_w_180 = np.polyadd(np.polymul(inum, rden), np.polymul(rnum, -iden))
166166
w_180 = np.roots(test_w_180)
167-
#print ('1:w_180', w_180)
168167

169168
# first remove imaginary and negative frequencies, epsw removes the
170169
# "0" frequency for type-2 systems
171170
w_180 = np.real(w_180[(np.imag(w_180) == 0) * (w_180 >= epsw)])
172-
#print ('2:w_180', w_180)
173171

174172
# evaluate response at remaining frequencies, to test for phase 180 vs 0
175173
with np.errstate(all='ignore'):
@@ -182,7 +180,6 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
182180

183181
# and sort
184182
w_180.sort()
185-
#print ('3:w_180', w_180)
186183

187184
# test magnitude is 1 for gain crossover/phase margins
188185
test_wc = np.polysub(np.polyadd(_polysqr(rnum), _polysqr(inum)),
@@ -203,32 +200,28 @@ def stability_margins(sysdata, returnall=False, epsw=0.0):
203200

204201
# find the solutions, for positive omega, and only real ones
205202
wstab = np.roots(test_wstab)
206-
#print('wstabr', wstab)
207203
wstab = np.real(wstab[(np.imag(wstab) == 0) *
208204
(np.real(wstab) >= 0)])
209-
#print('wstab', wstab)
210205

211206
# and find the value of the 2nd derivative there, needs to be positive
212207
wstabplus = np.polyval(np.polyder(test_wstab), wstab)
213-
#print('wstabplus', wstabplus)
214208
wstab = np.real(wstab[(np.imag(wstab) == 0) * (wstab > epsw) *
215209
(wstabplus > 0.)])
216-
#print('wstab', wstab)
217210
wstab.sort()
218211

219212
else:
220213
# a bit coarse, have the interpolated frd evaluated again
221214
def mod(w):
222215
"""to give the function to calculate |G(jw)| = 1"""
223-
return np.abs(sys.evalfr(w)[0][0]) - 1
216+
return np.abs(sys._evalfr(w)[0][0]) - 1
224217

225218
def arg(w):
226219
"""function to calculate the phase angle at -180 deg"""
227-
return np.angle(-sys.evalfr(w)[0][0])
220+
return np.angle(-sys._evalfr(w)[0][0])
228221

229222
def dstab(w):
230223
"""function to calculate the distance from -1 point"""
231-
return np.abs(sys.evalfr(w)[0][0] + 1.)
224+
return np.abs(sys._evalfr(w)[0][0] + 1.)
232225

233226
# Find all crossings, note that this depends on omega having
234227
# a correct range
@@ -239,37 +232,28 @@ def dstab(w):
239232

240233
# find the phase crossings ang(H(jw) == -180
241234
widx = np.where(np.diff(np.sign(arg(sys.omega))))[0]
242-
#print('widx (180)', widx, sys.omega[widx])
243-
#print('x', sys.evalfr(sys.omega[widx])[0][0])
244-
widx = widx[np.real(sys.evalfr(sys.omega[widx])[0][0]) <= 0]
245-
#print('widx (180,2)', widx)
235+
widx = widx[np.real(sys._evalfr(sys.omega[widx])[0][0]) <= 0]
246236
w_180 = np.array(
247237
[ sp.optimize.brentq(arg, sys.omega[i], sys.omega[i+1])
248238
for i in widx if i+1 < len(sys.omega) ])
249-
#print('x', sys.evalfr(w_180)[0][0])
250-
#print('w_180', w_180)
251239

252240
# find all stab margins?
253241
widx = np.where(np.diff(np.sign(np.diff(dstab(sys.omega)))))[0]
254-
#print('widx', widx)
255-
#print('wstabx', sys.omega[widx])
256242
wstab = np.array([ sp.optimize.minimize_scalar(
257243
dstab, bracket=(sys.omega[i], sys.omega[i+1])).x
258244
for i in widx if i+1 < len(sys.omega) and
259245
np.diff(np.diff(dstab(sys.omega[i-1:i+2])))[0] > 0 ])
260-
#print('wstabf0', wstab)
261246
wstab = wstab[(wstab >= sys.omega[0]) *
262247
(wstab <= sys.omega[-1])]
263-
#print ('wstabf', wstab)
264248

265249

266250
# margins, as iterables, converted frdata and xferfcn calculations to
267251
# vector for this
268252
with np.errstate(all='ignore'):
269-
gain_w_180 = np.abs(sys.evalfr(w_180)[0][0])
253+
gain_w_180 = np.abs(sys._evalfr(w_180)[0][0])
270254
GM = 1.0/gain_w_180
271-
SM = np.abs(sys.evalfr(wstab)[0][0]+1)
272-
PM = np.remainder(np.angle(sys.evalfr(wc)[0][0], deg=True), 360.0) - 180.0
255+
SM = np.abs(sys._evalfr(wstab)[0][0]+1)
256+
PM = np.remainder(np.angle(sys._evalfr(wc)[0][0], deg=True), 360.0) - 180.0
273257

274258
if returnall:
275259
return GM, PM, SM, w_180, wc, wstab
@@ -331,7 +315,7 @@ def phase_crossover_frequencies(sys):
331315

332316
# using real() to avoid rounding errors and results like 1+0j
333317
# it would be nice to have a vectorized version of self.evalfr here
334-
gain = np.real(np.asarray([tf.evalfr(f)[0][0] for f in realposfreq]))
318+
gain = np.real(np.asarray([tf._evalfr(f)[0][0] for f in realposfreq]))
335319

336320
return realposfreq, gain
337321

‎control/statesp.py‎

Lines changed: 16 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -61,8 +61,7 @@
6161
from numpy.linalg.linalg import LinAlgError
6262
import scipy as sp
6363
from scipy.signal import lti, cont2discrete
64-
# from exceptions import Exception
65-
import warnings
64+
from warnings import warn
6665
from .lti import LTI, timebase, timebaseEqual, isdtime
6766
from .xferfcn import _convertToTransferFunction
6867
from copy import deepcopy
@@ -380,20 +379,26 @@ def __rdiv__(self, other):
380379

381380
raise NotImplementedError("StateSpace.__rdiv__ is not implemented yet.")
382381

383-
# TODO: add discrete time check
384382
def evalfr(self, omega):
385383
"""Evaluate a SS system's transfer function at a single frequency.
386384
387-
self.evalfr(omega) returns the value of the transfer function matrix with
388-
input value s = i * omega.
385+
self._evalfr(omega) returns the value of the transfer function matrix
386+
with input value s = i * omega.
389387
390388
"""
389+
warn("StateSpace.evalfr(omega) will be depracted in a future "
390+
"release of python-control; use evalfr(sys, omega*1j) instead",
391+
PendingDeprecationWarning)
392+
return self._evalfr(omega)
393+
394+
def _evalfr(self, omega):
395+
"""Evaluate a SS system's transfer function at a single frequency"""
391396
# Figure out the point to evaluate the transfer function
392397
if isdtime(self, strict=True):
393398
dt = timebase(self)
394399
s = exp(1.j * omega * dt)
395400
if (omega * dt > math.pi):
396-
warnings.warn("evalfr: frequency evaluation above Nyquist frequency")
401+
warn("_evalfr: frequency evaluation above Nyquist frequency")
397402
else:
398403
s = omega * 1.j
399404

@@ -997,9 +1002,9 @@ def _mimo2siso(sys, input, output, warn_conversion=False):
9971002
#Convert sys to SISO if necessary
9981003
if sys.inputs > 1 or sys.outputs > 1:
9991004
if warn_conversion:
1000-
warnings.warn("Converting MIMO system to SISO system. "
1001-
"Only input {i} and output {o} are used."
1002-
.format(i=input, o=output))
1005+
warn("Converting MIMO system to SISO system. "
1006+
"Only input {i} and output {o} are used."
1007+
.format(i=input, o=output))
10031008
# $X = A*X + B*U
10041009
# Y = C*X + D*U
10051010
new_B = sys.B[:, input]
@@ -1047,9 +1052,8 @@ def _mimo2simo(sys, input, warn_conversion=False):
10471052
#Convert sys to SISO if necessary
10481053
if sys.inputs > 1:
10491054
if warn_conversion:
1050-
warnings.warn("Converting MIMO system to SIMO system. "
1051-
"Only input {i} is used."
1052-
.format(i=input))
1055+
warn("Converting MIMO system to SIMO system. "
1056+
"Only input {i} is used." .format(i=input))
10531057
# $X = A*X + B*U
10541058
# Y = C*X + D*U
10551059
new_B = sys.B[:, input]

‎control/tests/statesp_test.py‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
from control import matlab
1111
from control.statesp import StateSpace, _convertToStateSpace
1212
from control.xferfcn import TransferFunction
13+
from control.lti import evalfr
1314
from control.exception import slycot_check
1415

1516
class TestStateSpace(unittest.TestCase):
@@ -113,7 +114,17 @@ def testEvalFr(self):
113114
[-0.331544857768052 + 0.0576105032822757j,
114115
0.128919037199125 - 0.143824945295405j]]
115116

116-
np.testing.assert_almost_equal(sys.evalfr(1.), resp)
117+
# Correct versions of the call
118+
np.testing.assert_almost_equal(evalfr(sys, 1j), resp)
119+
np.testing.assert_almost_equal(sys._evalfr(1.), resp)
120+
121+
# Deprecated version of the call (should generate warning)
122+
import warnings
123+
with warnings.catch_warnings(record=True) as w:
124+
warnings.simplefilter("always")
125+
sys.evalfr(1.)
126+
assert len(w) == 1
127+
assert issubclass(w[-1].category, PendingDeprecationWarning)
117128

118129
@unittest.skipIf(not slycot_check(), "slycot not installed")
119130
def testFreqResp(self):

‎control/tests/xferfcn_test.py‎

Lines changed: 18 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import numpy as np
88
from control.statesp import StateSpace, _convertToStateSpace
99
from control.xferfcn import TransferFunction, _convertToTransferFunction
10+
from control.lti import evalfr
1011
from control.exception import slycot_check
1112
# from control.lti import isdtime
1213

@@ -319,22 +320,35 @@ def testDivSISO(self):
319320
np.testing.assert_array_equal(sys4.den, sys3.num)
320321

321322
# Tests for TransferFunction.evalfr.
322-
323323
def testEvalFrSISO(self):
324324
"""Evaluate the frequency response of a SISO system at one frequency."""
325325

326326
sys = TransferFunction([1., 3., 5], [1., 6., 2., -1])
327327

328-
np.testing.assert_array_almost_equal(sys.evalfr(1.),
328+
np.testing.assert_array_almost_equal(evalfr(sys, 1j),
329329
np.array([[-0.5 - 0.5j]]))
330-
np.testing.assert_array_almost_equal(sys.evalfr(32.),
330+
np.testing.assert_array_almost_equal(evalfr(sys, 32j),
331331
np.array([[0.00281959302585077 - 0.030628473607392j]]))
332332

333333
# Test call version as well
334334
np.testing.assert_almost_equal(sys(1.j), -0.5 - 0.5j)
335335
np.testing.assert_almost_equal(sys(32.j),
336336
0.00281959302585077 - 0.030628473607392j)
337337

338+
# Test internal version (with real argument)
339+
np.testing.assert_array_almost_equal(sys._evalfr(1.),
340+
np.array([[-0.5 - 0.5j]]))
341+
np.testing.assert_array_almost_equal(sys._evalfr(32.),
342+
np.array([[0.00281959302585077 - 0.030628473607392j]]))
343+
344+
# Deprecated version of the call (should generate warning)
345+
import warnings
346+
with warnings.catch_warnings(record=True) as w:
347+
warnings.simplefilter("always")
348+
sys.evalfr(1.)
349+
assert len(w) == 1
350+
assert issubclass(w[-1].category, PendingDeprecationWarning)
351+
338352
@unittest.skipIf(not slycot_check(), "slycot not installed")
339353
def testEvalFrMIMO(self):
340354
"""Evaluate the frequency response of a MIMO system at one frequency."""
@@ -348,7 +362,7 @@ def testEvalFrMIMO(self):
348362
[-0.083333333333333, -0.188235294117647 - 0.847058823529412j,
349363
-1. - 8.j]]
350364

351-
np.testing.assert_array_almost_equal(sys.evalfr(2.), resp)
365+
np.testing.assert_array_almost_equal(sys._evalfr(2.), resp)
352366

353367
# Test call version as well
354368
np.testing.assert_array_almost_equal(sys(2.j), resp)

‎control/xferfcn.py‎

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,7 @@
5959
import scipy as sp
6060
from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete
6161
from copy import deepcopy
62+
import warnings
6263
from warnings import warn
6364
from .lti import LTI, timebaseEqual, timebase, isdtime
6465

@@ -491,19 +492,24 @@ def __pow__(self, other):
491492
def evalfr(self, omega):
492493
"""Evaluate a transfer function at a single angular frequency.
493494
494-
self.evalfr(omega) returns the value of the
495-
transfer function matrix with
496-
input value s = i * omega.
495+
self._evalfr(omega) returns the value of the transfer function
496+
matrix with input value s = i * omega.
497497
498498
"""
499+
warn("TransferFunction.evalfr(omega) will be deprecated in a "
500+
"future release of python-control; use evalfr(sys, omega*1j) "
501+
"instead", PendingDeprecationWarning)
502+
return self._evalfr(omega)
499503

504+
def _evalfr(self, omega):
505+
"""Evaluate a transfer function at a single angular frequency."""
500506
# TODO: implement for discrete time systems
501507
if isdtime(self, strict=True):
502508
# Convert the frequency to discrete time
503509
dt = timebase(self)
504510
s = exp(1.j * omega * dt)
505511
if (omega * dt > pi):
506-
warn("evalfr: frequency evaluation above Nyquist frequency")
512+
warn("_evalfr: frequency evaluation above Nyquist frequency")
507513
else:
508514
s = 1.j * omega
509515

@@ -552,7 +558,7 @@ def freqresp(self, omega):
552558
dt = timebase(self)
553559
slist = np.array([exp(1.j * w * dt) for w in omega])
554560
if (max(omega) * dt > pi):
555-
warn("evalfr: frequency evaluation above Nyquist frequency")
561+
warn("freqresp: frequency evaluation above Nyquist frequency")
556562
else:
557563
slist = np.array([1j * w for w in omega])
558564

0 commit comments

Comments
 (0)