Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commitbb65f31

Browse files
authored
Merge pull request#179 from murrayrm/update_evalfr
Update evalfr() to be consistent with MATLAB usage
2 parents7795a13 +2915baf commitbb65f31

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
fork,winenumerate(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
defevalfr(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+
returnself._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+
defeval(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+
returnself._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
ifgetattr(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
fork,winenumerate(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
fork,winenumerate(omega):
453-
fresp[:, :,k]=sys.evalfr(w)
476+
fresp[:, :,k]=sys._evalfr(w)
454477

455478
returnFRD(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
withnp.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
defmod(w):
222215
"""to give the function to calculate |G(jw)| = 1"""
223-
returnnp.abs(sys.evalfr(w)[0][0])-1
216+
returnnp.abs(sys._evalfr(w)[0][0])-1
224217

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

229222
defdstab(w):
230223
"""function to calculate the distance from -1 point"""
231-
returnnp.abs(sys.evalfr(w)[0][0]+1.)
224+
returnnp.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
foriinwidxifi+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
foriinwidxifi+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
withnp.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
ifreturnall:
275259
returnGM,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]forfinrealposfreq]))
318+
gain=np.real(np.asarray([tf._evalfr(f)[0][0]forfinrealposfreq]))
335319

336320
returnrealposfreq,gain
337321

‎control/statesp.py

Lines changed: 16 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -61,8 +61,7 @@
6161
fromnumpy.linalg.linalgimportLinAlgError
6262
importscipyassp
6363
fromscipy.signalimportlti,cont2discrete
64-
# from exceptions import Exception
65-
importwarnings
64+
fromwarningsimportwarn
6665
from .ltiimportLTI,timebase,timebaseEqual,isdtime
6766
from .xferfcnimport_convertToTransferFunction
6867
fromcopyimportdeepcopy
@@ -380,20 +379,26 @@ def __rdiv__(self, other):
380379

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

383-
# TODO: add discrete time check
384382
defevalfr(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+
withinput 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+
returnself._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
ifisdtime(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
ifsys.inputs>1orsys.outputs>1:
9991004
ifwarn_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
ifsys.inputs>1:
10491054
ifwarn_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
fromcontrolimportmatlab
1111
fromcontrol.statespimportStateSpace,_convertToStateSpace
1212
fromcontrol.xferfcnimportTransferFunction
13+
fromcontrol.ltiimportevalfr
1314
fromcontrol.exceptionimportslycot_check
1415

1516
classTestStateSpace(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+
importwarnings
123+
withwarnings.catch_warnings(record=True)asw:
124+
warnings.simplefilter("always")
125+
sys.evalfr(1.)
126+
assertlen(w)==1
127+
assertissubclass(w[-1].category,PendingDeprecationWarning)
117128

118129
@unittest.skipIf(notslycot_check(),"slycot not installed")
119130
deftestFreqResp(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
importnumpyasnp
88
fromcontrol.statespimportStateSpace,_convertToStateSpace
99
fromcontrol.xferfcnimportTransferFunction,_convertToTransferFunction
10+
fromcontrol.ltiimportevalfr
1011
fromcontrol.exceptionimportslycot_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
deftestEvalFrSISO(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+
importwarnings
346+
withwarnings.catch_warnings(record=True)asw:
347+
warnings.simplefilter("always")
348+
sys.evalfr(1.)
349+
assertlen(w)==1
350+
assertissubclass(w[-1].category,PendingDeprecationWarning)
351+
338352
@unittest.skipIf(notslycot_check(),"slycot not installed")
339353
deftestEvalFrMIMO(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
importscipyassp
6060
fromscipy.signalimportlti,tf2zpk,zpk2tf,cont2discrete
6161
fromcopyimportdeepcopy
62+
importwarnings
6263
fromwarningsimportwarn
6364
from .ltiimportLTI,timebaseEqual,timebase,isdtime
6465

@@ -491,19 +492,24 @@ def __pow__(self, other):
491492
defevalfr(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+
returnself._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
ifisdtime(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)forwinomega])
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*wforwinomega])
558564

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp