Skip to content

Update pole/zero and root locus plots to use _map/_plot pattern #953

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 our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 16 commits into from
Jan 12, 2024
Merged
Show file tree
Hide file tree
Changes from all 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
2 changes: 1 addition & 1 deletion control/freqplot.py
Original file line number Diff line number Diff line change
Expand Up @@ -162,7 +162,7 @@ def bode_plot(
values with no plot.
rcParams : dict
Override the default parameters used for generating plots.
Default is set up config.default['freqplot.rcParams'].
Default is set by config.default['freqplot.rcParams'].
wrap_phase : bool or float
If wrap_phase is `False` (default), then the phase will be unwrapped
so that it is continuously increasing or decreasing. If wrap_phase is
Expand Down
86 changes: 52 additions & 34 deletions control/grid.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,21 @@
import numpy as np
from numpy import cos, sin, sqrt, linspace, pi, exp
# grid.py - code to add gridlines to root locus and pole-zero diagrams
#
# This code generates grids for pole-zero diagrams (including root locus
# diagrams). Rather than just draw a grid in place, it uses the AxisArtist
# package to generate a custom grid that will scale with the figure.
#

import matplotlib.pyplot as plt
from mpl_toolkits.axisartist import SubplotHost
from mpl_toolkits.axisartist.grid_helper_curvelinear \
import GridHelperCurveLinear
import mpl_toolkits.axisartist.angle_helper as angle_helper
import numpy as np
from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D
from mpl_toolkits.axisartist import SubplotHost
from mpl_toolkits.axisartist.grid_helper_curvelinear import \
GridHelperCurveLinear
from numpy import cos, exp, linspace, pi, sin, sqrt

from .iosys import isdtime


class FormatterDMS(object):
Expand Down Expand Up @@ -65,14 +74,15 @@ def __call__(self, transform_xy, x1, y1, x2, y2):
return lon_min, lon_max, lat_min, lat_max


def sgrid():
def sgrid(scaling=None):
# From matplotlib demos:
# https://matplotlib.org/gallery/axisartist/demo_curvelinear_grid.html
# https://matplotlib.org/gallery/axisartist/demo_floating_axis.html

# PolarAxes.PolarTransform takes radian. However, we want our coordinate
# system in degree
# system in degrees
tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()

# polar projection, which involves cycle, and also has limits in
# its coordinates, needs a special method to find the extremes
# (min, max of the coordinate within the view).
Expand All @@ -89,6 +99,7 @@ def sgrid():
tr, extreme_finder=extreme_finder, grid_locator1=grid_locator1,
tick_formatter1=tick_formatter1)

# Set up an axes with a specialized grid helper
fig = plt.gcf()
ax = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)

Expand All @@ -97,15 +108,20 @@ def sgrid():
ax.axis[:].major_ticklabels.set_visible(visible)
ax.axis[:].major_ticks.set_visible(False)
ax.axis[:].invert_ticklabel_direction()
ax.axis[:].major_ticklabels.set_color('gray')

# Set up internal tickmarks and labels along the real/imag axes
ax.axis["wnxneg"] = axis = ax.new_floating_axis(0, 180)
axis.set_ticklabel_direction("-")
axis.label.set_visible(False)

ax.axis["wnxpos"] = axis = ax.new_floating_axis(0, 0)
axis.label.set_visible(False)

ax.axis["wnypos"] = axis = ax.new_floating_axis(0, 90)
axis.label.set_visible(False)
axis.set_axis_direction("left")
axis.set_axis_direction("right")

ax.axis["wnyneg"] = axis = ax.new_floating_axis(0, 270)
axis.label.set_visible(False)
axis.set_axis_direction("left")
Expand All @@ -119,43 +135,41 @@ def sgrid():
ax.axis["bottom"].get_helper().nth_coord_ticks = 0

fig.add_subplot(ax)

# RECTANGULAR X Y AXES WITH SCALE
# par2 = ax.twiny()
# par2.axis["top"].toggle(all=False)
# par2.axis["right"].toggle(all=False)
# new_fixed_axis = par2.get_grid_helper().new_fixed_axis
# par2.axis["left"] = new_fixed_axis(loc="left",
# axes=par2,
# offset=(0, 0))
# par2.axis["bottom"] = new_fixed_axis(loc="bottom",
# axes=par2,
# offset=(0, 0))
# FINISH RECTANGULAR

ax.grid(True, zorder=0, linestyle='dotted')

_final_setup(ax)
_final_setup(ax, scaling=scaling)
return ax, fig


def _final_setup(ax):
# Utility function used by all grid code
def _final_setup(ax, scaling=None):
ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.axhline(y=0, color='black', lw=1)
ax.axvline(x=0, color='black', lw=1)
plt.axis('equal')
ax.axhline(y=0, color='black', lw=0.25)
ax.axvline(x=0, color='black', lw=0.25)

# Set up the scaling for the axes
scaling = 'equal' if scaling is None else scaling
plt.axis(scaling)

def nogrid():
f = plt.gcf()
ax = plt.axes()

_final_setup(ax)
return ax, f
# If not grid is given, at least separate stable/unstable regions
def nogrid(dt=None, ax=None, scaling=None):
fig = plt.gcf()
if ax is None:
ax = fig.gca()

# Draw the unit circle for discrete time systems
if isdtime(dt=dt, strict=True):
s = np.linspace(0, 2*pi, 100)
ax.plot(np.cos(s), np.sin(s), 'k--', lw=0.5, dashes=(5, 5))

_final_setup(ax, scaling=scaling)
return ax, fig

def zgrid(zetas=None, wns=None, ax=None):
# Grid for discrete time system (drawn, not rendered by AxisArtist)
# TODO (at some point): think about using customized grid generator?
def zgrid(zetas=None, wns=None, ax=None, scaling=None):
"""Draws discrete damping and frequency grid"""

fig = plt.gcf()
Expand Down Expand Up @@ -206,5 +220,9 @@ def zgrid(zetas=None, wns=None, ax=None):
ax.annotate(r"$\frac{"+num+r"\pi}{T}$", xy=(an_x, an_y),
xytext=(an_x, an_y), size=9)

_final_setup(ax)
# Set default axes to allow some room around the unit circle
ax.set_xlim([-1.1, 1.1])
ax.set_ylim([-1.1, 1.1])

_final_setup(ax, scaling=scaling)
return ax, fig
48 changes: 35 additions & 13 deletions control/iosys.py
Original file line number Diff line number Diff line change
Expand Up @@ -503,42 +503,64 @@ def common_timebase(dt1, dt2):
raise ValueError("Systems have incompatible timebases")

# Check to see if a system is a discrete time system
def isdtime(sys, strict=False):
def isdtime(sys=None, strict=False, dt=None):
"""
Check to see if a system is a discrete time system.

Parameters
----------
sys : I/O or LTI system
System to be checked
strict: bool (default = False)
If strict is True, make sure that timebase is not None
sys : I/O system, optional
System to be checked.
dt : None or number, optional
Timebase to be checked.
strict: bool, default=False
If strict is True, make sure that timebase is not None.
"""

# Check to see if this is a constant
# See if we were passed a timebase instead of a system
if sys is None:
if dt is None:
return True if not strict else False
else:
return dt > 0
elif dt is not None:
raise TypeError("passing both system and timebase not allowed")

# Check timebase of the system
if isinstance(sys, (int, float, complex, np.number)):
# OK as long as strict checking is off
# Constants OK as long as strict checking is off
return True if not strict else False
else:
return sys.isdtime(strict)


# Check to see if a system is a continuous time system
def isctime(sys, strict=False):
def isctime(sys=None, dt=None, strict=False):
"""
Check to see if a system is a continuous-time system.

Parameters
----------
sys : I/O or LTI system
System to be checked
sys : I/O system, optional
System to be checked.
dt : None or number, optional
Timebase to be checked.
strict: bool (default = False)
If strict is True, make sure that timebase is not None
If strict is True, make sure that timebase is not None.
"""

# Check to see if this is a constant
# See if we were passed a timebase instead of a system
if sys is None:
if dt is None:
return True if not strict else False
else:
return dt == 0
elif dt is not None:
raise TypeError("passing both system and timebase not allowed")

# Check timebase of the system
if isinstance(sys, (int, float, complex, np.number)):
# OK as long as strict checking is off
# Constants OK as long as strict checking is off
return True if not strict else False
else:
return sys.isctime(strict)
Expand Down
105 changes: 103 additions & 2 deletions control/matlab/wrappers.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,14 @@
from ..lti import LTI
from ..exception import ControlArgument

__all__ = ['bode', 'nyquist', 'ngrid', 'dcgain', 'connect']
__all__ = ['bode', 'nyquist', 'ngrid', 'rlocus', 'pzmap', 'dcgain', 'connect']

def bode(*args, **kwargs):
"""bode(syslist[, omega, dB, Hz, deg, ...])

Bode plot of the frequency response.

Plots a bode gain and phase diagram
Plots a bode gain and phase diagram.

Parameters
----------
Expand Down Expand Up @@ -195,6 +195,106 @@ def _parse_freqplot_args(*args):
return syslist, omega, plotstyle, other


# TODO: rewrite to call root_locus_map, without using legacy plot keyword
def rlocus(*args, **kwargs):
"""rlocus(sys[, klist, xlim, ylim, ...])

Root locus diagram.

Calculate the root locus by finding the roots of 1 + k * G(s) where G
is a linear system with transfer function num(s)/den(s) and each k is
an element of kvect.

Parameters
----------
sys : LTI object
Linear input/output systems (SISO only, for now).
kvect : array_like, optional
Gains to use in computing plot of closed-loop poles.
xlim : tuple or list, optional
Set limits of x axis, normally with tuple
(see :doc:`matplotlib:api/axes_api`).
ylim : tuple or list, optional
Set limits of y axis, normally with tuple
(see :doc:`matplotlib:api/axes_api`).

Returns
-------
roots : ndarray
Closed-loop root locations, arranged in which each row corresponds
to a gain in gains.
gains : ndarray
Gains used. Same as kvect keyword argument if provided.

Notes
-----
This function is a wrapper for :func:`~control.root_locus_plot`,
with legacy return arguments.

"""
from ..rlocus import root_locus_plot

# Use the plot keyword to get legacy behavior
kwargs = dict(kwargs) # make a copy since we modify this
if 'plot' not in kwargs:
kwargs['plot'] = True

# Turn off deprecation warning
with warnings.catch_warnings():
warnings.filterwarnings(
'ignore', message='.* return values of .* is deprecated',
category=DeprecationWarning)
retval = root_locus_plot(*args, **kwargs)

return retval


# TODO: rewrite to call pole_zero_map, without using legacy plot keyword
def pzmap(*args, **kwargs):
"""pzmap(sys[, grid, plot])

Plot a pole/zero map for a linear system.

Parameters
----------
sys: LTI (StateSpace or TransferFunction)
Linear system for which poles and zeros are computed.
plot: bool, optional
If ``True`` a graph is generated with Matplotlib,
otherwise the poles and zeros are only computed and returned.
grid: boolean (default = False)
If True plot omega-damping grid.

Returns
-------
poles: array
The system's poles.
zeros: array
The system's zeros.

Notes
-----
This function is a wrapper for :func:`~control.pole_zero_plot`,
with legacy return arguments.

"""
from ..pzmap import pole_zero_plot

# Use the plot keyword to get legacy behavior
kwargs = dict(kwargs) # make a copy since we modify this
if 'plot' not in kwargs:
kwargs['plot'] = True

# Turn off deprecation warning
with warnings.catch_warnings():
warnings.filterwarnings(
'ignore', message='.* return values of .* is deprecated',
category=DeprecationWarning)
retval = pole_zero_plot(*args, **kwargs)

return retval


from ..nichols import nichols_grid
def ngrid():
return nichols_grid()
Expand Down Expand Up @@ -254,6 +354,7 @@ def dcgain(*args):

from ..bdalg import connect as ct_connect
def connect(*args):

"""Index-based interconnection of an LTI system.

The system `sys` is a system typically constructed with `append`, with
Expand Down
Loading