Skip to content

Commit a11d3be

Browse files
murrayrmmurrayrm
authored
Merge pull request #924 from murrayrm/freq_plots-30Jun2023
2 parents 0a6146b + 8e84d00 commit a11d3be

39 files changed

+9899
-6771
lines changed

‎LICENSE

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
Copyright (c) 2009-2016 by California Institute of Technology
2+
Copyright (c) 2016-2023 by python-control developers
23
All rights reserved.
34

45
Redistribution and use in source and binary forms, with or without

‎control/config.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -326,13 +326,13 @@ def use_legacy_defaults(version):
326326
#
327327
# Use this function to handle a legacy keyword that has been renamed. This
328328
# function pops the old keyword off of the kwargs dictionary and issues a
329-
# warning. if both the old and new keyword are present, a ControlArgument
329+
# warning. If both the old and new keyword are present, a ControlArgument
330330
# exception is raised.
331331
#
332332
def _process_legacy_keyword(kwargs, oldkey, newkey, newval):
333333
if kwargs.get(oldkey) is not None:
334334
warnings.warn(
335-
f"keyworld '{oldkey}' is deprecated; use '{newkey}'",
335+
f"keyword '{oldkey}' is deprecated; use '{newkey}'",
336336
DeprecationWarning)
337337
if newval is not None:
338338
raise ControlArgument(

‎control/ctrlutil.py

Lines changed: 2 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -86,18 +86,9 @@ def unwrap(angle, period=2*math.pi):
8686
return angle
8787

8888
def issys(obj):
89-
"""Return True if an object is a Linear Time Invariant (LTI) system,
90-
otherwise False.
89+
"""Deprecated function to check if an object is an LTI system.
9190
92-
Examples
93-
--------
94-
>>> G = ct.tf([1], [1, 1])
95-
>>> ct.issys(G)
96-
True
97-
98-
>>> K = np.array([[1, 1]])
99-
>>> ct.issys(K)
100-
False
91+
Use isinstance(obj, ct.LTI)
10192
10293
"""
10394
warnings.warn("issys() is deprecated; use isinstance(obj, ct.LTI)",

‎control/descfcn.py

Lines changed: 195 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -18,9 +18,11 @@
1818
import scipy
1919
from warnings import warn
2020

21-
from .freqplot import nyquist_plot
21+
from .freqplot import nyquist_response
22+
from . import config
2223

2324
__all__ = ['describing_function', 'describing_function_plot',
25+
'describing_function_response', 'DescribingFunctionResponse',
2426
'DescribingFunctionNonlinearity', 'friction_backlash_nonlinearity',
2527
'relay_hysteresis_nonlinearity', 'saturation_nonlinearity']
2628

@@ -205,14 +207,74 @@ def describing_function(
205207
# Return the values in the same shape as they were requested
206208
return retdf
207209

210+
#
211+
# Describing function response/plot
212+
#
208213

209-
def describing_function_plot(
210-
H, F, A, omega=None, refine=True, label="%5.2g @ %-5.2g",
211-
warn=None, **kwargs):
212-
"""Plot a Nyquist plot with a describing function for a nonlinear system.
214+
# Simple class to store the describing function response
215+
class DescribingFunctionResponse:
216+
"""Results of describing function analysis.
217+
218+
Describing functions allow analysis of a linear I/O systems with a
219+
static nonlinear feedback function. The DescribingFunctionResponse
220+
class is used by the :func:`~control.describing_function_response`
221+
function to return the results of a describing function analysis. The
222+
response object can be used to obtain information about the describing
223+
function analysis or generate a Nyquist plot showing the frequency
224+
response of the linear systems and the describing function for the
225+
nonlinear element.
226+
227+
Attributes
228+
----------
229+
response : :class:`~control.FrequencyResponseData`
230+
Frequency response of the linear system component of the system.
231+
intersections : 1D array of 2-tuples or None
232+
A list of all amplitudes and frequencies in which
233+
:math:`H(j\\omega) N(a) = -1`, where :math:`N(a)` is the describing
234+
function associated with `F`, or `None` if there are no such
235+
points. Each pair represents a potential limit cycle for the
236+
closed loop system with amplitude given by the first value of the
237+
tuple and frequency given by the second value.
238+
N_vals : complex array
239+
Complex value of the describing function.
240+
positions : list of complex
241+
Location of the intersections in the complex plane.
242+
243+
"""
244+
def __init__(self, response, N_vals, positions, intersections):
245+
"""Create a describing function response data object."""
246+
self.response = response
247+
self.N_vals = N_vals
248+
self.positions = positions
249+
self.intersections = intersections
250+
251+
def plot(self, **kwargs):
252+
"""Plot the results of a describing function analysis.
253+
254+
See :func:`~control.describing_function_plot` for details.
255+
"""
256+
return describing_function_plot(self, **kwargs)
257+
258+
# Implement iter, getitem, len to allow recovering the intersections
259+
def __iter__(self):
260+
return iter(self.intersections)
261+
262+
def __getitem__(self, index):
263+
return list(self.__iter__())[index]
264+
265+
def __len__(self):
266+
return len(self.intersections)
213267

214-
This function generates a Nyquist plot for a closed loop system consisting
215-
of a linear system with a static nonlinear function in the feedback path.
268+
269+
# Compute the describing function response + intersections
270+
def describing_function_response(
271+
H, F, A, omega=None, refine=True, warn_nyquist=None,
272+
plot=False, check_kwargs=True, **kwargs):
273+
"""Compute the describing function response of a system.
274+
275+
This function uses describing function analysis to analyze a closed
276+
loop system consisting of a linear system with a static nonlinear
277+
function in the feedback path.
216278
217279
Parameters
218280
----------
@@ -226,53 +288,53 @@ def describing_function_plot(
226288
List of amplitudes to be used for the describing function plot.
227289
omega : list, optional
228290
List of frequencies to be used for the linear system Nyquist curve.
229-
label : str, optional
230-
Formatting string used to label intersection points on the Nyquist
231-
plot. Defaults to "%5.2g @ %-5.2g". Set to `None` to omit labels.
232-
warn : bool, optional
291+
warn_nyquist : bool, optional
233292
Set to True to turn on warnings generated by `nyquist_plot` or False
234293
to turn off warnings. If not set (or set to None), warnings are
235294
turned off if omega is specified, otherwise they are turned on.
236295
237296
Returns
238297
-------
239-
intersections : 1D array of 2-tuples or None
240-
A list of all amplitudes and frequencies in which :math:`H(j\\omega)
241-
N(a) = -1`, where :math:`N(a)` is the describing function associated
242-
with `F`, or `None` if there are no such points. Each pair represents
243-
a potential limit cycle for the closed loop system with amplitude
244-
given by the first value of the tuple and frequency given by the
245-
second value.
298+
response : :class:`~control.DescribingFunctionResponse` object
299+
Response object that contains the result of the describing function
300+
analysis. The following information can be retrieved from this
301+
object:
302+
response.intersections : 1D array of 2-tuples or None
303+
A list of all amplitudes and frequencies in which
304+
:math:`H(j\\omega) N(a) = -1`, where :math:`N(a)` is the describing
305+
function associated with `F`, or `None` if there are no such
306+
points. Each pair represents a potential limit cycle for the
307+
closed loop system with amplitude given by the first value of the
308+
tuple and frequency given by the second value.
246309
247310
Examples
248311
--------
249312
>>> H_simple = ct.tf([8], [1, 2, 2, 1])
250313
>>> F_saturation = ct.saturation_nonlinearity(1)
251314
>>> amp = np.linspace(1, 4, 10)
252-
>>> ct.describing_function_plot(H_simple, F_saturation, amp) # doctest: +SKIP
315+
>>> response = ct.describing_function_response(H_simple, F_saturation, amp)
316+
>>> response.intersections # doctest: +SKIP
253317
[(3.343844998258643, 1.4142293090899216)]
318+
>>> lines = response.plot()
254319
255320
"""
256321
# Decide whether to turn on warnings or not
257-
if warn is None:
322+
if warn_nyquist is None:
258323
# Turn warnings on unless omega was specified
259-
warn = omega is None
324+
warn_nyquist = omega is None
260325

261326
# Start by drawing a Nyquist curve
262-
count, contour = nyquist_plot(
263-
H, omega, plot=True, return_contour=True,
264-
warn_encirclements=warn, warn_nyquist=warn, **kwargs)
265-
H_omega, H_vals = contour.imag, H(contour)
327+
response = nyquist_response(
328+
H, omega, warn_encirclements=warn_nyquist, warn_nyquist=warn_nyquist,
329+
check_kwargs=check_kwargs, **kwargs)
330+
H_omega, H_vals = response.contour.imag, H(response.contour)
266331

267332
# Compute the describing function
268333
df = describing_function(F, A)
269334
N_vals = -1/df
270335

271-
# Now add the describing function curve to the plot
272-
plt.plot(N_vals.real, N_vals.imag)
273-
274336
# Look for intersection points
275-
intersections = []
337+
positions, intersections = [], []
276338
for i in range(N_vals.size - 1):
277339
for j in range(H_vals.size - 1):
278340
intersect = _find_intersection(
@@ -305,17 +367,114 @@ def _cost(x):
305367
else:
306368
a_final, omega_final = res.x[0], res.x[1]
307369

308-
# Add labels to the intersection points
309-
if isinstance(label, str):
310-
pos = H(1j * omega_final)
311-
plt.text(pos.real, pos.imag, label % (a_final, omega_final))
312-
elif label is not None or label is not False:
313-
raise ValueError("label must be formatting string or None")
370+
pos = H(1j * omega_final)
314371

315372
# Save the final estimate
373+
positions.append(pos)
316374
intersections.append((a_final, omega_final))
317375

318-
return intersections
376+
return DescribingFunctionResponse(
377+
response, N_vals, positions, intersections)
378+
379+
380+
def describing_function_plot(
381+
*sysdata, label="%5.2g @ %-5.2g", **kwargs):
382+
"""describing_function_plot(data, *args, **kwargs)
383+
384+
Plot a Nyquist plot with a describing function for a nonlinear system.
385+
386+
This function generates a Nyquist plot for a closed loop system
387+
consisting of a linear system with a static nonlinear function in the
388+
feedback path.
389+
390+
The function may be called in one of two forms:
391+
392+
describing_function_plot(response[, options])
393+
394+
describing_function_plot(H, F, A[, omega[, options]])
395+
396+
In the first form, the response should be generated using the
397+
:func:`~control.describing_function_response` function. In the second
398+
form, that function is called internally, with the listed arguments.
399+
400+
Parameters
401+
----------
402+
data : :class:`~control.DescribingFunctionData`
403+
A describing function response data object created by
404+
:func:`~control.describing_function_response`.
405+
H : LTI system
406+
Linear time-invariant (LTI) system (state space, transfer function, or
407+
FRD)
408+
F : static nonlinear function
409+
A static nonlinearity, either a scalar function or a single-input,
410+
single-output, static input/output system.
411+
A : list
412+
List of amplitudes to be used for the describing function plot.
413+
omega : list, optional
414+
List of frequencies to be used for the linear system Nyquist
415+
curve. If not specified (or None), frequencies are computed
416+
automatically based on the properties of the linear system.
417+
refine : bool, optional
418+
If True (default), refine the location of the intersection of the
419+
Nyquist curve for the linear system and the describing function to
420+
determine the intersection point
421+
label : str, optional
422+
Formatting string used to label intersection points on the Nyquist
423+
plot. Defaults to "%5.2g @ %-5.2g". Set to `None` to omit labels.
424+
425+
Returns
426+
-------
427+
lines : 1D array of Line2D
428+
Arrray of Line2D objects for each line in the plot. The first
429+
element of the array is a list of lines (typically only one) for
430+
the Nyquist plot of the linear I/O styem. The second element of
431+
the array is a list of lines (typically only one) for the
432+
describing function curve.
433+
434+
Examples
435+
--------
436+
>>> H_simple = ct.tf([8], [1, 2, 2, 1])
437+
>>> F_saturation = ct.saturation_nonlinearity(1)
438+
>>> amp = np.linspace(1, 4, 10)
439+
>>> lines = ct.describing_function_plot(H_simple, F_saturation, amp)
440+
441+
"""
442+
# Process keywords
443+
warn_nyquist = config._process_legacy_keyword(
444+
kwargs, 'warn', 'warn_nyquist', kwargs.pop('warn_nyquist', None))
445+
446+
if label not in (False, None) and not isinstance(label, str):
447+
raise ValueError("label must be formatting string, False, or None")
448+
449+
# Get the describing function response
450+
if len(sysdata) == 3:
451+
sysdata = sysdata + (None, ) # set omega to default value
452+
if len(sysdata) == 4:
453+
dfresp = describing_function_response(
454+
*sysdata, refine=kwargs.pop('refine', True),
455+
warn_nyquist=warn_nyquist)
456+
elif len(sysdata) == 1:
457+
dfresp = sysdata[0]
458+
else:
459+
raise TypeError("1, 3, or 4 position arguments required")
460+
461+
# Create a list of lines for the output
462+
out = np.empty(2, dtype=object)
463+
464+
# Plot the Nyquist response
465+
out[0] = dfresp.response.plot(**kwargs)[0]
466+
467+
# Add the describing function curve to the plot
468+
lines = plt.plot(dfresp.N_vals.real, dfresp.N_vals.imag)
469+
out[1] = lines
470+
471+
# Label the intersection points
472+
if label:
473+
for pos, (a, omega) in zip(dfresp.positions, dfresp.intersections):
474+
# Add labels to the intersection points
475+
plt.text(pos.real, pos.imag, label % (a, omega))
476+
477+
return out
319478

320479

321480
# Utility function to figure out whether two line segments intersection

0 commit comments

Comments
 (0)