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

'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

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 dictonary 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('bode', 'dB', kwargs, _bode_defaults, pop=True)

deg = config._get_param('bode', 'deg', kwargs, _bode_defaults, pop=True)

Hz = config._get_param('bode', 'Hz', kwargs, _bode_defaults, pop=True)

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 tuple

if not hasattr(syslist, '__iter__'):

syslist = (syslist,)

# Decide whether to go above Nyquist frequency

omega_range_given = True if omega is not None else False

if omega is None:

omega_num = config._get_param(

'freqplot', 'number_of_samples', omega_num)

if omega_limits is None:

# Select a default range if none is provided

omega = _default_frequency_range(syslist,

number_of_samples=omega_num)

else:

omega_range_given = True

omega_limits = np.asarray(omega_limits)

if len(omega_limits) != 2:

raise ValueError("len(omega_limits) must be 2")

if Hz:

omega_limits *= 2. * math.pi

omega = np.logspace(np.log10(omega_limits[0]),

np.log10(omega_limits[1]), num=omega_num,

endpoint=True)

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['fig']

ax_mag = fig.axes[0]

ax_phase = fig.axes[2]

sisotool = kwargs['sisotool']

del kwargs['fig']

del kwargs['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 = -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)

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, nyquistfrq))

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.mirror_style': '--',

'nyquist.arrows': 2,

'nyquist.arrow_size': 8,

'nyquist.indent_radius': 1e-1,

'nyquist.indent_direction': 'right',

omega_num=None, label_freq=0, color=None,

return_contour=False, warn_nyquist=True, *args, **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 = config._get_param('freqplot', 'number_of_samples', omega_num)

mirror_style = config._get_param(

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

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)

indent_direction = config._get_param(

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

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

if not hasattr(syslist, '__iter__'):

syslist = (syslist,)

# Decide whether to go above Nyquist frequency

omega_range_given = True if omega is not None else False

# Figure out the frequency limits

if omega is None:

if omega_limits is None:

# Select a default range if none is provided

omega = _default_frequency_range(

syslist, number_of_samples=omega_num)

# Replace first point with the origin

omega[0] = 0

else:

omega_range_given = True

omega_limits = np.asarray(omega_limits)

if len(omega_limits) != 2:

raise ValueError("len(omega_limits) must be 2")

omega = np.logspace(np.log10(omega_limits[0]),

np.log10(omega_limits[1]), num=omega_num,

endpoint=True)

# 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):

# Transform frequencies in 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")

# Transform frequencies to continuous domain

contour = np.exp(1j * omega * sys.dt)

else:

contour = 1j * omega_sys

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

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

sys.isctime() and indent_direction != 'none':

poles = sys.pole()

for i, s in enumerate(contour):

# Find the nearest pole

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

# See if we need to indent around it

if abs(s - p) < indent_radius:

if p.real < 0 or \

(p.real == 0 and indent_direction == 'right'):

# Indent to the right

contour[i] += \

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

elif p.real > 0 or \

(p.real == 0 and indent_direction == 'left'):

# Indent to the left

contour[i] -= \

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

else:

ValueError("unknown value for indent_direction")

# TODO: add code to indent around discrete poles on unit circle

# Compute the primary curve

resp = sys(contour)

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

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

count = int(np.round(np.sum(np.diff(phase)) / np.pi, 0))

counts.append(count)

contours.append(contour)

if plot:

# Parse the arrows keyword

if 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))

elif not arrows:

arrow_pos = []

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)

# Save the components of the response

x, y = resp.real, resp.imag

# Plot the primary curve

p = plt.plot(x, y, '-', color=color, *args, **kwargs)

c = p[0].get_color()

ax = plt.gca()

_add_arrows_to_line2D(

ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=1)

# Plot the mirror image

if mirror_style is not False:

p = plt.plot(x, -y, mirror_style, color=c, *args, **kwargs)

_add_arrows_to_line2D(

ax, p[0], arrow_pos, arrowstyle=arrow_style, dir=-1)

# Mark the -1 point

plt.plot([-1], [0], 'r+')

# Label the frequencies of the points

if label_freq:

ind = slice(None, None, label_freq)

for xpt, ypt, omegapt in zip(x[ind], y[ind], omega_sys[ind]):

# Convert to Hz

f = omegapt / (2 * np.pi)

# Factor out multiples of 1000 and limit the

# result to the range [-8, 8].

pow1000 = max(min(get_pow1000(f), 8), -8)

# Get the SI prefix.

prefix = gen_prefix(pow1000)

# Apply the text. (Use a space before the text to

# prevent overlap with the data.)

#

# np.round() is used because 0.99... appears

# instead of 1.0, and this would otherwise be

# truncated to 0.

plt.text(xpt, ypt, ' ' +

str(int(np.round(f / 1000 ** pow1000, 0))) + ' ' +

prefix + 'Hz')

if plot:

ax = plt.gca()

ax.set_xlabel("Real axis")

ax.set_ylabel("Imaginary axis")

ax.grid(color="lightgray")

# "Squeeze" the results

if len(syslist) == 1:

counts, contours = counts[0], contours[0]

# Return counts and (optionally) the contour we used

axes, line, arrow_locs=[0.2, 0.4, 0.6, 0.8],

arrowstyle='-|>', arrowsize=1, dir=1, transform=None):

"""

if not isinstance(line, mpl.lines.Line2D):

raise ValueError("expected a matplotlib.lines.Line2D object")

x, y = line.get_xdata(), line.get_ydata()

arrow_kw = {

"arrowstyle": arrowstyle,

}

color = line.get_color()

use_multicolor_lines = isinstance(color, np.ndarray)

if use_multicolor_lines:

raise NotImplementedError("multicolor lines not supported")

else:

arrow_kw['color'] = color

linewidth = line.get_linewidth()

if isinstance(linewidth, np.ndarray):

raise NotImplementedError("multiwidth lines not supported")

else:

arrow_kw['linewidth'] = linewidth

if transform is None:

transform = axes.transData

# Compute the arc length along the curve

s = np.cumsum(np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2))

arrows = []

for loc in arrow_locs:

n = np.searchsorted(s, s[-1] * loc)

# Figure out what direction to paint the arrow

if dir == 1:

arrow_tail = (x[n], y[n])

arrow_head = (np.mean(x[n:n + 2]), np.mean(y[n:n + 2]))

elif dir == -1:

# Orient the arrow in the other direction on the segment

arrow_tail = (x[n + 1], y[n + 1])

arrow_head = (np.mean(x[n:n + 2]), np.mean(y[n:n + 2]))

else:

raise ValueError("unknown value for keyword 'dir'")

p = mpl.patches.FancyArrowPatch(

arrow_tail, arrow_head, transform=transform, lw=0,

**arrow_kw)

axes.add_patch(p)

arrows.append(p)

"""Plot the "Gang of 4" transfer functions for a system

if not P.issiso() or not C.issiso():

# TODO: Add MIMO go4 plots.

raise ControlMIMONotImplemented(

"Gang of four is currently only implemented for SISO systems.")

# Get the default parameter values

dB = config._get_param('bode', 'dB', kwargs, _bode_defaults, pop=True)

Hz = config._get_param('bode', 'Hz', kwargs, _bode_defaults, pop=True)

grid = config._get_param('bode', 'grid', kwargs, _bode_defaults, pop=True)

# Compute the senstivity functions

L = P * C

S = feedback(1, L)

T = L * S

# Select a default range if none is provided

# TODO: This needs to be made more intelligent

if omega is None:

omega = _default_frequency_range((P, C, S))

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

# gangof4_plot will superimpose the data. See details in bode_plot.

plot_axes = {'t': None, 's': None, 'ps': None, 'cs': None}

for ax in plt.gcf().axes:

label = ax.get_label()

if label.startswith('control-gangof4-'):

key = label[len('control-gangof4-'):]

if key not in plot_axes:

raise RuntimeError(

"unknown gangof4 axis type '{}'".format(label))

plot_axes[key] = ax

# if any of the axes are missing, start from scratch

if any((ax is None for ax in plot_axes.values())):

plt.clf()

plot_axes = {'s': plt.subplot(221, label='control-gangof4-s'),

'ps': plt.subplot(222, label='control-gangof4-ps'),

'cs': plt.subplot(223, label='control-gangof4-cs'),

't': plt.subplot(224, label='control-gangof4-t')}

#

# Plot the four sensitivity functions

#

omega_plot = omega / (2. * math.pi) if Hz else omega

# TODO: Need to add in the mag = 1 lines

mag_tmp, phase_tmp, omega = S.frequency_response(omega)

mag = np.squeeze(mag_tmp)

if dB: