Add hysteresis nonlinearity

Author: baggepinnen

docs/src/lib/nonlinear.md

@@ -229,6 +229,39 @@ f2 = plot(step(feedback(duffing, C), 8), plotx=true, plot_title="Controlled wigg
 plot(f1,f2, size=(1300,800))
 ```
 
+### Hysteresis
+Here we demonstrate that we may use this simple framework to model also stateful nonlinearities, such as hysteresis. The `hysteresis` function internally creates a feedback interconnection between a fast first-order system and a `sign` or `tanh` nonlinearity to create a simple hysteresis loop. The width and amplitude of the loop can be adjusted through the parameters `width` and `amplitude`, respectively.
+```@example HYSTERESIS
+using ControlSystems, Plots
+import ControlSystemsBase: hysteresis
+
+amplitude = 0.7
+width = 1.5
+sys_hyst =  hysteresis(; width, amplitude)
+
+t = 0:0.01:20
+ufun(y,x,t) = y .= 5.0 .* sin(t) ./ (1+t/5) # A sine wave that sweeps back and forth with decreasing amplitude
+
+res = lsim(sys_hyst, ufun, t)
+
+p1 = plot(res.u[:], res.y[:], 
+    title="Hysteresis Loop",
+    xlabel="Input u(t)", 
+    ylabel="Output y(t)",
+    lw=2, color=:blue, lab="", framestyle=:zerolines)
+
+hline!([amplitude -amplitude], l=:dash, c=:red, label=["Amplitude" ""])
+vline!([width -width], l=:dash, c=:green, label=["Width" ""])
+
+# Plotting time series to see the "jumps" in the response
+p2 = plot(t, [res.u[:] res.y[:]], 
+    title="Time Domain Response",
+    label=["Input (Sine)" "Output (Hysteresis)"], 
+    xlabel="Time (s)",
+    lw=2)
+
+plot(p1, p2, layout=(1,2), size=(900, 400))
+```
 
 ## Limitations
 - Remember, this functionality is experimental and subject to breakage.
@@ -260,4 +293,5 @@ ControlSystemsBase.saturation
 ControlSystemsBase.ratelimit
 ControlSystemsBase.deadzone
 ControlSystemsBase.linearize
+ControlSystemsBase.hysteresis
 ```

lib/ControlSystemsBase/src/nonlinear_components.jl

@@ -131,4 +131,42 @@ See [`ControlSystemsBase.offset`](@ref) for handling operating points.
 """
 deadzone(args...) = nonlinearity(DeadZone(args...))
 deadzone(v::AbstractVector, args...) = nonlinearity(DeadZone.(v, args...))
-Base.show(io::IO, f::DeadZone) = f.u == -f.l ? print(io, "deadzone($(f.u))") : print(io, "deadzone($(f.l), $(f.u))")
+Base.show(io::IO, f::DeadZone) = f.u == -f.l ? print(io, "deadzone($(f.u))") : print(io, "deadzone($(f.l), $(f.u))")
+
+
+## Hysteresis ==================================================================
+
+"""
+    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┌───┼───┬─►  
+   │   │   ▲    
+   │   │   │    
+ ──┼───┼───┼───► u
+   │   │   │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
+    nl = nonlinearity(nl_func)
+    amplitude/width*(feedback(G, -nl) - 1)
+end
+
+

test/test_nonlinear_timeresp.jl

@@ -235,3 +235,25 @@ res = lsim(nl, (x,t)->2sin(t), 0:0.01:4pi)
 @test minimum(res.y) ≈ -0.6 rtol=1e-3
 @test maximum(res.y) ≈ 1.3 rtol=1e-3
 
+## Hysteresis
+using ControlSystemsBase: hysteresis
+
+# Construction
+@test hysteresis() isa HammersteinWienerSystem
+@test hysteresis(; amplitude=0.7, width=1.5) isa HammersteinWienerSystem
+
+# Simulate with sine input and verify output bounded by amplitude
+amplitude = 0.7
+width = 1.5
+sys_hyst = hysteresis(; width, amplitude)
+t = 0:0.001:10
+res = lsim(sys_hyst, (x,t) -> 5sin(t), t)
+@test maximum(res.y) ≈ amplitude rtol=0.1
+@test minimum(res.y) ≈ -amplitude rtol=0.1
+
+# Finite hardness (smooth tanh transition)
+sys_hyst_smooth = hysteresis(; width, amplitude, hardness=5.0)
+res_smooth = lsim(sys_hyst_smooth, (x,t) -> 5sin(t), t)
+@test maximum(res_smooth.y) ≈ amplitude rtol=0.15
+@test minimum(res_smooth.y) ≈ -amplitude rtol=0.15
+