Skip to content

Update Nyquist rescaling + other improvements #1155

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.

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

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
5 changes: 4 additions & 1 deletion control/ctrlplot.py
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,10 @@
# # Customize axes (curvilinear grids, shared axes, etc)
#
# # Plot the data
# lines = np.full(ax_array.shape, [])
# lines = np.empty(ax_array.shape, dtype=object)
# for i in range(ax_array.shape[0]):
# for j in range(ax_array.shape[1]):
# lines[i, j] = []
# line_labels = _process_line_labels(label, ntraces, nrows, ncols)
# color_offset, color_cycle = _get_color_offset(ax)
# for i, j in itertools.product(range(nrows), range(ncols)):
Expand Down
3 changes: 2 additions & 1 deletion control/descfcn.py
Original file line number Diff line number Diff line change
Expand Up @@ -521,7 +521,8 @@ def describing_function_plot(
# Plot the Nyquist response
cplt = dfresp.response.plot(**kwargs)
ax = cplt.axes[0, 0] # Get the axes where the plot was made
lines[0] = cplt.lines[0] # Return Nyquist lines for first system
lines[0] = np.concatenate( # Return Nyquist lines for first system
cplt.lines.flatten()).tolist()

# Add the describing function curve to the plot
lines[1] = ax.plot(dfresp.N_vals.real, dfresp.N_vals.imag)
Expand Down
77 changes: 52 additions & 25 deletions control/freqplot.py
Original file line number Diff line number Diff line change
Expand Up @@ -1100,13 +1100,14 @@ def gen_zero_centered_series(val_min, val_max, period):
_nyquist_defaults = {
'nyquist.primary_style': ['-', '-.'], # style for primary curve
'nyquist.mirror_style': ['--', ':'], # style for mirror curve
'nyquist.arrows': 2, # number of arrows around curve
'nyquist.arrows': 3, # number of arrows around curve
'nyquist.arrow_size': 8, # pixel size for arrows
'nyquist.encirclement_threshold': 0.05, # warning threshold
'nyquist.indent_radius': 1e-4, # indentation radius
'nyquist.indent_direction': 'right', # indentation direction
'nyquist.indent_points': 50, # number of points to insert
'nyquist.max_curve_magnitude': 20, # clip large values
'nyquist.indent_points': 200, # number of points to insert
'nyquist.max_curve_magnitude': 15, # rescale large values
'nyquist.blend_fraction': 0.15, # when to start scaling
'nyquist.max_curve_offset': 0.02, # offset of primary/mirror
'nyquist.start_marker': 'o', # marker at start of curve
'nyquist.start_marker_size': 4, # size of the marker
Expand Down Expand Up @@ -1638,6 +1639,10 @@ def nyquist_plot(
The matplotlib axes to draw the figure on. If not specified and
the current figure has a single axes, that axes is used.
Otherwise, a new figure is created.
blend_fraction : float, optional
For portions of the Nyquist curve that are scaled to have a maximum
magnitude of `max_curve_magnitude`, begin a smooth rescaling at
this fraction of `max_curve_magnitude`. Default value is 0.15.
encirclement_threshold : float, optional
Define the threshold for generating a warning if the number of net
encirclements is a non-integer value. Default value is 0.05 and can
Expand Down Expand Up @@ -1751,8 +1756,8 @@ def nyquist_plot(
to avoid poles, resulting in a scaling of the Nyquist plot, the line
styles are according to the settings of the `primary_style` and
`mirror_style` keywords. By default the scaled portions of the primary
curve use a dotted line style and the scaled portion of the mirror
image use a dashdot line style.
curve use a dashdot line style and the scaled portions of the mirror
image use a dotted line style.

Examples
--------
Expand Down Expand Up @@ -1784,6 +1789,8 @@ def nyquist_plot(
ax_user = ax
max_curve_magnitude = config._get_param(
'nyquist', 'max_curve_magnitude', kwargs, _nyquist_defaults, pop=True)
blend_fraction = config._get_param(
'nyquist', 'blend_fraction', kwargs, _nyquist_defaults, pop=True)
max_curve_offset = config._get_param(
'nyquist', 'max_curve_offset', kwargs, _nyquist_defaults, pop=True)
rcParams = config._get_param('ctrlplot', 'rcParams', kwargs, pop=True)
Expand Down Expand Up @@ -1878,10 +1885,16 @@ def _parse_linestyle(style_name, allow_false=False):
legend_loc, _, show_legend = _process_legend_keywords(
kwargs, None, 'upper right')

# Figure out where the blended curve should start
if blend_fraction < 0 or blend_fraction > 1:
raise ValueError("blend_fraction must be between 0 and 1")
blend_curve_start = (1 - blend_fraction) * max_curve_magnitude

# Create a list of lines for the output
out = np.empty(len(nyquist_responses), dtype=object)
for i in range(out.shape[0]):
out[i] = [] # unique list in each element
out = np.empty((len(nyquist_responses), 4), dtype=object)
for i in range(len(nyquist_responses)):
for j in range(4):
out[i, j] = [] # unique list in each element

for idx, response in enumerate(nyquist_responses):
resp = response.response
Expand All @@ -1892,20 +1905,31 @@ def _parse_linestyle(style_name, allow_false=False):

# Find the different portions of the curve (with scaled pts marked)
reg_mask = np.logical_or(
np.abs(resp) > max_curve_magnitude,
splane_contour.real != 0)
# reg_mask = np.logical_or(
# np.abs(resp.real) > max_curve_magnitude,
# np.abs(resp.imag) > max_curve_magnitude)
np.abs(resp) > blend_curve_start,
np.logical_not(np.isclose(splane_contour.real, 0)))

scale_mask = ~reg_mask \
& np.concatenate((~reg_mask[1:], ~reg_mask[-1:])) \
& np.concatenate((~reg_mask[0:1], ~reg_mask[:-1]))

# Rescale the points with large magnitude
rescale = np.logical_and(
reg_mask, abs(resp) > max_curve_magnitude)
resp[rescale] *= max_curve_magnitude / abs(resp[rescale])
rescale_idx = (np.abs(resp) > blend_curve_start)

if np.any(rescale_idx): # Only process if rescaling is needed
subset = resp[rescale_idx]
abs_subset = np.abs(subset)
unit_vectors = subset / abs_subset # Preserve phase/direction

if blend_curve_start == max_curve_magnitude:
# Clip at max_curve_magnitude
resp[rescale_idx] = max_curve_magnitude * unit_vectors
else:
# Logistic scaling
newmag = blend_curve_start + \
(max_curve_magnitude - blend_curve_start) * \
(abs_subset - blend_curve_start) / \
(abs_subset + max_curve_magnitude - 2 * blend_curve_start)
resp[rescale_idx] = newmag * unit_vectors

# Get the label to use for the line
label = response.sysname if line_labels is None else line_labels[idx]
Expand All @@ -1916,7 +1940,7 @@ def _parse_linestyle(style_name, allow_false=False):
p = ax.plot(
x_reg, y_reg, primary_style[0], color=color, label=label, **kwargs)
c = p[0].get_color()
out[idx] += p
out[idx, 0] += p

# Figure out how much to offset the curve: the offset goes from
# zero at the start of the scaled section to max_curve_offset as
Expand All @@ -1928,12 +1952,12 @@ def _parse_linestyle(style_name, allow_false=False):
x_scl = np.ma.masked_where(scale_mask, resp.real)
y_scl = np.ma.masked_where(scale_mask, resp.imag)
if x_scl.count() >= 1 and y_scl.count() >= 1:
out[idx] += ax.plot(
out[idx, 1] += ax.plot(
x_scl * (1 + curve_offset),
y_scl * (1 + curve_offset),
primary_style[1], color=c, **kwargs)
else:
out[idx] += [None]
out[idx, 1] += [None]

# Plot the primary curve (invisible) for setting arrows
x, y = resp.real.copy(), resp.imag.copy()
Expand All @@ -1948,15 +1972,15 @@ def _parse_linestyle(style_name, allow_false=False):
# Plot the mirror image
if mirror_style is not False:
# Plot the regular and scaled segments
out[idx] += ax.plot(
out[idx, 2] += ax.plot(
x_reg, -y_reg, mirror_style[0], color=c, **kwargs)
if x_scl.count() >= 1 and y_scl.count() >= 1:
out[idx] += ax.plot(
out[idx, 3] += ax.plot(
x_scl * (1 - curve_offset),
-y_scl * (1 - curve_offset),
mirror_style[1], color=c, **kwargs)
else:
out[idx] += [None]
out[idx, 3] += [None]

# Add the arrows (on top of an invisible contour)
x, y = resp.real.copy(), resp.imag.copy()
Expand All @@ -1966,12 +1990,15 @@ def _parse_linestyle(style_name, allow_false=False):
_add_arrows_to_line2D(
ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1)
else:
out[idx] += [None, None]
out[idx, 2] += [None]
out[idx, 3] += [None]

# Mark the start of the curve
if start_marker:
ax.plot(resp[0].real, resp[0].imag, start_marker,
color=c, markersize=start_marker_size)
segment = 0 if 0 in rescale_idx else 1 # regular vs scaled
out[idx, segment] += ax.plot(
resp[0].real, resp[0].imag, start_marker,
color=c, markersize=start_marker_size)

# Mark the -1 point
ax.plot([-1], [0], 'r+')
Expand Down
4 changes: 2 additions & 2 deletions control/tests/descfcn_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -190,11 +190,11 @@ def test_describing_function_plot():

cplt = response.plot()
assert len(plt.gcf().get_axes()) == 1 # make sure there is a plot
assert len(cplt.lines[0]) == 4 and len(cplt.lines[1]) == 1
assert len(cplt.lines[0]) == 5 and len(cplt.lines[1]) == 1

# Call plot directly
cplt = ct.describing_function_plot(H_larger, F_saturation, amp, omega)
assert len(cplt.lines[0]) == 4 and len(cplt.lines[1]) == 1
assert len(cplt.lines[0]) == 5 and len(cplt.lines[1]) == 1


def test_describing_function_exceptions():
Expand Down
76 changes: 70 additions & 6 deletions control/tests/nyquist_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -291,7 +291,7 @@ def test_nyquist_indent_default(indentsys):

def test_nyquist_indent_dont(indentsys):
# 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
# indent_radius is larger than 0.1 -> no extra quarter circle around origin
with pytest.warns() as record:
count, contour = ct.nyquist_response(
indentsys, omega=[0, 0.2, 0.3, 0.4], indent_radius=.1007,
Expand Down Expand Up @@ -406,15 +406,16 @@ def test_linestyle_checks():
# Set the line styles
cplt = ct.nyquist_plot(
sys, primary_style=[':', ':'], mirror_style=[':', ':'])
assert all([line.get_linestyle() == ':' for line in cplt.lines[0]])
assert all([lines[0].get_linestyle() == ':' for lines in cplt.lines[0, :]])

# Set the line colors
cplt = ct.nyquist_plot(sys, color='g')
assert all([line.get_color() == 'g' for line in cplt.lines[0]])
assert all([line.get_color() == 'g' for line in cplt.lines[0, 0]])

# Turn off the mirror image
cplt = ct.nyquist_plot(sys, mirror_style=False)
assert cplt.lines[0][2:] == [None, None]
assert cplt.lines[0, 2] == [None]
assert cplt.lines[0, 3] == [None]

with pytest.raises(ValueError, match="invalid 'primary_style'"):
ct.nyquist_plot(sys, primary_style=False)
Expand All @@ -428,6 +429,7 @@ def test_linestyle_checks():
ct.nyquist_plot(sys, primary_style=':', mirror_style='-.')

@pytest.mark.usefixtures("editsdefaults")
@pytest.mark.xfail(reason="updated code avoids warning")
def test_nyquist_legacy():
ct.use_legacy_defaults('0.9.1')

Expand Down Expand Up @@ -526,6 +528,42 @@ def test_no_indent_pole():
sys, warn_encirclements=False, indent_direction='none')


def test_nyquist_rescale():
sys = 2 * ct.tf([1], [1, 1]) * ct.tf([1], [1, 0])**2
sys.name = 'How example'

# Default case
resp = ct.nyquist_response(sys, indent_direction='left')
cplt = resp.plot(label='default [0.15]')
assert len(cplt.lines[0, 0]) == 2
assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)])

# Sharper corner
cplt = ct.nyquist_plot(
sys*4, indent_direction='left',
max_curve_magnitude=17, blend_fraction=0.05, label='fraction=0.05')
assert len(cplt.lines[0, 0]) == 2
assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)])

# More gradual corner
cplt = ct.nyquist_plot(
sys*0.25, indent_direction='left',
max_curve_magnitude=13, blend_fraction=0.25, label='fraction=0.25')
assert len(cplt.lines[0, 0]) == 2
assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)])

# No corner
cplt = ct.nyquist_plot(
sys*12, indent_direction='left',
max_curve_magnitude=19, blend_fraction=0, label='fraction=0')
assert len(cplt.lines[0, 0]) == 2
assert all([len(cplt.lines[0, i]) == 1 for i in range(1, 4)])

# Bad value
with pytest.raises(ValueError, match="blend_fraction must be between"):
ct.nyquist_plot(sys, indent_direction='left', blend_fraction=1.2)


if __name__ == "__main__":
#
# Interactive mode: generate plots for manual viewing
Expand Down Expand Up @@ -566,8 +604,8 @@ def test_no_indent_pole():
sys = 3 * (s+6)**2 / (s * (s**2 + 1e-4 * s + 1))
plt.figure()
ct.nyquist_plot(sys)
ct.nyquist_plot(sys, max_curve_magnitude=15)
ct.nyquist_plot(sys, indent_radius=1e-6, max_curve_magnitude=25)
ct.nyquist_plot(sys, max_curve_magnitude=10)
ct.nyquist_plot(sys, indent_radius=1e-6, max_curve_magnitude=20)

print("Unusual Nyquist plot")
sys = ct.tf([1], [1, 3, 2]) * ct.tf([1], [1, 0, 1])
Expand Down Expand Up @@ -595,3 +633,29 @@ def test_no_indent_pole():
plt.figure()
cplt = ct.nyquist_plot([sys, sys1, sys2])
cplt.set_plot_title("Mixed FRD, tf data")

plt.figure()
print("Jon How example")
test_nyquist_rescale()

#
# Save the figures in a PDF file for later comparisons
#
import subprocess
from matplotlib.backends.backend_pdf import PdfPages
from datetime import date

# Create the file to store figures
try:
git_info = subprocess.check_output(
['git', 'describe'], text=True).strip()
except subprocess.CalledProcessError:
git_info = 'UNKNOWN-REPO-INFO'
pdf = PdfPages(
f'nyquist_gallery-{git_info}-{date.today().isoformat()}.pdf')

# Go through each figure and save it
for fignum in plt.get_fignums():
pdf.savefig(plt.figure(fignum))

pdf.close()
Binary file modified doc/figures/freqplot-nyquist-custom.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified doc/figures/freqplot-nyquist-default.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
2 changes: 1 addition & 1 deletion doc/response.rst
Original file line number Diff line number Diff line change
Expand Up @@ -547,7 +547,7 @@ the computation of the Nyquist curve and the way the data are plotted:
sys = ct.tf([1, 0.2], [1, 0, 1]) * ct.tf([1], [1, 0])
nyqresp = ct.nyquist_response(sys)
nyqresp.plot(
max_curve_magnitude=6, max_curve_offset=1,
max_curve_magnitude=6, max_curve_offset=1, blend_fraction=0.05,
arrows=[0, 0.15, 0.3, 0.6, 0.7, 0.925],
title='Custom Nyquist plot')
print("Encirclements =", nyqresp.count)
Expand Down
Loading
Loading