'freqplot.feature_periphery_decades': 1,

'freqplot.number_of_samples': 1000,

'freqplot.dB': False, # Plot gain in dB

'freqplot.deg': True, # Plot phase in degrees

'freqplot.Hz': False, # Plot frequency in Hertz

'freqplot.grid': True, # Turn on grid for gain and phase

'freqplot.wrap_phase': False, # Wrap the phase plot at a given value

# deprecations

'deprecated.bode.dB': 'freqplot.dB',

'deprecated.bode.deg': 'freqplot.deg',

'deprecated.bode.Hz': 'freqplot.Hz',

'deprecated.bode.grid': 'freqplot.grid',

'deprecated.bode.wrap_phase': 'freqplot.wrap_phase',

plot=True, omega_limits=None, omega_num=None,

margins=None, method='best', *args, **kwargs):

"""Bode plot for a system

# Make a copy of the kwargs dictionary since we will modify it

kwargs = dict(kwargs)

# Check to see if legacy 'Plot' keyword was used

if 'Plot' in kwargs:

import warnings

warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",

FutureWarning)

# Map 'Plot' keyword to 'plot' keyword

plot = kwargs.pop('Plot')

# Get values for params (and pop from list to allow keyword use in plot)

dB = config._get_param(

'freqplot', 'dB', kwargs, _freqplot_defaults, pop=True)

deg = config._get_param(

'freqplot', 'deg', kwargs, _freqplot_defaults, pop=True)

Hz = config._get_param(

'freqplot', 'Hz', kwargs, _freqplot_defaults, pop=True)

grid = config._get_param(

'freqplot', 'grid', kwargs, _freqplot_defaults, pop=True)

plot = config._get_param('freqplot', 'plot', plot, True)

margins = config._get_param(

'freqplot', 'margins', margins, False)

wrap_phase = config._get_param(

'freqplot', 'wrap_phase', kwargs, _freqplot_defaults, pop=True)

initial_phase = config._get_param(

'freqplot', 'initial_phase', kwargs, None, pop=True)

omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)

# If argument was a singleton, turn it into a tuple

if not isinstance(syslist, (list, tuple)):

syslist = (syslist,)

omega, omega_range_given = _determine_omega_vector(

syslist, omega, omega_limits, omega_num, Hz=Hz)

if plot:

# Set up the axes with labels so that multiple calls to

# bode_plot will superimpose the data. This was implicit

# before matplotlib 2.1, but changed after that (See

# https://github.com/matplotlib/matplotlib/issues/9024).

# The code below should work on all cases.

# Get the current figure

if 'sisotool' in kwargs:

fig = kwargs.pop('fig')

ax_mag = fig.axes[0]

ax_phase = fig.axes[2]

sisotool = kwargs.pop('sisotool')

else:

fig = plt.gcf()

ax_mag = None

ax_phase = None

sisotool = False

# Get the current axes if they already exist

for ax in fig.axes:

if ax.get_label() == 'control-bode-magnitude':

ax_mag = ax

elif ax.get_label() == 'control-bode-phase':

ax_phase = ax

# If no axes present, create them from scratch

if ax_mag is None or ax_phase is None:

plt.clf()

ax_mag = plt.subplot(211, label='control-bode-magnitude')

ax_phase = plt.subplot(

212, label='control-bode-phase', sharex=ax_mag)

mags, phases, omegas, nyquistfrqs = [], [], [], []

for sys in syslist:

if not sys.issiso():

# TODO: Add MIMO bode plots.

raise ControlMIMONotImplemented(

"Bode is currently only implemented for SISO systems.")

else:

omega_sys = np.asarray(omega)

if sys.isdtime(strict=True):

nyquistfrq = math.pi / sys.dt

if not omega_range_given:

# limit up to and including nyquist frequency

omega_sys = np.hstack((

omega_sys[omega_sys < nyquistfrq], nyquistfrq))

else:

nyquistfrq = None

mag, phase, omega_sys = sys.frequency_response(omega_sys)

mag = np.atleast_1d(mag)

phase = np.atleast_1d(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_value = -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_value = initial_phase/180. * math.pi

else:

initial_phase_value = initial_phase

else:

raise ValueError("initial_phase must be a number.")

# Shift the phase if needed

if abs(phase[0] - initial_phase_value) > math.pi:

phase -= 2*math.pi * \

round((phase[0] - initial_phase_value) / (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)

omegas.append(omega_sys)

nyquistfrqs.append(nyquistfrq)

# Get the dimensions of the current axis, which we will divide up

# TODO: Not current implemented; just use subplot for now

if plot:

nyquistfrq_plot = None

if Hz:

omega_plot = omega_sys / (2. * math.pi)

if nyquistfrq:

nyquistfrq_plot = nyquistfrq / (2. * math.pi)

else:

omega_plot = omega_sys

if nyquistfrq:

nyquistfrq_plot = nyquistfrq

phase_plot = phase * 180. / math.pi if deg else phase

mag_plot = mag

if nyquistfrq_plot:

# append data for vertical nyquist freq indicator line.

# if this extra nyquist lime is is plotted in a single plot

# command then line order is preserved when

# creating a legend eg. legend(('sys1', 'sys2'))

omega_nyq_line = np.array(

(np.nan, nyquistfrq_plot, nyquistfrq_plot))

omega_plot = np.hstack((omega_plot, omega_nyq_line))

mag_nyq_line = np.array((

np.nan, 0.7*min(mag_plot), 1.3*max(mag_plot)))

mag_plot = np.hstack((mag_plot, mag_nyq_line))

phase_range = max(phase_plot) - min(phase_plot)

phase_nyq_line = np.array(

(np.nan,

min(phase_plot) - 0.2 * phase_range,

max(phase_plot) + 0.2 * phase_range))

phase_plot = np.hstack((phase_plot, phase_nyq_line))

#

# Magnitude plot

#

if dB:

ax_mag.semilogx(omega_plot, 20 * np.log10(mag_plot),

*args, **kwargs)

else:

ax_mag.loglog(omega_plot, mag_plot, *args, **kwargs)

# Add a grid to the plot + labeling

ax_mag.grid(grid and not margins, which='both')

ax_mag.set_ylabel("Magnitude (dB)" if dB else "Magnitude")

#

# Phase plot

#

# 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, method=method)

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.

else:

phase_limit = -180.

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

math.radians(phase_limit),

color='k', linestyle=':', zorder=-20)

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(

[Wcp, Wcp], [0., -1e5],

color='k', linestyle=':', zorder=-20)

else:

ax_mag.loglog(

[Wcp, Wcp], [1., 1e-8],

color='k', linestyle=':', zorder=-20)

if deg:

ax_phase.semilogx(

[Wcp, Wcp], [1e5, phase_limit + pm],

color='k', linestyle=':', zorder=-20)

ax_phase.semilogx(

[Wcp, Wcp], [phase_limit + pm, phase_limit],

color='k', zorder=-20)

else:

ax_phase.semilogx(

[Wcp, Wcp], [1e5, math.radians(phase_limit) +

math.radians(pm)],

color='k', linestyle=':', zorder=-20)

ax_phase.semilogx(

[Wcp, Wcp], [math.radians(phase_limit) +

math.radians(pm),

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(

[Wcg, Wcg], [-20.*np.log10(gm), -1e5],

color='k', linestyle=':', zorder=-20)

ax_mag.semilogx(

[Wcg, Wcg], [0, -20*np.log10(gm)],

color='k', zorder=-20)

else:

ax_mag.loglog(

[Wcg, Wcg], [1./gm, 1e-8], color='k',

linestyle=':', zorder=-20)

ax_mag.loglog(

[Wcg, Wcg], [1., 1./gm], color='k', zorder=-20)

if deg:

ax_phase.semilogx(

[Wcg, Wcg], [0, phase_limit],

color='k', linestyle=':', zorder=-20)

else:

ax_phase.semilogx(

[Wcg, Wcg], [0, math.radians(phase_limit)],

color='k', linestyle=':', zorder=-20)

ax_mag.set_ylim(mag_ylim)

ax_phase.set_ylim(phase_ylim)

if sisotool:

ax_mag.text(

0.04, 0.06,

'G.M.: %.2f %s\nFreq: %.2f %s' %

(20*np.log10(gm) if dB else gm,

'dB ' if dB else '',

Wcg, 'Hz' if Hz else 'rad/s'),

horizontalalignment='left',

verticalalignment='bottom',

transform=ax_mag.transAxes,

fontsize=8 if int(mpl.__version__[0]) == 1 else 6)

ax_phase.text(

0.04, 0.06,

'P.M.: %.2f %s\nFreq: %.2f %s' %

(pm if deg else math.radians(pm),

'deg' if deg else 'rad',

Wcp, 'Hz' if Hz else 'rad/s'),

horizontalalignment='left',

verticalalignment='bottom',

transform=ax_phase.transAxes,

fontsize=8 if int(mpl.__version__[0]) == 1 else 6)

else:

plt.suptitle(

"Gm = %.2f %s(at %.2f %s), "

"Pm = %.2f %s (at %.2f %s)" %

(20*np.log10(gm) if dB else gm,

'dB ' if dB else '',

Wcg, 'Hz' if Hz else 'rad/s',

pm if deg else math.radians(pm),

'deg' if deg else 'rad',

Wcp, 'Hz' if Hz else 'rad/s'))

# Add a grid to the plot + labeling

ax_phase.set_ylabel("Phase (deg)" if deg else "Phase (rad)")

def gen_zero_centered_series(val_min, val_max, period):

v1 = np.ceil(val_min / period - 0.2)

v2 = np.floor(val_max / period + 0.2)

return np.arange(v1, v2 + 1) * period

if deg:

ylim = ax_phase.get_ylim()

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], 45.))

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], 15.), minor=True)

else:

ylim = ax_phase.get_ylim()

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], math.pi / 4.))

ax_phase.set_yticks(gen_zero_centered_series(

ylim[0], ylim[1], math.pi / 12.), minor=True)

ax_phase.grid(grid and not margins, which='both')

# ax_mag.grid(which='minor', alpha=0.3)

# ax_mag.grid(which='major', alpha=0.9)

# ax_phase.grid(which='minor', alpha=0.3)

# ax_phase.grid(which='major', alpha=0.9)

# Label the frequency axis

ax_phase.set_xlabel("Frequency (Hz)" if Hz

else "Frequency (rad/sec)")

if len(syslist) == 1:

return mags[0], phases[0], omegas[0]

else:

'nyquist.primary_style': ['-', '-.'], # style for primary curve

'nyquist.mirror_style': ['--', ':'], # style for mirror curve

'nyquist.arrows': 2, # number of arrors 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.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 maker

syslist, omega=None, plot=True, omega_limits=None, omega_num=None,

label_freq=0, color=None, return_contour=False,

warn_encirclements=True, warn_nyquist=True, **kwargs):

"""Nyquist plot for a system

# Check to see if legacy 'Plot' keyword was used

if 'Plot' in kwargs:

warnings.warn("'Plot' keyword is deprecated in nyquist_plot; "

"use 'plot'", FutureWarning)

# Map 'Plot' keyword to 'plot' keyword

plot = kwargs.pop('Plot')

# Check to see if legacy 'labelFreq' keyword was used

if 'labelFreq' in kwargs:

warnings.warn("'labelFreq' keyword is deprecated in nyquist_plot; "

"use 'label_freq'", FutureWarning)

# Map 'labelFreq' keyword to 'label_freq' keyword

label_freq = kwargs.pop('labelFreq')

# Check to see if legacy 'arrow_width' or 'arrow_length' were used

if 'arrow_width' in kwargs or 'arrow_length' in kwargs:

warnings.warn(

"'arrow_width' and 'arrow_length' keywords are deprecated in "

"nyquist_plot; use `arrow_size` instead", FutureWarning)

kwargs['arrow_size'] = \

(kwargs.get('arrow_width', 0) + kwargs.get('arrow_length', 0)) / 2

kwargs.pop('arrow_width', False)

kwargs.pop('arrow_length', False)

# Get values for params (and pop from list to allow keyword use in plot)

omega_num_given = omega_num is not None

omega_num = config._get_param('freqplot', 'number_of_samples', omega_num)

arrows = config._get_param(

'nyquist', 'arrows', kwargs, _nyquist_defaults, pop=True)

arrow_size = config._get_param(

'nyquist', 'arrow_size', kwargs, _nyquist_defaults, pop=True)

arrow_style = config._get_param('nyquist', 'arrow_style', kwargs, None)

indent_radius = config._get_param(

'nyquist', 'indent_radius', kwargs, _nyquist_defaults, pop=True)

encirclement_threshold = config._get_param(

'nyquist', 'encirclement_threshold', kwargs,

_nyquist_defaults, pop=True)

indent_direction = config._get_param(

'nyquist', 'indent_direction', kwargs, _nyquist_defaults, pop=True)

indent_points = config._get_param(

'nyquist', 'indent_points', kwargs, _nyquist_defaults, pop=True)

max_curve_magnitude = config._get_param(

'nyquist', 'max_curve_magnitude', kwargs, _nyquist_defaults, pop=True)

max_curve_offset = config._get_param(

'nyquist', 'max_curve_offset', kwargs, _nyquist_defaults, pop=True)

start_marker = config._get_param(

'nyquist', 'start_marker', kwargs, _nyquist_defaults, pop=True)

start_marker_size = config._get_param(

'nyquist', 'start_marker_size', kwargs, _nyquist_defaults, pop=True)

# Set line styles for the curves

def _parse_linestyle(style_name, allow_false=False):

style = config._get_param(

'nyquist', style_name, kwargs, _nyquist_defaults, pop=True)

if isinstance(style, str):

# Only one style provided, use the default for the other

style = [style, _nyquist_defaults['nyquist.' + style_name][1]]

warnings.warn(

"use of a single string for linestyle will be deprecated "

" in a future release", PendingDeprecationWarning)

if (allow_false and style is False) or \

(isinstance(style, list) and len(style) == 2):

return style

else:

raise ValueError(f"invalid '{style_name}': {style}")

primary_style = _parse_linestyle('primary_style')

mirror_style = _parse_linestyle('mirror_style', allow_false=True)

# If argument was a singleton, turn it into a tuple

if not isinstance(syslist, (list, tuple)):

syslist = (syslist,)

# Determine the range of frequencies to use, based on args/features

omega, omega_range_given = _determine_omega_vector(

syslist, omega, omega_limits, omega_num, feature_periphery_decades=2)

# If omega was not specified explicitly, start at omega = 0

if not omega_range_given:

if omega_num_given:

# Just reset the starting point

omega[0] = 0.0

else:

# Insert points between the origin and the first frequency point

omega = np.concatenate((

np.linspace(0, omega[0], indent_points), omega[1:]))

# Go through each system and keep track of the results

counts, contours = [], []

for sys in syslist:

if not sys.issiso():

# TODO: Add MIMO nyquist plots.

raise ControlMIMONotImplemented(

"Nyquist plot currently only supports SISO systems.")

# Figure out the frequency range

omega_sys = np.asarray(omega)

# Determine the contour used to evaluate the Nyquist curve

if sys.isdtime(strict=True):

# Restrict frequencies for discrete-time systems

nyquistfrq = math.pi / sys.dt

if not omega_range_given:

# limit up to and including nyquist frequency

omega_sys = np.hstack((

omega_sys[omega_sys < nyquistfrq], nyquistfrq))

# Issue a warning if we are sampling above Nyquist

if np.any(omega_sys * sys.dt > np.pi) and warn_nyquist:

warnings.warn("evaluation above Nyquist frequency")

# do indentations in s-plane where it is more convenient

splane_contour = 1j * omega_sys

# Bend the contour around any poles on/near the imaginary axis

# TODO: smarter indent radius that depends on dcgain of system

# and timebase of discrete system.

if isinstance(sys, (StateSpace, TransferFunction)) \

and indent_direction != 'none':

if sys.isctime():

splane_poles = sys.poles()

splane_cl_poles = sys.feedback().poles()

else:

# map z-plane poles to s-plane, ignoring any at the origin

# because we don't need to indent for them

zplane_poles = sys.poles()

zplane_poles = zplane_poles[~np.isclose(abs(zplane_poles), 0.)]

splane_poles = np.log(zplane_poles) / sys.dt

zplane_cl_poles = sys.feedback().poles()

zplane_cl_poles = zplane_cl_poles[

~np.isclose(abs(zplane_poles), 0.)]

splane_cl_poles = np.log(zplane_cl_poles) / sys.dt

#

# Check to make sure indent radius is small enough

#

# If there is a closed loop pole that is near the imaginary access

# at a point that is near an open loop pole, it is possible that

# indentation might skip or create an extraneous encirclement.

# We check for that situation here and generate a warning if that

# could happen.

#

for p_cl in splane_cl_poles:

# See if any closed loop poles are near the imaginary axis

if abs(p_cl.real) <= indent_radius:

# See if any open loop poles are close to closed loop poles

p_ol = splane_poles[

(np.abs(splane_poles - p_cl)).argmin()]

if abs(p_ol - p_cl) <= indent_radius and \

warn_encirclements:

warnings.warn(

"indented contour may miss closed loop pole; "

"consider reducing indent_radius to be less than "

f"{abs(p_ol - p_cl):5.2g}", stacklevel=2)

#

# See if we should add some frequency points near imaginary poles

#

for p in splane_poles:

# See if we need to process this pole (skip if on the negative

# imaginary axis or not near imaginary axis + user override)

if p.imag < 0 or abs(p.real) > indent_radius or \

omega_range_given:

continue

# Find the frequencies before the pole frequency

below_points = np.argwhere(

splane_contour.imag - abs(p.imag) < -indent_radius)

if below_points.size > 0:

first_point = below_points[-1].item()

start_freq = p.imag - indent_radius

else:

# Add the points starting at the beginning of the contour

assert splane_contour[0] == 0

first_point = 0

start_freq = 0

# Find the frequencies after the pole frequency

above_points = np.argwhere(

splane_contour.imag - abs(p.imag) > indent_radius)

last_point = above_points[0].item()

# Add points for half/quarter circle around pole frequency

# (these will get indented left or right below)

splane_contour = np.concatenate((

splane_contour[0:first_point+1],

(1j * np.linspace(

start_freq, p.imag + indent_radius, indent_points)),

splane_contour[last_point:]))

# Indent points that are too close to a pole

for i, s in enumerate(splane_contour):

# Find the nearest pole

p = splane_poles[(np.abs(splane_poles - s)).argmin()]

# See if we need to indent around it

if abs(s - p) < indent_radius:

# Figure out how much to offset (simple trigonometry)

offset = np.sqrt(indent_radius ** 2 - (s - p).imag ** 2) \

- (s - p).real

# Figure out which way to offset the contour point

if p.real < 0 or (p.real == 0 and

indent_direction == 'right'):

# Indent to the right

splane_contour[i] += offset

elif p.real > 0 or (p.real == 0 and

indent_direction == 'left'):

# Indent to the left

splane_contour[i] -= offset

else:

raise ValueError("unknown value for indent_direction")

# change contour to z-plane if necessary

if sys.isctime():

contour = splane_contour

else:

contour = np.exp(splane_contour * sys.dt)

# Compute the primary curve

resp = sys(contour)

# Compute CW encirclements of -1 by integrating the (unwrapped) angle

phase = -unwrap(np.angle(resp + 1))

encirclements = np.sum(np.diff(phase)) / np.pi

count = int(np.round(encirclements, 0))

# Let the user know if the count might not make sense

if abs(encirclements - count) > encirclement_threshold and \

warn_encirclements:

warnings.warn(

"number of encirclements was a non-integer value; this can"

" happen is contour is not closed, possibly based on a"

" frequency range that does not include zero.")

#

# Make sure that the enciriclements match the Nyquist criterion

#

# If the user specifies the frequency points to use, it is possible

# to miss enciriclements, so we check here to make sure that the

# Nyquist criterion is actually satisfied.

#

if isinstance(sys, (StateSpace, TransferFunction)):

# Count the number of open/closed loop RHP poles

if sys.isctime():

if indent_direction == 'right':

P = (sys.poles().real > 0).sum()

else:

P = (sys.poles().real >= 0).sum()

Z = (sys.feedback().poles().real >= 0).sum()

else:

if indent_direction == 'right':

P = (np.abs(sys.poles()) > 1).sum()

else:

P = (np.abs(sys.poles()) >= 1).sum()

Z = (np.abs(sys.feedback().poles()) >= 1).sum()

# Check to make sure the results make sense; warn if not

if Z != count + P and warn_encirclements:

warnings.warn(

"number of encirclements does not match Nyquist criterion;"

" check frequency range and indent radius/direction",

UserWarning, stacklevel=2)

elif indent_direction == 'none' and any(sys.poles().real == 0) and \

warn_encirclements:

warnings.warn(

"system has pure imaginary poles but indentation is"

" turned off; results may be meaningless",

RuntimeWarning, stacklevel=2)

counts.append(count)

contours.append(contour)

if plot:

# Parse the arrows keyword

if not arrows:

arrow_pos = []

elif isinstance(arrows, int):

N = arrows

# Space arrows out, starting midway along each "region"

arrow_pos = np.linspace(0.5/N, 1 + 0.5/N, N, endpoint=False)

elif isinstance(arrows, (list, np.ndarray)):

arrow_pos = np.sort(np.atleast_1d(arrows))

else:

raise ValueError("unknown or unsupported arrow location")

# Set the arrow style

if arrow_style is None:

arrow_style = mpl.patches.ArrowStyle(

'simple', head_width=arrow_size, head_length=arrow_size)