Skip to content

Commit 2021a75

Browse files
murrayrmmurrayrm
committed
improvements to Nyquist contour tracing
* start tracing from omega = 0 * add support for indenting around imaginary poles * change return value to number of encirclements (+ contour, if desired) * update MATLAB interface to retain legacy format * new unit tests
1 parent 6de5e91 commit 2021a75

6 files changed

Lines changed: 349 additions & 52 deletions

File tree

‎control/freqplot.py‎

Lines changed: 129 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -46,11 +46,14 @@
4646
import matplotlib as mpl
4747
import matplotlib.pyplot as plt
4848
import numpy as np
49+
import warnings
4950

5051
from .ctrlutil import unwrap
5152
from .bdalg import feedback
5253
from .margins import stability_margins
5354
from .exception import ControlMIMONotImplemented
55+
from .statesp import StateSpace
56+
from .xferfcn import TransferFunction
5457
from . import config
5558

5659
__all__ = ['bode_plot', 'nyquist_plot', 'gangof4_plot',
@@ -195,7 +198,7 @@ def bode_plot(syslist, omega=None,
195198
if not hasattr(syslist, '__iter__'):
196199
syslist = (syslist,)
197200

198-
# decide whether to go above nyquist. freq
201+
# Decide whether to go above Nyquist frequency
199202
omega_range_given = True if omega is not None else False
200203

201204
if omega is None:
@@ -526,20 +529,29 @@ def gen_zero_centered_series(val_min, val_max, period):
526529
'nyquist.mirror_style': '--',
527530
'nyquist.arrows': 2,
528531
'nyquist.arrow_size': 8,
532+
'nyquist.indent_radius': 1e-1,
533+
'nyquist.indent_direction': 'right',
529534
}
530535

531536

532537
def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
533-
omega_num=None, label_freq=0, color=None, *args, **kwargs):
534-
"""
535-
Nyquist plot for a system
538+
omega_num=None, label_freq=0, color=None,
539+
return_contour=False, warn_nyquist=True, *args, **kwargs):
540+
"""Nyquist plot for a system
536541
537542
Plots a Nyquist plot for the system over a (optional) frequency range.
543+
The curve is computed by evaluating the Nyqist segment along the positive
544+
imaginary axis, with a mirror image generated to reflect the negative
545+
imaginary axis. Poles on or near the imaginary axis are avoided using a
546+
small indentation. The portion of the Nyquist contour at infinity is not
547+
explicitly computed (since it maps to a constant value for any system with
548+
a proper transfer function).
538549
539550
Parameters
540551
----------
541552
syslist : list of LTI
542-
List of linear input/output systems (single system is OK)
553+
List of linear input/output systems (single system is OK). Nyquist
554+
curves for each system are plotted on the same graph.
543555
544556
plot : boolean
545557
If True, plot magnitude
@@ -548,8 +560,8 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
548560
Set of frequencies to be evaluated, in rad/sec.
549561
550562
omega_limits : array_like of two values
551-
Limits to the range of frequencies. Ignored if omega
552-
is provided, and auto-generated if omitted.
563+
Limits to the range of frequencies. Ignored if omega is provided, and
564+
auto-generated if omitted.
553565
554566
omega_num : int
555567
Number of frequency samples to plot. Defaults to
@@ -563,6 +575,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
563575
omit completely. Default linestyle ('--') is determined by
564576
config.defaults['nyquist.mirror_style'].
565577
578+
return_contour : bool
579+
If 'True', return the contour used to evaluate the Nyquist plot.
580+
566581
label_freq : int
567582
Label every nth frequency on the plot. If not specified, no labels
568583
are generated.
@@ -584,6 +599,17 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
584599
arrow_style : matplotlib.patches.ArrowStyle
585600
Define style used for Nyquist curve arrows (overrides `arrow_size`).
586601
602+
indent_radius : float
603+
Amount to indent the Nyquist contour around poles that are at or near
604+
the imaginary axis.
605+
606+
indent_direction : str
607+
For poles on the imaginary axis, set the direction of indentation to
608+
be 'right' (default), 'left', or 'none'.
609+
610+
warn_nyquist : bool, optional
611+
If set to `False', turn off warnings about frequencies above Nyquist.
612+
587613
*args : :func:`matplotlib.pyplot.plot` positional properties, optional
588614
Additional arguments for `matplotlib` plots (color, linestyle, etc)
589615
@@ -592,67 +618,92 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
592618
593619
Returns
594620
-------
595-
real : ndarray (or list of ndarray if len(syslist) > 1))
596-
real part of the frequency response array
621+
count : int (or list of int if len(syslist) > 1)
622+
Number of encirclements of the point -1 by the Nyquist curve. If
623+
multiple systems are given, an array of counts is returned.
597624
598-
imag : ndarray (or list of ndarray if len(syslist) > 1))
599-
imaginary part of the frequency response array
625+
contour : ndarray (or list of ndarray if len(syslist) > 1)), optional
626+
The contour used to create the primary Nyquist curve segment. To
627+
obtain the Nyquist curve values, evaluate system(s) along contour.
600628
601-
omega : ndarray (or list of ndarray if len(syslist) > 1))
602-
frequencies in rad/s
629+
Notes
630+
-----
631+
1. If a discrete time model is given, the frequency response is computed
632+
along the upper branch of the unit circle, using the mapping ``z =
633+
exp(1j * omega * dt)`` where `omega` ranges from 0 to `pi/dt` and `dt`
634+
is the discrete timebase. If timebase not specified (``dt=True``),
635+
`dt` is set to 1.
636+
637+
2. If a continuous-time system contains poles on or near the imaginary
638+
axis, a small indentation will be used to avoid the pole. The radius
639+
of the indentation is given by `indent_radius` and it is taken the the
640+
right of stable poles and the left of unstable poles. If a pole is
641+
exactly on the imaginary axis, the `indent_direction` parameter can be
642+
used to set the direction of indentation. Setting `indent_direction`
643+
to `none` will turn off indentation. If `return_contour` is True, the
644+
exact contour used for evaluation is returned.
603645
604646
Examples
605647
--------
606648
>>> sys = ss([[1, -2], [3, -4]], [[5], [7]], [[6, 8]], [[9]])
607-
>>> real, imag, freq = nyquist_plot(sys)
649+
>>> count = nyquist_plot(sys)
608650
609651
"""
610652
# Check to see if legacy 'Plot' keyword was used
611653
if 'Plot' in kwargs:
612-
import warnings
613654
warnings.warn("'Plot' keyword is deprecated in nyquist_plot; "
614655
"use 'plot'", FutureWarning)
615656
# Map 'Plot' keyword to 'plot' keyword
616657
plot = kwargs.pop('Plot')
617658

618659
# Check to see if legacy 'labelFreq' keyword was used
619660
if 'labelFreq' in kwargs:
620-
import warnings
621661
warnings.warn("'labelFreq' keyword is deprecated in nyquist_plot; "
622662
"use 'label_freq'", FutureWarning)
623663
# Map 'labelFreq' keyword to 'label_freq' keyword
624664
label_freq = kwargs.pop('labelFreq')
625665

626666
# Check to see if legacy 'arrow_width' or 'arrow_length' were used
627667
if 'arrow_width' in kwargs or 'arrow_length' in kwargs:
628-
warn("'arrow_width' and 'arrow_length' keywords are deprecated in "
629-
"nyquist_plot; use `arrow_size` instead", FutureWarning)
668+
warnings.warn(
669+
"'arrow_width' and 'arrow_length' keywords are deprecated in "
670+
"nyquist_plot; use `arrow_size` instead", FutureWarning)
630671
kwargs['arrow_size'] = \
631-
(kwargs.get('arrow_width', 0) + kwargs.get('arrow_width', 0)) / 2
672+
(kwargs.get('arrow_width', 0) + kwargs.get('arrow_length', 0)) / 2
673+
kwargs.pop('arrow_width', False)
674+
kwargs.pop('arrow_length', False)
632675

633676
# Get values for params (and pop from list to allow keyword use in plot)
677+
omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)
634678
mirror_style = config._get_param(
635679
'nyquist', 'mirror_style', kwargs, _nyquist_defaults, pop=True)
636680
arrows = config._get_param(
637681
'nyquist', 'arrows', kwargs, _nyquist_defaults, pop=True)
638682
arrow_size = config._get_param(
639683
'nyquist', 'arrow_size', kwargs, _nyquist_defaults, pop=True)
640684
arrow_style = config._get_param('nyquist', 'arrow_style', kwargs, None)
685+
indent_radius = config._get_param(
686+
'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True)
687+
indent_direction = config._get_param(
688+
'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True)
641689

642690
# If argument was a singleton, turn it into a list
643-
if not hasattr(syslist, '__iter__'):
691+
isscalar = not hasattr(syslist, '__iter__')
692+
if isscalar:
644693
syslist = (syslist,)
645694

646-
# decide whether to go above nyquist. freq
695+
# Decide whether to go above Nyquist frequency
647696
omega_range_given = True if omega is not None else False
648697

698+
# Figure out the frequency limits
649699
if omega is None:
650-
omega_num = config._get_param(
651-
'freqplot', 'number_of_samples', omega_num)
652700
if omega_limits is None:
653701
# Select a default range if none is provided
654-
omega = _default_frequency_range(syslist,
655-
number_of_samples=omega_num)
702+
omega = _default_frequency_range(
703+
syslist, number_of_samples=omega_num)
704+
705+
# Replace first point with the origin
706+
omega[0] = 0
656707
else:
657708
omega_range_given = True
658709
omega_limits = np.asarray(omega_limits)
@@ -662,30 +713,69 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
662713
np.log10(omega_limits[1]), num=omega_num,
663714
endpoint=True)
664715

665-
xs, ys, omegas = [], [], []
716+
# Go through each system and keep track of the results
717+
counts, contours = [], []
666718
for sys in syslist:
667719
if not sys.issiso():
668720
# TODO: Add MIMO nyquist plots.
669721
raise ControlMIMONotImplemented(
670722
"Nyquist plot currently only supports SISO systems.")
671723

724+
# Figure out the frequency range
672725
omega_sys = np.asarray(omega)
726+
727+
# Determine the contour used to evaluate the Nyquist curve
673728
if sys.isdtime(strict=True):
729+
# Transform frequencies in for discrete-time systems
674730
nyquistfrq = math.pi / sys.dt
675731
if not omega_range_given:
676732
# limit up to and including nyquist frequency
677733
omega_sys = np.hstack((
678734
omega_sys[omega_sys < nyquistfrq], nyquistfrq))
679735

680-
mag, phase, omega_sys = sys.frequency_response(omega_sys)
736+
# Issue a warning if we are sampling above Nyquist
737+
if np.any(omega_sys * sys.dt > np.pi) and warn_nyquist:
738+
warnings.warn("evaluation above Nyquist frequency")
739+
740+
# Transform frequencies to continuous domain
741+
contour = np.exp(1j * omega * sys.dt)
742+
else:
743+
contour = 1j * omega_sys
744+
745+
# Bend the contour around any poles on/near the imaginary access
746+
if isinstance(sys, (StateSpace, TransferFunction)) and \
747+
sys.isctime() and indent_direction != 'none':
748+
poles = sys.pole()
749+
for i, s in enumerate(contour):
750+
# Find the nearest pole
751+
p = poles[(np.abs(poles - s)).argmin()]
752+
753+
# See if we need to indent around it
754+
if abs(s - p) < indent_radius:
755+
if p.real < 0 or \
756+
(p.real == 0 and indent_direction == 'right'):
757+
# Indent to the right
758+
contour[i] += \
759+
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
760+
elif p.real > 0 or \
761+
(p.real == 0 and indent_direction == 'left'):
762+
# Indent to the left
763+
contour[i] -= \
764+
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
765+
else:
766+
ValueError("unknown value for indent_direction")
767+
768+
# TODO: add code to indent around discrete poles on unit circle
681769

682770
# Compute the primary curve
683-
x = mag * np.cos(phase)
684-
y = mag * np.sin(phase)
771+
resp = sys(contour)
685772

686-
xs.append(x)
687-
ys.append(y)
688-
omegas.append(omega_sys)
773+
# Compute CW encirclements of -1 by integrating the (unwrapped) angle
774+
phase = -unwrap(np.angle(resp + 1))
775+
count = int(np.round(np.sum(np.diff(phase)) / np.pi, 0))
776+
777+
counts.append(count)
778+
contours.append(contour)
689779

690780
if plot:
691781
# Parse the arrows keyword
@@ -705,6 +795,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
705795
arrow_style = mpl.patches.ArrowStyle(
706796
'simple', head_width=arrow_size, head_length=arrow_size)
707797

798+
# Save the components of the response
799+
x, y = resp.real, resp.imag
800+
708801
# Plot the primary curve
709802
p = plt.plot(x, y, '-', color=color, *args, **kwargs)
710803
c = p[0].get_color()
@@ -751,10 +844,11 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
751844
ax.set_ylabel("Imaginary axis")
752845
ax.grid(color="lightgray")
753846

754-
if len(syslist) == 1:
755-
return xs[0], ys[0], omegas[0]
847+
# Return counts and (optionally) the contour we used
848+
if return_contour:
849+
return (counts[0], contours[0]) if isscalar else (counts, contours)
756850
else:
757-
return xs, ys, omegas
851+
return counts[0] if isscalar else counts
758852

759853

760854
# Internal function to add arrows to a curve

‎control/matlab/wrappers.py‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -105,7 +105,13 @@ def nyquist(*args, **kwargs):
105105
kwargs.update(other)
106106

107107
# Call the nyquist command
108-
return nyquist_plot(syslist, omega, *args, **kwargs)
108+
kwargs['return_contour'] = True
109+
_, contour = nyquist_plot(syslist, omega, *args, **kwargs)
110+
111+
# Create the MATLAB output arguments
112+
freqresp = syslist(contour)
113+
real, imag = freqresp.real, freqresp.imag
114+
return real, imag, contour.imag
109115

110116

111117
def _parse_freqplot_args(*args):

‎control/tests/freqresp_test.py‎

Lines changed: 12 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -75,11 +75,18 @@ def test_nyquist_basic(ss_siso):
7575
tf_siso = tf(ss_siso)
7676
nyquist_plot(ss_siso)
7777
nyquist_plot(tf_siso)
78-
assert len(nyquist_plot(tf_siso, plot=False, omega_num=20)[0] == 20)
79-
omega = nyquist_plot(tf_siso, plot=False, omega_limits=(1, 100))[2]
80-
assert_allclose(omega[0], 1)
81-
assert_allclose(omega[-1], 100)
82-
assert len(nyquist_plot(tf_siso, plot=False, omega=np.logspace(-1, 1, 10))[0])==10
78+
count, contour = nyquist_plot(
79+
tf_siso, plot=False, return_contour=True, omega_num=20)
80+
assert len(contour) == 20
81+
82+
count, contour = nyquist_plot(
83+
tf_siso, plot=False, omega_limits=(1, 100), return_contour=True)
84+
assert_allclose(contour[0], 1j)
85+
assert_allclose(contour[-1], 100j)
86+
87+
count, contour = nyquist_plot(
88+
tf_siso, plot=False, omega=np.logspace(-1, 1, 10), return_contour=True)
89+
assert len(contour) == 10
8390

8491

8592
@pytest.mark.filterwarnings("ignore:.*non-positive left xlim:UserWarning")

0 commit comments

Comments
 (0)