diff --git a/docs/Project.toml b/docs/Project.toml
index f3773a5..05667b3 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -15,7 +15,9 @@ OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
@@ -35,5 +37,7 @@ OptimizationNLopt = "0.1, 0.2, 0.3"
 OptimizationOptimJL = "0.1, 0.2, 0.3, 0.4"
 Plots = "1"
 RecursiveArrayTools = "2, 3, 4.0"
+SciMLLogging = "2"
 SciMLSensitivity = "7"
+StochasticDiffEq = "6, 7"
 Zygote = "0.6, 0.7"
diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md
index 4b94011..e1c7eca 100644
--- a/docs/src/getting_started.md
+++ b/docs/src/getting_started.md
@@ -33,6 +33,7 @@ by defining the equation as a function with a single parameter `p=[a]`:
 ```@example ode
 using DifferentialEquations, RecursiveArrayTools, Plots, DiffEqParamEstim
 using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO
+using SciMLLogging: None
 
 function f(du, u, p, t)
     a = p[]
@@ -79,7 +80,7 @@ function:
 ```@example ode
 cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t, data),
     Optimization.AutoForwardDiff(),
-    maxiters = 10000, verbose = false)
+    maxiters = 10000, verbose = None())
 ```
 
 This objective function internally is calling the ODE solver to get solutions
@@ -88,7 +89,8 @@ Note that we set `maxiters` in a way that causes the differential equation solve
 error more quickly when in bad regions of the parameter space, speeding up the
 process. If the integrator stops early (due to divergence), then those parameters
 are given an infinite loss, and thus this is a quick way to avoid bad parameters.
-We set `verbose=false` because this divergence can get noisy. The `Optimization.AutoForwardDiff()`
+We set `verbose = None()` (from [SciMLLogging.jl](https://docs.sciml.ai/SciMLLogging/stable/))
+to silence the solver because this divergence can get noisy. The `Optimization.AutoForwardDiff()`
 is a choice of automatic differentiation, i.e., how the gradients are calculated.
 For more information on this choice, see
 [the automatic differentiation choice API](https://docs.sciml.ai/Optimization/stable/API/optimization_function/#Automatic-Differentiation-Construction-Choice-Recommendations).
@@ -171,7 +173,7 @@ We can build an objective function and solve the multiple parameter version just
 ```@example ode
 cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t, data),
     Optimization.AutoForwardDiff(),
-    maxiters = 10000, verbose = false)
+    maxiters = 10000, verbose = None())
 optprob = Optimization.OptimizationProblem(cost_function, [1.3, 0.8, 2.8, 1.2])
 result_bfgs = solve(optprob, BFGS())
 ```
@@ -189,7 +191,7 @@ cost_function = build_loss_objective(prob, Tsit5(),
     L2Loss(t, data, differ_weight = 0.3,
         data_weight = 0.7),
     Optimization.AutoForwardDiff(),
-    maxiters = 10000, verbose = false)
+    maxiters = 10000, verbose = None())
 optprob = Optimization.OptimizationProblem(cost_function, [1.3, 0.8, 2.8, 1.2])
 result_bfgs = solve(optprob, BFGS())
 ```
diff --git a/docs/src/tutorials/ensemble.md b/docs/src/tutorials/ensemble.md
index ffa83e8..469d5f8 100644
--- a/docs/src/tutorials/ensemble.md
+++ b/docs/src/tutorials/ensemble.md
@@ -89,7 +89,7 @@ losses = [L2Loss(data_times, data[:, :, i]) for i in 1:N]
 So `losses[i]` is a function which computes the loss of a solution against the data of the ith trajectory. So to build our true loss function, we sum the losses:
 
 ```@example ensemble
-loss(sim) = sum(losses[i](sim[i]) for i in 1:N)
+loss(sim) = sum(losses[i](sim.u[i]) for i in 1:N)
 ```
 
 As a double check, make sure that `loss(sim)` outputs zero (since we generated the data from sim). Now we generate data with other parameters:
diff --git a/docs/src/tutorials/generalized_likelihood.md b/docs/src/tutorials/generalized_likelihood.md
index d3fb91d..34af2b9 100644
--- a/docs/src/tutorials/generalized_likelihood.md
+++ b/docs/src/tutorials/generalized_likelihood.md
@@ -6,6 +6,7 @@ but in this case, fit just two parameters.
 
 ```@example likelihood
 using DifferentialEquations, DiffEqParamEstim, Optimization, OptimizationBBO
+using SciMLLogging: None
 f1 = function (du, u, p, t)
     du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
     du[2] = -3.0 * u[2] + u[1] * u[2]
@@ -69,7 +70,7 @@ corresponding to that distribution fit:
 
 ```@example likelihood
 obj = build_loss_objective(prob1, Tsit5(), LogLikeLoss(t, distributions),
-    maxiters = 10000, verbose = false)
+    maxiters = 10000, verbose = None())
 ```
 
 First, let's use the objective function to plot the likelihood landscape:
diff --git a/docs/src/tutorials/stochastic_evaluations.md b/docs/src/tutorials/stochastic_evaluations.md
index d8f8300..247bf98 100644
--- a/docs/src/tutorials/stochastic_evaluations.md
+++ b/docs/src/tutorials/stochastic_evaluations.md
@@ -9,6 +9,8 @@ Let's use the same Lotka-Volterra equation as before, but this time add noise:
 ```@example sde
 using DifferentialEquations, DiffEqParamEstim, Plots, Optimization, ForwardDiff,
       OptimizationOptimJL
+using StochasticDiffEq: SRIW1, SOSRI
+using SciMLLogging: None
 
 pf_func = function (du, u, p, t)
     du[1] = p[1] * u[1] - p[2] * u[1] * u[2]
@@ -53,7 +55,7 @@ We use Optim.jl for optimization below
 ```@example sde
 obj = build_loss_objective(monte_prob, SOSRI(), L2Loss(t, aggregate_data),
     Optimization.AutoForwardDiff(),
-    maxiters = 10000, verbose = false, trajectories = 1000)
+    maxiters = 10000, verbose = None(), trajectories = 1000)
 optprob = Optimization.OptimizationProblem(obj, [1.0, 0.5])
 result = solve(optprob, Optim.BFGS())
 ```
@@ -70,7 +72,7 @@ Instead, when we use `L2Loss` with first differencing enabled, we get much more
 obj = build_loss_objective(monte_prob, SRIW1(),
     L2Loss(t, aggregate_data, differ_weight = 1.0,
         data_weight = 0.5), Optimization.AutoForwardDiff(),
-    verbose = false, trajectories = 1000, maxiters = 1000)
+    verbose = None(), trajectories = 1000, maxiters = 1000)
 optprob = Optimization.OptimizationProblem(obj, [1.0, 0.5])
 result = solve(optprob, Optim.BFGS())
 result.original
diff --git a/test/dde_tests.jl b/test/dde_tests.jl
index f37efe9..952fa56 100644
--- a/test/dde_tests.jl
+++ b/test/dde_tests.jl
@@ -19,7 +19,7 @@ t = collect(range(0, stop = 10, length = 30))
 randomized = VectorOfArray([(sol(t[i]) + 0.03randn(2)) for i in 1:length(t)])
 data = convert(Array, randomized)
 
-using DiffEqParamEstim, OptimizationNLopt
+using DiffEqParamEstim, Optimization, OptimizationOptimJL
 
 function f_lotka2(du, u, h, p, t)
     du[1] = 0.5 * u[1] - p[1] * u[1] * u[2]
@@ -32,13 +32,12 @@ p = [0.5]
 prob_opt = DDEProblem(f_lotka2, u0, h, tspan, p, constant_lags = [0.5])
 cost_function = build_loss_objective(
     prob_opt, MethodOfSteps(Tsit5()),
-    L2Loss(t, data), Optimization.AutoZygote(),
+    L2Loss(t, data), Optimization.AutoForwardDiff(),
     abstol = 1.0e-8,
     reltol = 1.0e-8
 )
 
 optprob = Optimization.OptimizationProblem(cost_function, [1.0], lb = [0.0], ub = [1.0])
-opt = Opt(:GN_ESCH, 1)
 res = solve(optprob, BFGS())
 
 @test res.u[1] ≈ 0.5 atol = 5.0e-3
diff --git a/test/likelihood.jl b/test/likelihood.jl
index 455e789..5b98a46 100644
--- a/test/likelihood.jl
+++ b/test/likelihood.jl
@@ -1,5 +1,5 @@
 using OrdinaryDiffEq, DiffEqParamEstim, Distributions, Test, RecursiveArrayTools,
-    OptimizationBBO
+    Optimization, OptimizationBBO
 # ,BlackBoxOptim,
 
 pf_func = function (du, u, p, t)
diff --git a/test/multiple_shooting_objective_test.jl b/test/multiple_shooting_objective_test.jl
index 0fde9a2..aeae56e 100644
--- a/test/multiple_shooting_objective_test.jl
+++ b/test/multiple_shooting_objective_test.jl
@@ -55,4 +55,4 @@ optprob = Optimization.OptimizationProblem(
     ub = last.(bound)
 )
 result = solve(optprob, BFGS(), maxiters = 500)
-@test result.u[(end - 1):end] ≈ [1.5, 1.0] atol = 2.0e-1 broken = true
+@test result.u[(end - 1):end] ≈ [1.5, 1.0] atol = 2.0e-1
diff --git a/test/out_of_place_odes.jl b/test/out_of_place_odes.jl
index 6fbafe8..8926d35 100644
--- a/test/out_of_place_odes.jl
+++ b/test/out_of_place_odes.jl
@@ -1,4 +1,5 @@
-using OrdinaryDiffEq, Test, SciMLSensitivity, Optimization, OptimizationOptimJL
+using OrdinaryDiffEq, Test, SciMLSensitivity, Optimization, OptimizationOptimJL,
+    DiffEqParamEstim
 
 function LotkaVolterraTest_not_inplace(u, a, t)
     b, c, d = 1.0, 3.0, 1.0
diff --git a/test/test_on_monte.jl b/test/test_on_monte.jl
index fac35c0..64f1997 100644
--- a/test/test_on_monte.jl
+++ b/test/test_on_monte.jl
@@ -1,5 +1,5 @@
 using DiffEqParamEstim, OrdinaryDiffEq, StochasticDiffEq,
-    DiffEqBase, RecursiveArrayTools, OptimizationOptimJL, Zygote
+    DiffEqBase, RecursiveArrayTools, Optimization, OptimizationOptimJL, Zygote
 using Test
 
 pf_func = function (du, u, p, t)