16081041622

Committed 04 Jul 2025 09:18PM UTC coverage: 94.733% (-0.01%) from 94.745%

Build # 16081041622

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1155 from murrayrm/fix_nyquist_rescaling-24Mar2025

Update Nyquist rescaling + other improvements

Run Details

9946 of 10499 relevant lines covered (94.73%)

8.28 hits per line

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

97.99

control/descfcn.py

# descfcn.py - describing function analysis
# RMM, 23 Jan 2021

"""This module contains functions for performing closed loop analysis of
systems with memoryless nonlinearities using describing function analysis.

"""

import math
from warnings import warn

import numpy as np
import scipy

from . import config
from .ctrlplot import ControlPlot
from .freqplot import nyquist_response

__all__ = ['describing_function', 'describing_function_plot',
           'describing_function_response', 'DescribingFunctionResponse',
           'DescribingFunctionNonlinearity', 'friction_backlash_nonlinearity',
           'relay_hysteresis_nonlinearity', 'saturation_nonlinearity']

# Class for nonlinearities with a built-in describing function
class DescribingFunctionNonlinearity():
    """Base class for nonlinear systems with a describing function.

    This class is intended to be used as a base class for nonlinear functions
    that have an analytically defined describing function.  Subclasses should
    override the `__call__` and `describing_function` methods and (optionally)
    the `_isstatic` method (should be False if `__call__` updates the
    instance state).

    """
    def __init__(self):
        """Initialize a describing function nonlinearity (optional)."""
        pass

    def __call__(self, A):
        """Evaluate the nonlinearity at a (scalar) input value."""
        raise NotImplementedError(
            "__call__() not implemented for this function (internal error)")

    def describing_function(self, A):
        """Return the describing function for a nonlinearity.

        This method is used to allow analytical representations of the
        describing function for a nonlinearity.  It turns the (complex) value
        of the describing function for sinusoidal input of amplitude `A`.

        Parameters
        ----------
        A : float
            Amplitude of the sinusoidal input to the nonlinearity.

        Returns
        -------
        float
            Value of the describing function at the given amplitude.

        """
        raise NotImplementedError(
            "describing function not implemented for this function")

    def _isstatic(self):
        """Return True if the function has no internal state (memoryless).

        This internal function is used to optimize numerical computation of
        the describing function.  It can be set to True if the instance
        maintains no internal memory of the instance state.  Assumed False by
        default.

        """
        return False

    # Utility function used to compute common describing functions
    def _f(self, x):
        return math.copysign(1, x) if abs(x) > 1 else \
            (math.asin(x) + x * math.sqrt(1 - x**2)) * 2 / math.pi


def describing_function(
        F, A, num_points=100, zero_check=True, try_method=True):
    """Numerically compute describing function of a nonlinear function.

    The describing function of a nonlinearity is given by magnitude and phase
    of the first harmonic of the function when evaluated along a sinusoidal
    input :math:`A \\sin \\omega t`.  This function returns the magnitude and
    phase of the describing function at amplitude :math:`A`.

    Parameters
    ----------
    F : callable
        The function F() should accept a scalar number as an argument and
        return a scalar number.  For compatibility with (static) nonlinear
        input/output systems, the output can also return a 1D array with a
        single element.

        If the function is an object with a method `describing_function`
        then this method will be used to computing the describing function
        instead of a nonlinear computation.  Some common nonlinearities
        use the `DescribingFunctionNonlinearity` class,
        which provides this functionality.

    A : array_like
        The amplitude(s) at which the describing function should be calculated.

    num_points : int, optional
        Number of points to use in computing describing function (default =
        100).

    zero_check : bool, optional
        If True (default) then `A` is zero, the function will be evaluated
        and checked to make sure it is zero.  If not, a `TypeError` exception
        is raised.  If zero_check is False, no check is made on the value of
        the function at zero.

    try_method : bool, optional
        If True (default), check the `F` argument to see if it is an object
        with a `describing_function` method and use this to compute the
        describing function.  More information in the `describing_function`
        method for the `DescribingFunctionNonlinearity` class.

    Returns
    -------
    df : ndarray of complex
        The (complex) value of the describing function at the given amplitudes.

    Raises
    ------
    TypeError
        If A[i] < 0 or if A[i] = 0 and the function F(0) is non-zero.

    Examples
    --------
    >>> F = lambda x: np.exp(-x)      # Basic diode description
    >>> A = np.logspace(-1, 1, 20)    # Amplitudes from 0.1 to 10.0
    >>> df_values = ct.describing_function(F, A)
    >>> len(df_values)
    20

    """
    # If there is an analytical solution, trying using that first
    if try_method and hasattr(F, 'describing_function'):
        try:
            return np.vectorize(F.describing_function, otypes=[complex])(A)
        except NotImplementedError:
            # Drop through and do the numerical computation
            pass

    #
    # The describing function of a nonlinear function F() can be computed by
    # evaluating the nonlinearity over a sinusoid.  The Fourier series for a
    # nonlinear function evaluated on a sinusoid can be written as
    #
    # F(A\sin\omega t) = \sum_{k=1}^\infty M_k(A) \sin(k\omega t + \phi_k(A))
    #
    # The describing function is given by the complex number
    #
    #    N(A) = M_1(A) e^{j \phi_1(A)} / A
    #
    # To compute this, we compute F(A \sin\theta) for \theta between 0 and 2
    # \pi, use the identities
    #
    #   \sin(\theta + \phi) = \sin\theta \cos\phi + \cos\theta \sin\phi
    #   \int_0^{2\pi} \sin^2 \theta d\theta = \pi
    #   \int_0^{2\pi} \cos^2 \theta d\theta = \pi
    #
    # and then integrate the product against \sin\theta and \cos\theta to obtain
    #
    #   \int_0^{2\pi} F(A\sin\theta) \sin\theta d\theta = M_1 \pi \cos\phi
    #   \int_0^{2\pi} F(A\sin\theta) \cos\theta d\theta = M_1 \pi \sin\phi
    #
    # From these we can compute M1 and \phi.
    #

    # Evaluate over a full range of angles (leave off endpoint a la DFT)
    theta, dtheta = np.linspace(
        0, 2*np.pi, num_points, endpoint=False, retstep=True)
    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)

    # See if this is a static nonlinearity (assume not, just in case)
    if not hasattr(F, '_isstatic') or not F._isstatic():
        # Initialize any internal state by going through an initial cycle
        for x in np.atleast_1d(A).min() * sin_theta:
            F(x)                # ignore the result

    # Go through all of the amplitudes we were given
    retdf = np.empty(np.shape(A), dtype=complex)
    df = retdf                  # Access to the return array
    df.shape = (-1, )           # as a 1D array
    for i, a in enumerate(np.atleast_1d(A)):
        # Make sure we got a valid argument
        if a == 0:
            # Check to make sure the function has zero output with zero input
            if zero_check and np.squeeze(F(0.)) != 0:
                raise ValueError("function must evaluate to zero at zero")
            df[i] = 1.
            continue
        elif a < 0:
            raise ValueError("cannot evaluate describing function for A < 0")

        # Save the scaling factor to make the formulas simpler
        scale = dtheta / np.pi / a

        # Evaluate the function along a sinusoid
        F_eval = np.array([F(x) for x in a*sin_theta]).squeeze()

        # Compute the projections onto sine and cosine
        df_real = (F_eval @ sin_theta) * scale     # = M_1 \cos\phi / a
        df_imag = (F_eval @ cos_theta) * scale     # = M_1 \sin\phi / a

        df[i] = df_real + 1j * df_imag

    # Return the values in the same shape as they were requested
    return retdf

#
# Describing function response/plot
#

# Simple class to store the describing function response
class DescribingFunctionResponse:
    """Results of describing function analysis.

    Describing functions allow analysis of a linear I/O systems with a
    nonlinear feedback function.  The DescribingFunctionResponse class
    is used by the `describing_function_response` function to return
    the results of a describing function analysis.  The response
    object can be used to obtain information about the describing
    function analysis or generate a Nyquist plot showing the frequency
    response of the linear systems and the describing function for the
    nonlinear element.

    Parameters
    ----------
    response : `FrequencyResponseData`
        Frequency response of the linear system component of the system.
    intersections : 1D array of 2-tuples or None
        A list of all amplitudes and frequencies in which
        :math:`H(j\\omega) N(A) = -1`, where :math:`N(A)` is the describing
        function associated with `F`, or None if there are no such
        points.  Each pair represents a potential limit cycle for the
        closed loop system with amplitude given by the first value of the
        tuple and frequency given by the second value.
    N_vals : complex array
        Complex value of the describing function, indexed by amplitude.
    positions : list of complex
        Location of the intersections in the complex plane.

    """
    def __init__(self, response, N_vals, positions, intersections):
        """Create a describing function response data object."""
        self.response = response
        self.N_vals = N_vals
        self.positions = positions
        self.intersections = intersections

    def plot(self, **kwargs):
        """Plot the results of a describing function analysis.

        See `describing_function_plot` for details.
        """
        return describing_function_plot(self, **kwargs)

    # Implement iter, getitem, len to allow recovering the intersections
    def __iter__(self):
        return iter(self.intersections)

    def __getitem__(self, index):
        return list(self.__iter__())[index]

    def __len__(self):
        return len(self.intersections)


# Compute the describing function response + intersections
def describing_function_response(
        H, F, A, omega=None, refine=True, warn_nyquist=None,
        _check_kwargs=True, **kwargs):
    """Compute the describing function response of a system.

    This function uses describing function analysis to analyze a closed
    loop system consisting of a linear system with a nonlinear function in
    the feedback path.

    Parameters
    ----------
    H : LTI system
        Linear time-invariant (LTI) system (state space, transfer function,
        or FRD).
    F : nonlinear function
        Feedback nonlinearity, either a scalar function or a single-input,
        single-output, static input/output system.
    A : list
        List of amplitudes to be used for the describing function plot.
    omega : list, optional
        List of frequencies to be used for the linear system Nyquist curve.
    warn_nyquist : bool, optional
        Set to True to turn on warnings generated by `nyquist_plot` or
        False to turn off warnings.  If not set (or set to None),
        warnings are turned off if omega is specified, otherwise they are
        turned on.
    refine : bool, optional
        If True, `scipy.optimize.minimize` to refine the estimate
        of the intersection of the frequency response and the describing
        function.

    Returns
    -------
    response : `DescribingFunctionResponse` object
        Response object that contains the result of the describing function
        analysis.  The results can plotted using the
        `~DescribingFunctionResponse.plot` method.
    response.intersections : 1D ndarray of 2-tuples or None
        A list of all amplitudes and frequencies in which
        :math:`H(j\\omega) N(a) = -1`, where :math:`N(a)` is the describing
        function associated with `F`, or None if there are no such
        points.  Each pair represents a potential limit cycle for the
        closed loop system with amplitude given by the first value of the
        tuple and frequency given by the second value.
    response.Nvals : complex ndarray
        Complex value of the describing function, indexed by amplitude.

    See Also
    --------
    DescribingFunctionResponse, describing_function_plot

    Examples
    --------
    >>> H_simple = ct.tf([8], [1, 2, 2, 1])
    >>> F_saturation = ct.saturation_nonlinearity(1)
    >>> amp = np.linspace(1, 4, 10)
    >>> response = ct.describing_function_response(H_simple, F_saturation, amp)
    >>> response.intersections  # doctest: +SKIP
    [(3.343844998258643, 1.4142293090899216)]
    >>> cplt = response.plot()

    """
    # Decide whether to turn on warnings or not
    if warn_nyquist is None:
        # Turn warnings on unless omega was specified
        warn_nyquist = omega is None

    # Start by drawing a Nyquist curve
    response = nyquist_response(
        H, omega, warn_encirclements=warn_nyquist, warn_nyquist=warn_nyquist,
        _check_kwargs=_check_kwargs, **kwargs)
    H_omega, H_vals = response.contour.imag, H(response.contour)

    # Compute the describing function
    df = describing_function(F, A)
    N_vals = -1/df

    # Look for intersection points
    positions, intersections = [], []
    for i in range(N_vals.size - 1):
        for j in range(H_vals.size - 1):
            intersect = _find_intersection(
                N_vals[i], N_vals[i+1], H_vals[j], H_vals[j+1])
            if intersect == None:
                continue

            # Found an intersection, compute a and omega
            s_amp, s_omega = intersect
            a_guess = (1 - s_amp) * A[i] + s_amp * A[i+1]
            omega_guess = (1 - s_omega) * H_omega[j] + s_omega * H_omega[j+1]

            # Refine the coarse estimate to get better intersection point
            a_final, omega_final = a_guess, omega_guess
            if refine:
                # Refine the answer to get more accuracy
                def _cost(x):
                    # If arguments are invalid, return a "large" value
                    # Note: imposing bounds messed up the optimization (?)
                    if x[0] < 0 or x[1] < 0:
                        return 1
                    return abs(1 + H(1j * x[1]) *
                               describing_function(F, x[0]))**2
                res = scipy.optimize.minimize(
                    _cost, [a_guess, omega_guess])
                # bounds=[(A[i], A[i+1]), (H_omega[j], H_omega[j+1])])

                if not res.success:
                    warn("not able to refine result; returning estimate")
                else:
                    a_final, omega_final = res.x[0], res.x[1]

            pos = H(1j * omega_final)

            # Save the final estimate
            positions.append(pos)
            intersections.append((a_final, omega_final))

    return DescribingFunctionResponse(
        response, N_vals, positions, intersections)


def describing_function_plot(
        *sysdata, point_label="%5.2g @ %-5.2g", label=None, **kwargs):
    """describing_function_plot(data, *args, **kwargs)

    Nyquist plot with describing function for a nonlinear system.

    This function generates a Nyquist plot for a closed loop system
    consisting of a linear system with a nonlinearity in the feedback path.

    The function may be called in one of two forms:

        describing_function_plot(response[, options])

        describing_function_plot(H, F, A[, omega[, options]])

    In the first form, the response should be generated using the
    `describing_function_response` function.  In the second
    form, that function is called internally, with the listed arguments.

    Parameters
    ----------
    data : `DescribingFunctionResponse`
        A describing function response data object created by
        `describing_function_response`.
    H : LTI system
        Linear time-invariant (LTI) system (state space, transfer function,
        or FRD).
    F : nonlinear function
        Nonlinearity in the feedback path, either a scalar function or a
        single-input, single-output, static input/output system.
    A : list
        List of amplitudes to be used for the describing function plot.
    omega : list, optional
        List of frequencies to be used for the linear system Nyquist
        curve. If not specified (or None), frequencies are computed
        automatically based on the properties of the linear system.
    refine : bool, optional
        If True (default), refine the location of the intersection of the
        Nyquist curve for the linear system and the describing function to
        determine the intersection point.
    label : str or array_like of str, optional
        If present, replace automatically generated label with the given label.
    point_label : str, optional
        Formatting string used to label intersection points on the Nyquist
        plot.  Defaults to "%5.2g @ %-5.2g".  Set to None to omit labels.
    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.
    title : str, optional
        Set the title of the plot.  Defaults to plot type and system name(s).
    warn_nyquist : bool, optional
        Set to True to turn on warnings generated by `nyquist_plot` or
        False to turn off warnings.  If not set (or set to None),
        warnings are turned off if omega is specified, otherwise they are
        turned on.
    **kwargs : `matplotlib.pyplot.plot` keyword properties, optional
        Additional keywords passed to `matplotlib` to specify line properties
        for Nyquist curve.

    Returns
    -------
    cplt : `ControlPlot` object
        Object containing the data that were plotted.  See `ControlPlot`
        for more detailed information.
    cplt.lines : array of `matplotlib.lines.Line2D`
        Array containing information on each line in the plot.  The first
        element of the array is a list of lines (typically only one) for
        the Nyquist plot of the linear I/O system.  The second element of
        the array is a list of lines (typically only one) for the
        describing function curve.
    cplt.axes : 2D array of `matplotlib.axes.Axes`
        Axes for each subplot.
    cplt.figure : `matplotlib.figure.Figure`
        Figure containing the plot.

    See Also
    --------
    DescribingFunctionResponse, describing_function_response

    Examples
    --------
    >>> H_simple = ct.tf([8], [1, 2, 2, 1])
    >>> F_saturation = ct.saturation_nonlinearity(1)
    >>> amp = np.linspace(1, 4, 10)
    >>> cplt = ct.describing_function_plot(H_simple, F_saturation, amp)

    """
    # Process keywords
    warn_nyquist = config._process_legacy_keyword(
        kwargs, 'warn', 'warn_nyquist', kwargs.pop('warn_nyquist', None))
    point_label = config._process_legacy_keyword(
        kwargs, 'label', 'point_label', point_label)

    # TODO: update to be consistent with ctrlplot use of `label`
    if label not in (False, None) and not isinstance(label, str):
        raise ValueError("label must be formatting string, False, or None")

    # Get the describing function response
    if len(sysdata) == 3:
        sysdata = sysdata + (None, )    # set omega to default value
    if len(sysdata) == 4:
        dfresp = describing_function_response(
            *sysdata, refine=kwargs.pop('refine', True),
            warn_nyquist=warn_nyquist)
    elif len(sysdata) == 1:
        if not isinstance(sysdata[0], DescribingFunctionResponse):
            raise TypeError("data must be DescribingFunctionResponse")
        else:
            dfresp = sysdata[0]
    else:
        raise TypeError("1, 3, or 4 position arguments required")

    # Don't allow legend keyword arguments
    for kw in ['legend_loc', 'legend_map', 'show_legend']:
        if kw in kwargs:
            raise TypeError(f"unexpected keyword argument '{kw}'")

    # Create a list of lines for the output
    lines = np.empty(2, dtype=object)

    # Plot the Nyquist response
    cplt = dfresp.response.plot(**kwargs)
    ax = cplt.axes[0, 0]        # Get the axes where the plot was made
    lines[0] = np.concatenate(  # Return Nyquist lines for first system
        cplt.lines.flatten()).tolist()

    # Add the describing function curve to the plot
    lines[1] = ax.plot(dfresp.N_vals.real, dfresp.N_vals.imag)

    # Label the intersection points
    if point_label:
        for pos, (a, omega) in zip(dfresp.positions, dfresp.intersections):
            # Add labels to the intersection points
            ax.text(pos.real, pos.imag, point_label % (a, omega))

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


# Utility function to figure out whether two line segments intersection
def _find_intersection(L1a, L1b, L2a, L2b):
    # Compute the tangents for the segments
    L1t = L1b - L1a
    L2t = L2b - L2a

    # Set up components of the solution: b = M s
    b = L1a - L2a
    detM = L1t.imag * L2t.real - L1t.real * L2t.imag
    if abs(detM) < 1e-8:        # TODO: fix magic number
        return None

    # Solve for the intersection points on each line segment
    s1 = (L2t.imag * b.real - L2t.real * b.imag) / detM
    if s1 < 0 or s1 > 1:
        return None

    s2 = (L1t.imag * b.real - L1t.real * b.imag) / detM
    if s2 < 0 or s2 > 1:
        return None

    # Debugging test
    # np.testing.assert_almost_equal(L1a + s1 * L1t, L2a + s2 * L2t)

    # Intersection is within segments; return proportional distance
    return (s1, s2)


# Saturation nonlinearity
class saturation_nonlinearity(DescribingFunctionNonlinearity):
    """Saturation nonlinearity for describing function analysis.

    This class creates a nonlinear function representing a saturation with
    given upper and lower bounds, including the describing function for the
    nonlinearity.  The following call creates a nonlinear function suitable
    for describing function analysis:

        F = saturation_nonlinearity(ub[, lb])

    By default, the lower bound is set to the negative of the upper bound.
    Asymmetric saturation functions can be created, but note that these
    functions will not have zero bias and hence care must be taken in using
    the nonlinearity for analysis.

    Parameters
    ----------
    lb, ub : float
        Upper and lower saturation bounds.

    Examples
    --------
    >>> nl = ct.saturation_nonlinearity(5)
    >>> nl(1)
    np.int64(1)
    >>> nl(10)
    np.int64(5)
    >>> nl(-10)
    np.int64(-5)

    """
    def __init__(self, ub=1, lb=None):
        # Create the describing function nonlinearity object
        super(saturation_nonlinearity, self).__init__()

        # Process arguments
        if lb == None:
            # Only received one argument; assume symmetric around zero
            lb, ub = -abs(ub), abs(ub)

        # Make sure the bounds are sensible
        if lb > 0 or ub < 0 or lb + ub != 0:
            warn("asymmetric saturation; ignoring non-zero bias term")

        self.lb = lb
        self.ub = ub

    def __call__(self, x):
        return np.clip(x, self.lb, self.ub)

    def _isstatic(self):
        return True

    def describing_function(self, A):
        """Return the describing function for a saturation nonlinearity.

        Parameters
        ----------
        A : float
            Amplitude of the sinusoidal input to the nonlinearity.

        Returns
        -------
        float
            Value of the describing function at the given amplitude.

        """
        # Check to make sure the amplitude is positive
        if A < 0:
            raise ValueError("cannot evaluate describing function for A < 0")

        if self.lb <= A and A <= self.ub:
            return 1.
        else:
            alpha, beta = math.asin(self.ub/A), math.asin(-self.lb/A)
            return (math.sin(alpha + beta) * math.cos(alpha - beta) +
                    (alpha + beta)) / math.pi


# Relay with hysteresis (FBS2e, Example 10.12)
class relay_hysteresis_nonlinearity(DescribingFunctionNonlinearity):
    """Relay w/ hysteresis for describing function analysis.

    This class creates a nonlinear function representing a a relay with
    symmetric upper and lower bounds of magnitude `b` and a hysteretic region
    of width `c` (using the notation from [FBS2e](https://fbsbook.org),
    Example 10.12, including the describing function for the nonlinearity.
    The following call creates a nonlinear function suitable for describing
    function analysis::

        F = relay_hysteresis_nonlinearity(b, c)

    The output of this function is b if x > c and -b if x < -c.  For -c <=
    x <= c, the value depends on the branch of the hysteresis loop (as
    illustrated in Figure 10.20 of FBS2e).

    Parameters
    ----------
    b : float
        Hysteresis bound.
    c : float
        Width of hysteresis region.

    Examples
    --------
    >>> nl = ct.relay_hysteresis_nonlinearity(1, 2)
    >>> nl(0)
    -1
    >>> nl(1)  # not enough for switching on
    -1
    >>> nl(5)
    1
    >>> nl(-1)  # not enough for switching off
    1
    >>> nl(-5)
    -1

    """
    def __init__(self, b, c):
        # Create the describing function nonlinearity object
        super(relay_hysteresis_nonlinearity, self).__init__()

        # Initialize the state to bottom branch
        self.branch = -1        # lower branch
        self.b = b              # relay output value
        self.c = c              # size of hysteresis region

    def __call__(self, x):
        if x > self.c:
            y = self.b
            self.branch = 1
        elif x < -self.c:
            y = -self.b
            self.branch = -1
        elif self.branch == -1:
            y = -self.b
        elif self.branch == 1:
            y = self.b
        return y

    def _isstatic(self):
        return False

    def describing_function(self, A):
        """Return the describing function for a hysteresis nonlinearity.

        Parameters
        ----------
        A : float
            Amplitude of the sinusoidal input to the nonlinearity.

        Returns
        -------
        float
            Value of the describing function at the given amplitude.

        """
        # Check to make sure the amplitude is positive
        if A < 0:
            raise ValueError("cannot evaluate describing function for A < 0")

        if A < self.c:
            return np.nan

        df_real = 4 * self.b * math.sqrt(1 - (self.c/A)**2) / (A * math.pi)
        df_imag = -4 * self.b * self.c / (math.pi * A**2)
        return df_real + 1j * df_imag


# Friction-dominated backlash nonlinearity (#48 in Gelb and Vander Velde, 1968)
class friction_backlash_nonlinearity(DescribingFunctionNonlinearity):
    """Backlash nonlinearity for describing function analysis.

    This class creates a nonlinear function representing a friction-dominated
    backlash nonlinearity ,including the describing function for the
    nonlinearity.  The following call creates a nonlinear function suitable
    for describing function analysis::

        F = friction_backlash_nonlinearity(b)

    This function maintains an internal state representing the 'center' of
    a mechanism with backlash.  If the new input is within b/2 of the
    current center, the output is unchanged.  Otherwise, the output is
    given by the input shifted by b/2.

    Parameters
    ----------
    b : float
        Backlash amount.

    Examples
    --------
    >>> nl = ct.friction_backlash_nonlinearity(2)  # backlash of +/- 1
    >>> nl(0)
    0
    >>> nl(1)  # not enough to overcome backlash
    0
    >>> nl(2)
    1.0
    >>> nl(1)
    1.0
    >>> nl(0)  # not enough to overcome backlash
    1.0
    >>> nl(-1)
    0.0

    """

    def __init__(self, b):
        # Create the describing function nonlinearity object
        super(friction_backlash_nonlinearity, self).__init__()

        self.b = b              # backlash distance
        self.center = 0         # current center position

    def __call__(self, x):
        # If we are outside the backlash, move and shift the center
        if x - self.center > self.b/2:
            self.center = x - self.b/2
        elif x - self.center < -self.b/2:
            self.center = x + self.b/2
        return self.center

    def _isstatic(self):
        return False

    def describing_function(self, A):
        """Return the describing function for a backlash nonlinearity.

        Parameters
        ----------
        A : float
            Amplitude of the sinusoidal input to the nonlinearity.

        Returns
        -------
        float
            Value of the describing function at the given amplitude.

        """
        # Check to make sure the amplitude is positive
        if A < 0:
            raise ValueError("cannot evaluate describing function for A < 0")

        if A <= self.b/2:
            return 0

        df_real = (1 + self._f(1 - self.b/A)) / 2
        df_imag = -(2 * self.b/A - (self.b/A)**2) / math.pi
        return df_real + 1j * df_imag

python-control / python-control / 16081041622

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