Describing function analysis

Richard M. Murray, 27 Jan 2021

This Jupyter notebook shows how to use the descfcn module of the Python Control Toolbox to perform describing function analysis of a nonlinear system. A brief introduction to describing functions can be found in Feedback Systems, Section 10.5 (Generalized Notions of Gain and Phase).

[1]:
import control as ct
import numpy as np
import matplotlib.pyplot as plt
import math

Built-in describing functions

The Python Control Toobox has a number of built-in functions that provide describing functions for some standard nonlinearities.

Saturation nonlinearity

A saturation nonlinearity can be obtained using the ct.saturation_nonlinearity function. This function takes the saturation level as an argument.

[2]:
saturation=ct.saturation_nonlinearity(0.75)
x = np.linspace(-2, 2, 50)
plt.plot(x, saturation(x))
plt.xlabel("Input, x")
plt.ylabel("Output, y = sat(x)")
plt.title("Input/output map for a saturation nonlinearity");
_images/describing_functions_4_0.png
[3]:
amp_range = np.linspace(0, 2, 50)
plt.plot(amp_range, ct.describing_function(saturation, amp_range))
plt.xlabel("Amplitude A")
plt.ylabel("Describing function, N(A)")
plt.title("Describing function for a saturation nonlinearity");
/Users/murray/anaconda3/envs/python3.10-slycot/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)
_images/describing_functions_5_1.png

Backlash nonlinearity

A friction-dominated backlash nonlinearity can be obtained using the ct.friction_backlash_nonlinearity function. This function takes as is argument the size of the backlash region.

[4]:
backlash = ct.friction_backlash_nonlinearity(0.5)
theta = np.linspace(0, 2*np.pi, 50)
x = np.sin(theta)
plt.plot(x, [backlash(z) for z in x])
plt.xlabel("Input, x")
plt.ylabel("Output, y = backlash(x)")
plt.title("Input/output map for a friction-dominated backlash nonlinearity");
_images/describing_functions_7_0.png
[5]:
amp_range = np.linspace(0, 2, 50)
N_a = ct.describing_function(backlash, amp_range)

plt.figure()
plt.plot(amp_range, abs(N_a))
plt.xlabel("Amplitude A")
plt.ylabel("Amplitude of describing function, N(A)")
plt.title("Describing function for a backlash nonlinearity")

plt.figure()
plt.plot(amp_range, np.angle(N_a))
plt.xlabel("Amplitude A")
plt.ylabel("Phase of describing function, N(A)")
plt.title("Describing function for a backlash nonlinearity");
_images/describing_functions_8_0.png
_images/describing_functions_8_1.png

User-defined, static nonlinearities

In addition to pre-defined nonlinearies, it is possible to computing describing functions for static nonlinearities. The describing function for any suitable nonlinear function can be computed numerically using the describing_function function.

[6]:
# Define a saturation nonlinearity as a simple function
def my_saturation(x):
    if abs(x) >= 1:
        return math.copysign(1, x)
    else:
        return x

amp_range = np.linspace(0, 2, 50)
plt.plot(amp_range, ct.describing_function(my_saturation, amp_range).real)
plt.xlabel("Amplitude A")
plt.ylabel("Describing function, N(A)")
plt.title("Describing function for a saturation nonlinearity");
_images/describing_functions_10_0.png

Stability analysis using describing functions

Describing functions can be used to assess stability of closed loop systems consisting of a linear system and a static nonlinear using a Nyquist plot.

Limit cycle position for a third order system with saturation nonlinearity

Consider a nonlinear feedback system consisting of a third-order linear system with transfer function H(s) and a saturation nonlinearity having describing function N(a). Stability can be assessed by looking for points at which

H(j\omega) N(a) = -1

The describing_function_plot function plots H(j\omega) and -1/N(a) and prints out the the amplitudes and frequencies corresponding to intersections of these curves.

[7]:
# Linear dynamics
H_simple = ct.tf([8], [1, 2, 2, 1])
omega = np.logspace(-3, 3, 500)

# Nonlinearity
F_saturation = ct.saturation_nonlinearity(1)
amp = np.linspace(00, 5, 50)

# Describing function plot (return value = amp, freq)
ct.describing_function_plot(H_simple, F_saturation, amp, omega)
[7]:
[(3.343977839541308, 1.4142156916816762)]
_images/describing_functions_13_1.png

The intersection occurs at amplitude 3.3 and frequency 1.4 rad/sec (= 0.2 Hz) and thus we predict a limit cycle with amplitude 3.3 and period of approximately 5 seconds.

[8]:
# Create an I/O system simulation to see what happens
io_saturation = ct.NonlinearIOSystem(
    None,
    lambda t, x, u, params: F_saturation(u),
    inputs=1, outputs=1
)

sys = ct.feedback(ct.ss(H_simple), io_saturation)
T = np.linspace(0, 30, 200)
t, y = ct.input_output_response(sys, T, 0.1, 0)
plt.plot(t, y);
_images/describing_functions_15_0.png

Limit cycle prediction with for a time-delay system with backlash

This example demonstrates a more complicated interaction between a (non-static) nonlinearity and a higher order transfer function, resulting in multiple intersection points.

[9]:
# Linear dynamics
H_simple = ct.tf([1], [1, 2, 2, 1])
H_multiple = H_simple * ct.tf(*ct.pade(5, 4)) * 4
omega = np.logspace(-3, 3, 500)

# Nonlinearity
F_backlash = ct.friction_backlash_nonlinearity(1)
amp = np.linspace(0.6, 5, 50)

# Describing function plot
ct.describing_function_plot(H_multiple, F_backlash, amp, omega, mirror_style=False)
[9]:
[(0.6260158833534124, 0.3102619497970334),
 (0.8741930326860968, 1.2156410944770426)]
_images/describing_functions_17_1.png
[ ]: