Skip to content

Commit c383014

Browse files
murrayrmmurrayrm
committed
add frequency_limit and line style processing + unit tests, fixes
1 parent de5ea5a commit c383014

File tree

4 files changed

+179
-75
lines changed

4 files changed

+179
-75
lines changed

‎control/freqplot.py

+98-64
Original file line numberDiff line numberDiff line change
@@ -7,14 +7,6 @@
77
# Nyquist plots and other frequency response plots. The code for Nichols
88
# charts is in nichols.py. The code for pole-zero diagrams is in pzmap.py
99
# and rlocus.py.
10-
#
11-
# Functionality to add/check (Jul 2023, working list)
12-
# [?] Allow line colors/styles to be set in plot() command (also time plots)
13-
# [ ] Get sisotool working in iPython and document how to make it work
14-
# [ ] Allow use of subplot labels instead of output/input subtitles
15-
# [i] Allow frequency range to be overridden in bode_plot
16-
# [i] Unit tests for discrete time systems with different sample times
17-
# [ ] Add unit tests for ct.config.defaults['freqplot_number_of_samples']
1810

1911
import numpy as np
2012
import matplotlib as mpl
@@ -103,7 +95,7 @@ def bode_plot(
10395
overlay_outputs=None, overlay_inputs=None, phase_label=None,
10496
magnitude_label=None, display_margins=None,
10597
margins_method='best', legend_map=None, legend_loc=None,
106-
sharex=None, sharey=None, title=None, relabel=True, **kwargs):
98+
sharex=None, sharey=None, title=None, **kwargs):
10799
"""Bode plot for a system.
108100
109101
Plot the magnitude and phase of the frequency response over a
@@ -116,7 +108,8 @@ def bode_plot(
116108
single system or frequency response can also be passed.
117109
omega : array_like, optoinal
118110
List of frequencies in rad/sec over to plot over. If not specified,
119-
this will be determined from the proporties of the systems.
111+
this will be determined from the proporties of the systems. Ignored
112+
if `data` is not a list of systems.
120113
*fmt : :func:`matplotlib.pyplot.plot` format string, optional
121114
Passed to `matplotlib` as the format string for all lines in the plot.
122115
The `omega` parameter must be present (use omega=None if needed).
@@ -147,16 +140,6 @@ def bode_plot(
147140
148141
Other Parameters
149142
----------------
150-
plot : bool, optional
151-
(legacy) If given, `bode_plot` returns the legacy return values
152-
of magnitude, phase, and frequency. If False, just return the
153-
values with no plot.
154-
omega_limits : array_like of two values
155-
Limits of the to generate frequency vector. If Hz=True the limits
156-
are in Hz otherwise in rad/s.
157-
omega_num : int
158-
Number of samples to plot. Defaults to
159-
config.defaults['freqplot.number_of_samples'].
160143
grid : bool
161144
If True, plot grid lines on gain and phase plots. Default is set by
162145
`config.defaults['freqplot.grid']`.
@@ -166,6 +149,20 @@ def bode_plot(
166149
value specified. Units are in either degrees or radians, depending on
167150
the `deg` parameter. Default is -180 if wrap_phase is False, 0 if
168151
wrap_phase is True.
152+
omega_limits : array_like of two values
153+
Set limits for plotted frequency range. If Hz=True the limits
154+
are in Hz otherwise in rad/s.
155+
omega_num : int
156+
Number of samples to use for the frequeny range. Defaults to
157+
config.defaults['freqplot.number_of_samples']. Ignore if data is
158+
not a list of systems.
159+
plot : bool, optional
160+
(legacy) If given, `bode_plot` returns the legacy return values
161+
of magnitude, phase, and frequency. If False, just return the
162+
values with no plot.
163+
rcParams : dict
164+
Override the default parameters used for generating plots.
165+
Default is set up config.default['freqplot.rcParams'].
169166
wrap_phase : bool or float
170167
If wrap_phase is `False` (default), then the phase will be unwrapped
171168
so that it is continuously increasing or decreasing. If wrap_phase is
@@ -220,7 +217,6 @@ def bode_plot(
220217
'freqplot', 'wrap_phase', kwargs, _freqplot_defaults, pop=True)
221218
initial_phase = config._get_param(
222219
'freqplot', 'initial_phase', kwargs, None, pop=True)
223-
omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)
224220
freqplot_rcParams = config._get_param(
225221
'freqplot', 'rcParams', kwargs, _freqplot_defaults, pop=True)
226222

@@ -262,7 +258,7 @@ def bode_plot(
262258
data = [data]
263259

264260
#
265-
# Pre-process the data to be plotted (unwrap phase)
261+
# Pre-process the data to be plotted (unwrap phase, limit frequencies)
266262
#
267263
# To maintain compatibility with legacy uses of bode_plot(), we do some
268264
# initial processing on the data, specifically phase unwrapping and
@@ -277,6 +273,17 @@ def bode_plot(
277273
data = frequency_response(
278274
data, omega=omega, omega_limits=omega_limits,
279275
omega_num=omega_num, Hz=Hz)
276+
else:
277+
# Generate warnings if frequency keywords were given
278+
if omega_num is not None:
279+
warnings.warn("`omega_num` ignored when passed response data")
280+
elif omega is not None:
281+
warnings.warn("`omega` ignored when passed response data")
282+
283+
# Check to make sure omega_limits is sensible
284+
if omega_limits is not None and \
285+
(len(omega_limits) != 2 or omega_limits[1] <= omega_limits[0]):
286+
raise ValueError(f"invalid limits: {omega_limits=}")
280287

281288
# If plot_phase is not specified, check the data first, otherwise true
282289
if plot_phase is None:
@@ -288,7 +295,6 @@ def bode_plot(
288295

289296
mag_data, phase_data, omega_data = [], [], []
290297
for response in data:
291-
phase = response.phase.copy()
292298
noutputs, ninputs = response.noutputs, response.ninputs
293299

294300
if initial_phase is None:
@@ -306,9 +312,9 @@ def bode_plot(
306312
raise ValueError("initial_phase must be a number.")
307313

308314
# Reshape the phase to allow standard indexing
309-
phase = phase.reshape((noutputs, ninputs, -1))
315+
phase = response.phase.copy().reshape((noutputs, ninputs, -1))
310316

311-
# Shift and wrap
317+
# Shift and wrap the phase
312318
for i, j in itertools.product(range(noutputs), range(ninputs)):
313319
# Shift the phase if needed
314320
if abs(phase[i, j, 0] - initial_phase_value) > math.pi:
@@ -641,7 +647,7 @@ def _make_line_label(response, output_index, input_index):
641647
# Get the (pre-processed) data in fully indexed form
642648
mag = mag_data[index].reshape((noutputs, ninputs, -1))
643649
phase = phase_data[index].reshape((noutputs, ninputs, -1))
644-
omega_sys, sysname = response.omega, response.sysname
650+
omega_sys, sysname = omega_data[index], response.sysname
645651

646652
for i, j in itertools.product(range(noutputs), range(ninputs)):
647653
# Get the axes to use for magnitude and phase
@@ -831,21 +837,17 @@ def _make_line_label(response, output_index, input_index):
831837

832838
for i in range(noutputs):
833839
for j in range(ninputs):
840+
# Utility function to generate phase labels
834841
def gen_zero_centered_series(val_min, val_max, period):
835842
v1 = np.ceil(val_min / period - 0.2)
836843
v2 = np.floor(val_max / period + 0.2)
837844
return np.arange(v1, v2 + 1) * period
838845

839-
# TODO: put Nyquist lines here?
840-
841-
# TODO: what is going on here
842-
# TODO: fix to use less dense labels, when needed
843-
# TODO: make sure turning sharey on and off makes labels come/go
846+
# Label the phase axes using multiples of 45 degrees
844847
if plot_phase:
845848
ax_phase = ax_array[phase_map[i, j]]
846849

847850
# Set the labels
848-
# TODO: tighten up code
849851
if deg:
850852
ylim = ax_phase.get_ylim()
851853
num = np.floor((ylim[1] - ylim[0]) / 45)
@@ -877,6 +879,11 @@ def gen_zero_centered_series(val_min, val_max, period):
877879
if share_frequency in [True, 'all', 'col']:
878880
ax_array[i, j].tick_params(labelbottom=False)
879881

882+
# If specific omega_limits were given, use them
883+
if omega_limits is not None:
884+
for i, j in itertools.product(range(nrows), range(ncols)):
885+
ax_array[i, j].set_xlim(omega_limits)
886+
880887
#
881888
# Update the plot title (= figure suptitle)
882889
#
@@ -895,7 +902,7 @@ def gen_zero_centered_series(val_min, val_max, period):
895902
else:
896903
title = data[0].title
897904

898-
if fig is not None and title is not None:
905+
if fig is not None and isinstance(title, str):
899906
# Get the current title, if it exists
900907
old_title = None if fig._suptitle is None else fig._suptitle._text
901908
new_title = title
@@ -1294,11 +1301,11 @@ def nyquist_response(
12941301
# Determine the contour used to evaluate the Nyquist curve
12951302
if sys.isdtime(strict=True):
12961303
# Restrict frequencies for discrete-time systems
1297-
nyquistfrq = math.pi / sys.dt
1304+
nyq_freq = math.pi / sys.dt
12981305
if not omega_range_given:
12991306
# limit up to and including Nyquist frequency
13001307
omega_sys = np.hstack((
1301-
omega_sys[omega_sys < nyquistfrq], nyquistfrq))
1308+
omega_sys[omega_sys < nyq_freq], nyq_freq))
13021309

13031310
# Issue a warning if we are sampling above Nyquist
13041311
if np.any(omega_sys * sys.dt > np.pi) and warn_nyquist:
@@ -1817,7 +1824,7 @@ def _parse_linestyle(style_name, allow_false=False):
18171824
x, y = resp.real.copy(), resp.imag.copy()
18181825
x[reg_mask] *= (1 + curve_offset[reg_mask])
18191826
y[reg_mask] *= (1 + curve_offset[reg_mask])
1820-
p = plt.plot(x, y, linestyle='None', color=c, **kwargs)
1827+
p = plt.plot(x, y, linestyle='None', color=c)
18211828

18221829
# Add arrows
18231830
ax = plt.gca()
@@ -2210,10 +2217,6 @@ def singular_values_plot(
22102217
Hz : bool
22112218
If True, plot frequency in Hz (omega must be provided in rad/sec).
22122219
Default value (False) set by config.defaults['freqplot.Hz'].
2213-
plot : bool, optional
2214-
(legacy) If given, `singular_values_plot` returns the legacy return
2215-
values of magnitude, phase, and frequency. If False, just return
2216-
the values with no plot.
22172220
legend_loc : str, optional
22182221
For plots with multiple lines, a legend will be included in the
22192222
given location. Default is 'center right'. Use False to supress.
@@ -2233,6 +2236,26 @@ def singular_values_plot(
22332236
omega : ndarray (or list of ndarray if len(data) > 1))
22342237
If plot=False, frequency in rad/sec (deprecated).
22352238
2239+
Other Parameters
2240+
----------------
2241+
grid : bool
2242+
If True, plot grid lines on gain and phase plots. Default is set by
2243+
`config.defaults['freqplot.grid']`.
2244+
omega_limits : array_like of two values
2245+
Set limits for plotted frequency range. If Hz=True the limits
2246+
are in Hz otherwise in rad/s.
2247+
omega_num : int
2248+
Number of samples to use for the frequeny range. Defaults to
2249+
config.defaults['freqplot.number_of_samples']. Ignore if data is
2250+
not a list of systems.
2251+
plot : bool, optional
2252+
(legacy) If given, `singular_values_plot` returns the legacy return
2253+
values of magnitude, phase, and frequency. If False, just return
2254+
the values with no plot.
2255+
rcParams : dict
2256+
Override the default parameters used for generating plots.
2257+
Default is set up config.default['freqplot.rcParams'].
2258+
22362259
"""
22372260
# Keyword processing
22382261
dB = config._get_param(
@@ -2241,7 +2264,6 @@ def singular_values_plot(
22412264
'freqplot', 'Hz', kwargs, _freqplot_defaults, pop=True)
22422265
grid = config._get_param(
22432266
'freqplot', 'grid', kwargs, _freqplot_defaults, pop=True)
2244-
omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)
22452267
freqplot_rcParams = config._get_param(
22462268
'freqplot', 'rcParams', kwargs, _freqplot_defaults, pop=True)
22472269

@@ -2255,6 +2277,17 @@ def singular_values_plot(
22552277
data, omega=omega, omega_limits=omega_limits,
22562278
omega_num=omega_num)
22572279
else:
2280+
# Generate warnings if frequency keywords were given
2281+
if omega_num is not None:
2282+
warnings.warn("`omega_num` ignored when passed response data")
2283+
elif omega is not None:
2284+
warnings.warn("`omega` ignored when passed response data")
2285+
2286+
# Check to make sure omega_limits is sensible
2287+
if omega_limits is not None and \
2288+
(len(omega_limits) != 2 or omega_limits[1] <= omega_limits[0]):
2289+
raise ValueError(f"invalid limits: {omega_limits=}")
2290+
22582291
responses = data
22592292

22602293
# Process (legacy) plot keyword
@@ -2308,48 +2341,49 @@ def singular_values_plot(
23082341
# Create a list of lines for the output
23092342
out = np.empty(len(data), dtype=object)
23102343

2344+
# Plot the singular values for each response
23112345
for idx_sys, response in enumerate(responses):
23122346
sigma = sigmas[idx_sys].transpose() # frequency first for plotting
2313-
omega_sys = omegas[idx_sys]
2347+
omega = omegas[idx_sys] / (2 * math.pi) if Hz else omegas[idx_sys]
2348+
23142349
if response.isdtime(strict=True):
2315-
nyquistfrq = math.pi / response.dt
2350+
nyq_freq = (0.5/response.dt) if Hz else (math.pi/response.dt)
23162351
else:
2317-
nyquistfrq = None
2352+
nyq_freq = None
23182353

2319-
color = color_cycle[(idx_sys + color_offset) % len(color_cycle)]
2320-
color = kwargs.pop('color', color)
2321-
2322-
# TODO: copy from above
2323-
nyquistfrq_plot = None
2324-
if Hz:
2325-
omega_plot = omega_sys / (2. * math.pi)
2326-
if nyquistfrq:
2327-
nyquistfrq_plot = nyquistfrq / (2. * math.pi)
2354+
# See if the color was specified, otherwise rotate
2355+
if kwargs.get('color', None) or any(
2356+
[isinstance(arg, str) and
2357+
any([c in arg for c in "bgrcmykw#"]) for arg in fmt]):
2358+
color_arg = {} # color set by *fmt, **kwargs
23282359
else:
2329-
omega_plot = omega_sys
2330-
if nyquistfrq:
2331-
nyquistfrq_plot = nyquistfrq
2332-
sigma_plot = sigma
2360+
color_arg = {'color': color_cycle[
2361+
(idx_sys + color_offset) % len(color_cycle)]}
23332362

23342363
# Decide on the system name
23352364
sysname = response.sysname if response.sysname is not None \
23362365
else f"Unknown-{idx_sys}"
23372366

2367+
# Plot the data
23382368
if dB:
23392369
with plt.rc_context(freqplot_rcParams):
23402370
out[idx_sys] = ax_sigma.semilogx(
2341-
omega_plot, 20 * np.log10(sigma_plot), *fmt, color=color,
2342-
label=sysname, **kwargs)
2371+
omega, 20 * np.log10(sigma), *fmt,
2372+
label=sysname, **color_arg, **kwargs)
23432373
else:
23442374
with plt.rc_context(freqplot_rcParams):
23452375
out[idx_sys] = ax_sigma.loglog(
2346-
omega_plot, sigma_plot, color=color, label=sysname,
2347-
*fmt, **kwargs)
2376+
omega, sigma, label=sysname, *fmt, **color_arg, **kwargs)
23482377

2349-
if nyquistfrq_plot is not None:
2378+
# Plot the Nyquist frequency
2379+
if nyq_freq is not None:
23502380
ax_sigma.axvline(
2351-
nyquistfrq_plot, color=color, linestyle='--',
2352-
label='_nyq_freq_' + sysname)
2381+
nyq_freq, linestyle='--', label='_nyq_freq_' + sysname,
2382+
**color_arg)
2383+
2384+
# If specific omega_limits were given, use them
2385+
if omega_limits is not None:
2386+
ax_sigma.set_xlim(omega_limits)
23532387

23542388
# Add a grid to the plot + labeling
23552389
if grid:

‎control/sisotool.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -149,7 +149,7 @@ def _SisotoolUpdate(sys, fig, K, bode_plot_params, tvect=None):
149149

150150
# Update the bodeplot
151151
bode_plot_params['data'] = frequency_response(sys_loop*K.real)
152-
bode_plot(**bode_plot_params)
152+
bode_plot(**bode_plot_params, title=False)
153153

154154
# Set the titles and labels
155155
ax_mag.set_title('Bode magnitude',fontsize = title_font_size)

0 commit comments

Comments
 (0)