Skip to content

Commit 5125ff7

Browse files
murrayrmmurrayrm
authored
implementation of initial_phase, wrap_phase keywords for bode() (#494)
* implementation of initial_phase, wrap_phase keywords for bode() * add additional documentation on initial_phase if deg=False * clear figure in unit tests to avoid slow plotting problem
1 parent 383e7a4 commit 5125ff7

File tree

2 files changed

+139
-13
lines changed

2 files changed

+139
-13
lines changed

‎control/freqplot.py‎

Lines changed: 73 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,7 @@
7878
'bode.deg': True, # Plot phase in degrees
7979
'bode.Hz': False, # Plot frequency in Hertz
8080
'bode.grid': True, # Turn on grid for gain and phase
81+
'bode.wrap_phase': False, # Wrap the phase plot at a given value
8182
}
8283

8384

@@ -131,7 +132,19 @@ def bode_plot(syslist, omega=None,
131132
grid : bool
132133
If True, plot grid lines on gain and phase plots. Default is set by
133134
`config.defaults['bode.grid']`.
134-
135+
initial_phase : float
136+
Set the reference phase to use for the lowest frequency. If set, the
137+
initial phase of the Bode plot will be set to the value closest to the
138+
value specified. Units are in either degrees or radians, depending on
139+
the `deg` parameter. Default is -180 if wrap_phase is False, 0 if
140+
wrap_phase is True.
141+
wrap_phase : bool or float
142+
If wrap_phase is `False`, then the phase will be unwrapped so that it
143+
is continuously increasing or decreasing. If wrap_phase is `True` the
144+
phase will be restricted to the range [-180, 180) (or [:math:`-\\pi`,
145+
:math:`\\pi`) radians). If `wrap_phase` is specified as a float, the
146+
phase will be offset by 360 degrees if it falls below the specified
147+
value. Default to `False`, set by config.defaults['bode.wrap_phase'].
135148
136149
The default values for Bode plot configuration parameters can be reset
137150
using the `config.defaults` dictionary, with module name 'bode'.
@@ -171,6 +184,10 @@ def bode_plot(syslist, omega=None,
171184
grid = config._get_param('bode', 'grid', kwargs, _bode_defaults, pop=True)
172185
plot = config._get_param('bode', 'grid', plot, True)
173186
margins = config._get_param('bode', 'margins', margins, False)
187+
wrap_phase = config._get_param(
188+
'bode', 'wrap_phase', kwargs, _bode_defaults, pop=True)
189+
initial_phase = config._get_param(
190+
'bode', 'initial_phase', kwargs, None, pop=True)
174191

175192
# If argument was a singleton, turn it into a list
176193
if not getattr(syslist, '__iter__', False):
@@ -209,11 +226,47 @@ def bode_plot(syslist, omega=None,
209226
# TODO: What distance to the Nyquist frequency is appropriate?
210227
else:
211228
nyquistfrq = None
229+
212230
# Get the magnitude and phase of the system
213231
mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
214232
mag = np.atleast_1d(np.squeeze(mag_tmp))
215233
phase = np.atleast_1d(np.squeeze(phase_tmp))
216-
phase = unwrap(phase)
234+
235+
#
236+
# Post-process the phase to handle initial value and wrapping
237+
#
238+
239+
if initial_phase is None:
240+
# Start phase in the range 0 to -360 w/ initial phase = -180
241+
# If wrap_phase is true, use 0 instead (phase \in (-pi, pi])
242+
initial_phase = -math.pi if wrap_phase is not True else 0
243+
elif isinstance(initial_phase, (int, float)):
244+
# Allow the user to override the default calculation
245+
if deg:
246+
initial_phase = initial_phase/180. * math.pi
247+
else:
248+
raise ValueError("initial_phase must be a number.")
249+
250+
# Shift the phase if needed
251+
if abs(phase[0] - initial_phase) > math.pi:
252+
phase -= 2*math.pi * \
253+
round((phase[0] - initial_phase) / (2*math.pi))
254+
255+
# Phase wrapping
256+
if wrap_phase is False:
257+
phase = unwrap(phase) # unwrap the phase
258+
elif wrap_phase is True:
259+
pass # default calculation OK
260+
elif isinstance(wrap_phase, (int, float)):
261+
phase = unwrap(phase) # unwrap the phase first
262+
if deg:
263+
wrap_phase *= math.pi/180.
264+
265+
# Shift the phase if it is below the wrap_phase
266+
phase += 2*math.pi * np.maximum(
267+
0, np.ceil((wrap_phase - phase)/(2*math.pi)))
268+
else:
269+
raise ValueError("wrap_phase must be bool or float.")
217270

218271
mags.append(mag)
219272
phases.append(phase)
@@ -270,7 +323,9 @@ def bode_plot(syslist, omega=None,
270323
label='control-bode-phase',
271324
sharex=ax_mag)
272325

326+
#
273327
# Magnitude plot
328+
#
274329
if dB:
275330
pltline = ax_mag.semilogx(omega_plot, 20 * np.log10(mag),
276331
*args, **kwargs)
@@ -285,19 +340,22 @@ def bode_plot(syslist, omega=None,
285340
ax_mag.grid(grid and not margins, which='both')
286341
ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")
287342

343+
#
288344
# Phase plot
289-
if deg:
290-
phase_plot = phase * 180. / math.pi
291-
else:
292-
phase_plot = phase
345+
#
346+
phase_plot = phase * 180. / math.pi if deg else phase
347+
348+
# Plot the data
293349
ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs)
294350

295351
# Show the phase and gain margins in the plot
296352
if margins:
353+
# Compute stability margins for the system
297354
margin = stability_margins(sys)
298-
gm, pm, Wcg, Wcp = \
299-
margin[0], margin[1], margin[3], margin[4]
300-
# TODO: add some documentation describing why this is here
355+
gm, pm, Wcg, Wcp = (margin[i] for i in (0, 1, 3, 4))
356+
357+
# Figure out sign of the phase at the first gain crossing
358+
# (needed if phase_wrap is True)
301359
phase_at_cp = phases[0][(np.abs(omegas[0] - Wcp)).argmin()]
302360
if phase_at_cp >= 0.:
303361
phase_limit = 180.
@@ -307,6 +365,7 @@ def bode_plot(syslist, omega=None,
307365
if Hz:
308366
Wcg, Wcp = Wcg/(2*math.pi), Wcp/(2*math.pi)
309367

368+
# Draw lines at gain and phase limits
310369
ax_mag.axhline(y=0 if dB else 1, color='k', linestyle=':',
311370
zorder=-20)
312371
ax_phase.axhline(y=phase_limit if deg else
@@ -315,6 +374,7 @@ def bode_plot(syslist, omega=None,
315374
mag_ylim = ax_mag.get_ylim()
316375
phase_ylim = ax_phase.get_ylim()
317376

377+
# Annotate the phase margin (if it exists)
318378
if pm != float('inf') and Wcp != float('nan'):
319379
if dB:
320380
ax_mag.semilogx(
@@ -327,7 +387,7 @@ def bode_plot(syslist, omega=None,
327387

328388
if deg:
329389
ax_phase.semilogx(
330-
[Wcp, Wcp], [1e5, phase_limit+pm],
390+
[Wcp, Wcp], [1e5, phase_limit + pm],
331391
color='k', linestyle=':', zorder=-20)
332392
ax_phase.semilogx(
333393
[Wcp, Wcp], [phase_limit + pm, phase_limit],
@@ -343,6 +403,7 @@ def bode_plot(syslist, omega=None,
343403
math.radians(phase_limit)],
344404
color='k', zorder=-20)
345405

406+
# Annotate the gain margin (if it exists)
346407
if gm != float('inf') and Wcg != float('nan'):
347408
if dB:
348409
ax_mag.semilogx(
@@ -360,11 +421,11 @@ def bode_plot(syslist, omega=None,
360421

361422
if deg:
362423
ax_phase.semilogx(
363-
[Wcg, Wcg], [1e-8, phase_limit],
424+
[Wcg, Wcg], [0, phase_limit],
364425
color='k', linestyle=':', zorder=-20)
365426
else:
366427
ax_phase.semilogx(
367-
[Wcg, Wcg], [1e-8, math.radians(phase_limit)],
428+
[Wcg, Wcg], [0, math.radians(phase_limit)],
368429
color='k', linestyle=':', zorder=-20)
369430

370431
ax_mag.set_ylim(mag_ylim)

‎control/tests/freqresp_test.py‎

Lines changed: 66 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
import matplotlib.pyplot as plt
1010
import numpy as np
1111
from numpy.testing import assert_allclose
12+
import math
1213
import pytest
1314

1415
import control as ctrl
@@ -173,7 +174,7 @@ def test_bode_margin(dB, maginfty1, maginfty2, gminv,
173174
rtol=1e-5)
174175

175176
phase_to_infinity = (np.array([Wcg, Wcg]),
176-
np.array([1e-8, p0]))
177+
np.array([0, p0]))
177178
assert_allclose(phase_to_infinity, allaxes[1].lines[4].get_data(),
178179
rtol=1e-5)
179180

@@ -271,3 +272,67 @@ def test_options(editsdefaults):
271272
# Make sure we got the right number of points
272273
assert numpoints1 != numpoints3
273274
assert numpoints3 == 13
275+
276+
@pytest.mark.parametrize(
277+
"TF, initial_phase, default_phase, expected_phase",
278+
[pytest.param(ctrl.tf([1], [1, 0]),
279+
None, -math.pi/2, -math.pi/2, id="order1, default"),
280+
pytest.param(ctrl.tf([1], [1, 0]),
281+
180, -math.pi/2, 3*math.pi/2, id="order1, 180"),
282+
pytest.param(ctrl.tf([1], [1, 0, 0]),
283+
None, -math.pi, -math.pi, id="order2, default"),
284+
pytest.param(ctrl.tf([1], [1, 0, 0]),
285+
180, -math.pi, math.pi, id="order2, 180"),
286+
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
287+
None, -3*math.pi/2, -3*math.pi/2, id="order2, default"),
288+
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
289+
180, -3*math.pi/2, math.pi/2, id="order2, 180"),
290+
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0]),
291+
None, 0, 0, id="order4, default"),
292+
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0]),
293+
180, 0, 0, id="order4, 180"),
294+
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0]),
295+
-360, 0, -2*math.pi, id="order4, -360"),
296+
])
297+
def test_initial_phase(TF, initial_phase, default_phase, expected_phase):
298+
# Check initial phase of standard transfer functions
299+
mag, phase, omega = ctrl.bode(TF)
300+
assert(abs(phase[0] - default_phase) < 0.1)
301+
302+
# Now reset the initial phase to +180 and see if things work
303+
mag, phase, omega = ctrl.bode(TF, initial_phase=initial_phase)
304+
assert(abs(phase[0] - expected_phase) < 0.1)
305+
306+
# Make sure everything works in rad/sec as well
307+
if initial_phase:
308+
plt.clf() # clear previous figure (speeds things up)
309+
mag, phase, omega = ctrl.bode(
310+
TF, initial_phase=initial_phase/180. * math.pi, deg=False)
311+
assert(abs(phase[0] - expected_phase) < 0.1)
312+
313+
314+
@pytest.mark.parametrize(
315+
"TF, wrap_phase, min_phase, max_phase",
316+
[pytest.param(ctrl.tf([1], [1, 0]),
317+
None, -math.pi/2, 0, id="order1, default"),
318+
pytest.param(ctrl.tf([1], [1, 0]),
319+
True, -math.pi, math.pi, id="order1, True"),
320+
pytest.param(ctrl.tf([1], [1, 0]),
321+
-270, -3*math.pi/2, math.pi/2, id="order1, -270"),
322+
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
323+
None, -3*math.pi/2, 0, id="order3, default"),
324+
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
325+
True, -math.pi, math.pi, id="order3, True"),
326+
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
327+
-270, -3*math.pi/2, math.pi/2, id="order3, -270"),
328+
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
329+
True, -3*math.pi/2, 0, id="order5, default"),
330+
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
331+
True, -math.pi, math.pi, id="order5, True"),
332+
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
333+
-270, -3*math.pi/2, math.pi/2, id="order5, -270"),
334+
])
335+
def test_phase_wrap(TF, wrap_phase, min_phase, max_phase):
336+
mag, phase, omega = ctrl.bode(TF, wrap_phase=wrap_phase)
337+
assert(min(phase) >= min_phase)
338+
assert(max(phase) <= max_phase)

0 commit comments

Comments
 (0)