add some features for plotting describing functions

docs/src/lib/nonlinear.md

@@ -263,6 +263,45 @@ p2 = plot(t, [res.u[:] res.y[:]],
 plot(p1, p2, layout=(1,2), size=(900, 400))
 ```
 
+### Describing Function Analysis
+
+The [describing function](https://en.wikipedia.org/wiki/Describing_function) is a quasi-linearization technique for predicting limit cycles in nonlinear feedback systems. Given a nonlinearity ``f(x)`` and an input ``A\sin(\theta)``, the describing function ``N(A)`` is the complex gain of the fundamental Fourier component of the output.
+
+A limit cycle at amplitude ``A`` and frequency ``\omega`` is predicted when the Nyquist curve of the linear part ``L(j\omega)`` intersects ``-1/N(A)``.
+
+Analytical formulas are provided for [`Saturation`](@ref ControlSystemsBase.saturation), [`DeadZone`](@ref ControlSystemsBase.deadzone), and [`Hysteresis`](@ref ControlSystemsBase.Hysteresis). For other nonlinearities, a generic numerical method is used.
+
+Please note that the describing function is an approximation that only considers the fundamental harmonic of the output, and thus may not always give accurate predictions, especially for strongly nonlinear systems or systems that are not of lowpass character.
+
+```@example DF
+using ControlSystemsBase
+using ControlSystemsBase: describing_function, Saturation, DeadZone, Hysteresis
+# Analytical describing function for saturation
+A = 2.0
+describing_function(Saturation(1.0), A)
+```
+
+```@example DF
+# Numerical describing function for an arbitrary nonlinearity
+describing_function(x -> clamp(x, -1, 1), A)
+```
+
+```@example DF
+# Describing function for dead-zone
+describing_function(DeadZone(0.5), A)
+```
+
+To perform a graphical limit-cycle analysis, use [`ControlSystemsBase.describing_function_plot`](@ref) which overlays ``-1/N(A)`` on the Nyquist plot:
+```@example DF
+using Plots
+using ControlSystemsBase: describing_function_plot
+s = tf("s")
+G = 10 / (s^3 + 2s^2 + s + 1)
+
+describing_function_plot(G, Saturation(1.0); A_range=0.01:0.01:100)
+```
+The intersection of the Nyquist curve and the ``-1/N(A)`` curve indicates a potential limit cycle. The amplitude of the limit cycle can be estimated from the corresponding ``A`` value.
+
 ## Limitations
 - Remember, this functionality is experimental and subject to breakage.
 - Currently only `Continuous` systems supported.
@@ -294,4 +333,7 @@ ControlSystemsBase.ratelimit
 ControlSystemsBase.deadzone
 ControlSystemsBase.linearize
 ControlSystemsBase.hysteresis
+ControlSystemsBase.Hysteresis
+ControlSystemsBase.describing_function
+ControlSystemsBase.describing_function_plot
 ```

example/describing_function.jl

@@ -1,101 +1,66 @@
 #=
 This example demonstrates how to compute and plot describing functions in the Nyquist plane
+using the built-in `describing_function` from ControlSystems.jl.
 =#
 using ControlSystems
-using QuadGK
-
-"""
-    compute_describing_function(f, A)
-
-Numerically computes the describing function N(A) for a nonlinearity f(x)
-at input amplitude A. Returns a Complex number.
-"""
-function compute_describing_function(f, A)
-    # Fundamental sine component coefficient (b1)
-    # b1 = (1/π) * ∫ [f(A*sin(θ)) * sin(θ)] dθ from 0 to 2π
-    b1_integral, _ = quadgk(θ -> f(A * sin(θ)) * sin(θ), 0, 2π, atol=1e-6, rtol=1e-6)
-    b1 = b1_integral / π
-
-    # Fundamental cosine component coefficient (a1)
-    # a1 = (1/π) * ∫ [f(A*sin(θ)) * cos(θ)] dθ from 0 to 2π
-    a1_integral, _ = quadgk(θ -> f(A * sin(θ)) * cos(θ), 0, 2π, atol=1e-6, rtol=1e-6)
-    a1 = a1_integral / π
-
-    # N(A) = (b1 + j*a1) / A
-    return Complex(b1, a1) / A
-end
-
-"""
-    plot_negative_inverse_df(f; A_range=range(0.1, 10, length=100), label="Nonlinearity")
-Computes and plots -1/N(A) for a given function f over a range of amplitudes.
-"""
-function plot_negative_inverse_df(f; A_range=range(0.1, 10, length=100), label="Nonlinearity", fig=nothing)
-    # Compute the negative inverse values
-    vals = similar(A_range, ComplexF64)
-    Threads.@threads for i = eachindex(A_range)
-        A = A_range[i]
-        vals[i] = -1.0 / compute_describing_function(f, A)
-    end
-
-    # Extract real and imaginary parts for plotting
-    re_vals = real.(vals)
-    im_vals = imag.(vals)
-
-    # Create plot or add to existing one
-    if fig === nothing
-        fig = plot(title="Negative Inverse Describing Function Analysis", 
-                   xlabel="Real", ylabel="Imaginary", grid=true, zeroline=true)
-    end
-
-    # Plot the curve with arrows to show increasing amplitude direction
-    plot!(fig, re_vals, im_vals, label="-1/N(A): $label", lw=2.5, arrow=true, hover=A_range)
-
-    # Mark the start (low A) and end (high A)
-    scatter!(fig, [re_vals[1]], [im_vals[1]], label="Low A", markershape=:circle)
-    scatter!(fig, [re_vals[end]], [im_vals[end]], label="High A", markershape=:square)
-
-    return fig
-end
+
+using ControlSystemsBase: describing_function, describing_function_plot, Saturation, DeadZone, Hysteresis, nonlinearity, saturation, hysteresis
+
 
 # ==============================================================================
 ## Examples of usage
 # ==============================================================================
 # 1. Sign (Relay) function: f(x) = sgn(x)
 relay(x) = sign(x)
-println("Relay N(2): ", compute_describing_function(relay, 2.0)) 
+println("Relay N(2): ", describing_function(relay, 2.0))
 # Theoretical check: 4/(π*A) = 4/(2π) ≈ 0.6366
 
-# 2. Saturation function
-saturation(x, limit=1.0) = clamp(x, -limit, limit)
-println("Saturation N(2): ", compute_describing_function(x -> saturation(x, 1.0), 2.0))
+# 2. Saturation function — analytical
+println("Saturation N(2) (analytical): ", describing_function(Saturation(1.0), 2.0))
+# Compare with numerical
+println("Saturation N(2) (numerical):  ", describing_function(x -> clamp(x, -1, 1), 2.0))
 
-# 3. Dead-zone function
-deadzone(x, d=0.1) = abs(x) < d ? 0.0 : x - sign(x)*d
-println("Dead-zone N(2): ", compute_describing_function(x -> deadzone(x, 0.1), 2.0))
+# 3. Dead-zone function — analytical
+println("Dead-zone N(2) (analytical): ", describing_function(DeadZone(0.1), 2.0))
+# Compare with numerical
+dz = DeadZone(0.1)
+println("Dead-zone N(2) (numerical):  ", describing_function(x -> dz(x), 2.0))
 
 
 # ==============================================================================
 ## Nyquist Plot with Describing Function Overlay
 # ==============================================================================
 
-using ControlSystems, Plots
-
-nonlin = saturation
+using Plots
 
 s = tf("s")
 G = 10 / (s^3 + 2s^2 + s + 1)
-# Create the linear Nyquist plot first
-p_nyq = nyquistplot(G)
-# Overlay the describing function
-plot_negative_inverse_df(nonlin, label="deadzone"; fig=p_nyq,
-    A_range=0.01:0.01:100)
-# The intersection between the Nyquist plot and -1/N(A) indicates potential limit cycles. The amplitude of the limit cycle can be estimated from the corresponding A value, use the plotly() backend to display the amplitude on hover.
+
+# Create Nyquist plot with -1/N(A) overlay for saturation
+describing_function_plot(G, Saturation(1.0); A_range=0.01:0.01:100)
+# The intersection between the Nyquist plot and -1/N(A) indicates potential limit cycles.
+# The amplitude of the limit cycle can be estimated from the corresponding A value;
+# use the plotly() backend to display the amplitude on hover.
+
 
 # ==============================================================================
 ## Simulation
 # We can simulate the closed-loop system with the nonlinearity to verify the limit cycle prediction.
 # ==============================================================================
 
-nl = nonlinearity(nonlin)
-sys_cl = feedback(G*nl)
-plot(step(sys_cl, 200))
+nl = saturation(1.0)
+sys_cl = feedback(G, nl)
+plot(impulse(sys_cl, 100))
+
+
+# ==============================================================================
+## Hysteresis describing function
+# ==============================================================================
+G = 2 / ((s+1)*(s+2)) * pid(1,1,0)
+
+H = Hysteresis(3.0, 1.5, Inf)
+describing_function_plot(G, H; A_range=1.5:0.01:100)
+
+h = hysteresis(amplitude=3.0, width=1.5, hardness=Inf)
+sys_cl = feedback(G*h)
+plot(impulse(sys_cl, 50))