import scipy

from warnings import warn

from .freqplot import nyquist_plot

from .freqplot import nyquist_response

from . import config

__all__ = ['describing_function', 'describing_function_plot',

'describing_function_response', 'DescribingFunctionResponse',

'DescribingFunctionNonlinearity', 'friction_backlash_nonlinearity',

'relay_hysteresis_nonlinearity', 'saturation_nonlinearity']

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

return retdf

#

# Describing function response/plot

#

def describing_function_plot(

H, F, A, omega=None, refine=True, label="%5.2g @ %-5.2g",

warn=None, **kwargs):

"""Plot a Nyquist plot with a describing function for a nonlinear system.

# 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

static nonlinear feedback function. The DescribingFunctionResponse

class is used by the :func:`~control.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.

Attributes

----------

response : :class:`~control.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.

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 :func:`~control.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)

This function generates a Nyquist plot for a closed loop system consisting

of a linear system with a static nonlinear function in the feedback path.

# Compute the describing function response + intersections

def describing_function_response(

H, F, A, omega=None, refine=True, warn_nyquist=None,

plot=False, 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 static nonlinear

function in the feedback path.

Parameters

----------

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.

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.

warn : bool, optional

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.

Returns

-------

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.

response : :class:`~control.DescribingFunctionResponse` object

Response object that contains the result of the describing function

analysis. The following information can be retrieved from this

object:

response.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.

Examples

--------

>>> H_simple = ct.tf([8], [1, 2, 2, 1])

>>> F_saturation = ct.saturation_nonlinearity(1)

>>> amp = np.linspace(1, 4, 10)

>>> ct.describing_function_plot(H_simple, F_saturation, amp) # doctest: +SKIP

>>> response = ct.describing_function_response(H_simple, F_saturation, amp)

>>> response.intersections # doctest: +SKIP

[(3.343844998258643, 1.4142293090899216)]

>>> lines = response.plot()

"""

# Decide whether to turn on warnings or not

if warn is None:

if warn_nyquist is None:

# Turn warnings on unless omega was specified

warn = omega is None

warn_nyquist = omega is None

# Start by drawing a Nyquist curve

count, contour = nyquist_plot(

H, omega, plot=True, return_contour=True,

warn_encirclements=warn, warn_nyquist=warn, **kwargs)

H_omega, H_vals = contour.imag, H(contour)

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

# Now add the describing function curve to the plot

plt.plot(N_vals.real, N_vals.imag)

# Look for intersection points

intersections = []

positions, intersections = [], []

for i in range(N_vals.size - 1):

for j in range(H_vals.size - 1):

intersect = _find_intersection(

else:

a_final, omega_final = res.x[0], res.x[1]

# Add labels to the intersection points

if isinstance(label, str):

pos = H(1j * omega_final)

plt.text(pos.real, pos.imag, label % (a_final, omega_final))

elif label is not None or label is not False:

raise ValueError("label must be formatting string or None")

pos = H(1j * omega_final)

# Save the final estimate

positions.append(pos)

intersections.append((a_final, omega_final))

return intersections

return DescribingFunctionResponse(

response, N_vals, positions, intersections)

def describing_function_plot(

*sysdata, label="%5.2g @ %-5.2g", **kwargs):

"""describing_function_plot(data, *args, **kwargs)

Plot a Nyquist plot with a describing function for a nonlinear system.

This function generates a Nyquist plot for a closed loop system

consisting of a linear system with a static nonlinear function 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

:func:`~control.describing_function_response` function. In the second

form, that function is called internally, with the listed arguments.

Parameters

----------

data : :class:`~control.DescribingFunctionData`

A describing function response data object created by

:func:`~control.describing_function_response`.

H : LTI system

Linear time-invariant (LTI) system (state space, transfer function, or

FRD)

F : static nonlinear function

A static 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. 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, optional

Formatting string used to label intersection points on the Nyquist

plot. Defaults to "%5.2g @ %-5.2g". Set to `None` to omit labels.

Returns

-------

lines : 1D array of Line2D

Arrray of Line2D objects for 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 styem. The second element of

the array is a list of lines (typically only one) for the

describing function curve.

Examples

--------

>>> H_simple = ct.tf([8], [1, 2, 2, 1])

>>> F_saturation = ct.saturation_nonlinearity(1)

>>> amp = np.linspace(1, 4, 10)

>>> lines = 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))

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:

dfresp = sysdata[0]

else:

raise TypeError("1, 3, or 4 position arguments required")

# Create a list of lines for the output

out = np.empty(2, dtype=object)

# Plot the Nyquist response

out[0] = dfresp.response.plot(**kwargs)[0]

# Add the describing function curve to the plot

lines = plt.plot(dfresp.N_vals.real, dfresp.N_vals.imag)

out[1] = lines

# Label the intersection points

if label:

for pos, (a, omega) in zip(dfresp.positions, dfresp.intersections):

# Add labels to the intersection points

plt.text(pos.real, pos.imag, label % (a, omega))

return out

# Utility function to figure out whether two line segments intersection