Skip to content

Commit 402b45f

Browse files
murrayrmmurrayrm
committed
restore [wmin, wmax] functionality + documentation updates
1 parent eeebec3 commit 402b45f

File tree

9 files changed

+282
-119
lines changed

9 files changed

+282
-119
lines changed

‎control/freqplot.py

Lines changed: 115 additions & 80 deletions
Large diffs are not rendered by default.

‎control/lti.py

Lines changed: 10 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -386,16 +386,18 @@ def frequency_response(
386386
sysdata : LTI system or list of LTI systems
387387
Linear system(s) for which frequency response is computed.
388388
omega : float or 1D array_like, optional
389-
A list of frequencies in radians/sec at which the system should be
390-
evaluated. The list can be either a Python list or a numpy array
391-
and will be sorted before evaluation. If None (default), a common
392-
set of frequencies that works across all given systems is computed.
389+
Frequencies in radians/sec at which the system should be
390+
evaluated. Can be a single frequency or array of frequencies, which
391+
will be sorted before evaluation. If None (default), a common set
392+
of frequencies that works across all given systems is computed.
393393
omega_limits : array_like of two values, optional
394-
Limits to the range of frequencies, in rad/sec. Ignored if
395-
omega is provided, and auto-generated if omitted.
394+
Limits to the range of frequencies, in rad/sec. Specifying
395+
``omega`` as a list of two elements is equivalent to providing
396+
``omega_limits``. Ignored if omega is provided.
396397
omega_num : int, optional
397-
Number of frequency samples to plot. Defaults to
398-
config.defaults['freqplot.number_of_samples'].
398+
Number of frequency samples at which to compute the response.
399+
Defaults to config.defaults['freqplot.number_of_samples']. Ignored
400+
if omega is provided.
399401
400402
Returns
401403
-------

‎control/tests/freqplot_test.py

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -181,6 +181,12 @@ def test_basic_freq_plots(savefigs=False):
181181
if savefigs:
182182
plt.savefig('freqplot-siso_bode-default.png')
183183

184+
plt.figure()
185+
omega = np.logspace(-2, 2, 500)
186+
ct.frequency_response([sys1, sys2], omega).plot(initial_phase=0)
187+
if savefigs:
188+
plt.savefig('freqplot-siso_bode-omega.png')
189+
184190
# Basic MIMO Bode plot
185191
plt.figure()
186192
sys_mimo = ct.tf(
@@ -213,6 +219,24 @@ def test_basic_freq_plots(savefigs=False):
213219
if savefigs:
214220
plt.savefig('freqplot-siso_nichols-default.png')
215221

222+
# Nyquist plot - default settings
223+
plt.figure()
224+
sys = ct.tf([1, 0.2], [1, 1, 3, 1, 1], name='sys')
225+
ct.nyquist(sys)
226+
if savefigs:
227+
plt.savefig('freqplot-nyquist-default.png')
228+
229+
# Nyquist plot - custom settings
230+
plt.figure()
231+
sys = ct.tf([1, 0.2], [1, 0, 1]) * ct.tf([1], [1, 0])
232+
nyqresp = ct.nyquist_response(sys)
233+
nyqresp.plot(
234+
max_curve_magnitude=6, max_curve_offset=1,
235+
arrows=[0, 0.15, 0.3, 0.6, 0.7, 0.925], label='sys')
236+
print("Encirclements =", nyqresp.count)
237+
if savefigs:
238+
plt.savefig('freqplot-nyquist-custom.png')
239+
216240

217241
def test_gangof4_plots(savefigs=False):
218242
proc = ct.tf([1], [1, 1, 1], name="process")

‎control/tests/freqresp_test.py

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -709,3 +709,25 @@ def test_singular_values_plot_mpl_superimpose_nyq(ss_mimo_ct, ss_mimo_dt):
709709
assert(len(nyquist_line[0]) == 2)
710710
assert(nyquist_line[0][0] == nyquist_line[0][1])
711711
assert(nyquist_line[0][0] == np.pi/sys_dt.dt)
712+
713+
714+
def test_freqresp_omega_limits():
715+
sys = ctrl.rss(4, 1, 1)
716+
717+
# Generate a standard frequency response (no limits specified)
718+
resp0 = ctrl.frequency_response(sys)
719+
720+
# Regenerate the response using omega_limits
721+
resp1 = ctrl.frequency_response(
722+
sys, omega_limits=[resp0.omega[0], resp0.omega[-1]])
723+
np.testing.assert_equal(resp0.omega, resp1.omega)
724+
725+
# Regenerate the response using omega as a list of two elements
726+
resp2 = ctrl.frequency_response(sys, [resp0.omega[0], resp0.omega[-1]])
727+
np.testing.assert_equal(resp0.omega, resp2.omega)
728+
assert resp2.omega.size > 100
729+
730+
# Make sure that generating response using array does the right thing
731+
resp3 = ctrl.frequency_response(
732+
sys, np.array([resp0.omega[0], resp0.omega[-1]]))
733+
np.testing.assert_equal(resp3.omega, [resp0.omega[0], resp0.omega[-1]])

‎control/tests/nyquist_test.py

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -438,6 +438,33 @@ def test_discrete_nyquist():
438438
ct.nyquist_response(sys)
439439

440440

441+
def test_freqresp_omega_limits():
442+
sys = ct.rss(4, 1, 1)
443+
444+
# Generate a standard frequency response (no limits specified)
445+
resp0 = ct.nyquist_response(sys)
446+
assert resp0.contour.size > 2
447+
448+
# Regenerate the response using omega_limits
449+
resp1 = ct.nyquist_response(
450+
sys, omega_limits=[resp0.contour[1].imag, resp0.contour[-1].imag])
451+
assert resp1.contour.size > 2
452+
assert np.isclose(resp1.contour[0], resp0.contour[1])
453+
assert np.isclose(resp1.contour[-1], resp0.contour[-1])
454+
455+
# Regenerate the response using omega as a list of two elements
456+
resp2 = ct.nyquist_response(
457+
sys, [resp0.contour[1].imag, resp0.contour[-1].imag])
458+
np.testing.assert_equal(resp1.contour, resp2.contour)
459+
460+
# Make sure that generating response using array does the right thing
461+
resp3 = ct.nyquist_response(
462+
sys, np.array([resp0.contour[1].imag, resp0.contour[-1].imag]))
463+
np.testing.assert_equal(
464+
resp3.contour,
465+
np.array([resp0.contour[1], resp0.contour[-1]]))
466+
467+
441468
if __name__ == "__main__":
442469
#
443470
# Interactive mode: generate plots for manual viewing

‎doc/freqplot-nyquist-custom.png

42.7 KB
Loading

‎doc/freqplot-nyquist-default.png

40.8 KB
Loading

‎doc/freqplot-siso_bode-omega.png

43.8 KB
Loading

‎doc/plotting.rst

Lines changed: 84 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -107,13 +107,13 @@ each other). The following plot shows the use of `plot_inputs='overlay'`
107107
as well as the ability to reposition the legends using the `legend_map`
108108
keyword::
109109

110-
timepts = np.linspace(0, 10, 100)
111-
U = np.vstack([np.sin(timepts), np.cos(2*timepts)])
112-
ct.input_output_response(sys_mimo, timepts, U).plot(
113-
plot_inputs='overlay',
114-
legend_map=np.array([['lower right'], ['lower right']]),
115-
title="I/O response for 2x2 MIMO system " +
116-
"[plot_inputs='overlay', legend_map]")
110+
timepts = np.linspace(0, 10, 100)
111+
U = np.vstack([np.sin(timepts), np.cos(2*timepts)])
112+
ct.input_output_response(sys_mimo, timepts, U).plot(
113+
plot_inputs='overlay',
114+
legend_map=np.array([['lower right'], ['lower right']]),
115+
title="I/O response for 2x2 MIMO system " +
116+
"[plot_inputs='overlay', legend_map]")
117117

118118
.. image:: timeplot-mimo_ioresp-ov_lm.png
119119

@@ -122,17 +122,17 @@ instead of plotting the outputs on the top and inputs on the bottom, the
122122
inputs are plotted on the left and outputs on the right, as shown in the
123123
following figure::
124124

125-
U1 = np.vstack([np.sin(timepts), np.cos(2*timepts)])
126-
resp1 = ct.input_output_response(sys_mimo, timepts, U1)
125+
U1 = np.vstack([np.sin(timepts), np.cos(2*timepts)])
126+
resp1 = ct.input_output_response(sys_mimo, timepts, U1)
127127

128-
U2 = np.vstack([np.cos(2*timepts), np.sin(timepts)])
129-
resp2 = ct.input_output_response(sys_mimo, timepts, U2)
128+
U2 = np.vstack([np.cos(2*timepts), np.sin(timepts)])
129+
resp2 = ct.input_output_response(sys_mimo, timepts, U2)
130130

131-
ct.combine_time_responses(
132-
[resp1, resp2], trace_labels=["Scenario #1", "Scenario #2"]).plot(
133-
transpose=True,
134-
title="I/O responses for 2x2 MIMO system, multiple traces "
135-
"[transpose]")
131+
ct.combine_time_responses(
132+
[resp1, resp2], trace_labels=["Scenario #1", "Scenario #2"]).plot(
133+
transpose=True,
134+
title="I/O responses for 2x2 MIMO system, multiple traces "
135+
"[transpose]")
136136

137137
.. image:: timeplot-mimo_ioresp-mt_tr.png
138138

@@ -146,11 +146,11 @@ Additional customization is possible using the `input_props`,
146146
`output_props`, and `trace_props` keywords to set complementary line colors
147147
and styles for various signals and traces::
148148

149-
out = ct.step_response(sys_mimo).plot(
150-
plot_inputs='overlay', overlay_signals=True, overlay_traces=True,
151-
output_props=[{'color': c} for c in ['blue', 'orange']],
152-
input_props=[{'color': c} for c in ['red', 'green']],
153-
trace_props=[{'linestyle': s} for s in ['-', '--']])
149+
out = ct.step_response(sys_mimo).plot(
150+
plot_inputs='overlay', overlay_signals=True, overlay_traces=True,
151+
output_props=[{'color': c} for c in ['blue', 'orange']],
152+
input_props=[{'color': c} for c in ['red', 'green']],
153+
trace_props=[{'linestyle': s} for s in ['-', '--']])
154154

155155
.. image:: timeplot-mimo_step-linestyle.png
156156

@@ -196,7 +196,7 @@ overlaying the inputs or outputs::
196196

197197
.. image:: freqplot-mimo_bode-magonly.png
198198

199-
The :func:`~ct.singular_values_response` function can be used to
199+
The :func:`~control.singular_values_response` function can be used to
200200
generate Bode plots that show the singular values of a transfer
201201
function::
202202

@@ -213,16 +213,69 @@ plot, use `plot_type='nichols'`::
213213
.. image:: freqplot-siso_nichols-default.png
214214

215215
Another response function that can be used to generate Bode plots is
216-
the :func:`~ct.gangof4` function, which computes the four primary
216+
the :func:`~control.gangof4` function, which computes the four primary
217217
sensitivity functions for a feedback control system in standard form::
218218

219-
proc = ct.tf([1], [1, 1, 1], name="process")
220-
ctrl = ct.tf([100], [1, 5], name="control")
221-
response = rect.gangof4_response(proc, ctrl)
222-
ct.bode_plot(response) # or response.plot()
219+
proc = ct.tf([1], [1, 1, 1], name="process")
220+
ctrl = ct.tf([100], [1, 5], name="control")
221+
response = rect.gangof4_response(proc, ctrl)
222+
ct.bode_plot(response) # or response.plot()
223223

224224
.. image:: freqplot-gangof4.png
225225

226+
Nyquist analysys can be done using the :func:`~control.nyquist_response`
227+
function, which evaluates an LTI system along the Nyquist contour, and
228+
the :func:`~control.nyquist_plot` function, which generates a Nyquist plot::
229+
230+
sys = ct.tf([1, 0.2], [1, 1, 3, 1, 1], name='sys')
231+
nyquist_plot(sys)
232+
233+
.. image:: freqplot-nyquist-default.png
234+
235+
The :func:`~control.nyquist_response` function can be used to compute
236+
the number of encirclement of the -1 point and can return the Nyquist
237+
contour that was used to generate the Nyquist curve.
238+
239+
By default, the Nyquist response will generate small semicircles around
240+
poles that are on the imaginary axis. In addition, portions of the Nyquist
241+
curve that far from the origin are scaled to a maximum value, with the line
242+
style is changed to reflect the scaling, and it is possible to offset the
243+
scaled portions to separate out the portions of the Nyquist curve at
244+
<math>\infty</math>. A number of keyword parameters for both are available
245+
for :func:`~control.nyquist_response`and :func:`~control.nyquist_plot` to
246+
tune the computation of the Nyquist curve and the way the data are
247+
plotted::
248+
249+
sys = ct.tf([1, 0.2], [1, 0, 1]) * ct.tf([1], [1, 0])
250+
nyqresp = ct.nyquist_response(sys)
251+
nyqresp.plot(
252+
max_curve_magnitude=6, max_curve_offset=1,
253+
arrows=[0, 0.15, 0.3, 0.6, 0.7, 0.925], label='sys')
254+
print("Encirclements =", nyqresp.count)
255+
256+
.. image:: freqplot-nyquist-custom.png
257+
258+
All frequency domain plotting functions will automatically compute the
259+
range of frequencies to plot based on the poles and zeros of the frequency
260+
response. Frequency points can be explicitly specified by including an
261+
array of frequencies as a second argument (after the list of systems)::
262+
263+
sys1 = ct.tf([1], [1, 2, 1], name='sys1')
264+
sys2 = ct.tf([1, 0.2], [1, 1, 3, 1, 1], name='sys2')
265+
omega = np.logspace(-2, 2, 500)
266+
ct.frequency_response([sys1, sys2], omega).plot(initial_phase=0)
267+
268+
.. image:: freqplot-siso_bode-omega.png
269+
270+
Alternatively. frequency ranges can be specified by passing a list of the
271+
form ``[wmin, wmax]``, where ``wmin`` and ``wmax`` are the minimum and
272+
maximum frequencies in the (log-spaced) frequency range::
273+
274+
response = ct.frequency_response([sys1, sys2], [1e-2, 1e2])
275+
276+
The number of (log-spaced) points in the frequency can be specified using
277+
the ``omega_num`` keyword parameter.
278+
226279

227280
Pole/zero data
228281
==============
@@ -288,7 +341,7 @@ The default method for generating a phase plane plot is to provide a
288341
2D dynamical system along with a range of coordinates and time limit::
289342

290343
sys = ct.nlsys(
291-
lambda t, x, u, params: np.array([[0, 1], [-1, -1]]) @ x,
344+
lambda t, x, u, params: np.array([[0, 1], [-1, -1]]) @ x,
292345
states=['position', 'velocity'], inputs=0, name='damped oscillator')
293346
axis_limits = [-1, 1, -1, 1]
294347
T = 8
@@ -310,15 +363,15 @@ an inverted pendulum system, which is created using a mesh grid::
310363
m, l, b, g = params['m'], params['l'], params['b'], params['g']
311364
return [x[1], -b/m * x[1] + (g * l / m) * np.sin(x[0]) + u[0]/m]
312365
invpend = ct.nlsys(invpend_update, states=2, inputs=1, name='invpend')
313-
366+
314367
ct.phase_plane_plot(
315368
invpend, [-2*pi, 2*pi, -2, 2], 5,
316369
gridtype='meshgrid', gridspec=[5, 8], arrows=3,
317370
plot_equilpoints={'gridspec': [12, 9]},
318371
params={'m': 1, 'l': 1, 'b': 0.2, 'g': 1})
319372
plt.xlabel(r"$\theta$ [rad]")
320373
plt.ylabel(r"$\dot\theta$ [rad/sec]")
321-
374+
322375
.. image:: phaseplot-invpend-meshgrid.png
323376

324377
This figure shows several features of more complex phase plane plots:
@@ -341,7 +394,7 @@ are part of the :mod:`~control.phaseplot` (pp) module::
341394
-x[0] + x[1] * (1 - x[0]**2 - x[1]**2)]
342395
oscillator = ct.nlsys(
343396
oscillator_update, states=2, inputs=0, name='nonlinear oscillator')
344-
397+
345398
ct.phase_plane_plot(oscillator, [-1.5, 1.5, -1.5, 1.5], 0.9)
346399
pp.streamlines(
347400
oscillator, np.array([[0, 0]]), 1.5,

0 commit comments

Comments
 (0)