11994777565

Committed 24 Nov 2024 09:02AM UTC coverage: 94.693% (+0.001%) from 94.692%

Build # 11994777565

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1065 from murrayrm/rlocus_singleton-21Nov2024

allow root locus maps with only 1 point

Run Details

9172 of 9686 relevant lines covered (94.69%)

8.27 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.72

control/rlocus.py

# rlocus.py - code for computing a root locus plot
# Code contributed by Ryan Krauss, 2010
#
# RMM, 17 June 2010: modified to be a standalone piece of code
#   * Added BSD copyright info to file (per Ryan)
#   * Added code to convert (num, den) to poly1d's if they aren't already.
#     This allows Ryan's code to run on a standard signal.ltisys object
#     or a control.TransferFunction object.
#   * Added some comments to make sure I understand the code
#
# RMM, 2 April 2011: modified to work with new LTI structure (see ChangeLog)
#   * Not tested: should still work on signal.ltisys objects
#
# Sawyer B. Fuller (minster@uw.edu) 21 May 2020:
#   * added compatibility with discrete-time systems.
#

import warnings
from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal  # signal processing toolbox
from numpy import array, imag, poly1d, real, vstack, zeros_like

from . import config
from .ctrlplot import ControlPlot
from .exception import ControlMIMONotImplemented
from .iosys import isdtime
from .lti import LTI
from .xferfcn import _convert_to_transfer_function

__all__ = ['root_locus_map', 'root_locus_plot', 'root_locus', 'rlocus']

# Default values for module parameters
_rlocus_defaults = {
    'rlocus.grid': True,
}


# Root locus map
def root_locus_map(sysdata, gains=None):
    """Compute the root locus map for an LTI system.

    Calculate the root locus by finding the roots of 1 + k * G(s) where G
    is a linear system and k varies over a range of gains.

    Parameters
    ----------
    sysdata : LTI system or list of LTI systems
        Linear input/output systems (SISO only, for now).
    gains : array_like, optional
        Gains to use in computing plot of closed-loop poles.  If not given,
        gains are chosen to include the main features of the root locus map.

    Returns
    -------
    rldata : PoleZeroData or list of PoleZeroData
        Root locus data object(s).  The loci of the root locus diagram are
        available in the array `rldata.loci`, indexed by the gain index and
        the locus index, and the gains are in the array `rldata.gains`.

    Notes
    -----
    For backward compatibility, the `rldata` return object can be
    assigned to the tuple `roots, gains`.

    """
    from .pzmap import PoleZeroData, PoleZeroList

    # Convert the first argument to a list
    syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]

    responses = []
    for idx, sys in enumerate(syslist):
        if not sys.issiso():
            raise ControlMIMONotImplemented(
                "sys must be single-input single-output (SISO)")

        # Convert numerator and denominator to polynomials if they aren't
        nump, denp = _systopoly1d(sys[0, 0])

        if gains is None:
            kvect, root_array, _, _ = _default_gains(nump, denp, None, None)
        else:
            kvect = np.atleast_1d(gains)
            root_array = _RLFindRoots(nump, denp, kvect)
            root_array = _RLSortRoots(root_array)

        responses.append(PoleZeroData(
            sys.poles(), sys.zeros(), kvect, root_array,
            dt=sys.dt, sysname=sys.name, sys=sys))

    if isinstance(sysdata, (list, tuple)):
        return PoleZeroList(responses)
    else:
        return responses[0]


def root_locus_plot(
        sysdata, gains=None, grid=None, plot=None, **kwargs):

    """Root locus plot.

    Calculate the root locus by finding the roots of 1 + k * G(s) where G
    is a linear system and k varies over a range of gains.

    Parameters
    ----------
    sysdata : PoleZeroMap or LTI object or list
        Linear input/output systems (SISO only, for now).
    gains : array_like, optional
        Gains to use in computing plot of closed-loop poles.  If not given,
        gains are chosen to include the main features of the root locus map.
    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`).
    plot : bool, optional
        (legacy) If given, `root_locus_plot` returns the legacy return values
        of roots and gains.  If False, just return the values with no plot.
    grid : bool or str, optional
        If `True` plot omega-damping grid, if `False` show imaginary axis
        for continuous time systems, unit circle for discrete time systems.
        If `empty`, do not draw any additonal lines.  Default value is set
        by config.default['rlocus.grid'].
    initial_gain : float, optional
        Mark the point on the root locus diagram corresponding to the
        given gain.
    color : matplotlib color spec, optional
        Specify the color of the markers and lines.

    Returns
    -------
    cplt : :class:`ControlPlot` object
        Object containing the data that were plotted:

          * cplt.lines: Array of :class:`matplotlib.lines.Line2D` objects
            for each set of markers in the plot. The shape of the array is
            given by (nsys, 3) where nsys is the number of systems or
            responses passed to the function.  The second index specifies
            the object type:

              - lines[idx, 0]: poles
              - lines[idx, 1]: zeros
              - lines[idx, 2]: loci

          * cplt.axes: 2D array of :class:`matplotlib.axes.Axes` for the plot.

          * cplt.figure: :class:`matplotlib.figure.Figure` containing the plot.

        See :class:`ControlPlot` for more detailed information.

    roots, gains : ndarray
        (legacy) If the `plot` keyword is given, returns the closed-loop
        root locations, arranged such that each row corresponds to a gain,
        and the array of gains (same as `gains` keyword argument if provided).

    Other Parameters
    ----------------
    ax : matplotlib.axes.Axes, optional
        The matplotlib axes to draw the figure on.  If not specified and
        the current figure has a single axes, that axes is used.
        Otherwise, a new figure is created.
    label : str or array_like of str, optional
        If present, replace automatically generated label(s) with the given
        label(s).  If sysdata is a list, strings should be specified for each
        system.
    legend_loc : int or str, optional
        Include a legend in the given location. Default is 'center right',
        with no legend for a single response.  Use False to suppress legend.
    show_legend : bool, optional
        Force legend to be shown if ``True`` or hidden if ``False``.  If
        ``None``, then show legend when there is more than one line on the
        plot or ``legend_loc`` has been specified.
    title : str, optional
        Set the title of the plot.  Defaults to plot type and system name(s).

    Notes
    -----
    The root_locus_plot function calls matplotlib.pyplot.axis('equal'), which
    means that trying to reset the axis limits may not behave as expected.
    To change the axis limits, use matplotlib.pyplot.gca().axis('auto') and
    then set the axis limits to the desired values.

    """
    from .pzmap import pole_zero_plot

    # Legacy parameters
    for oldkey in ['kvect', 'k']:
        gains = config._process_legacy_keyword(kwargs, oldkey, 'gains', gains)

    if isinstance(sysdata, list) and all(
            [isinstance(sys, LTI) for sys in sysdata]) or \
            isinstance(sysdata, LTI):
        responses = root_locus_map(sysdata, gains=gains)
    else:
        responses = sysdata

    #
    # Process `plot` keyword
    #
    # See bode_plot for a description of how this keyword is handled to
    # support legacy implementatoins of root_locus.
    #
    if plot is not None:
        warnings.warn(
            "root_locus() return value of roots, gains is deprecated; "
            "use root_locus_map()", FutureWarning)

    if plot is False:
        return responses.loci, responses.gains

    # Plot the root loci
    cplt = responses.plot(grid=grid, **kwargs)

    # Legacy processing: return locations of poles and zeros as a tuple
    if plot is True:
        return responses.loci, responses.gains

    return ControlPlot(cplt.lines, cplt.axes, cplt.figure)


def _default_gains(num, den, xlim, ylim):
    """Unsupervised gains calculation for root locus plot.

    References
    ----------
    Ogata, K. (2002). Modern control engineering (4th ed.). Upper
    Saddle River, NJ : New Delhi: Prentice Hall..

    """
    # Compute the break points on the real axis for the root locus plot
    k_break, real_break = _break_points(num, den)

    # Decide on the maximum gain to use and create the gain vector
    kmax = _k_max(num, den, real_break, k_break)
    kvect = np.hstack((np.linspace(0, kmax, 50), np.real(k_break)))
    kvect.sort()

    # Find the roots for all of the gains and sort them
    root_array = _RLFindRoots(num, den, kvect)
    root_array = _RLSortRoots(root_array)

    # Keep track of the open loop poles and zeros
    open_loop_poles = den.roots
    open_loop_zeros = num.roots

    # ???
    if open_loop_zeros.size != 0 and \
       open_loop_zeros.size < open_loop_poles.size:
        open_loop_zeros_xl = np.append(
            open_loop_zeros,
            np.ones(open_loop_poles.size - open_loop_zeros.size)
            * open_loop_zeros[-1])
        root_array_xl = np.append(root_array, open_loop_zeros_xl)
    else:
        root_array_xl = root_array
    singular_points = np.concatenate((num.roots, den.roots), axis=0)
    important_points = np.concatenate((singular_points, real_break), axis=0)
    important_points = np.concatenate((important_points, np.zeros(2)), axis=0)
    root_array_xl = np.append(root_array_xl, important_points)

    false_gain = float(den.coeffs[0]) / float(num.coeffs[0])
    if false_gain < 0 and not den.order > num.order:
        # TODO: make error message more understandable
        raise ValueError("Not implemented support for 0 degrees root locus "
                         "with equal order of numerator and denominator.")

    if xlim is None and false_gain > 0:
        x_tolerance = 0.05 * (np.max(np.real(root_array_xl))
                              - np.min(np.real(root_array_xl)))
        xlim = _ax_lim(root_array_xl)
    elif xlim is None and false_gain < 0:
        axmin = np.min(np.real(important_points)) \
            - (np.max(np.real(important_points))
               - np.min(np.real(important_points)))
        axmin = np.min(np.array([axmin, np.min(np.real(root_array_xl))]))
        axmax = np.max(np.real(important_points)) \
            + np.max(np.real(important_points)) \
            - np.min(np.real(important_points))
        axmax = np.max(np.array([axmax, np.max(np.real(root_array_xl))]))
        xlim = [axmin, axmax]
        x_tolerance = 0.05 * (axmax - axmin)
    else:
        x_tolerance = 0.05 * (xlim[1] - xlim[0])

    if ylim is None:
        y_tolerance = 0.05 * (np.max(np.imag(root_array_xl))
                              - np.min(np.imag(root_array_xl)))
        ylim = _ax_lim(root_array_xl * 1j)
    else:
        y_tolerance = 0.05 * (ylim[1] - ylim[0])

    # Figure out which points are spaced too far apart
    if x_tolerance == 0:
        # Root locus is on imaginary axis (rare), use just y distance
        tolerance = y_tolerance
    elif y_tolerance == 0:
        # Root locus is on imaginary axis (common), use just x distance
        tolerance = x_tolerance
    else:
        tolerance = np.min([x_tolerance, y_tolerance])
    indexes_too_far = _indexes_filt(root_array, tolerance)

    # Add more points into the root locus for points that are too far apart
    while len(indexes_too_far) > 0 and kvect.size < 5000:
        for counter, index in enumerate(indexes_too_far):
            index = index + counter*3
            new_gains = np.linspace(kvect[index], kvect[index + 1], 5)
            new_points = _RLFindRoots(num, den, new_gains[1:4])
            kvect = np.insert(kvect, index + 1, new_gains[1:4])
            root_array = np.insert(root_array, index + 1, new_points, axis=0)

        root_array = _RLSortRoots(root_array)
        indexes_too_far = _indexes_filt(root_array, tolerance)

    new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4)))
    new_points = _RLFindRoots(num, den, new_gains[1:4])
    kvect = np.append(kvect, new_gains[1:4])
    root_array = np.concatenate((root_array, new_points), axis=0)
    root_array = _RLSortRoots(root_array)
    return kvect, root_array, xlim, ylim


def _indexes_filt(root_array, tolerance):
    """Calculate the distance between points and return the indices.

    Filter the indexes so only the resolution of points within the xlim and
    ylim is improved when zoom is used.

    """
    distance_points = np.abs(np.diff(root_array, axis=0))
    indexes_too_far = list(np.unique(np.where(distance_points > tolerance)[0]))
    indexes_too_far.sort()
    return indexes_too_far


def _break_points(num, den):
    """Extract break points over real axis and gains given these locations"""
    # type: (np.poly1d, np.poly1d) -> (np.array, np.array)
    dnum = num.deriv(m=1)
    dden = den.deriv(m=1)
    polynom = den * dnum - num * dden
    real_break_pts = polynom.r
    # don't care about infinite break points
    real_break_pts = real_break_pts[num(real_break_pts) != 0]
    k_break = -den(real_break_pts) / num(real_break_pts)
    idx = k_break >= 0   # only positives gains
    k_break = k_break[idx]
    real_break_pts = real_break_pts[idx]
    if len(k_break) == 0:
        k_break = [0]
        real_break_pts = den.roots
    return k_break, real_break_pts


def _ax_lim(root_array):
    """Utility to get the axis limits"""
    axmin = np.min(np.real(root_array))
    axmax = np.max(np.real(root_array))
    if axmax != axmin:
        deltax = (axmax - axmin) * 0.02
    else:
        deltax = np.max([1., axmax / 2])
    axlim = [axmin - deltax, axmax + deltax]
    return axlim


def _k_max(num, den, real_break_points, k_break_points):
    """"Calculate the maximum gain for the root locus shown in the figure."""
    asymp_number = den.order - num.order
    singular_points = np.concatenate((num.roots, den.roots), axis=0)
    important_points = np.concatenate(
        (singular_points, real_break_points), axis=0)
    false_gain = den.coeffs[0] / num.coeffs[0]

    if asymp_number > 0:
        asymp_center = (np.sum(den.roots) - np.sum(num.roots))/asymp_number
        distance_max = 4 * np.max(np.abs(important_points - asymp_center))
        asymp_angles = (2 * np.arange(0, asymp_number) - 1) \
            * np.pi / asymp_number
        if false_gain > 0:
            # farthest points over asymptotes
            farthest_points = asymp_center \
                + distance_max * np.exp(asymp_angles * 1j)
        else:
            asymp_angles = asymp_angles + np.pi
            # farthest points over asymptotes
            farthest_points = asymp_center \
                + distance_max * np.exp(asymp_angles * 1j)
        kmax_asymp = np.real(np.abs(den(farthest_points)
                                    / num(farthest_points)))
    else:
        kmax_asymp = np.abs([np.abs(den.coeffs[0])
                             / np.abs(num.coeffs[0]) * 3])

    kmax = np.max(np.concatenate((np.real(kmax_asymp),
                                  np.real(k_break_points)), axis=0))
    if np.abs(false_gain) > kmax:
        kmax = np.abs(false_gain)
    return kmax


def _systopoly1d(sys):
    """Extract numerator and denominator polynomails for a system"""
    # Allow inputs from the signal processing toolbox
    if (isinstance(sys, scipy.signal.lti)):
        nump = sys.num
        denp = sys.den

    else:
        # Convert to a transfer function, if needed
        sys = _convert_to_transfer_function(sys)

        # Make sure we have a SISO system
        if not sys.issiso():
            raise ControlMIMONotImplemented()

        # Start by extracting the numerator and denominator from system object
        nump = sys.num[0][0]
        denp = sys.den[0][0]

    # Check to see if num, den are already polynomials; otherwise convert
    if (not isinstance(nump, poly1d)):
        nump = poly1d(nump)

    if (not isinstance(denp, poly1d)):
        denp = poly1d(denp)

    return (nump, denp)


def _RLFindRoots(nump, denp, kvect):
    """Find the roots for the root locus."""
    # Convert numerator and denominator to polynomials if they aren't
    roots = []
    for k in np.atleast_1d(kvect):
        curpoly = denp + k * nump
        curroots = curpoly.r
        if len(curroots) < denp.order:
            # if I have fewer poles than open loop, it is because i have
            # one at infinity
            curroots = np.append(curroots, np.inf)

        curroots.sort()
        roots.append(curroots)

    return vstack(roots)


def _RLSortRoots(roots):
    """Sort the roots from _RLFindRoots, so that the root
    locus doesn't show weird pseudo-branches as roots jump from
    one branch to another."""

    sorted = zeros_like(roots)
    for n, row in enumerate(roots):
        if n == 0:
            sorted[n, :] = row
        else:
            # sort the current row by finding the element with the
            # smallest absolute distance to each root in the
            # previous row
            available = list(range(len(prevrow)))
            for elem in row:
                evect = elem - prevrow[available]
                ind1 = abs(evect).argmin()
                ind = available.pop(ind1)
                sorted[n, ind] = elem
        prevrow = sorted[n, :]
    return sorted


# Alternative ways to call these functions
root_locus = root_locus_plot
rlocus = root_locus_plot

1	# rlocus.py - code for computing a root locus plot
2	# Code contributed by Ryan Krauss, 2010
3	#
4	# RMM, 17 June 2010: modified to be a standalone piece of code
5	# * Added BSD copyright info to file (per Ryan)
6	# * Added code to convert (num, den) to poly1d's if they aren't already.
7	# This allows Ryan's code to run on a standard signal.ltisys object
8	# or a control.TransferFunction object.
9	# * Added some comments to make sure I understand the code
10	#
11	# RMM, 2 April 2011: modified to work with new LTI structure (see ChangeLog)
12	# * Not tested: should still work on signal.ltisys objects
13	#
14	# Sawyer B. Fuller (minster@uw.edu) 21 May 2020:
15	# * added compatibility with discrete-time systems.
16	#
17
18	import warnings	9✔
19	from functools import partial	9✔
20
21	import matplotlib.pyplot as plt	9✔
22	import numpy as np	9✔
23	import scipy.signal # signal processing toolbox	9✔
24	from numpy import array, imag, poly1d, real, vstack, zeros_like	9✔
25
26	from . import config	9✔
27	from .ctrlplot import ControlPlot	9✔
28	from .exception import ControlMIMONotImplemented	9✔
29	from .iosys import isdtime	9✔
30	from .lti import LTI	9✔
31	from .xferfcn import _convert_to_transfer_function	9✔
32
33	__all__ = ['root_locus_map', 'root_locus_plot', 'root_locus', 'rlocus']	9✔
34
35	# Default values for module parameters
36	_rlocus_defaults = {	9✔
37	'rlocus.grid': True,
38	}
39
40
41	# Root locus map
42	def root_locus_map(sysdata, gains=None):	9✔
43	"""Compute the root locus map for an LTI system.
44
45	Calculate the root locus by finding the roots of 1 + k * G(s) where G
46	is a linear system and k varies over a range of gains.
47
48	Parameters
49	----------
50	sysdata : LTI system or list of LTI systems
51	Linear input/output systems (SISO only, for now).
52	gains : array_like, optional
53	Gains to use in computing plot of closed-loop poles. If not given,
54	gains are chosen to include the main features of the root locus map.
55
56	Returns
57	-------
58	rldata : PoleZeroData or list of PoleZeroData
59	Root locus data object(s). The loci of the root locus diagram are
60	available in the array `rldata.loci`, indexed by the gain index and
61	the locus index, and the gains are in the array `rldata.gains`.
62
63	Notes
64	-----
65	For backward compatibility, the `rldata` return object can be
66	assigned to the tuple `roots, gains`.
67
68	"""
69	from .pzmap import PoleZeroData, PoleZeroList	9✔
70
71	# Convert the first argument to a list
72	syslist = sysdata if isinstance(sysdata, (list, tuple)) else [sysdata]	9✔
73
74	responses = []	9✔
75	for idx, sys in enumerate(syslist):	9✔
76	if not sys.issiso():	9✔
77	raise ControlMIMONotImplemented(	×
78	"sys must be single-input single-output (SISO)")
79
80	# Convert numerator and denominator to polynomials if they aren't
81	nump, denp = _systopoly1d(sys[0, 0])	9✔
82
83	if gains is None:	9✔
84	kvect, root_array, _, _ = _default_gains(nump, denp, None, None)	9✔
85	else:
86	kvect = np.atleast_1d(gains)	9✔
87	root_array = _RLFindRoots(nump, denp, kvect)	9✔
88	root_array = _RLSortRoots(root_array)	9✔
89
90	responses.append(PoleZeroData(	9✔
91	sys.poles(), sys.zeros(), kvect, root_array,
92	dt=sys.dt, sysname=sys.name, sys=sys))
93
94	if isinstance(sysdata, (list, tuple)):	9✔
95	return PoleZeroList(responses)	9✔
96	else:
97	return responses[0]	9✔
98
99
100	def root_locus_plot(	9✔
101	sysdata, gains=None, grid=None, plot=None, **kwargs):
102
103	"""Root locus plot.
104
105	Calculate the root locus by finding the roots of 1 + k * G(s) where G
106	is a linear system and k varies over a range of gains.
107
108	Parameters
109	----------
110	sysdata : PoleZeroMap or LTI object or list
111	Linear input/output systems (SISO only, for now).
112	gains : array_like, optional
113	Gains to use in computing plot of closed-loop poles. If not given,
114	gains are chosen to include the main features of the root locus map.
115	xlim : tuple or list, optional
116	Set limits of x axis, normally with tuple
117	(see :doc:`matplotlib:api/axes_api`).
118	ylim : tuple or list, optional
119	Set limits of y axis, normally with tuple
120	(see :doc:`matplotlib:api/axes_api`).
121	plot : bool, optional
122	(legacy) If given, `root_locus_plot` returns the legacy return values
123	of roots and gains. If False, just return the values with no plot.
124	grid : bool or str, optional
125	If `True` plot omega-damping grid, if `False` show imaginary axis
126	for continuous time systems, unit circle for discrete time systems.
127	If `empty`, do not draw any additonal lines. Default value is set
128	by config.default['rlocus.grid'].
129	initial_gain : float, optional
130	Mark the point on the root locus diagram corresponding to the
131	given gain.
132	color : matplotlib color spec, optional
133	Specify the color of the markers and lines.
134
135	Returns
136	-------
137	cplt : :class:`ControlPlot` object
138	Object containing the data that were plotted:
139
140	* cplt.lines: Array of :class:`matplotlib.lines.Line2D` objects
141	for each set of markers in the plot. The shape of the array is
142	given by (nsys, 3) where nsys is the number of systems or
143	responses passed to the function. The second index specifies
144	the object type:
145
146	- lines[idx, 0]: poles
147	- lines[idx, 1]: zeros
148	- lines[idx, 2]: loci
149
150	* cplt.axes: 2D array of :class:`matplotlib.axes.Axes` for the plot.
151
152	* cplt.figure: :class:`matplotlib.figure.Figure` containing the plot.
153
154	See :class:`ControlPlot` for more detailed information.
155
156	roots, gains : ndarray
157	(legacy) If the `plot` keyword is given, returns the closed-loop
158	root locations, arranged such that each row corresponds to a gain,
159	and the array of gains (same as `gains` keyword argument if provided).
160
161	Other Parameters
162	----------------
163	ax : matplotlib.axes.Axes, optional
164	The matplotlib axes to draw the figure on. If not specified and
165	the current figure has a single axes, that axes is used.
166	Otherwise, a new figure is created.
167	label : str or array_like of str, optional
168	If present, replace automatically generated label(s) with the given
169	label(s). If sysdata is a list, strings should be specified for each
170	system.
171	legend_loc : int or str, optional
172	Include a legend in the given location. Default is 'center right',
173	with no legend for a single response. Use False to suppress legend.
174	show_legend : bool, optional
175	Force legend to be shown if ``True`` or hidden if ``False``. If
176	``None``, then show legend when there is more than one line on the
177	plot or ``legend_loc`` has been specified.
178	title : str, optional
179	Set the title of the plot. Defaults to plot type and system name(s).
180
181	Notes
182	-----
183	The root_locus_plot function calls matplotlib.pyplot.axis('equal'), which
184	means that trying to reset the axis limits may not behave as expected.
185	To change the axis limits, use matplotlib.pyplot.gca().axis('auto') and
186	then set the axis limits to the desired values.
187
188	"""
189	from .pzmap import pole_zero_plot	9✔
190
191	# Legacy parameters
192	for oldkey in ['kvect', 'k']:	9✔
193	gains = config._process_legacy_keyword(kwargs, oldkey, 'gains', gains)	9✔
194
195	if isinstance(sysdata, list) and all(	9✔
196	[isinstance(sys, LTI) for sys in sysdata]) or \
197	isinstance(sysdata, LTI):
198	responses = root_locus_map(sysdata, gains=gains)	9✔
199	else:
200	responses = sysdata	9✔
201
202	#
203	# Process `plot` keyword
204	#
205	# See bode_plot for a description of how this keyword is handled to
206	# support legacy implementatoins of root_locus.
207	#
208	if plot is not None:	9✔
209	warnings.warn(	9✔
210	"root_locus() return value of roots, gains is deprecated; "
211	"use root_locus_map()", FutureWarning)
212
213	if plot is False:	9✔
214	return responses.loci, responses.gains	9✔
215
216	# Plot the root loci
217	cplt = responses.plot(grid=grid, **kwargs)	9✔
218
219	# Legacy processing: return locations of poles and zeros as a tuple
220	if plot is True:	9✔
221	return responses.loci, responses.gains	9✔
222
223	return ControlPlot(cplt.lines, cplt.axes, cplt.figure)	9✔
224
225
226	def _default_gains(num, den, xlim, ylim):	9✔
227	"""Unsupervised gains calculation for root locus plot.
228
229	References
230	----------
231	Ogata, K. (2002). Modern control engineering (4th ed.). Upper
232	Saddle River, NJ : New Delhi: Prentice Hall..
233
234	"""
235	# Compute the break points on the real axis for the root locus plot
236	k_break, real_break = _break_points(num, den)	9✔
237
238	# Decide on the maximum gain to use and create the gain vector
239	kmax = _k_max(num, den, real_break, k_break)	9✔
240	kvect = np.hstack((np.linspace(0, kmax, 50), np.real(k_break)))	9✔
241	kvect.sort()	9✔
242
243	# Find the roots for all of the gains and sort them
244	root_array = _RLFindRoots(num, den, kvect)	9✔
245	root_array = _RLSortRoots(root_array)	9✔
246
247	# Keep track of the open loop poles and zeros
248	open_loop_poles = den.roots	9✔
249	open_loop_zeros = num.roots	9✔
250
251	# ???
252	if open_loop_zeros.size != 0 and \	9✔
253	open_loop_zeros.size < open_loop_poles.size:
254	open_loop_zeros_xl = np.append(	9✔
255	open_loop_zeros,
256	np.ones(open_loop_poles.size - open_loop_zeros.size)
257	* open_loop_zeros[-1])
258	root_array_xl = np.append(root_array, open_loop_zeros_xl)	9✔
259	else:
260	root_array_xl = root_array	9✔
261	singular_points = np.concatenate((num.roots, den.roots), axis=0)	9✔
262	important_points = np.concatenate((singular_points, real_break), axis=0)	9✔
263	important_points = np.concatenate((important_points, np.zeros(2)), axis=0)	9✔
264	root_array_xl = np.append(root_array_xl, important_points)	9✔
265
266	false_gain = float(den.coeffs[0]) / float(num.coeffs[0])	9✔
267	if false_gain < 0 and not den.order > num.order:	9✔
268	# TODO: make error message more understandable
269	raise ValueError("Not implemented support for 0 degrees root locus "	9✔
270	"with equal order of numerator and denominator.")
271
272	if xlim is None and false_gain > 0:	9✔
273	x_tolerance = 0.05 * (np.max(np.real(root_array_xl))	9✔
274	- np.min(np.real(root_array_xl)))
275	xlim = _ax_lim(root_array_xl)	9✔
276	elif xlim is None and false_gain < 0:	9✔
277	axmin = np.min(np.real(important_points)) \	9✔
278	- (np.max(np.real(important_points))
279	- np.min(np.real(important_points)))
280	axmin = np.min(np.array([axmin, np.min(np.real(root_array_xl))]))	9✔
281	axmax = np.max(np.real(important_points)) \	9✔
282	+ np.max(np.real(important_points)) \
283	- np.min(np.real(important_points))
284	axmax = np.max(np.array([axmax, np.max(np.real(root_array_xl))]))	9✔
285	xlim = [axmin, axmax]	9✔
286	x_tolerance = 0.05 * (axmax - axmin)	9✔
287	else:
288	x_tolerance = 0.05 * (xlim[1] - xlim[0])	×
289
290	if ylim is None:	9✔
291	y_tolerance = 0.05 * (np.max(np.imag(root_array_xl))	9✔
292	- np.min(np.imag(root_array_xl)))
293	ylim = _ax_lim(root_array_xl * 1j)	9✔
294	else:
295	y_tolerance = 0.05 * (ylim[1] - ylim[0])	×
296
297	# Figure out which points are spaced too far apart
298	if x_tolerance == 0:	9✔
299	# Root locus is on imaginary axis (rare), use just y distance
300	tolerance = y_tolerance	×
301	elif y_tolerance == 0:	9✔
302	# Root locus is on imaginary axis (common), use just x distance
303	tolerance = x_tolerance	9✔
304	else:
305	tolerance = np.min([x_tolerance, y_tolerance])	9✔
306	indexes_too_far = _indexes_filt(root_array, tolerance)	9✔
307
308	# Add more points into the root locus for points that are too far apart
309	while len(indexes_too_far) > 0 and kvect.size < 5000:	9✔
310	for counter, index in enumerate(indexes_too_far):	9✔
311	index = index + counter*3	9✔
312	new_gains = np.linspace(kvect[index], kvect[index + 1], 5)	9✔
313	new_points = _RLFindRoots(num, den, new_gains[1:4])	9✔
314	kvect = np.insert(kvect, index + 1, new_gains[1:4])	9✔
315	root_array = np.insert(root_array, index + 1, new_points, axis=0)	9✔
316
317	root_array = _RLSortRoots(root_array)	9✔
318	indexes_too_far = _indexes_filt(root_array, tolerance)	9✔
319
320	new_gains = kvect[-1] * np.hstack((np.logspace(0, 3, 4)))	9✔
321	new_points = _RLFindRoots(num, den, new_gains[1:4])	9✔
322	kvect = np.append(kvect, new_gains[1:4])	9✔
323	root_array = np.concatenate((root_array, new_points), axis=0)	9✔
324	root_array = _RLSortRoots(root_array)	9✔
325	return kvect, root_array, xlim, ylim	9✔
326
327
328	def _indexes_filt(root_array, tolerance):	9✔
329	"""Calculate the distance between points and return the indices.
330
331	Filter the indexes so only the resolution of points within the xlim and
332	ylim is improved when zoom is used.
333
334	"""
335	distance_points = np.abs(np.diff(root_array, axis=0))	9✔
336	indexes_too_far = list(np.unique(np.where(distance_points > tolerance)[0]))	9✔
337	indexes_too_far.sort()	9✔
338	return indexes_too_far	9✔
339
340
341	def _break_points(num, den):	9✔
342	"""Extract break points over real axis and gains given these locations"""
343	# type: (np.poly1d, np.poly1d) -> (np.array, np.array)
344	dnum = num.deriv(m=1)	9✔
345	dden = den.deriv(m=1)	9✔
346	polynom = den * dnum - num * dden	9✔
347	real_break_pts = polynom.r	9✔
348	# don't care about infinite break points
349	real_break_pts = real_break_pts[num(real_break_pts) != 0]	9✔
350	k_break = -den(real_break_pts) / num(real_break_pts)	9✔
351	idx = k_break >= 0 # only positives gains	9✔
352	k_break = k_break[idx]	9✔
353	real_break_pts = real_break_pts[idx]	9✔
354	if len(k_break) == 0:	9✔
355	k_break = [0]	9✔
356	real_break_pts = den.roots	9✔
357	return k_break, real_break_pts	9✔
358
359
360	def _ax_lim(root_array):	9✔
361	"""Utility to get the axis limits"""
362	axmin = np.min(np.real(root_array))	9✔
363	axmax = np.max(np.real(root_array))	9✔
364	if axmax != axmin:	9✔
365	deltax = (axmax - axmin) * 0.02	9✔
366	else:
367	deltax = np.max([1., axmax / 2])	9✔
368	axlim = [axmin - deltax, axmax + deltax]	9✔
369	return axlim	9✔
370
371
372	def _k_max(num, den, real_break_points, k_break_points):	9✔
373	""""Calculate the maximum gain for the root locus shown in the figure."""
374	asymp_number = den.order - num.order	9✔
375	singular_points = np.concatenate((num.roots, den.roots), axis=0)	9✔
376	important_points = np.concatenate(	9✔
377	(singular_points, real_break_points), axis=0)
378	false_gain = den.coeffs[0] / num.coeffs[0]	9✔
379
380	if asymp_number > 0:	9✔
381	asymp_center = (np.sum(den.roots) - np.sum(num.roots))/asymp_number	9✔
382	distance_max = 4 * np.max(np.abs(important_points - asymp_center))	9✔
383	asymp_angles = (2 * np.arange(0, asymp_number) - 1) \	9✔
384	* np.pi / asymp_number
385	if false_gain > 0:	9✔
386	# farthest points over asymptotes
387	farthest_points = asymp_center \	9✔
388	+ distance_max * np.exp(asymp_angles * 1j)
389	else:
390	asymp_angles = asymp_angles + np.pi	9✔
391	# farthest points over asymptotes
392	farthest_points = asymp_center \	9✔
393	+ distance_max * np.exp(asymp_angles * 1j)
394	kmax_asymp = np.real(np.abs(den(farthest_points)	9✔
395	/ num(farthest_points)))
396	else:
397	kmax_asymp = np.abs([np.abs(den.coeffs[0])	9✔
398	/ np.abs(num.coeffs[0]) * 3])
399
400	kmax = np.max(np.concatenate((np.real(kmax_asymp),	9✔
401	np.real(k_break_points)), axis=0))
402	if np.abs(false_gain) > kmax:	9✔
403	kmax = np.abs(false_gain)	9✔
404	return kmax	9✔
405
406
407	def _systopoly1d(sys):	9✔
408	"""Extract numerator and denominator polynomails for a system"""
409	# Allow inputs from the signal processing toolbox
410	if (isinstance(sys, scipy.signal.lti)):	9✔
411	nump = sys.num	×
412	denp = sys.den	×
413
414	else:
415	# Convert to a transfer function, if needed
416	sys = _convert_to_transfer_function(sys)	9✔
417
418	# Make sure we have a SISO system
419	if not sys.issiso():	9✔
420	raise ControlMIMONotImplemented()	×
421
422	# Start by extracting the numerator and denominator from system object
423	nump = sys.num[0][0]	9✔
424	denp = sys.den[0][0]	9✔
425
426	# Check to see if num, den are already polynomials; otherwise convert
427	if (not isinstance(nump, poly1d)):	9✔
428	nump = poly1d(nump)	9✔
429
430	if (not isinstance(denp, poly1d)):	9✔
431	denp = poly1d(denp)	9✔
432
433	return (nump, denp)	9✔
434
435
436	def _RLFindRoots(nump, denp, kvect):	9✔
437	"""Find the roots for the root locus."""
438	# Convert numerator and denominator to polynomials if they aren't
439	roots = []	9✔
440	for k in np.atleast_1d(kvect):	9✔
441	curpoly = denp + k * nump	9✔
442	curroots = curpoly.r	9✔
443	if len(curroots) < denp.order:	9✔
444	# if I have fewer poles than open loop, it is because i have
445	# one at infinity
446	curroots = np.append(curroots, np.inf)	×
447
448	curroots.sort()	9✔
449	roots.append(curroots)	9✔
450
451	return vstack(roots)	9✔
452
453
454	def _RLSortRoots(roots):	9✔
455	"""Sort the roots from _RLFindRoots, so that the root
456	locus doesn't show weird pseudo-branches as roots jump from
457	one branch to another."""
458
459	sorted = zeros_like(roots)	9✔
460	for n, row in enumerate(roots):	9✔
461	if n == 0:	9✔
462	sorted[n, :] = row	9✔
463	else:
464	# sort the current row by finding the element with the
465	# smallest absolute distance to each root in the
466	# previous row
467	available = list(range(len(prevrow)))	9✔
468	for elem in row:	9✔
469	evect = elem - prevrow[available]	9✔
470	ind1 = abs(evect).argmin()	9✔
471	ind = available.pop(ind1)	9✔
472	sorted[n, ind] = elem	9✔
473	prevrow = sorted[n, :]	9✔
474	return sorted	9✔
475
476
477	# Alternative ways to call these functions
478	root_locus = root_locus_plot	9✔
479	rlocus = root_locus_plot	9✔

python-control / python-control / 11994777565

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous