From 487dcef563059915bb0dd200a6d4dec5ea7c6e91 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 13 Mar 2026 10:57:59 +0100 Subject: [PATCH 1/4] add some features for plotting describing functions --- docs/src/lib/nonlinear.md | 42 ++++ example/describing_function.jl | 111 +++------ .../src/nonlinear_components.jl | 225 ++++++++++++++++-- 3 files changed, 289 insertions(+), 89 deletions(-) diff --git a/docs/src/lib/nonlinear.md b/docs/src/lib/nonlinear.md index 46a80bd5b..1cc9bddec 100644 --- a/docs/src/lib/nonlinear.md +++ b/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 +describing_function +ControlSystemsBase.describing_function_plot ``` \ No newline at end of file diff --git a/example/describing_function.jl b/example/describing_function.jl index a36aa16c2..b4201a3e1 100644 --- a/example/describing_function.jl +++ b/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)) \ No newline at end of file +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)) diff --git a/lib/ControlSystemsBase/src/nonlinear_components.jl b/lib/ControlSystemsBase/src/nonlinear_components.jl index 8194d43fa..5ea5f1352 100644 --- a/lib/ControlSystemsBase/src/nonlinear_components.jl +++ b/lib/ControlSystemsBase/src/nonlinear_components.jl @@ -136,37 +136,230 @@ Base.show(io::IO, f::DeadZone) = f.u == -f.l ? print(io, "deadzone($(f.u))") : p ## Hysteresis ================================================================== +""" + Hysteresis{T} <: Function + +A named struct representing the internal nonlinearity used inside [`hysteresis`](@ref). +Stores the `amplitude`, `width`, and `hardness` parameters so that analytical +[`describing_function`](@ref) dispatch is possible. +""" +struct Hysteresis{T} <: Function + amplitude::T + width::T + hardness::T +end + +function (h::Hysteresis)(y) + if isfinite(h.hardness) + h.width * tanh(h.hardness * y) + else + h.width * sign(y) + end +end + +Base.show(io::IO, h::Hysteresis) = print(io, "Hysteresis(amplitude=$(h.amplitude), width=$(h.width), hardness=$(h.hardness))") + """ hysteresis(; amplitude, width, Tf, hardness) Create a hysteresis nonlinearity. The signal switches between `±amplitude` when the input crosses `±width`. `Tf` controls the time constant of the internal state that tracks the hysteresis, and `hardness` controls how sharp the transition is between the two states. ``` - - y▲ - │ -amp┌───┼───┬─► - │ │ ▲ - │ │ │ + + y▲ + │ +amp┌───┼───┬─► + │ │ ▲ + │ │ │ ──┼───┼───┼───► u - │ │ │w - ▼ │ │ - ◄─┴───┼───┘ - │ + │ │ │w + ▼ │ │ + ◄─┴───┼───┘ + │ ``` $nonlinear_warning """ function hysteresis(; amplitude=1.0, width=1.0, Tf=0.001, hardness=20.0) T = promote_type(typeof(amplitude), typeof(width), typeof(Tf), typeof(hardness)) - G = tf(1, [Tf, 1]) - if isfinite(hardness) - nl_func = y -> width * tanh(hardness*y) - else - nl_func = y -> width * sign(y) - end + G = tf(1, [Tf, 1]) + nl_func = Hysteresis(T(amplitude), T(width), T(hardness)) nl = nonlinearity(nl_func) amplitude/width*(feedback(G, -nl) - 1) end +## Describing Function ========================================================= + +""" + describing_function(f, A; N=1000) + +Compute the describing function ``N(A)`` for a static nonlinearity `f` at input +amplitude `A`. The describing function is the complex ratio of the fundamental +Fourier component of the output to the input amplitude. + +The generic method uses numerical integration (trapezoidal rule with `N` points, +which converges exponentially for smooth periodic integrands). + +Analytical overloads are provided for [`Saturation`](@ref ControlSystemsBase.saturation), +[`DeadZone`](@ref ControlSystemsBase.deadzone), and +[`Hysteresis`](@ref ControlSystemsBase.Hysteresis). + +# Arguments +- `f`: A scalar nonlinear function ``f: \\mathbb{R} \\to \\mathbb{R}``, or a + [`HammersteinWienerSystem`](@ref) with a single nonlinearity. +- `A`: Input sinusoidal amplitude (positive real). +- `N`: Number of quadrature points for the numerical method. + +# Examples +```julia +julia> using ControlSystemsBase + +julia> describing_function(Saturation(1.0), 2.0) # Analytical +0.6089977810442294 + +julia> describing_function(x -> clamp(x, -1, 1), 2.0) # Numerical +0.6089971740034709 + 2.263467191454538e-17im +``` +""" +function describing_function(f, A; N=1000) + A > 0 || throw(ArgumentError("Amplitude A must be positive, got $A")) + dθ = 2π / N + b1 = 0.0 + a1 = 0.0 + for i in 0:N-1 + θ = i * dθ + sinθ = sin(θ) + cosθ = cos(θ) + val = f(A * sinθ) + b1 += val * sinθ + a1 += val * cosθ + end + Complex(b1, a1) * dθ / (π * A) +end + +""" + describing_function(s::Saturation, A) + +Analytical describing function for symmetric saturation with limit `d`. +For `A ≤ d`, returns `1`. For `A > d`: +```math +N(A) = \\frac{2}{\\pi}\\left[\\arcsin\\frac{d}{A} + \\frac{d}{A}\\sqrt{1 - \\frac{d^2}{A^2}}\\right] +``` +Falls back to numerical computation for asymmetric saturation (`l ≠ -u`). +""" +function describing_function(s::Saturation, A; kwargs...) + if s.l != -s.u + return describing_function(x -> s(x), A; kwargs...) + end + d = s.u + A ≤ d && return complex(float(one(d))) + r = d / A + complex((2/π) * (asin(r) + r * sqrt(1 - r^2))) +end + +""" + describing_function(dz::DeadZone, A) + +Analytical describing function for symmetric dead-zone with threshold `d`. +For `A ≤ d`, returns `0`. For `A > d`: +```math +N(A) = 1 - \\frac{2}{\\pi}\\left[\\arcsin\\frac{d}{A} + \\frac{d}{A}\\sqrt{1 - \\frac{d^2}{A^2}}\\right] +``` +Falls back to numerical computation for asymmetric dead-zone (`l ≠ -u`). +""" +function describing_function(dz::DeadZone, A; kwargs...) + if dz.l != -dz.u + return describing_function(x -> dz(x), A; kwargs...) + end + d = dz.u + A ≤ d && return complex(zero(float(d))) + r = d / A + complex(1 - (2/π) * (asin(r) + r * sqrt(1 - r^2))) +end + +""" + describing_function(h::Hysteresis, A) + +Analytical describing function for ideal relay with hysteresis. +For `A ≤ width`, returns `0`. For `A > width`: +```math +N(A) = \\frac{4M}{\\pi A}\\left[\\sqrt{1 - \\frac{b^2}{A^2}} - j\\frac{b}{A}\\right] +``` +where `M = amplitude` and `b = width`. +""" +function describing_function(h::Hysteresis, A; kwargs...) + M = h.amplitude + b = h.width + A ≤ b && return complex(zero(float(M))) + r = b / A + (4M / (π * A)) * complex(sqrt(1 - r^2), -r) + # describing_function(x->h(x), A; kwargs...) +end + +""" + describing_function(sys::HammersteinWienerSystem, A; kwargs...) + +Compute the describing function for the nonlinearity contained in a +[`HammersteinWienerSystem`](@ref). The system must contain exactly one nonlinearity. +""" +function describing_function(sys::HammersteinWienerSystem, A; kwargs...) + length(sys.f) == 1 || error("describing_function requires a single nonlinearity, got $(length(sys.f))") + isempty(sys.A) && sys.D == [0 1; 1 0] || error("Computing the describing function of a dynamic nonlinearity or a generic nonlinear system is not supported.") + describing_function(sys.f[1], A; kwargs...) +end + +""" + describing_function_plot(L, f; A_range, kwargs...) + describing_function_plot(sys::HammersteinWienerSystem; A_range, kwargs...) + +Create a Nyquist plot of the linear system `L` with the curve ``-1/N(A)`` overlaid, +where ``N(A)`` is the describing function of the nonlinearity `f`. + +If a [`HammersteinWienerSystem`](@ref) is passed, the linear part and nonlinearity +are extracted automatically. + +Intersections of the Nyquist curve and the ``-1/N(A)`` curve indicate potential +limit cycles. + +# Arguments +- `L`: A linear system (transfer function or state-space). +- `f`: A static nonlinearity (callable or struct like [`Saturation`](@ref ControlSystemsBase.saturation)). +- `A_range`: Range of amplitudes for the describing function (default `range(0.1, 10, length=200)`). +- `kwargs...`: Additional keyword arguments passed to `nyquistplot`. + +# Example +```julia +using ControlSystems, Plots +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) +``` +""" +function describing_function_plot(L, f; A_range=range(0.1, 10, length=200), N=1000, kwargs...) + neg_inv = map(A_range) do a + n = describing_function(f, a; N) + iszero(n) ? complex(NaN) : -1 / n + end + fig = nyquistplot(L; kwargs...) + RecipesBase.plot!(fig, real.(neg_inv), imag.(neg_inv), lab="-1/N(A)", lw=2.5, hover = A_range) + fig +end + +function describing_function_plot(sys::HammersteinWienerSystem; kwargs...) + length(sys.f) == 1 || error("describing_function_plot requires a single nonlinearity, got $(length(sys.f))") + # Close the outer (upper) loop with unity negative feedback to get the + # linear transfer from the nonlinearity output (Δu) back to its input (Δy). + # This is L_eff = -L_cl where L_cl is the positive-feedback loop transfer. + # The condition for limit cycles is L_eff(jω) = -1/N_f(A), where N_f is + # the DF of the actual static nonlinearity f (not a composite LFT system). + L = sminreal(-lft(sys.P.P, ss(-1.0), :u)) + f = sys.f[1] + f isa Hysteresis && @warn "This method of computing the describing function for a hysteresis element is known to be inaccurate. Hysteresis is a dynamic nonlinearity, and when approximated by an LFT system the linear part does not attenuate high frequencies as assumed by the describing function theory. For more accurate limit cycle predictions, consider using the `Hysteresis` struct (capital `H`) and call `describing_function_plot(P, H)` rather than forming the composite system by composition with `hysteresis(...)` (lowercase `h`). Computing the describing function with an instance of the `Hysteresis` struct uses the analytical formula for the describing function of an ideal hysteresis element." + # Wrap in a closure to bypass analytical DF dispatch (e.g., Hysteresis + # has an analytical DF for the *composite* hysteresis behavior, but here + # we need the DF of the raw static function since the internal linear + # dynamics are already included in L). + describing_function_plot(L, x -> f(x); kwargs...) +end + From 0052d7f502b6ed0e879d5885eb3aab74511a81bf Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 13 Mar 2026 11:52:59 +0100 Subject: [PATCH 2/4] qualify --- docs/src/lib/nonlinear.md | 2 +- lib/ControlSystemsBase/src/nonlinear_components.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/lib/nonlinear.md b/docs/src/lib/nonlinear.md index 1cc9bddec..604f73167 100644 --- a/docs/src/lib/nonlinear.md +++ b/docs/src/lib/nonlinear.md @@ -334,6 +334,6 @@ ControlSystemsBase.deadzone ControlSystemsBase.linearize ControlSystemsBase.hysteresis ControlSystemsBase.Hysteresis -describing_function +ControlSystemsBase.describing_function ControlSystemsBase.describing_function_plot ``` \ No newline at end of file diff --git a/lib/ControlSystemsBase/src/nonlinear_components.jl b/lib/ControlSystemsBase/src/nonlinear_components.jl index 5ea5f1352..3d11bb721 100644 --- a/lib/ControlSystemsBase/src/nonlinear_components.jl +++ b/lib/ControlSystemsBase/src/nonlinear_components.jl @@ -137,7 +137,7 @@ Base.show(io::IO, f::DeadZone) = f.u == -f.l ? print(io, "deadzone($(f.u))") : p ## Hysteresis ================================================================== """ - Hysteresis{T} <: Function + Hysteresis(amplitude, width, hardness) A named struct representing the internal nonlinearity used inside [`hysteresis`](@ref). Stores the `amplitude`, `width`, and `hardness` parameters so that analytical From c2b77e67408807fc1f3ecbc25198c355f1ed4c14 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Sun, 15 Mar 2026 12:54:11 +0100 Subject: [PATCH 3/4] add tests --- .../test/test_hammerstein_wiener.jl | 59 +++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl b/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl index 035e74056..87d369332 100644 --- a/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl +++ b/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl @@ -126,3 +126,62 @@ G = ss(1) A, B = ControlSystemsBase.linearize((x,u)->x.^2 + sin.(u), [1.0], [2.0]) @test A ≈ [2.0;;] @test B ≈ [cos(2);;] + +## Test nonlinear_components coverage ========================================== +using ControlSystemsBase: Saturation, DeadZone, Offset, Hysteresis, describing_function, deadzone + +@testset "nonlinear_components" begin + + # Show methods + @test sprint(show, Saturation(1.0)) == "saturation(1.0)" + @test sprint(show, Saturation(-2.0, 1.0)) == "saturation(-2.0, 1.0)" + @test sprint(show, Offset(1.5)) == "offset(1.5)" + @test sprint(show, DeadZone(1.0)) == "deadzone(1.0)" + @test sprint(show, DeadZone(-2.0, 1.0)) == "deadzone(-2.0, 1.0)" + @test sprint(show, Hysteresis(1.0, 0.5, 20.0)) == "Hysteresis(amplitude=1.0, width=0.5, hardness=20.0)" + + # Hysteresis callable + h = Hysteresis(1.0, 0.5, 20.0) + @test isfinite(h(10.0)) + @test h(10.0) ≈ 0.5 * tanh(20.0 * 10.0) + h_inf = Hysteresis(1.0, 0.5, Inf) + @test h_inf(10.0) == 0.5 + @test h_inf(-10.0) == -0.5 + + # Vector constructors (MIMO) + @test saturation([1.0, 2.0]) isa HammersteinWienerSystem + @test deadzone([0.5, 1.0]) isa HammersteinWienerSystem + @test offset([1.0, 2.0]) isa HammersteinWienerSystem + + # describing_function — numerical + df_num = describing_function(x -> clamp(x, -1, 1), 2.0) + df_ana = describing_function(Saturation(1.0), 2.0) + @test df_num ≈ df_ana rtol=1e-3 + @test_throws ArgumentError describing_function(abs, -1.0) + + # describing_function — Saturation analytical + @test describing_function(Saturation(2.0), 1.0) == 1.0 + 0im # A ≤ d + @test describing_function(Saturation(1.0), 2.0) ≈ 0.6089977810442294 + @test describing_function(Saturation(-0.5, 1.0), 2.0) isa Complex # asymmetric fallback + + # describing_function — DeadZone analytical + @test describing_function(DeadZone(2.0), 1.0) == 0.0 + 0im # A ≤ d + df_sat = describing_function(Saturation(1.0), 2.0) + df_dz = describing_function(DeadZone(1.0), 2.0) + @test df_sat + df_dz ≈ 1.0 + 0im # complementary + @test describing_function(DeadZone(-0.5, 1.0), 2.0) isa Complex # asymmetric fallback + + # describing_function — Hysteresis analytical + @test describing_function(Hysteresis(1.0, 0.5, 20.0), 0.3) == 0.0 + 0im # A ≤ width + h = Hysteresis(1.0, 0.5, 20.0) + A = 2.0 + r = 0.5 / A + expected = (4 * 1.0 / (π * A)) * complex(sqrt(1 - r^2), -r) + @test describing_function(h, A) ≈ expected + + # describing_function — HammersteinWienerSystem + nl = saturation(1.0) + @test describing_function(nl, 2.0) ≈ describing_function(Saturation(1.0), 2.0) + @test_throws ErrorException describing_function(nonlinearity(abs2) + nonlinearity(abs), 1.0) + +end From 6327e5f5e46b88fbc8016532166d9167339da793 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 6 Jul 2026 07:14:13 +0000 Subject: [PATCH 4/4] Address review findings for describing-function functionality Correctness: - describing_function(sys::HammersteinWienerSystem) now computes the DF of the static callable when the nonlinearity is a Hysteresis, instead of dispatching to the analytical relay-with-hysteresis formula which describes the composite element realized by hysteresis(), not the static function the system applies. - describing_function_plot(sys::HammersteinWienerSystem) errors for systems whose external part is not SISO, and the docstring now states the open-loop unity-negative-feedback assumption (do not pass feedback(G*nl)). - Add promoting outer constructor Hysteresis(amplitude, width, hardness) so mixed numeric types work as documented. - Analytical describing_function overloads validate A > 0 like the generic numerical method, instead of silently returning a plausible value. - Fix docstring example output (the analytical Saturation DF is complex). - Docstring no longer claims exponential convergence of the trapezoidal rule unconditionally; document increasing N for discontinuous nonlinearities. - Restore direction-of-increasing-amplitude cues on the -1/N(A) curve (arrow and labeled endpoint markers) that the deleted example code had; hover only works with the plotly backend. Cleanup: - DeadZone DF computed via the exact identity N_dz(A) = 1 - N_sat(A), removing the duplicated formula and handling the asymmetric case exactly. - Name the numerical quadrature core describing_function_numerical and call it directly where analytical dispatch must be bypassed. - Remove commented-out dead code in the Hysteresis DF method. - Clarify the dual semantics of the Hysteresis struct (static callable vs analytical DF token) in its docstring. Co-Authored-By: Claude Fable 5 --- example/describing_function.jl | 3 + .../src/nonlinear_components.jl | 93 ++++++++++++++----- .../test/test_hammerstein_wiener.jl | 23 +++++ 3 files changed, 94 insertions(+), 25 deletions(-) diff --git a/example/describing_function.jl b/example/describing_function.jl index b4201a3e1..0eca14f51 100644 --- a/example/describing_function.jl +++ b/example/describing_function.jl @@ -14,6 +14,9 @@ using ControlSystemsBase: describing_function, describing_function_plot, Saturat relay(x) = sign(x) println("Relay N(2): ", describing_function(relay, 2.0)) # Theoretical check: 4/(π*A) = 4/(2π) ≈ 0.6366 +# The numerical method uses a fixed-grid trapezoidal rule; for discontinuous +# nonlinearities like the relay it converges slowly, increase the keyword +# argument N for higher accuracy. # 2. Saturation function — analytical println("Saturation N(2) (analytical): ", describing_function(Saturation(1.0), 2.0)) diff --git a/lib/ControlSystemsBase/src/nonlinear_components.jl b/lib/ControlSystemsBase/src/nonlinear_components.jl index 3d11bb721..b726bf4d0 100644 --- a/lib/ControlSystemsBase/src/nonlinear_components.jl +++ b/lib/ControlSystemsBase/src/nonlinear_components.jl @@ -142,12 +142,18 @@ Base.show(io::IO, f::DeadZone) = f.u == -f.l ? print(io, "deadzone($(f.u))") : p A named struct representing the internal nonlinearity used inside [`hysteresis`](@ref). Stores the `amplitude`, `width`, and `hardness` parameters so that analytical [`describing_function`](@ref) dispatch is possible. + +Note: as a callable, `Hysteresis` evaluates the *static* internal function +`width * tanh(hardness * y)`, while the analytical [`describing_function`](@ref) +overload describes the complete relay-with-hysteresis element realized by +[`hysteresis`](@ref) (which embeds this static function in a feedback loop). """ struct Hysteresis{T} <: Function amplitude::T width::T hardness::T end +Hysteresis(amplitude, width, hardness) = Hysteresis(promote(amplitude, width, hardness)...) function (h::Hysteresis)(y) if isfinite(h.hardness) @@ -198,8 +204,10 @@ Compute the describing function ``N(A)`` for a static nonlinearity `f` at input amplitude `A`. The describing function is the complex ratio of the fundamental Fourier component of the output to the input amplitude. -The generic method uses numerical integration (trapezoidal rule with `N` points, -which converges exponentially for smooth periodic integrands). +The generic method uses numerical integration (trapezoidal rule with `N` points). +The trapezoidal rule converges very rapidly for smooth periodic integrands, but only +slowly when `f` is discontinuous (e.g., a relay); increase `N` if high accuracy is +required for such nonlinearities. Analytical overloads are provided for [`Saturation`](@ref ControlSystemsBase.saturation), [`DeadZone`](@ref ControlSystemsBase.deadzone), and @@ -216,13 +224,18 @@ Analytical overloads are provided for [`Saturation`](@ref ControlSystemsBase.sat julia> using ControlSystemsBase julia> describing_function(Saturation(1.0), 2.0) # Analytical -0.6089977810442294 +0.6089977810442294 + 0.0im julia> describing_function(x -> clamp(x, -1, 1), 2.0) # Numerical 0.6089971740034709 + 2.263467191454538e-17im ``` """ -function describing_function(f, A; N=1000) +describing_function(f, A; kwargs...) = describing_function_numerical(f, A; kwargs...) + +# Numerical quadrature core. Call this directly instead of `describing_function` +# where analytical dispatch on the nonlinearity type must be bypassed, e.g., to +# compute the describing function of the static callable of a `Hysteresis`. +function describing_function_numerical(f, A; N=1000) A > 0 || throw(ArgumentError("Amplitude A must be positive, got $A")) dθ = 2π / N b1 = 0.0 @@ -249,8 +262,9 @@ N(A) = \\frac{2}{\\pi}\\left[\\arcsin\\frac{d}{A} + \\frac{d}{A}\\sqrt{1 - \\fra Falls back to numerical computation for asymmetric saturation (`l ≠ -u`). """ function describing_function(s::Saturation, A; kwargs...) + A > 0 || throw(ArgumentError("Amplitude A must be positive, got $A")) if s.l != -s.u - return describing_function(x -> s(x), A; kwargs...) + return describing_function_numerical(s, A; kwargs...) end d = s.u A ≤ d && return complex(float(one(d))) @@ -261,40 +275,41 @@ end """ describing_function(dz::DeadZone, A) -Analytical describing function for symmetric dead-zone with threshold `d`. -For `A ≤ d`, returns `0`. For `A > d`: +Analytical describing function for dead-zone. For symmetric dead-zone with threshold +`d` and `A ≤ d`, returns `0`. For `A > d`: ```math N(A) = 1 - \\frac{2}{\\pi}\\left[\\arcsin\\frac{d}{A} + \\frac{d}{A}\\sqrt{1 - \\frac{d^2}{A^2}}\\right] ``` -Falls back to numerical computation for asymmetric dead-zone (`l ≠ -u`). +Computed through the identity `dz(x) = x - clamp(x, l, u)`, which by linearity of the +describing function gives ``N_{dz}(A) = 1 - N_{sat}(A)``. This holds for asymmetric +dead-zone (`l ≠ -u`) as well, in which case the saturation part is computed numerically. """ function describing_function(dz::DeadZone, A; kwargs...) - if dz.l != -dz.u - return describing_function(x -> dz(x), A; kwargs...) - end - d = dz.u - A ≤ d && return complex(zero(float(d))) - r = d / A - complex(1 - (2/π) * (asin(r) + r * sqrt(1 - r^2))) + 1 - describing_function(Saturation(dz.l, dz.u), A; kwargs...) end """ describing_function(h::Hysteresis, A) -Analytical describing function for ideal relay with hysteresis. +Analytical describing function for an ideal relay with hysteresis, i.e., the complete +hysteresis element realized by [`hysteresis`](@ref) (`hardness` is ignored). For `A ≤ width`, returns `0`. For `A > width`: ```math N(A) = \\frac{4M}{\\pi A}\\left[\\sqrt{1 - \\frac{b^2}{A^2}} - j\\frac{b}{A}\\right] ``` where `M = amplitude` and `b = width`. + +Note: this is *not* the describing function of the static callable +`h(y) = width * tanh(hardness * y)`; to compute that, bypass analytical dispatch +with a closure: `describing_function(y -> h(y), A)`. """ function describing_function(h::Hysteresis, A; kwargs...) + A > 0 || throw(ArgumentError("Amplitude A must be positive, got $A")) M = h.amplitude b = h.width A ≤ b && return complex(zero(float(M))) r = b / A (4M / (π * A)) * complex(sqrt(1 - r^2), -r) - # describing_function(x->h(x), A; kwargs...) end """ @@ -302,11 +317,24 @@ end Compute the describing function for the nonlinearity contained in a [`HammersteinWienerSystem`](@ref). The system must contain exactly one nonlinearity. + +The result is the describing function of the static function actually applied by the +system. In particular, for `nonlinearity(Hysteresis(...))` this is the numerically +computed describing function of the static callable `width*tanh(hardness*y)`, not the +analytical relay-with-hysteresis formula of `describing_function(::Hysteresis, A)`. """ function describing_function(sys::HammersteinWienerSystem, A; kwargs...) length(sys.f) == 1 || error("describing_function requires a single nonlinearity, got $(length(sys.f))") isempty(sys.A) && sys.D == [0 1; 1 0] || error("Computing the describing function of a dynamic nonlinearity or a generic nonlinear system is not supported.") - describing_function(sys.f[1], A; kwargs...) + f = sys.f[1] + if f isa Hysteresis + # The system applies only the static callable of the Hysteresis struct; the + # analytical overload describes the composite relay-with-hysteresis element + # and would give the wrong answer here. + describing_function_numerical(f, A; kwargs...) + else + describing_function(f, A; kwargs...) + end end """ @@ -317,10 +345,16 @@ Create a Nyquist plot of the linear system `L` with the curve ``-1/N(A)`` overla where ``N(A)`` is the describing function of the nonlinearity `f`. If a [`HammersteinWienerSystem`](@ref) is passed, the linear part and nonlinearity -are extracted automatically. +are extracted automatically. The system must be an *open-loop* interconnection with a +single nonlinearity and a SISO external part, e.g., `G*saturation(1.0)`; the loop is +closed internally with unity negative feedback around the external channel. Do not pass +an already closed-loop system such as `feedback(G*nl)` — the loop would then be closed +twice and the plotted loop transfer would be wrong. Intersections of the Nyquist curve and the ``-1/N(A)`` curve indicate potential -limit cycles. +limit cycles. The ends of the ``-1/N(A)`` curve are marked to indicate the direction +of increasing amplitude; with the plotly backend, hovering over the curve shows the +amplitude ``A``. # Arguments - `L`: A linear system (transfer function or state-space). @@ -342,12 +376,21 @@ function describing_function_plot(L, f; A_range=range(0.1, 10, length=200), N=10 iszero(n) ? complex(NaN) : -1 / n end fig = nyquistplot(L; kwargs...) - RecipesBase.plot!(fig, real.(neg_inv), imag.(neg_inv), lab="-1/N(A)", lw=2.5, hover = A_range) + RecipesBase.plot!(fig, real.(neg_inv), imag.(neg_inv), lab="-1/N(A)", lw=2.5, arrow=true, hover = A_range) + # Mark the curve endpoints to indicate the direction of increasing amplitude; + # the arrow and hover information are not supported by all plotting backends. + i1 = findfirst(isfinite, neg_inv) + i2 = findlast(isfinite, neg_inv) + if i1 !== nothing && i1 != i2 + RecipesBase.plot!(fig, [real(neg_inv[i1])], [imag(neg_inv[i1])], seriestype=:scatter, markershape=:circle, lab="A = $(round(A_range[i1], sigdigits=3))") + RecipesBase.plot!(fig, [real(neg_inv[i2])], [imag(neg_inv[i2])], seriestype=:scatter, markershape=:square, lab="A = $(round(A_range[i2], sigdigits=3))") + end fig end function describing_function_plot(sys::HammersteinWienerSystem; kwargs...) length(sys.f) == 1 || error("describing_function_plot requires a single nonlinearity, got $(length(sys.f))") + sys.P.nu1 == 1 && sys.P.ny1 == 1 || error("describing_function_plot requires a SISO external part, got $(sys.P.ny1) external outputs and $(sys.P.nu1) external inputs.") # Close the outer (upper) loop with unity negative feedback to get the # linear transfer from the nonlinearity output (Δu) back to its input (Δy). # This is L_eff = -L_cl where L_cl is the positive-feedback loop transfer. @@ -356,10 +399,10 @@ function describing_function_plot(sys::HammersteinWienerSystem; kwargs...) L = sminreal(-lft(sys.P.P, ss(-1.0), :u)) f = sys.f[1] f isa Hysteresis && @warn "This method of computing the describing function for a hysteresis element is known to be inaccurate. Hysteresis is a dynamic nonlinearity, and when approximated by an LFT system the linear part does not attenuate high frequencies as assumed by the describing function theory. For more accurate limit cycle predictions, consider using the `Hysteresis` struct (capital `H`) and call `describing_function_plot(P, H)` rather than forming the composite system by composition with `hysteresis(...)` (lowercase `h`). Computing the describing function with an instance of the `Hysteresis` struct uses the analytical formula for the describing function of an ideal hysteresis element." - # Wrap in a closure to bypass analytical DF dispatch (e.g., Hysteresis - # has an analytical DF for the *composite* hysteresis behavior, but here - # we need the DF of the raw static function since the internal linear - # dynamics are already included in L). + # Wrap in a closure to bypass analytical DF dispatch and force + # describing_function_numerical (e.g., Hysteresis has an analytical DF for + # the *composite* hysteresis element, but here we need the DF of the raw + # static function since the internal linear dynamics are already included in L). describing_function_plot(L, x -> f(x); kwargs...) end diff --git a/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl b/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl index 87d369332..2c26e3b2c 100644 --- a/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl +++ b/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl @@ -184,4 +184,27 @@ using ControlSystemsBase: Saturation, DeadZone, Offset, Hysteresis, describing_f @test describing_function(nl, 2.0) ≈ describing_function(Saturation(1.0), 2.0) @test_throws ErrorException describing_function(nonlinearity(abs2) + nonlinearity(abs), 1.0) + # Promoting Hysteresis constructor + @test Hysteresis(3, 1.5, Inf) === Hysteresis(3.0, 1.5, Inf) + @test Hysteresis(1, 1, 20.0) === Hysteresis(1.0, 1.0, 20.0) + + # Amplitude validation in analytical overloads + @test_throws ArgumentError describing_function(Saturation(1.0), -2.0) + @test_throws ArgumentError describing_function(DeadZone(1.0), 0.0) + @test_throws ArgumentError describing_function(Hysteresis(1.0, 0.5, 20.0), -1.0) + + # DeadZone through the identity N_dz = 1 - N_sat, also for asymmetric dead-zone + @test describing_function(DeadZone(-0.5, 1.0), 2.0) ≈ 1 - describing_function(Saturation(-0.5, 1.0), 2.0) + + # A HammersteinWienerSystem containing a Hysteresis applies only the static + # callable, so the DF must be that of the callable, not the analytical + # relay-with-hysteresis formula + df_hw = describing_function(nonlinearity(h), 2.0) + @test df_hw ≈ describing_function(x -> h(x), 2.0) + @test abs(imag(df_hw)) < 1e-10 # a static odd nonlinearity has a real DF + @test !(df_hw ≈ describing_function(h, 2.0)) + + # describing_function_plot requires a SISO external part + @test_throws ErrorException ControlSystemsBase.describing_function_plot([tf(1); tf(1)] * nl) + end