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/3] 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/3] 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/3] 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