# 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