diff --git a/docs/src/lib/nonlinear.md b/docs/src/lib/nonlinear.md index 46a80bd5b..604f73167 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 +ControlSystemsBase.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..0eca14f51 100644 --- a/example/describing_function.jl +++ b/example/describing_function.jl @@ -1,101 +1,69 @@ #= 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 +# 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 -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..b726bf4d0 100644 --- a/lib/ControlSystemsBase/src/nonlinear_components.jl +++ b/lib/ControlSystemsBase/src/nonlinear_components.jl @@ -136,37 +136,273 @@ Base.show(io::IO, f::DeadZone) = f.u == -f.l ? print(io, "deadzone($(f.u))") : p ## Hysteresis ================================================================== +""" + 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 +[`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) + 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). +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 +[`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 + 0.0im + +julia> describing_function(x -> clamp(x, -1, 1), 2.0) # Numerical +0.6089971740034709 + 2.263467191454538e-17im +``` +""" +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 + 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...) + A > 0 || throw(ArgumentError("Amplitude A must be positive, got $A")) + if s.l != -s.u + return describing_function_numerical(s, 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 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] +``` +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...) + 1 - describing_function(Saturation(dz.l, dz.u), A; kwargs...) +end + +""" + describing_function(h::Hysteresis, A) + +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) +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. + +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.") + 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 + +""" + 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. 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. 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). +- `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, 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. + # 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 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 035e74056..2c26e3b2c 100644 --- a/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl +++ b/lib/ControlSystemsBase/test/test_hammerstein_wiener.jl @@ -126,3 +126,85 @@ 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) + + # 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