Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

implementation of initial_phase, wrap_phase keywords for bode()#494

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 ourterms of service andprivacy statement. We’ll occasionally send you account related emails.

Already on GitHub?Sign in to your account

Merged
Merged
Show file tree
Hide file tree
Changes fromall 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
85 changes: 73 additions & 12 deletionscontrol/freqplot.py
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -78,6 +78,7 @@
'bode.deg': True, # Plot phase in degrees
'bode.Hz': False, # Plot frequency in Hertz
'bode.grid': True, # Turn on grid for gain and phase
'bode.wrap_phase': False, # Wrap the phase plot at a given value
}


Expand DownExpand Up@@ -131,7 +132,19 @@ def bode_plot(syslist, omega=None,
grid : bool
If True, plot grid lines on gain and phase plots. Default is set by
`config.defaults['bode.grid']`.

initial_phase : float
Set the reference phase to use for the lowest frequency. If set, the
initial phase of the Bode plot will be set to the value closest to the
value specified. Units are in either degrees or radians, depending on
the `deg` parameter. Default is -180 if wrap_phase is False, 0 if
wrap_phase is True.
wrap_phase : bool or float
If wrap_phase is `False`, then the phase will be unwrapped so that it
is continuously increasing or decreasing. If wrap_phase is `True` the
phase will be restricted to the range [-180, 180) (or [:math:`-\\pi`,
:math:`\\pi`) radians). If `wrap_phase` is specified as a float, the
phase will be offset by 360 degrees if it falls below the specified
value. Default to `False`, set by config.defaults['bode.wrap_phase'].

The default values for Bode plot configuration parameters can be reset
using the `config.defaults` dictionary, with module name 'bode'.
Expand DownExpand Up@@ -171,6 +184,10 @@ def bode_plot(syslist, omega=None,
grid = config._get_param('bode', 'grid', kwargs, _bode_defaults, pop=True)
plot = config._get_param('bode', 'grid', plot, True)
margins = config._get_param('bode', 'margins', margins, False)
wrap_phase = config._get_param(
'bode', 'wrap_phase', kwargs, _bode_defaults, pop=True)
initial_phase = config._get_param(
'bode', 'initial_phase', kwargs, None, pop=True)

# If argument was a singleton, turn it into a list
if not getattr(syslist, '__iter__', False):
Expand DownExpand Up@@ -209,11 +226,47 @@ def bode_plot(syslist, omega=None,
# TODO: What distance to the Nyquist frequency is appropriate?
else:
nyquistfrq = None

# Get the magnitude and phase of the system
mag_tmp, phase_tmp, omega_sys = sys.freqresp(omega_sys)
mag = np.atleast_1d(np.squeeze(mag_tmp))
phase = np.atleast_1d(np.squeeze(phase_tmp))
phase = unwrap(phase)

#
# Post-process the phase to handle initial value and wrapping
#

if initial_phase is None:
# Start phase in the range 0 to -360 w/ initial phase = -180
# If wrap_phase is true, use 0 instead (phase \in (-pi, pi])
initial_phase = -math.pi if wrap_phase is not True else 0
elif isinstance(initial_phase, (int, float)):
# Allow the user to override the default calculation
if deg:
initial_phase = initial_phase/180. * math.pi
else:
raise ValueError("initial_phase must be a number.")

# Shift the phase if needed
if abs(phase[0] - initial_phase) > math.pi:
phase -= 2*math.pi * \
round((phase[0] - initial_phase) / (2*math.pi))

# Phase wrapping
if wrap_phase is False:
phase = unwrap(phase) # unwrap the phase
elif wrap_phase is True:
pass # default calculation OK
elif isinstance(wrap_phase, (int, float)):
phase = unwrap(phase) # unwrap the phase first
if deg:
wrap_phase *= math.pi/180.

# Shift the phase if it is below the wrap_phase
phase += 2*math.pi * np.maximum(
0, np.ceil((wrap_phase - phase)/(2*math.pi)))
else:
raise ValueError("wrap_phase must be bool or float.")

mags.append(mag)
phases.append(phase)
Expand DownExpand Up@@ -270,7 +323,9 @@ def bode_plot(syslist, omega=None,
label='control-bode-phase',
sharex=ax_mag)

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

#
# Phase plot
if deg:
phase_plot = phase * 180. / math.pi
else:
phase_plot = phase
#
phase_plot = phase * 180. / math.pi if deg else phase

# Plot the data
ax_phase.semilogx(omega_plot, phase_plot, *args, **kwargs)

# Show the phase and gain margins in the plot
if margins:
# Compute stability margins for the system
margin = stability_margins(sys)
gm, pm, Wcg, Wcp = \
margin[0], margin[1], margin[3], margin[4]
# TODO: add some documentation describing why this is here
gm, pm, Wcg, Wcp = (margin[i] for i in (0, 1, 3, 4))

# Figure out sign of the phase at the first gain crossing
# (needed if phase_wrap is True)
phase_at_cp = phases[0][(np.abs(omegas[0] - Wcp)).argmin()]
if phase_at_cp >= 0.:
phase_limit = 180.
Expand All@@ -307,6 +365,7 @@ def bode_plot(syslist, omega=None,
if Hz:
Wcg, Wcp = Wcg/(2*math.pi), Wcp/(2*math.pi)

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

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

if deg:
ax_phase.semilogx(
[Wcp, Wcp], [1e5, phase_limit+pm],
[Wcp, Wcp], [1e5, phase_limit +pm],
color='k', linestyle=':', zorder=-20)
ax_phase.semilogx(
[Wcp, Wcp], [phase_limit + pm, phase_limit],
Expand All@@ -343,6 +403,7 @@ def bode_plot(syslist, omega=None,
math.radians(phase_limit)],
color='k', zorder=-20)

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

if deg:
ax_phase.semilogx(
[Wcg, Wcg], [1e-8, phase_limit],
[Wcg, Wcg], [0, phase_limit],
color='k', linestyle=':', zorder=-20)
else:
ax_phase.semilogx(
[Wcg, Wcg], [1e-8, math.radians(phase_limit)],
[Wcg, Wcg], [0, math.radians(phase_limit)],
color='k', linestyle=':', zorder=-20)

ax_mag.set_ylim(mag_ylim)
Expand Down
67 changes: 66 additions & 1 deletioncontrol/tests/freqresp_test.py
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -9,6 +9,7 @@
import matplotlib.pyplot as plt
import numpy as np
from numpy.testing import assert_allclose
import math
import pytest

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

phase_to_infinity = (np.array([Wcg, Wcg]),
np.array([1e-8, p0]))
np.array([0, p0]))
assert_allclose(phase_to_infinity, allaxes[1].lines[4].get_data(),
rtol=1e-5)

Expand DownExpand Up@@ -271,3 +272,67 @@ def test_options(editsdefaults):
# Make sure we got the right number of points
assert numpoints1 != numpoints3
assert numpoints3 == 13

@pytest.mark.parametrize(
"TF, initial_phase, default_phase, expected_phase",
[pytest.param(ctrl.tf([1], [1, 0]),
None, -math.pi/2, -math.pi/2, id="order1, default"),
pytest.param(ctrl.tf([1], [1, 0]),
180, -math.pi/2, 3*math.pi/2, id="order1, 180"),
pytest.param(ctrl.tf([1], [1, 0, 0]),
None, -math.pi, -math.pi, id="order2, default"),
pytest.param(ctrl.tf([1], [1, 0, 0]),
180, -math.pi, math.pi, id="order2, 180"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
None, -3*math.pi/2, -3*math.pi/2, id="order2, default"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
180, -3*math.pi/2, math.pi/2, id="order2, 180"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0]),
None, 0, 0, id="order4, default"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0]),
180, 0, 0, id="order4, 180"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0]),
-360, 0, -2*math.pi, id="order4, -360"),
])
def test_initial_phase(TF, initial_phase, default_phase, expected_phase):
# Check initial phase of standard transfer functions
mag, phase, omega = ctrl.bode(TF)
assert(abs(phase[0] - default_phase) < 0.1)

# Now reset the initial phase to +180 and see if things work
mag, phase, omega = ctrl.bode(TF, initial_phase=initial_phase)
assert(abs(phase[0] - expected_phase) < 0.1)

# Make sure everything works in rad/sec as well
if initial_phase:
plt.clf() # clear previous figure (speeds things up)
mag, phase, omega = ctrl.bode(
TF, initial_phase=initial_phase/180. * math.pi, deg=False)
assert(abs(phase[0] - expected_phase) < 0.1)


@pytest.mark.parametrize(
"TF, wrap_phase, min_phase, max_phase",
[pytest.param(ctrl.tf([1], [1, 0]),
None, -math.pi/2, 0, id="order1, default"),
pytest.param(ctrl.tf([1], [1, 0]),
True, -math.pi, math.pi, id="order1, True"),
pytest.param(ctrl.tf([1], [1, 0]),
-270, -3*math.pi/2, math.pi/2, id="order1, -270"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
None, -3*math.pi/2, 0, id="order3, default"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
True, -math.pi, math.pi, id="order3, True"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0]),
-270, -3*math.pi/2, math.pi/2, id="order3, -270"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
True, -3*math.pi/2, 0, id="order5, default"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
True, -math.pi, math.pi, id="order5, True"),
pytest.param(ctrl.tf([1], [1, 0, 0, 0, 0, 0]),
-270, -3*math.pi/2, math.pi/2, id="order5, -270"),
])
def test_phase_wrap(TF, wrap_phase, min_phase, max_phase):
mag, phase, omega = ctrl.bode(TF, wrap_phase=wrap_phase)
assert(min(phase) >= min_phase)
assert(max(phase) <= max_phase)
[8]ページ先頭
©2009-2025 Movatter.jp