Skip to content

REBASE: use __call__ instead of evalfr in lti system classes #499

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

REBASE: use __call__ instead of evalfr in lti system classes #499

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions control/bench/time_freqresp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
sys_tf = tf(sys)
w = logspace(-1,1,50)
ntimes = 1000
time_ss = timeit("sys.freqresp(w)", setup="from __main__ import sys, w", number=ntimes)
time_tf = timeit("sys_tf.freqresp(w)", setup="from __main__ import sys_tf, w", number=ntimes)
time_ss = timeit("sys.freqquency_response(w)", setup="from __main__ import sys, w", number=ntimes)
time_tf = timeit("sys_tf.frequency_response(w)", setup="from __main__ import sys_tf, w", number=ntimes)
print("State-space model on %d states: %f" % (nstates, time_ss))
print("Transfer-function model on %d states: %f" % (nstates, time_tf))
192 changes: 95 additions & 97 deletions control/frdata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
from warnings import warn
import numpy as np
from numpy import angle, array, empty, ones, \
real, imag, absolute, eye, linalg, where, dot
real, imag, absolute, eye, linalg, where, dot, sort
from scipy.interpolate import splprep, splev
from .lti import LTI

Expand Down Expand Up @@ -100,24 +100,18 @@ def __init__(self, *args, **kwargs):
object, other than an FRD, call FRD(sys, omega)

"""
# TODO: discrete-time FRD systems?
smooth = kwargs.get('smooth', False)

if len(args) == 2:
if not isinstance(args[0], FRD) and isinstance(args[0], LTI):
# not an FRD, but still a system, second argument should be
# the frequency range
otherlti = args[0]
self.omega = array(args[1], dtype=float)
self.omega.sort()
self.omega = sort(np.asarray(args[1], dtype=float))
numfreq = len(self.omega)

# calculate frequency response at my points
self.fresp = empty(
(otherlti.outputs, otherlti.inputs, numfreq),
dtype=complex)
for k, w in enumerate(self.omega):
self.fresp[:, :, k] = otherlti._evalfr(w)

self.fresp = otherlti(1j * self.omega, squeeze=False)
else:
# The user provided a response and a freq vector
self.fresp = array(args[0], dtype=complex)
Expand All @@ -141,7 +135,7 @@ def __init__(self, *args, **kwargs):
self.omega = args[0].omega
self.fresp = args[0].fresp
else:
raise ValueError("Needs 1 or 2 arguments; receivd %i." % len(args))
raise ValueError("Needs 1 or 2 arguments; received %i." % len(args))

# create interpolation functions
if smooth:
Expand All @@ -163,17 +157,17 @@ def __str__(self):
mimo = self.inputs > 1 or self.outputs > 1
outstr = ['Frequency response data']

mt, pt, wt = self.freqresp(self.omega)
for i in range(self.inputs):
for j in range(self.outputs):
if mimo:
outstr.append("Input %i to output %i:" % (i + 1, j + 1))
outstr.append('Freq [rad/s] Response')
outstr.append('------------ ---------------------')
outstr.extend(
['%12.3f %10.4g%+10.4gj' % (w, m, p)
for m, p, w in zip(real(self.fresp[j, i, :]),
imag(self.fresp[j, i, :]), wt)])
['%12.3f %10.4g%+10.4gj' % (w, re, im)
for w, re, im in zip(self.omega,
real(self.fresp[j, i, :]),
imag(self.fresp[j, i, :]))])

return '\n'.join(outstr)

Expand Down Expand Up @@ -342,110 +336,115 @@ 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
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) # pragma: no coverage
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.
# update Sawyer B. Fuller 2020.08.14: __call__ added to provide a uniform
# interface to systems in general and the lti.frequency_response method
def eval(self, omega, squeeze=True):
"""Evaluate a transfer function at angular 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.

Parameters
----------
omega: float or array_like
Frequencies in radians per second
squeeze: bool, optional (default=True)
If True and `sys` is single input single output (SISO), returns a
1D array or scalar depending on the length of `omega`

Returns
-------
fresp : (self.outputs, self.inputs, len(x)) or len(x) complex ndarray
The frequency response of the system. Array is ``len(x)`` if and
only if system is SISO and ``squeeze=True``.

"""
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)
else:
out = empty((self.outputs, self.inputs), dtype=complex)
omega_array = np.array(omega, ndmin=1) # array-like version of omega
if any(omega_array.imag > 0):
raise ValueError("FRD.eval can only accept real-valued omega")

if self.ifunc is None:
try:
out = self.fresp[:, :, where(self.omega == omega)[0][0]]
except Exception:
elements = np.isin(self.omega, omega) # binary array
if sum(elements) < len(omega_array):
raise ValueError(
"Frequency %f not in frequency list, try an interpolating"
" FRD if you want additional points" % omega)
else:
if getattr(omega, '__iter__', False):
for i in range(self.outputs):
for j in range(self.inputs):
for k, w in enumerate(omega):
frraw = splev(w, self.ifunc[i, j], der=0)
out[i, j, k] = frraw[0] + 1.0j * frraw[1]
"not all frequencies omega are in frequency list of FRD "
"system. Try an interpolating FRD for additional points.")
else:
for i in range(self.outputs):
for j in range(self.inputs):
frraw = splev(omega, self.ifunc[i, j], der=0)
out[i, j] = frraw[0] + 1.0j * frraw[1]

out = self.fresp[:, :, elements]
else:
out = empty((self.outputs, self.inputs, len(omega_array)),
dtype=complex)
for i in range(self.outputs):
for j in range(self.inputs):
for k, w in enumerate(omega_array):
frraw = splev(w, self.ifunc[i, j], der=0)
out[i, j, k] = frraw[0] + 1.0j * frraw[1]
if not hasattr(omega, '__len__'):
# omega is a scalar, squeeze down array along last dim
out = np.squeeze(out, axis=2)
if squeeze and self.issiso():
out = out[0][0]
return out

# Method for generating the frequency response of the system
def freqresp(self, omega):
"""Evaluate the frequency response at a list of angular frequencies.
def __call__(self, s, squeeze=True):
"""Evaluate system's transfer function at complex frequencies.

Reports the value of the magnitude, phase, and angular frequency of
the requency response evaluated at omega, where omega is a list of
angular frequencies, and is a sorted version of the input omega.
Returns the complex frequency response `sys(s)`.

To evaluate at a frequency omega in radians per second, enter
``x = omega * 1j`` or use ``sys.eval(omega)``

Parameters
----------
omega : array_like
A list of frequencies in radians/sec at which the system should be
evaluated. The list can be either a python list or a numpy array
and will be sorted before evaluation.
s: complex scalar or array_like
Complex frequencies
squeeze: bool, optional (default=True)
If True and `sys` is single input single output (SISO), returns a
1D array or scalar depending on the length of x`.

Returns
-------
mag : (self.outputs, self.inputs, len(omega)) ndarray
The magnitude (absolute value, not dB or log10) of the system
frequency response.
phase : (self.outputs, self.inputs, len(omega)) ndarray
The wrapped phase in radians of the system frequency response.
omega : ndarray or list or tuple
The list of sorted frequencies at which the response was
evaluated.
"""
# Preallocate outputs.
numfreq = len(omega)
mag = empty((self.outputs, self.inputs, numfreq))
phase = empty((self.outputs, self.inputs, numfreq))
fresp : (self.outputs, self.inputs, len(s)) or len(s) complex ndarray
The frequency response of the system. Array is ``len(s)`` if and
only if system is SISO and ``squeeze=True``.

omega.sort()
Raises
------
ValueError
If `s` is not purely imaginary, because
:class:`FrequencyDomainData` systems are only defined at imaginary
frequency values.

"""
if any(abs(np.array(s, ndmin=1).real) > 0):
raise ValueError("__call__: FRD systems can only accept"
"purely imaginary frequencies")
# need to preserve array or scalar status
if hasattr(s, '__len__'):
return self.eval(np.asarray(s).imag, squeeze=squeeze)
else:
return self.eval(complex(s).imag, squeeze=squeeze)

for k, w in enumerate(omega):
fresp = self._evalfr(w)
mag[:, :, k] = abs(fresp)
phase[:, :, k] = angle(fresp)
def freqresp(self, omega):
"""(deprecated) Evaluate transfer function at complex frequencies.

return mag, phase, omega
.. deprecated::0.9.0
Method has been given the more pythonic name
:meth:`FrequencyResponseData.frequency_response`. Or use
:func:`freqresp` in the MATLAB compatibility module.
"""
warn("FrequencyResponseData.freqresp(omega) will be removed in a "
"future release of python-control; use "
"FrequencyResponseData.frequency_response(omega), or "
"freqresp(sys, omega) in the MATLAB compatibility module "
"instead", DeprecationWarning)
return self.frequency_response(omega)

def feedback(self, other=1, sign=-1):
"""Feedback interconnection between two FRD objects."""
Expand Down Expand Up @@ -515,11 +514,10 @@ def _convertToFRD(sys, omega, inputs=1, outputs=1):
"Frequency ranges of FRD do not match, conversion not implemented")

elif isinstance(sys, LTI):
omega.sort()
fresp = empty((sys.outputs, sys.inputs, len(omega)), dtype=complex)
for k, w in enumerate(omega):
fresp[:, :, k] = sys._evalfr(w)

omega = np.sort(omega)
fresp = sys(1j * omega)
if len(fresp.shape) == 1:
fresp = fresp[np.newaxis, np.newaxis, :]
return FRD(fresp, omega, smooth=True)

elif isinstance(sys, (int, float, complex, np.number)):
Expand Down
24 changes: 12 additions & 12 deletions control/freqplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,14 +151,14 @@ def bode_plot(syslist, omega=None,

Notes
-----
1. Alternatively, you may use the lower-level method
``(mag, phase, freq) = sys.freqresp(freq)`` to generate the frequency
response for a system, but it returns a MIMO response.
1. Alternatively, you may use the lower-level methods
:meth:`LTI.frequency_response` or ``sys(s)`` or ``sys(z)`` or to
generate the frequency response for a single system.

2. If a discrete time model is given, the frequency response is plotted
along the upper branch of the unit circle, using the mapping z = exp(j
\\omega dt) where omega ranges from 0 to pi/dt and dt is the discrete
timebase. If not timebase is specified (dt = True), dt is set to 1.
along the upper branch of the unit circle, using the mapping ``z = exp(1j
* omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt` is the discrete
timebase. If timebase not specified (``dt=True``), `dt` is set to 1.

Examples
--------
Expand Down Expand Up @@ -228,7 +228,7 @@ def bode_plot(syslist, omega=None,
nyquistfrq = None

# Get the magnitude and phase of the system
mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
mag_tmp, phase_tmp, omega_sys = sys.frequency_response(omega_sys)
mag = np.atleast_1d(np.squeeze(mag_tmp))
phase = np.atleast_1d(np.squeeze(phase_tmp))

Expand Down Expand Up @@ -588,7 +588,7 @@ def nyquist_plot(syslist, omega=None, plot=True, label_freq=0,
"Nyquist is currently only implemented for SISO systems.")
else:
# Get the magnitude and phase of the system
mag_tmp, phase_tmp, omega = sys.freqresp(omega)
mag_tmp, phase_tmp, omega = sys.frequency_response(omega)
mag = np.squeeze(mag_tmp)
phase = np.squeeze(phase_tmp)

Expand Down Expand Up @@ -718,7 +718,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
omega_plot = omega / (2. * math.pi) if Hz else omega

# TODO: Need to add in the mag = 1 lines
mag_tmp, phase_tmp, omega = S.freqresp(omega)
mag_tmp, phase_tmp, omega = S.frequency_response(omega)
mag = np.squeeze(mag_tmp)
if dB:
plot_axes['s'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
Expand All @@ -728,7 +728,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
plot_axes['s'].tick_params(labelbottom=False)
plot_axes['s'].grid(grid, which='both')

mag_tmp, phase_tmp, omega = (P * S).freqresp(omega)
mag_tmp, phase_tmp, omega = (P * S).frequency_response(omega)
mag = np.squeeze(mag_tmp)
if dB:
plot_axes['ps'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
Expand All @@ -738,7 +738,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
plot_axes['ps'].set_ylabel("$|PS|$" + " (dB)" if dB else "")
plot_axes['ps'].grid(grid, which='both')

mag_tmp, phase_tmp, omega = (C * S).freqresp(omega)
mag_tmp, phase_tmp, omega = (C * S).frequency_response(omega)
mag = np.squeeze(mag_tmp)
if dB:
plot_axes['cs'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
Expand All @@ -749,7 +749,7 @@ def gangof4_plot(P, C, omega=None, **kwargs):
plot_axes['cs'].set_ylabel("$|CS|$" + " (dB)" if dB else "")
plot_axes['cs'].grid(grid, which='both')

mag_tmp, phase_tmp, omega = T.freqresp(omega)
mag_tmp, phase_tmp, omega = T.frequency_response(omega)
mag = np.squeeze(mag_tmp)
if dB:
plot_axes['t'].semilogx(omega_plot, 20 * np.log10(mag), **kwargs)
Expand Down
Loading