diff --git a/docs/Project.toml b/docs/Project.toml
index edc58b352..3b92b3023 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,9 +12,8 @@ ImplicitDifferentiation = "57b37032-215b-411a-8a7c-41a003a55207"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
-Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationGCMAES = "6f0a0517-dbc2-4a7a-8a20-99ae7f27e911"
-OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
+OptimizationIpopt = "43fad042-7963-4b32-ab19-e2a4f9a67124"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TrajectoryLimiters = "9792e600-fe43-4e4e-833b-462f466b8006"
diff --git a/docs/src/examples/automatic_differentiation.md b/docs/src/examples/automatic_differentiation.md
index f34531a19..4011c1cbd 100644
--- a/docs/src/examples/automatic_differentiation.md
+++ b/docs/src/examples/automatic_differentiation.md
@@ -108,8 +108,8 @@ We start by defining a helper function `plot_optimized` that will evaluate the p
 The constraint function `constraints` enforces the peak of the sensitivity function to be below `Msc`. Finally, we use [Optimization.jl](https://github.com/SciML/Optimization.jl) to optimize the cost function and tell it to use ForwardDiff.jl to compute the gradient of the cost function. The optimizer we use in this example is `Ipopt`.
 
 ```@example autodiff
-using Optimization, Statistics, LinearAlgebra
-using Ipopt, OptimizationMOI; MOI = OptimizationMOI.MOI
+using Statistics, LinearAlgebra
+using OptimizationIpopt
 
 function plot_optimized(P, params, res, systems)
     fig = plot(layout=(1,3), size=(1200,400), bottommargin=2Plots.mm)
@@ -170,15 +170,18 @@ Msc = 1.3        # Constraint on Ms
 
 params  = [kp, ki, kd, 0.01] # Initial guess for parameters
 
-solver = Ipopt.Optimizer()
-MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), 0)
-MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), 200)
-MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_tol"), 1e-1)
-MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_constr_viol_tol"), 1e-2)
-MOI.set(solver, MOI.RawOptimizerAttribute("acceptable_iter"), 5)
-MOI.set(solver, MOI.RawOptimizerAttribute("hessian_approximation"), "limited-memory")
+solver = IpoptOptimizer(;
+    acceptable_tol = 1e-1,
+    acceptable_iter = 5,
+    hessian_approximation = "limited-memory",
+    additional_options = Dict(
+        "print_level" => 0,
+        "max_iter" => 200,
+        "acceptable_constr_viol_tol" => 1e-2,
+    ),
+)
 
-fopt = OptimizationFunction(cost, Optimization.AutoForwardDiff(); cons=constraints)
+fopt = OptimizationFunction(cost, OptimizationIpopt.AutoForwardDiff(); cons=constraints)
 
 prob = OptimizationProblem(fopt, params, (P, systemspid);
     lb    = fill(-10.0, length(params)),
diff --git a/docs/src/lib/nonlinear.md b/docs/src/lib/nonlinear.md
index 4caaced24..46a80bd5b 100644
--- a/docs/src/lib/nonlinear.md
+++ b/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
 ```
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/Project.toml b/lib/ControlSystemsBase/Project.toml
index 8bc7fed63..7b1e0f60f 100644
--- a/lib/ControlSystemsBase/Project.toml
+++ b/lib/ControlSystemsBase/Project.toml
@@ -36,7 +36,7 @@ ComponentArrays = "0.15"
 DSP = "0.6.1, 0.7, 0.8"
 ForwardDiff = "0.10, 1.0"
 Hungarian = "0.7.0"
-ImplicitDifferentiation = "0.9"
+ImplicitDifferentiation = "0.8, 0.9"
 LinearAlgebra = "<0.0.1, 1"
 MacroTools = "0.5"
 Makie = "0.22, 0.23, 0.24"
diff --git a/lib/ControlSystemsBase/src/nonlinear_components.jl b/lib/ControlSystemsBase/src/nonlinear_components.jl
index 06d26f211..8194d43fa 100644
--- a/lib/ControlSystemsBase/src/nonlinear_components.jl
+++ b/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))")
\ No newline at end of file
+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
+
+
diff --git a/test/test_nonlinear_timeresp.jl b/test/test_nonlinear_timeresp.jl
index d34108e8e..633998693 100644
--- a/test/test_nonlinear_timeresp.jl
+++ b/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
+