Skip to content

in nyquist plots, draw contour at specified magnitude if threshold is exceeded #671

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 6 commits into from
Closed

in nyquist plots, draw contour at specified magnitude if threshold is exceeded #671

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5a5ca77
nyquist plot above a certain radius is also plotted at specfied infin…
sawyerbfuller Nov 21, 2021
c008083
improve terminology and docstring
sawyerbfuller Nov 21, 2021
5b047f7
fix error when no large contours
sawyerbfuller Nov 21, 2021
372efb4
bugfixes
sawyerbfuller Nov 21, 2021
23b3e8c
new indent code that adds a consistent number of points at every inde…
sawyerbfuller Nov 22, 2021
5e9b179
rename keyword to plot_maximum_magintude and option to turn it off by…
sawyerbfuller Nov 23, 2021
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
136 changes: 96 additions & 40 deletions control/freqplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -326,8 +326,8 @@ def bode_plot(syslist, omega=None,

if nyquistfrq_plot:
# append data for vertical nyquist freq indicator line.
# if this extra nyquist lime is is plotted in a single plot
# command then line order is preserved when
# by including this in the same plot command, line order
# is preserved when
# creating a legend eg. legend(('sys1', 'sys2'))
omega_nyq_line = np.array(
(np.nan, nyquistfrq_plot, nyquistfrq_plot))
Expand Down Expand Up @@ -522,21 +522,23 @@ def gen_zero_centered_series(val_min, val_max, period):
'nyquist.mirror_style': '--',
'nyquist.arrows': 2,
'nyquist.arrow_size': 8,
'nyquist.indent_radius': 1e-1,
'nyquist.indent_radius': 1e-6,
'nyquist.plot_maximum_magnitude': 5,
'nyquist.indent_direction': 'right',
}


def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
omega_num=None, label_freq=0, color=None,
return_contour=False, warn_nyquist=True, *args, **kwargs):
return_contour=False, warn_nyquist=True,
*args, **kwargs):
"""Nyquist plot for a system

Plots a Nyquist plot for the system over a (optional) frequency range.
The curve is computed by evaluating the Nyqist segment along the positive
imaginary axis, with a mirror image generated to reflect the negative
imaginary axis. Poles on or near the imaginary axis are avoided using a
small indentation. The portion of the Nyquist contour at infinity is not
small indentation. If portions of the Nyquist plot The portion of the Nyquist contour at infinity is not
Comment on lines -539 to +541
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some unintended text insertion here, I guess.

explicitly computed (since it maps to a constant value for any system with
a proper transfer function).

Expand All @@ -550,7 +552,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
If True, plot magnitude

omega : array_like
Set of frequencies to be evaluated, in rad/sec.
Set of frequencies to be evaluated, in rad/sec. Remark: if omega is
specified, you may need to specify specify a larger `indent_radius`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

double specify

than the default to get reasonable contours.

omega_limits : array_like of two values
Limits to the range of frequencies. Ignored if omega is provided, and
Expand Down Expand Up @@ -592,6 +596,11 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
arrow_style : matplotlib.patches.ArrowStyle
Define style used for Nyquist curve arrows (overrides `arrow_size`).

plot_maximum_magnitude : float
If this value is not 0, an additional curve showing the path of
the Nyquist plot when it leaves the plot window is drawn at a radius
equal to this value.

indent_radius : float
Amount to indent the Nyquist contour around poles that are at or near
the imaginary axis.
Expand Down Expand Up @@ -679,11 +688,14 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True)
indent_direction = config._get_param(
'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True)
plot_maximum_magnitude = config._get_param(
'nyquist', 'plot_maximum_magnitude', kwargs, _nyquist_defaults, pop=True)

# If argument was a singleton, turn it into a list
if not hasattr(syslist, '__iter__'):
syslist = (syslist,)

omega_given = True if omega is not None else False
omega, omega_range_given = _determine_omega_vector(
syslist, omega, omega_limits, omega_num)
if not omega_range_given:
Expand All @@ -692,13 +704,13 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,

# Go through each system and keep track of the results
counts, contours = [], []
for sys in syslist:
for isys, sys in enumerate(syslist):
if not sys.issiso():
# TODO: Add MIMO nyquist plots.
raise ControlMIMONotImplemented(
"Nyquist plot currently only supports SISO systems.")

# Figure out the frequency range
# local copy of freq range that may be truncated above nyquist freq
omega_sys = np.asarray(omega)

# Determine the contour used to evaluate the Nyquist curve
Expand All @@ -717,46 +729,64 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
# do indentations in s-plane where it is more convenient
splane_contour = 1j * omega_sys

# Bend the contour around any poles on/near the imaginary axis
# TODO: smarter indent radius that depends on dcgain of system
# and timebase of discrete system.
# indent the contour around any poles on/near the imaginary axis
if isinstance(sys, (StateSpace, TransferFunction)) \
and indent_direction != 'none':
if sys.isctime():
splane_poles = sys.pole()
else:
# map z-plane poles to s-plane, ignoring any at the origin
# because we don't need to indent for them
# map z-plane poles to s-plane, ignoring any at z-plane origin
# that don't map to finite s-plane location, because we don't
# need to indent for them
zplane_poles = sys.pole()
zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)]
splane_poles = np.log(zplane_poles)/sys.dt

if splane_contour[1].imag > indent_radius \
and np.any(np.isclose(abs(splane_poles), 0)) \
and not omega_range_given:
# add some points for quarter circle around poles at origin
splane_contour = np.concatenate(
(1j * np.linspace(0., indent_radius, 50),
splane_contour[1:]))
for i, s in enumerate(splane_contour):
# Find the nearest pole
p = splane_poles[(np.abs(splane_poles - s)).argmin()]
# See if we need to indent around it
if abs(s - p) < indent_radius:
if p.real < 0 or (np.isclose(p.real, 0) \
and indent_direction == 'right'):
# Indent to the right
splane_contour[i] += \
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
elif p.real > 0 or (np.isclose(p.real, 0) \
and indent_direction == 'left'):
# Indent to the left
splane_contour[i] -= \
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
else:
ValueError("unknown value for indent_direction")

# change contour to z-plane if necessary
if omega_given:
# indent given contour
for i, s in enumerate(splane_contour):
# Find the nearest pole
p = splane_poles[(np.abs(splane_poles - s)).argmin()]
# See if we need to indent around it
if abs(s - p) < indent_radius:
if p.real < 0 or (np.isclose(p.real, 0) \
and indent_direction == 'right'):
# Indent to the right
splane_contour[i] += \
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
elif p.real > 0 or (np.isclose(p.real, 0) \
and indent_direction == 'left'):
# Indent to the left
splane_contour[i] -= \
np.sqrt(indent_radius ** 2 - (s-p).imag ** 2)
else:
ValueError("unknown value for indent_direction")
else:
# find poles that are near imag axis and replace contour
# near there with indented contour
if indent_direction == 'right':
indent = indent_radius * \
np.exp(1j * np.linspace(-np.pi/2, np.pi/2, 51))
elif indent_direction == 'left':
indent = -indent_radius * \
np.exp(1j * np.linspace(np.pi/2, -np.pi/2, 51))
else:
ValueError("unknown value for indent_direction")
for p in splane_poles[np.isclose(splane_poles.real, 0) \
& (splane_poles.imag >= 0)]:
start = np.searchsorted(
splane_contour.imag, p.imag - indent_radius)
end = np.searchsorted(
splane_contour.imag, p.imag + indent_radius)
# only keep indent parts that are > 0 and within contour
indent_mask = ((indent + p).imag >= 0) \
& ((indent + p).imag <= splane_contour[-1].imag)
splane_contour = np.concatenate((
splane_contour[:start],
indent[indent_mask] + p,
splane_contour[end:]))

# map contour to z-plane if necessary
if sys.isctime():
contour = splane_contour
else:
Expand All @@ -765,6 +795,15 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
# Compute the primary curve
resp = sys(contour)

if plot_maximum_magnitude != 0:
# compute large radius contour where it exceeds
# fill with nans that don't plot
large_mag_contour = np.full_like(contour, np.nan + 1j * np.nan)
# where too big, use angle of resp but specifified mag
too_big = abs(resp) > plot_maximum_magnitude
large_mag_contour[too_big] = plot_maximum_magnitude *\
(1+isys/20) * np.exp(1j * np.angle(resp[too_big]))

# Compute CW encirclements of -1 by integrating the (unwrapped) angle
phase = -unwrap(np.angle(resp + 1))
count = int(np.round(np.sum(np.diff(phase)) / np.pi, 0))
Expand Down Expand Up @@ -792,19 +831,31 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,

# Save the components of the response
x, y = resp.real, resp.imag
if plot_maximum_magnitude != 0:
x_lg, y_lg = large_mag_contour.real, large_mag_contour.imag

# Plot the primary curve
p = plt.plot(x, y, '-', color=color, *args, **kwargs)
c = p[0].get_color()
ax = plt.gca()
_add_arrows_to_line2D(
ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1)
if plot_maximum_magnitude != 0:
# plot large magnitude contour
p = plt.plot(x_lg, y_lg, color='red', *args, **kwargs)
_add_arrows_to_line2D(
ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1)

# Plot the mirror image
if mirror_style is not False:
p = plt.plot(x, -y, mirror_style, color=c, *args, **kwargs)
_add_arrows_to_line2D(
ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1)
if plot_maximum_magnitude != 0:
p = plt.plot(x_lg, -y_lg, mirror_style,
color='red', *args, **kwargs)
_add_arrows_to_line2D(
ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1)

# Mark the -1 point
plt.plot([-1], [0], 'r+')
Expand Down Expand Up @@ -838,6 +889,9 @@ def nyquist_plot(syslist, omega=None, plot=True, omega_limits=None,
ax.set_xlabel("Real axis")
ax.set_ylabel("Imaginary axis")
ax.grid(color="lightgray")
if plot_maximum_magnitude != 0:
ax.set_xlim(-plot_maximum_magnitude*1.1,plot_maximum_magnitude*1.1)
ax.set_ylim(-plot_maximum_magnitude*1.1,plot_maximum_magnitude*1.1)

# "Squeeze" the results
if len(syslist) == 1:
Expand Down Expand Up @@ -873,7 +927,9 @@ def _add_arrows_to_line2D(
if not isinstance(line, mpl.lines.Line2D):
raise ValueError("expected a matplotlib.lines.Line2D object")
x, y = line.get_xdata(), line.get_ydata()

x, y = x[np.isfinite(x)], y[np.isfinite(x)]
if len(x) in (0, 1):
return []
arrow_kw = {
"arrowstyle": arrowstyle,
}
Expand Down
34 changes: 25 additions & 9 deletions control/tests/nyquist_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,8 @@ def test_nyquist_basic():

# Make sure that we can turn off frequency modification
count, contour_indented = ct.nyquist_plot(
sys, np.linspace(1e-4, 1e2, 100), return_contour=True)
sys, np.linspace(1e-4, 1e2, 100), indent_radius=0.01,
return_contour=True)
assert not all(contour_indented.real == 0)
count, contour = ct.nyquist_plot(
sys, np.linspace(1e-4, 1e2, 100), return_contour=True,
Expand All @@ -87,14 +88,17 @@ def test_nyquist_basic():
assert _Z(sys) == count + _P(sys)

# Nyquist plot with poles on imaginary axis, omega specified
# when specifying omega, usually need to specify larger indent radius
sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1])
count = ct.nyquist_plot(sys, np.linspace(1e-3, 1e1, 1000))
count = ct.nyquist_plot(sys, np.linspace(1e-3, 1e1, 1000),
indent_radius=0.1)
assert _Z(sys) == count + _P(sys)

# Nyquist plot with poles on imaginary axis, omega specified, with contour
sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1])
count, contour = ct.nyquist_plot(
sys, np.linspace(1e-3, 1e1, 1000), return_contour=True)
sys, np.linspace(1e-3, 1e1, 1000), indent_radius=0.1,
return_contour=True)
assert _Z(sys) == count + _P(sys)

# Nyquist plot with poles on imaginary axis, return contour
Expand Down Expand Up @@ -137,12 +141,23 @@ def test_nyquist_fbs_examples():
count = ct.nyquist_plot(sys)
assert _Z(sys) == count + _P(sys)

# when omega is specified, no points are added at indentations, leading
# to incorrect encirclement count
plt.figure()
plt.title("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2) [zoom]")
count = ct.nyquist_plot(sys, omega=np.linspace(1.5, 1e3, 1000))
# assert _Z(sys) == count + _P(sys)
assert count == -1

# when only the range of the frequency is provided, this permits
# extra points to be added at indentations, leading to correct count
plt.figure()
plt.title("Figure 10.10: L(s) = 3 (s+6)^2 / (s (s+1)^2) [zoom]")
count = ct.nyquist_plot(sys, omega_limits=[1.5, 1e3])
# Frequency limits for zoom give incorrect encirclement count
# assert _Z(sys) == count + _P(sys)
assert count == -1
assert count == 0



@pytest.mark.parametrize("arrows", [
Expand Down Expand Up @@ -193,13 +208,10 @@ def test_nyquist_indent():
plt.title("Pole at origin; indent_radius=default")
assert _Z(sys) == count + _P(sys)

# first value of default omega vector was 0.1, replaced by 0. for contour
# indent_radius is larger than 0.1 -> no extra quater circle around origin
# confirm that first element of indent is properly indented
count, contour = ct.nyquist_plot(sys, plot=False, indent_radius=.1007,
return_contour=True)
np.testing.assert_allclose(contour[0], .1007+0.j)
# second value of omega_vector is larger than indent_radius: not indented
assert np.all(contour.real[2:] == 0.)

plt.figure();
count, contour = ct.nyquist_plot(sys, indent_radius=0.01,
Expand All @@ -208,7 +220,7 @@ def test_nyquist_indent():
assert _Z(sys) == count + _P(sys)
# indent radius is smaller than the start of the default omega vector
# check that a quarter circle around the pole at origin has been added.
np.testing.assert_allclose(contour[:50].real**2 + contour[:50].imag**2,
np.testing.assert_allclose(contour[:25].real**2 + contour[:25].imag**2,
0.01**2)

plt.figure();
Expand All @@ -226,6 +238,10 @@ def test_nyquist_indent():
plt.title("Imaginary poles; encirclements = %d" % count)
assert _Z(sys) == count + _P(sys)

# confirm we can turn off maximum_magitude plot
plt.figure();
count = ct.nyquist_plot(sys, plot_maximum_magnitude=0)

# Imaginary poles with indentation to the left
plt.figure();
count = ct.nyquist_plot(sys, indent_direction='left', label_freq=300)
Expand Down