diff --git a/Project.toml b/Project.toml
index 0d2ddfcf0..4f19994aa 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "Optimization"
 uuid = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
-version = "3.12.1"
+version = "3.13.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
diff --git a/docs/pages.jl b/docs/pages.jl
index 1d7d2971c..0f7fe2519 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -31,5 +31,6 @@ pages = ["index.md",
         "Optim.jl" => "optimization_packages/optim.md",
         "Optimisers.jl" => "optimization_packages/optimisers.md",
         "QuadDIRECT.jl" => "optimization_packages/quaddirect.md",
+        "SpeedMapping.jl" => "optimization_packages/speedmapping.md",
     ],
 ]
diff --git a/docs/src/examples/rosenbrock.md b/docs/src/examples/rosenbrock.md
index 240124b57..9f74da162 100644
--- a/docs/src/examples/rosenbrock.md
+++ b/docs/src/examples/rosenbrock.md
@@ -69,7 +69,7 @@ sol = solve(prob, IPNewton())
 prob = OptimizationProblem(optf, x0, _p, lcons = [0.5], ucons = [0.5],
                            lb = [-500.0, -500.0], ub = [50.0, 50.0])
 sol = solve(prob, IPNewton()) # Notice now that x[1]^2 + x[2]^2 ≈ 0.5:
-# cons(sol.minimizer, _p) = 0.49999999999999994
+# cons(sol.u, _p) = 0.49999999999999994
 
 function con_c(res, x, p)
     res .= [x[1]^2 + x[2]^2]
@@ -77,7 +77,7 @@ end
 
 optf = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(); cons = con_c)
 prob = OptimizationProblem(optf, x0, _p, lcons = [-Inf], ucons = [0.25^2])
-sol = solve(prob, IPNewton()) # -Inf < cons_circ(sol.minimizer, _p) = 0.25^2
+sol = solve(prob, IPNewton()) # -Inf < cons_circ(sol.u, _p) = 0.25^2
 
 function con2_c(res, x, p)
     res .= [x[1]^2 + x[2]^2, x[2] * sin(x[1]) - x[1]]
diff --git a/lib/OptimizationFlux/src/OptimizationFlux.jl b/lib/OptimizationFlux/src/OptimizationFlux.jl
index 4bb9d440a..fb5f0deff 100644
--- a/lib/OptimizationFlux/src/OptimizationFlux.jl
+++ b/lib/OptimizationFlux/src/OptimizationFlux.jl
@@ -4,43 +4,99 @@ using Reexport, Printf, ProgressLogging
 @reexport using Flux, Optimization
 using Optimization.SciMLBase
 
-function SciMLBase.__solve(prob::OptimizationProblem, opt::Flux.Optimise.AbstractOptimiser,
-                           data = Optimization.DEFAULT_DATA;
-                           maxiters::Number = 0, callback = (args...) -> (false),
-                           progress = false, save_best = true, kwargs...)
-    if data != Optimization.DEFAULT_DATA
-        maxiters = length(data)
+struct FluxOptimizationCache{F <: OptimizationFunction, RC, O, D} <:
+       SciMLBase.AbstractOptimizationCache
+    f::F
+    reinit_cache::RC
+    opt::O
+    data::D
+    solver_args::NamedTuple
+end
+
+function Base.getproperty(cache::FluxOptimizationCache, x::Symbol)
+    if x in fieldnames(Optimization.ReInitCache)
+        return getfield(cache.reinit_cache, x)
+    end
+    return getfield(cache, x)
+end
+
+function FluxOptimizationCache(prob::OptimizationProblem, opt, data; kwargs...)
+    reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit`
+    f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype)
+    return FluxOptimizationCache(f, reinit_cache, opt, data, NamedTuple(kwargs))
+end
+
+SciMLBase.supports_opt_cache_interface(opt::Flux.Optimise.AbstractOptimiser) = true
+SciMLBase.has_reinit(cache::FluxOptimizationCache) = true
+function SciMLBase.reinit!(cache::FluxOptimizationCache; p = missing, u0 = missing)
+    if p === missing && u0 === missing
+        p, u0 = cache.p, cache.u0
+    else # at least one of them has a value
+        if p === missing
+            p = cache.p
+        end
+        if u0 === missing
+            u0 = cache.u0
+        end
+        if (eltype(p) <: Pair && !isempty(p)) || (eltype(u0) <: Pair && !isempty(u0)) # one is a non-empty symbolic map
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :ps) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `p` keyword argument as a vector of values, paying attention to parameter order."))
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :states) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
+            p, u0 = SciMLBase.process_p_u0_symbolic(cache, p, u0)
+        end
+    end
+
+    cache.reinit_cache.p = p
+    cache.reinit_cache.u0 = u0
+
+    return cache
+end
+
+function SciMLBase.__init(prob::OptimizationProblem, opt::Flux.Optimise.AbstractOptimiser,
+                          data = Optimization.DEFAULT_DATA;
+                          maxiters::Number = 0, callback = (args...) -> (false),
+                          progress = false, save_best = true, kwargs...)
+    return FluxOptimizationCache(prob, opt, data; maxiters, callback, progress, save_best,
+                                 kwargs...)
+end
+
+function SciMLBase.__solve(cache::FluxOptimizationCache)
+    if cache.data != Optimization.DEFAULT_DATA
+        maxiters = length(cache.data)
+        data = cache.data
     else
-        maxiters = Optimization._check_and_convert_maxiters(maxiters)
-        data = Optimization.take(data, maxiters)
+        maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters)
+        data = Optimization.take(cache.data, maxiters)
     end
 
     # Flux is silly and doesn't have an abstract type on its optimizers, so assume
     # this is a Flux optimizer
-    θ = copy(prob.u0)
+    θ = copy(cache.u0)
     G = copy(θ)
+    opt = deepcopy(cache.opt)
 
     local x, min_err, min_θ
-    min_err = typemax(eltype(prob.u0)) #dummy variables
+    min_err = typemax(eltype(cache.u0)) #dummy variables
     min_opt = 1
-    min_θ = prob.u0
-
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
+    min_θ = cache.u0
 
     t0 = time()
-    Optimization.@withprogress progress name="Training" begin for (i, d) in enumerate(data)
-        f.grad(G, θ, d...)
-        x = f.f(θ, prob.p, d...)
-        cb_call = callback(θ, x...)
+    Optimization.@withprogress cache.solver_args.progress name="Training" begin for (i, d) in enumerate(data)
+        cache.f.grad(G, θ, d...)
+        x = cache.f.f(θ, cache.p, d...)
+        cb_call = cache.solver_args.callback(θ, x...)
         if !(typeof(cb_call) <: Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
         elseif cb_call
             break
         end
         msg = @sprintf("loss: %.3g", x[1])
-        progress && ProgressLogging.@logprogress msg i/maxiters
+        cache.solver_args.progress && ProgressLogging.@logprogress msg i/maxiters
 
-        if save_best
+        if cache.solver_args.save_best
             if first(x) < first(min_err)  #found a better solution
                 min_opt = opt
                 min_err = x
@@ -50,7 +106,7 @@ function SciMLBase.__solve(prob::OptimizationProblem, opt::Flux.Optimise.Abstrac
                 opt = min_opt
                 x = min_err
                 θ = min_θ
-                callback(θ, x...)
+                cache.solver_args.callback(θ, x...)
                 break
             end
         end
@@ -59,8 +115,7 @@ function SciMLBase.__solve(prob::OptimizationProblem, opt::Flux.Optimise.Abstrac
 
     t1 = time()
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt, θ,
-                             x[1], solve_time = t1 - t0)
+    SciMLBase.build_solution(cache, opt, θ, x[1], solve_time = t1 - t0)
     # here should be build_solution to create the output message
 end
 
diff --git a/lib/OptimizationFlux/test/runtests.jl b/lib/OptimizationFlux/test/runtests.jl
index 9f0f9c3ce..58a6b6451 100644
--- a/lib/OptimizationFlux/test/runtests.jl
+++ b/lib/OptimizationFlux/test/runtests.jl
@@ -11,10 +11,27 @@ using Test
 
     prob = OptimizationProblem(optprob, x0, _p)
 
-    sol = Optimization.solve(prob, Flux.ADAM(0.1), maxiters = 1000)
+    sol = Optimization.solve(prob, Flux.Adam(0.1), maxiters = 1000)
     @test 10 * sol.objective < l1
 
     prob = OptimizationProblem(optprob, x0, _p)
-    sol = solve(prob, Flux.ADAM(), maxiters = 1000, progress = false)
+    sol = solve(prob, Flux.Adam(), maxiters = 1000, progress = false)
     @test 10 * sol.objective < l1
+
+    @testset "cache" begin
+        objective(x, p) = (p[1] - x[1])^2
+        x0 = zeros(1)
+        p = [1.0]
+
+        prob = OptimizationProblem(OptimizationFunction(objective,
+                                                        Optimization.AutoForwardDiff()), x0,
+                                   p)
+        cache = Optimization.init(prob, Flux.Adam(0.1), maxiters = 1000)
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[1.0] atol=1e-3
+
+        cache = Optimization.reinit!(cache; p = [2.0])
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[2.0] atol=1e-3
+    end
 end
diff --git a/lib/OptimizationGCMAES/src/OptimizationGCMAES.jl b/lib/OptimizationGCMAES/src/OptimizationGCMAES.jl
index 01171eb6f..f02671d45 100644
--- a/lib/OptimizationGCMAES/src/OptimizationGCMAES.jl
+++ b/lib/OptimizationGCMAES/src/OptimizationGCMAES.jl
@@ -11,7 +11,65 @@ struct GCMAESOpt end
 SciMLBase.requiresbounds(::GCMAESOpt) = true
 SciMLBase.allowsbounds(::GCMAESOpt) = true
 
-function __map_optimizer_args(prob::OptimizationProblem, opt::GCMAESOpt;
+struct GCMAESOptimizationCache{F <: OptimizationFunction, RC, LB, UB, S, O, P, S0} <:
+       SciMLBase.AbstractOptimizationCache
+    f::F
+    reinit_cache::RC
+    lb::LB
+    ub::UB
+    sense::S
+    opt::O
+    progress::P
+    sigma0::S0
+    solver_args::NamedTuple
+end
+
+function Base.getproperty(cache::GCMAESOptimizationCache, x::Symbol)
+    if x in fieldnames(Optimization.ReInitCache)
+        return getfield(cache.reinit_cache, x)
+    end
+    return getfield(cache, x)
+end
+
+function GCMAESOptimizationCache(prob::OptimizationProblem, opt; progress, sigma0,
+                                 kwargs...)
+    reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit`
+    f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype)
+    return GCMAESOptimizationCache(f, reinit_cache, prob.lb, prob.ub, prob.sense, opt,
+                                   progress, sigma0,
+                                   NamedTuple(kwargs))
+end
+
+SciMLBase.supports_opt_cache_interface(opt::GCMAESOpt) = true
+SciMLBase.has_reinit(cache::GCMAESOptimizationCache) = true
+function SciMLBase.reinit!(cache::GCMAESOptimizationCache; p = missing, u0 = missing)
+    if p === missing && u0 === missing
+        p, u0 = cache.p, cache.u0
+    else # at least one of them has a value
+        if p === missing
+            p = cache.p
+        end
+        if u0 === missing
+            u0 = cache.u0
+        end
+        if (eltype(p) <: Pair && !isempty(p)) || (eltype(u0) <: Pair && !isempty(u0)) # one is a non-empty symbolic map
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :ps) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `p` keyword argument as a vector of values, paying attention to parameter order."))
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :states) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
+            p, u0 = SciMLBase.process_p_u0_symbolic(cache, p, u0)
+        end
+    end
+
+    cache.reinit_cache.p = p
+    cache.reinit_cache.u0 = u0
+
+    return cache
+end
+
+function __map_optimizer_args(cache::GCMAESOptimizationCache, opt::GCMAESOpt;
                               callback = nothing,
                               maxiters::Union{Number, Nothing} = nothing,
                               maxtime::Union{Number, Nothing} = nothing,
@@ -40,50 +98,56 @@ function __map_optimizer_args(prob::OptimizationProblem, opt::GCMAESOpt;
     return mapped_args
 end
 
-function SciMLBase.__solve(prob::OptimizationProblem, opt::GCMAESOpt;
-                           maxiters::Union{Number, Nothing} = nothing,
-                           maxtime::Union{Number, Nothing} = nothing,
-                           abstol::Union{Number, Nothing} = nothing,
-                           reltol::Union{Number, Nothing} = nothing,
-                           progress = false,
-                           σ0 = 0.2,
-                           kwargs...)
-    local x
-    local G = similar(prob.u0)
-
+function SciMLBase.__init(prob::OptimizationProblem, opt::GCMAESOpt;
+                          maxiters::Union{Number, Nothing} = nothing,
+                          maxtime::Union{Number, Nothing} = nothing,
+                          abstol::Union{Number, Nothing} = nothing,
+                          reltol::Union{Number, Nothing} = nothing,
+                          progress = false,
+                          σ0 = 0.2,
+                          kwargs...)
     maxiters = Optimization._check_and_convert_maxiters(maxiters)
     maxtime = Optimization._check_and_convert_maxtime(maxtime)
+    return GCMAESOptimizationCache(prob, opt; maxiters, maxtime, abstol, reltol, progress,
+                                   sigma0 = σ0, kwargs...)
+end
 
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
+function SciMLBase.__solve(cache::GCMAESOptimizationCache)
+    local x
+    local G = similar(cache.u0)
 
     _loss = function (θ)
-        x = f.f(θ, prob.p)
+        x = cache.f.f(θ, cache.p)
         return x[1]
     end
 
-    if !isnothing(f.grad)
+    if !isnothing(cache.f.grad)
         g = function (θ)
-            f.grad(G, θ)
+            cache.f.grad(G, θ)
             return G
         end
     end
 
-    opt_args = __map_optimizer_args(prob, opt, maxiters = maxiters, maxtime = maxtime,
-                                    abstol = abstol, reltol = reltol; kwargs...)
+    opt_args = __map_optimizer_args(cache, cache.opt, maxiters = cache.solver_args.maxiters,
+                                    maxtime = cache.solver_args.maxtime,
+                                    abstol = cache.solver_args.abstol,
+                                    reltol = cache.solver_args.reltol; cache.solver_args...)
 
     t0 = time()
-    if prob.sense === Optimization.MaxSense
-        opt_xmin, opt_fmin, opt_ret = GCMAES.maximize(isnothing(f.grad) ? _loss :
-                                                      (_loss, g), prob.u0, σ0, prob.lb,
-                                                      prob.ub; opt_args...)
+    if cache.sense === Optimization.MaxSense
+        opt_xmin, opt_fmin, opt_ret = GCMAES.maximize(isnothing(cache.f.grad) ? _loss :
+                                                      (_loss, g), cache.u0,
+                                                      cache.sigma0, cache.lb,
+                                                      cache.ub; opt_args...)
     else
-        opt_xmin, opt_fmin, opt_ret = GCMAES.minimize(isnothing(f.grad) ? _loss :
-                                                      (_loss, g), prob.u0, σ0, prob.lb,
-                                                      prob.ub; opt_args...)
+        opt_xmin, opt_fmin, opt_ret = GCMAES.minimize(isnothing(cache.f.grad) ? _loss :
+                                                      (_loss, g), cache.u0,
+                                                      cache.sigma0, cache.lb,
+                                                      cache.ub; opt_args...)
     end
     t1 = time()
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt,
+    SciMLBase.build_solution(cache, cache.opt,
                              opt_xmin, opt_fmin; retcode = Symbol(Bool(opt_ret)),
                              solve_time = t1 - t0)
 end
diff --git a/lib/OptimizationGCMAES/test/runtests.jl b/lib/OptimizationGCMAES/test/runtests.jl
index 506f28631..2c4c6d58d 100644
--- a/lib/OptimizationGCMAES/test/runtests.jl
+++ b/lib/OptimizationGCMAES/test/runtests.jl
@@ -18,4 +18,19 @@ using Test
                                             ub = [1.0, 1.0])
     sol = solve(prob, GCMAESOpt(), maxiters = 1000)
     @test 10 * sol.objective < l1
+
+    @testset "cache" begin
+        objective(x, p) = (p[1] - x[1])^2
+        x0 = zeros(1)
+        p = [1.0]
+
+        prob = OptimizationProblem(objective, x0, p, lb = [-10.0], ub = [10.0])
+        cache = Optimization.init(prob, GCMAESOpt())
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[1.0] atol=1e-3
+
+        cache = Optimization.reinit!(cache; p = [2.0])
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[2.0] atol=1e-3
+    end
 end
diff --git a/lib/OptimizationMOI/Project.toml b/lib/OptimizationMOI/Project.toml
index 891713082..b48fa9b6d 100644
--- a/lib/OptimizationMOI/Project.toml
+++ b/lib/OptimizationMOI/Project.toml
@@ -1,7 +1,7 @@
 name = "OptimizationMOI"
 uuid = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 authors = ["Vaibhav Dixit <vaibhavyashdixit@gmail.com> and contributors"]
-version = "0.1.8"
+version = "0.1.11"
 
 [deps]
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
@@ -10,8 +10,10 @@ Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
 
 [compat]
+Ipopt_jll = "=300.1400.400"
 MathOptInterface = "1"
 Juniper = "0.9"
 Optimization = "3.9"
diff --git a/lib/OptimizationMOI/test/runtests.jl b/lib/OptimizationMOI/test/runtests.jl
index cd3fffdea..ab437d821 100644
--- a/lib/OptimizationMOI/test/runtests.jl
+++ b/lib/OptimizationMOI/test/runtests.jl
@@ -22,9 +22,9 @@ function _test_sparse_derivatives_hs071(backend, optimizer)
     sol = solve(prob, optimizer)
     @test isapprox(sol.objective, 17.014017145179164; atol = 1e-6)
     x = [1.0, 4.7429996418092970, 3.8211499817883077, 1.3794082897556983]
-    @test isapprox(sol.minimizer, x; atol = 1e-6)
-    @test prod(sol.minimizer) >= 25.0 - 1e-6
-    @test isapprox(sum(sol.minimizer .^ 2), 40.0; atol = 1e-6)
+    @test isapprox(sol.u, x; atol = 1e-6)
+    @test prod(sol.u) >= 25.0 - 1e-6
+    @test isapprox(sum(sol.u .^ 2), 40.0; atol = 1e-6)
     return
 end
 
diff --git a/lib/OptimizationNLopt/Project.toml b/lib/OptimizationNLopt/Project.toml
index eec0472aa..b6eb4404d 100644
--- a/lib/OptimizationNLopt/Project.toml
+++ b/lib/OptimizationNLopt/Project.toml
@@ -1,7 +1,7 @@
 name = "OptimizationNLopt"
 uuid = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 authors = ["Vaibhav Dixit <vaibhavyashdixit@gmail.com> and contributors"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
diff --git a/lib/OptimizationNLopt/src/OptimizationNLopt.jl b/lib/OptimizationNLopt/src/OptimizationNLopt.jl
index 73f80349b..354e07239 100644
--- a/lib/OptimizationNLopt/src/OptimizationNLopt.jl
+++ b/lib/OptimizationNLopt/src/OptimizationNLopt.jl
@@ -8,7 +8,64 @@ using Optimization.SciMLBase
 
 SciMLBase.allowsbounds(opt::Union{NLopt.Algorithm, NLopt.Opt}) = true
 
-function __map_optimizer_args!(prob::OptimizationProblem, opt::NLopt.Opt;
+struct NLoptOptimizationCache{F <: OptimizationFunction, RC, LB, UB, S, O, P} <:
+       SciMLBase.AbstractOptimizationCache
+    f::F
+    reinit_cache::RC
+    lb::LB
+    ub::UB
+    sense::S
+    opt::O
+    progress::P
+    solver_args::NamedTuple
+end
+
+function Base.getproperty(cache::NLoptOptimizationCache, x::Symbol)
+    if x in fieldnames(Optimization.ReInitCache)
+        return getfield(cache.reinit_cache, x)
+    end
+    return getfield(cache, x)
+end
+
+function NLoptOptimizationCache(prob::OptimizationProblem, opt; progress, kwargs...)
+    reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit`
+    f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype)
+
+    return NLoptOptimizationCache(f, reinit_cache, prob.lb, prob.ub, prob.sense, opt,
+                                  progress,
+                                  NamedTuple(kwargs))
+end
+
+SciMLBase.supports_opt_cache_interface(opt::Union{NLopt.Algorithm, NLopt.Opt}) = true
+SciMLBase.has_reinit(cache::NLoptOptimizationCache) = true
+function SciMLBase.reinit!(cache::NLoptOptimizationCache; p = missing, u0 = missing)
+    if p === missing && u0 === missing
+        p, u0 = cache.p, cache.u0
+    else # at least one of them has a value
+        if p === missing
+            p = cache.p
+        end
+        if u0 === missing
+            u0 = cache.u0
+        end
+        if (eltype(p) <: Pair && !isempty(p)) || (eltype(u0) <: Pair && !isempty(u0)) # one is a non-empty symbolic map
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :ps) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `p` keyword argument as a vector of values, paying attention to parameter order."))
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :states) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
+            p, u0 = SciMLBase.process_p_u0_symbolic(cache, p, u0)
+        end
+    end
+
+    cache.reinit_cache.p = p
+    cache.reinit_cache.u0 = u0
+
+    return cache
+end
+
+function __map_optimizer_args!(cache::NLoptOptimizationCache, opt::NLopt.Opt;
                                callback = nothing,
                                maxiters::Union{Number, Nothing} = nothing,
                                maxtime::Union{Number, Nothing} = nothing,
@@ -21,12 +78,12 @@ function __map_optimizer_args!(prob::OptimizationProblem, opt::NLopt.Opt;
                                kwargs...)
     if local_method !== nothing
         if isa(local_method, NLopt.Opt)
-            if ndims(local_method) != length(prob.u0)
+            if ndims(local_method) != length(cache.u0)
                 error("Passed local NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
             end
             local_meth = local_method
         else
-            local_meth = NLopt.Opt(local_method, length(prob.u0))
+            local_meth = NLopt.Opt(local_method, length(cache.u0))
         end
 
         if !isnothing(local_options)
@@ -51,12 +108,12 @@ function __map_optimizer_args!(prob::OptimizationProblem, opt::NLopt.Opt;
         eval(Meta.parse("NLopt." * string(j.first) * "!"))(opt, j.second)
     end
 
-    if prob.ub !== nothing
-        NLopt.upper_bounds!(opt, prob.ub)
+    if cache.ub !== nothing
+        opt.upper_bounds = cache.ub
     end
 
-    if prob.lb !== nothing
-        NLopt.lower_bounds!(opt, prob.lb)
+    if cache.lb !== nothing
+        opt.lower_bounds = cache.lb
     end
 
     if !(isnothing(maxiters))
@@ -83,6 +140,7 @@ function __nlopt_status_to_ReturnCode(status::Symbol)
                              NLopt.STOPVAL_REACHED,
                              NLopt.FTOL_REACHED,
                              NLopt.XTOL_REACHED,
+                             NLopt.ROUNDOFF_LIMITED,
                          ])
         return ReturnCode.Success
     elseif status == Symbol(NLopt.MAXEVAL_REACHED)
@@ -93,7 +151,6 @@ function __nlopt_status_to_ReturnCode(status::Symbol)
                                  NLopt.OUT_OF_MEMORY,
                                  NLopt.INVALID_ARGS,
                                  NLopt.FAILURE,
-                                 NLopt.ROUNDOFF_LIMITED,
                                  NLopt.FORCED_STOP,
                              ])
         return ReturnCode.Failure
@@ -102,65 +159,80 @@ function __nlopt_status_to_ReturnCode(status::Symbol)
     end
 end
 
-function SciMLBase.__solve(prob::OptimizationProblem,
-                           opt::Union{NLopt.Algorithm, NLopt.Opt};
-                           maxiters::Union{Number, Nothing} = nothing,
-                           maxtime::Union{Number, Nothing} = nothing,
-                           local_method::Union{NLopt.Algorithm, NLopt.Opt, Nothing} = nothing,
-                           local_maxiters::Union{Number, Nothing} = nothing,
-                           local_maxtime::Union{Number, Nothing} = nothing,
-                           local_options::Union{NamedTuple, Nothing} = nothing,
-                           abstol::Union{Number, Nothing} = nothing,
-                           reltol::Union{Number, Nothing} = nothing,
-                           progress = false,
-                           callback = (args...) -> (false),
-                           kwargs...)
-    local x
-
+function SciMLBase.__init(prob::OptimizationProblem, opt::Union{NLopt.Algorithm, NLopt.Opt};
+                          maxiters::Union{Number, Nothing} = nothing,
+                          maxtime::Union{Number, Nothing} = nothing,
+                          local_method::Union{NLopt.Algorithm, NLopt.Opt, Nothing} = nothing,
+                          local_maxiters::Union{Number, Nothing} = nothing,
+                          local_maxtime::Union{Number, Nothing} = nothing,
+                          local_options::Union{NamedTuple, Nothing} = nothing,
+                          abstol::Union{Number, Nothing} = nothing,
+                          reltol::Union{Number, Nothing} = nothing,
+                          progress = false,
+                          callback = (args...) -> (false),
+                          kwargs...)
     maxiters = Optimization._check_and_convert_maxiters(maxiters)
     maxtime = Optimization._check_and_convert_maxtime(maxtime)
     local_maxiters = Optimization._check_and_convert_maxiters(local_maxiters)
     local_maxtime = Optimization._check_and_convert_maxtime(local_maxtime)
 
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
+    return NLoptOptimizationCache(prob, opt; maxiters, maxtime, local_method,
+                                  local_maxiters, local_maxtime, local_options, abstol,
+                                  reltol, progress, callback)
+end
+
+function SciMLBase.__solve(cache::NLoptOptimizationCache)
+    local x
 
     _loss = function (θ)
-        x = f.f(θ, prob.p)
-        callback(θ, x...)
+        x = cache.f.f(θ, cache.p)
+        cache.solver_args.callback(θ, x...)
         return x[1]
     end
 
     fg! = function (θ, G)
         if length(G) > 0
-            f.grad(G, θ)
+            cache.f.grad(G, θ)
         end
 
         return _loss(θ)
     end
 
-    if isa(opt, NLopt.Opt)
-        if ndims(opt) != length(prob.u0)
+    opt_setup = if isa(cache.opt, NLopt.Opt)
+        if ndims(cache.opt) != length(cache.u0)
             error("Passed NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
         end
-        opt_setup = opt
+        cache.opt
     else
-        opt_setup = NLopt.Opt(opt, length(prob.u0))
+        NLopt.Opt(cache.opt, length(cache.u0))
     end
 
-    prob.sense === Optimization.MaxSense ? NLopt.max_objective!(opt_setup, fg!) :
-    NLopt.min_objective!(opt_setup, fg!)
+    if cache.sense === Optimization.MaxSense
+        NLopt.max_objective!(opt_setup, fg!)
+    else
+        NLopt.min_objective!(opt_setup, fg!)
+    end
 
-    __map_optimizer_args!(prob, opt_setup, maxiters = maxiters, maxtime = maxtime,
-                          abstol = abstol, reltol = reltol, local_method = local_method,
-                          local_maxiters = local_maxiters, local_options = local_options;
-                          kwargs...)
+    __map_optimizer_args!(cache, opt_setup, maxiters = cache.solver_args.maxiters,
+                          maxtime = cache.solver_args.maxtime,
+                          abstol = cache.solver_args.abstol,
+                          reltol = cache.solver_args.reltol,
+                          local_method = cache.solver_args.local_method,
+                          local_maxiters = cache.solver_args.local_maxiters,
+                          local_options = cache.solver_args.local_options;
+                          cache.solver_args...)
 
     t0 = time()
-    (minf, minx, ret) = NLopt.optimize(opt_setup, prob.u0)
+    (minf, minx, ret) = NLopt.optimize(opt_setup, cache.u0)
     t1 = time()
-
     retcode = __nlopt_status_to_ReturnCode(ret)
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt, minx,
+
+    if retcode == ReturnCode.Failure
+        @warn "NLopt failed to converge: $(ret)"
+        minx = fill(NaN, length(cache.u0))
+        minf = NaN
+    end
+    SciMLBase.build_solution(cache, cache.opt, minx,
                              minf; original = opt_setup, retcode = retcode,
                              solve_time = t1 - t0)
 end
diff --git a/lib/OptimizationNLopt/test/runtests.jl b/lib/OptimizationNLopt/test/runtests.jl
index beef90b15..7158659f4 100644
--- a/lib/OptimizationNLopt/test/runtests.jl
+++ b/lib/OptimizationNLopt/test/runtests.jl
@@ -10,37 +10,65 @@ using Test
     optprob = OptimizationFunction((x, p) -> -rosenbrock(x, p), Optimization.AutoZygote())
     prob = OptimizationProblem(optprob, x0, _p; sense = Optimization.MaxSense)
     sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2))
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
     optprob = OptimizationFunction(rosenbrock, Optimization.AutoZygote())
     prob = OptimizationProblem(optprob, x0, _p)
 
     sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2))
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
+    prob = OptimizationProblem(optprob, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8])
+
     sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
-    prob = OptimizationProblem(optprob, x0, lb = [-1.0, -1.0], ub = [0.8, 0.8])
     sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
     sol = solve(prob, NLopt.Opt(:G_MLSL_LDS, 2), local_method = NLopt.Opt(:LD_LBFGS, 2),
                 maxiters = 10000)
+    @test sol.retcode == ReturnCode.MaxIters
     @test 10 * sol.objective < l1
 
-    prob = OptimizationProblem(optprob, x0)
+    prob = OptimizationProblem(optprob, x0, _p)
     sol = solve(prob, NLopt.LN_BOBYQA())
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
     sol = solve(prob, NLopt.LD_LBFGS())
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
-    prob = OptimizationProblem(optprob, x0, lb = [-1.0, -1.0], ub = [0.8, 0.8])
+    prob = OptimizationProblem(optprob, x0, _p, lb = [-1.0, -1.0], ub = [0.8, 0.8])
     sol = solve(prob, NLopt.LD_LBFGS())
+    @test sol.retcode == ReturnCode.Success
     @test 10 * sol.objective < l1
 
     sol = solve(prob, NLopt.G_MLSL_LDS(), local_method = NLopt.LD_LBFGS(),
                 local_maxiters = 10000, maxiters = 10000, population = 10)
+    @test sol.retcode == ReturnCode.MaxIters
     @test 10 * sol.objective < l1
+
+    @testset "cache" begin
+        objective(x, p) = (p[1] - x[1])^2
+        x0 = zeros(1)
+        p = [1.0]
+
+        optf = OptimizationFunction(objective, Optimization.AutoZygote())
+        prob = OptimizationProblem(optf, x0, p)
+        cache = Optimization.init(prob, NLopt.Opt(:LD_LBFGS, 1))
+        sol = Optimization.solve!(cache)
+        @test sol.retcode == ReturnCode.Success
+        @test sol.u≈[1.0] atol=1e-3
+
+        cache = Optimization.reinit!(cache; p = [2.0])
+        sol = Optimization.solve!(cache)
+        @test sol.retcode == ReturnCode.Success
+        @test sol.u≈[2.0] atol=1e-3
+    end
 end
diff --git a/lib/OptimizationOptimJL/src/OptimizationOptimJL.jl b/lib/OptimizationOptimJL/src/OptimizationOptimJL.jl
index 958c70eb8..56dddcfe0 100644
--- a/lib/OptimizationOptimJL/src/OptimizationOptimJL.jl
+++ b/lib/OptimizationOptimJL/src/OptimizationOptimJL.jl
@@ -13,7 +13,78 @@ SciMLBase.allowsbounds(opt::Optim.SimulatedAnnealing) = false
 SciMLBase.requiresbounds(opt::Optim.Fminbox) = true
 SciMLBase.requiresbounds(opt::Optim.SAMIN) = true
 
-function __map_optimizer_args(prob::OptimizationProblem,
+struct OptimJLOptimizationCache{F, RC, LB, UB, LC, UC, S, O, D, P, C} <:
+       SciMLBase.AbstractOptimizationCache
+    f::F
+    reinit_cache::RC
+    lb::LB
+    ub::UB
+    lcons::LC
+    ucons::UC
+    sense::S
+    opt::O
+    data::D
+    progress::P
+    callback::C
+    solver_args::NamedTuple
+end
+
+function Base.getproperty(cache::OptimJLOptimizationCache, x::Symbol)
+    if x in fieldnames(Optimization.ReInitCache)
+        return getfield(cache.reinit_cache, x)
+    end
+    return getfield(cache, x)
+end
+
+function OptimJLOptimizationCache(prob::OptimizationProblem, opt, data; progress, callback,
+                                  kwargs...)
+    reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit`
+    num_cons = prob.ucons === nothing ? 0 : length(prob.ucons)
+    f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype, num_cons)
+
+    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing &&
+        error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")
+
+    opt isa Optim.ConstrainedOptimizer && f.cons_j === nothing &&
+        error("This optimizer requires derivative definitions for nonlinear constraints. If the problem does not have nonlinear constraints, choose a different optimizer. Otherwise define the derivative for cons using OptimizationFunction either directly or automatically generate them with one of the autodiff backends")
+
+    return OptimJLOptimizationCache(f, reinit_cache, prob.lb, prob.ub, prob.lcons,
+                                    prob.ucons, prob.sense,
+                                    opt, data, progress, callback, NamedTuple(kwargs))
+end
+
+SciMLBase.supports_opt_cache_interface(opt::Optim.AbstractOptimizer) = true
+SciMLBase.supports_opt_cache_interface(opt::Union{Optim.Fminbox, Optim.SAMIN}) = true
+SciMLBase.supports_opt_cache_interface(opt::Optim.ConstrainedOptimizer) = true
+SciMLBase.has_reinit(cache::OptimJLOptimizationCache) = true
+function SciMLBase.reinit!(cache::OptimJLOptimizationCache; p = missing, u0 = missing)
+    if p === missing && u0 === missing
+        p, u0 = cache.p, cache.u0
+    else # at least one of them has a value
+        if p === missing
+            p = cache.p
+        end
+        if u0 === missing
+            u0 = cache.u0
+        end
+        if (eltype(p) <: Pair && !isempty(p)) || (eltype(u0) <: Pair && !isempty(u0)) # one is a non-empty symbolic map
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :ps) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `p` keyword argument as a vector of values, paying attention to parameter order."))
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :states) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
+            p, u0 = SciMLBase.process_p_u0_symbolic(cache, p, u0)
+        end
+    end
+
+    cache.reinit_cache.p = p
+    cache.reinit_cache.u0 = u0
+
+    return cache
+end
+
+function __map_optimizer_args(cache::OptimJLOptimizationCache,
                               opt::Union{Optim.AbstractOptimizer, Optim.Fminbox,
                                          Optim.SAMIN, Optim.ConstrainedOptimizer};
                               callback = nothing,
@@ -47,10 +118,15 @@ function __map_optimizer_args(prob::OptimizationProblem,
     return Optim.Options(; mapped_args...)
 end
 
-function SciMLBase.__solve(prob::OptimizationProblem,
-                           opt::Optim.AbstractOptimizer,
-                           data = Optimization.DEFAULT_DATA;
-                           kwargs...)
+function SciMLBase.__init(prob::OptimizationProblem, opt::Optim.AbstractOptimizer,
+                          data = Optimization.DEFAULT_DATA;
+                          callback = (args...) -> (false),
+                          maxiters::Union{Number, Nothing} = nothing,
+                          maxtime::Union{Number, Nothing} = nothing,
+                          abstol::Union{Number, Nothing} = nothing,
+                          reltol::Union{Number, Nothing} = nothing,
+                          progress = false,
+                          kwargs...)
     if !isnothing(prob.lb) || !isnothing(prob.ub)
         if !(opt isa Union{Optim.Fminbox, Optim.SAMIN, Optim.AbstractConstrainedOptimizer})
             if opt isa Optim.ParticleSwarm
@@ -64,34 +140,44 @@ function SciMLBase.__solve(prob::OptimizationProblem,
         end
     end
 
-    return ___solve(prob, opt, data; kwargs...)
+    maxiters = if data != Optimization.DEFAULT_DATA
+        length(data)
+    else
+        maxiters
+    end
+
+    maxiters = Optimization._check_and_convert_maxiters(maxiters)
+    maxtime = Optimization._check_and_convert_maxtime(maxtime)
+    return OptimJLOptimizationCache(prob, opt, data; callback, maxiters, maxtime, abstol,
+                                    reltol, progress,
+                                    kwargs...)
 end
 
-function ___solve(prob::OptimizationProblem, opt::Optim.AbstractOptimizer,
-                  data = Optimization.DEFAULT_DATA;
-                  callback = (args...) -> (false),
-                  maxiters::Union{Number, Nothing} = nothing,
-                  maxtime::Union{Number, Nothing} = nothing,
-                  abstol::Union{Number, Nothing} = nothing,
-                  reltol::Union{Number, Nothing} = nothing,
-                  progress = false,
-                  kwargs...)
+function SciMLBase.__solve(cache::OptimJLOptimizationCache{F, RC, LB, UB, LC, UC, S, O, D, P
+                                                           }) where {
+                                                                     F,
+                                                                     RC,
+                                                                     LB,
+                                                                     UB, LC, UC,
+                                                                     S,
+                                                                     O <:
+                                                                     Optim.AbstractOptimizer,
+                                                                     D,
+                                                                     P
+                                                                     }
     local x, cur, state
 
-    if data != Optimization.DEFAULT_DATA
-        maxiters = length(data)
-    end
-
-    cur, state = iterate(data)
+    cur, state = iterate(cache.data)
 
     function _cb(trace)
-        cb_call = opt isa Optim.NelderMead ?
-                  callback(decompose_trace(trace).metadata["centroid"], x...) :
-                  callback(decompose_trace(trace).metadata["x"], x...)
+        cb_call = cache.opt isa Optim.NelderMead ?
+                  cache.callback(decompose_trace(trace).metadata["centroid"],
+                                 x...) :
+                  cache.callback(decompose_trace(trace).metadata["x"], x...)
         if !(typeof(cb_call) <: Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process.")
         end
-        nx_itr = iterate(data, state)
+        nx_itr = iterate(cache.data, state)
         if isnothing(nx_itr)
             true
         else
@@ -100,103 +186,102 @@ function ___solve(prob::OptimizationProblem, opt::Optim.AbstractOptimizer,
         end
     end
 
-    maxiters = Optimization._check_and_convert_maxiters(maxiters)
-    maxtime = Optimization._check_and_convert_maxtime(maxtime)
-
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
-
-    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing &&
-        error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")
-
     _loss = function (θ)
-        x = f.f(θ, prob.p, cur...)
+        x = cache.f.f(θ, cache.p, cur...)
         __x = first(x)
-        return prob.sense === Optimization.MaxSense ? -__x : __x
+        return cache.sense === Optimization.MaxSense ? -__x : __x
     end
 
     fg! = function (G, θ)
         if G !== nothing
-            f.grad(G, θ, cur...)
-            if prob.sense === Optimization.MaxSense
+            cache.f.grad(G, θ, cur...)
+            if cache.sense === Optimization.MaxSense
                 G .*= false
             end
         end
         return _loss(θ)
     end
 
-    if opt isa Optim.KrylovTrustRegion
+    if cache.opt isa Optim.KrylovTrustRegion
         hv = function (H, θ, v)
-            f.hv(H, θ, v, cur...)
-            if prob.sense === Optimization.MaxSense
+            cache.f.hv(H, θ, v, cur...)
+            if cache.sense === Optimization.MaxSense
                 H .*= false
             end
         end
-        optim_f = Optim.TwiceDifferentiableHV(_loss, fg!, hv, prob.u0)
+        optim_f = Optim.TwiceDifferentiableHV(_loss, fg!, hv, cache.u0)
     else
         gg = function (G, θ)
-            f.grad(G, θ, cur...)
-            if prob.sense === Optimization.MaxSense
+            cache.f.grad(G, θ, cur...)
+            if cache.sense === Optimization.MaxSense
                 G .*= false
             end
         end
 
         hh = function (H, θ)
-            f.hess(H, θ, cur...)
-            if prob.sense === Optimization.MaxSense
+            cache.f.hess(H, θ, cur...)
+            if cache.sense === Optimization.MaxSense
                 H .*= false
             end
         end
-        F = eltype(prob.u0)
-        optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, prob.u0, real(zero(F)),
-                                            Optim.NLSolversBase.alloc_DF(prob.u0,
-                                                                         real(zero(F))),
-                                            isnothing(f.hess_prototype) ?
-                                            Optim.NLSolversBase.alloc_H(prob.u0,
-                                                                        real(zero(F))) :
-                                            convert.(F, f.hess_prototype))
+        u0_type = eltype(cache.u0)
+        optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, cache.u0,
+                                            real(zero(u0_type)),
+                                            Optim.NLSolversBase.alloc_DF(cache.u0,
+                                                                         real(zero(u0_type))),
+                                            isnothing(cache.f.hess_prototype) ?
+                                            Optim.NLSolversBase.alloc_H(cache.u0,
+                                                                        real(zero(u0_type))) :
+                                            convert.(u0_type, cache.f.hess_prototype))
     end
 
-    opt_args = __map_optimizer_args(prob, opt, callback = _cb, maxiters = maxiters,
-                                    maxtime = maxtime, abstol = abstol, reltol = reltol;
-                                    kwargs...)
+    opt_args = __map_optimizer_args(cache, cache.opt, callback = _cb,
+                                    maxiters = cache.solver_args.maxiters,
+                                    maxtime = cache.solver_args.maxtime,
+                                    abstol = cache.solver_args.abstol,
+                                    reltol = cache.solver_args.reltol;
+                                    cache.solver_args...)
 
     t0 = time()
-    opt_res = Optim.optimize(optim_f, prob.u0, opt, opt_args)
+    opt_res = Optim.optimize(optim_f, cache.u0, cache.opt, opt_args)
     t1 = time()
     opt_ret = Symbol(Optim.converged(opt_res))
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt,
+    SciMLBase.build_solution(cache, cache.opt,
                              opt_res.minimizer,
-                             prob.sense === Optimization.MaxSense ? -opt_res.minimum :
+                             cache.sense === Optimization.MaxSense ? -opt_res.minimum :
                              opt_res.minimum; original = opt_res, retcode = opt_ret,
                              solve_time = t1 - t0)
 end
 
-function ___solve(prob::OptimizationProblem, opt::Union{Optim.Fminbox, Optim.SAMIN},
-                  data = Optimization.DEFAULT_DATA;
-                  callback = (args...) -> (false),
-                  maxiters::Union{Number, Nothing} = nothing,
-                  maxtime::Union{Number, Nothing} = nothing,
-                  abstol::Union{Number, Nothing} = nothing,
-                  reltol::Union{Number, Nothing} = nothing,
-                  progress = false,
-                  kwargs...)
+function SciMLBase.__solve(cache::OptimJLOptimizationCache{F, RC, LB, UB, LC, UC, S, O, D, P
+                                                           }) where {
+                                                                     F,
+                                                                     RC,
+                                                                     LB,
+                                                                     UB, LC, UC,
+                                                                     S,
+                                                                     O <:
+                                                                     Union{
+                                                                           Optim.Fminbox,
+                                                                           Optim.SAMIN
+                                                                           },
+                                                                     D,
+                                                                     P
+                                                                     }
     local x, cur, state
 
-    if data != Optimization.DEFAULT_DATA
-        maxiters = length(data)
-    end
-
-    cur, state = iterate(data)
+    cur, state = iterate(cache.data)
 
     function _cb(trace)
-        cb_call = !(opt isa Optim.SAMIN) && opt.method == Optim.NelderMead() ?
-                  callback(decompose_trace(trace).metadata["centroid"], x...) :
-                  callback(decompose_trace(trace).metadata["x"], x...)
+        cb_call = !(cache.opt isa Optim.SAMIN) && cache.opt.method == Optim.NelderMead() ?
+                  cache.callback(decompose_trace(trace).metadata["centroid"],
+                                 x...) :
+                  cache.callback(decompose_trace(trace).metadata["x"], x...)
         if !(typeof(cb_call) <: Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process.")
         end
-        nx_itr = iterate(data, state)
+        nx_itr = iterate(cache.data, state)
         if isnothing(nx_itr)
             true
         else
@@ -205,23 +290,15 @@ function ___solve(prob::OptimizationProblem, opt::Union{Optim.Fminbox, Optim.SAM
         end
     end
 
-    maxiters = Optimization._check_and_convert_maxiters(maxiters)
-    maxtime = Optimization._check_and_convert_maxtime(maxtime)
-
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
-
-    !(opt isa Optim.ZerothOrderOptimizer) && f.grad === nothing &&
-        error("Use OptimizationFunction to pass the derivatives or automatically generate them with one of the autodiff backends")
-
     _loss = function (θ)
-        x = f.f(θ, prob.p, cur...)
+        x = cache.f.f(θ, cache.p, cur...)
         __x = first(x)
-        return prob.sense === Optimization.MaxSense ? -__x : __x
+        return cache.sense === Optimization.MaxSense ? -__x : __x
     end
     fg! = function (G, θ)
         if G !== nothing
-            f.grad(G, θ, cur...)
-            if prob.sense === Optimization.MaxSense
+            cache.f.grad(G, θ, cur...)
+            if cache.sense === Optimization.MaxSense
                 G .*= false
             end
         end
@@ -229,50 +306,52 @@ function ___solve(prob::OptimizationProblem, opt::Union{Optim.Fminbox, Optim.SAM
     end
 
     gg = function (G, θ)
-        f.grad(G, θ, cur...)
-        if prob.sense === Optimization.MaxSense
+        cache.f.grad(G, θ, cur...)
+        if cache.sense === Optimization.MaxSense
             G .*= false
         end
     end
-    optim_f = Optim.OnceDifferentiable(_loss, gg, fg!, prob.u0)
+    optim_f = Optim.OnceDifferentiable(_loss, gg, fg!, cache.u0)
 
-    opt_args = __map_optimizer_args(prob, opt, callback = _cb, maxiters = maxiters,
-                                    maxtime = maxtime, abstol = abstol, reltol = reltol;
-                                    kwargs...)
+    opt_args = __map_optimizer_args(cache, cache.opt, callback = _cb,
+                                    maxiters = cache.solver_args.maxiters,
+                                    maxtime = cache.solver_args.maxtime,
+                                    abstol = cache.solver_args.abstol,
+                                    reltol = cache.solver_args.reltol;
+                                    cache.solver_args...)
 
     t0 = time()
-    opt_res = Optim.optimize(optim_f, prob.lb, prob.ub, prob.u0, opt, opt_args)
+    opt_res = Optim.optimize(optim_f, cache.lb, cache.ub, cache.u0, cache.opt, opt_args)
     t1 = time()
     opt_ret = Symbol(Optim.converged(opt_res))
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt,
+    SciMLBase.build_solution(cache, cache.opt,
                              opt_res.minimizer, opt_res.minimum;
                              original = opt_res, retcode = opt_ret, solve_time = t1 - t0)
 end
 
-function ___solve(prob::OptimizationProblem, opt::Optim.ConstrainedOptimizer,
-                  data = Optimization.DEFAULT_DATA;
-                  callback = (args...) -> (false),
-                  maxiters::Union{Number, Nothing} = nothing,
-                  maxtime::Union{Number, Nothing} = nothing,
-                  abstol::Union{Number, Nothing} = nothing,
-                  reltol::Union{Number, Nothing} = nothing,
-                  progress = false,
-                  kwargs...)
+function SciMLBase.__solve(cache::OptimJLOptimizationCache{F, RC, LB, UB, LC, UC, S, O, D, P
+                                                           }) where {
+                                                                     F,
+                                                                     RC,
+                                                                     LB,
+                                                                     UB, LC, UC,
+                                                                     S,
+                                                                     O <:
+                                                                     Optim.ConstrainedOptimizer,
+                                                                     D,
+                                                                     P
+                                                                     }
     local x, cur, state
 
-    if data != Optimization.DEFAULT_DATA
-        maxiters = length(data)
-    end
-
-    cur, state = iterate(data)
+    cur, state = iterate(cache.data)
 
     function _cb(trace)
-        cb_call = callback(decompose_trace(trace).metadata["x"], x...)
+        cb_call = cache.callback(decompose_trace(trace).metadata["x"], x...)
         if !(typeof(cb_call) <: Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process.")
         end
-        nx_itr = iterate(data, state)
+        nx_itr = iterate(cache.data, state)
         if isnothing(nx_itr)
             true
         else
@@ -281,74 +360,70 @@ function ___solve(prob::OptimizationProblem, opt::Optim.ConstrainedOptimizer,
         end
     end
 
-    maxiters = Optimization._check_and_convert_maxiters(maxiters)
-    maxtime = Optimization._check_and_convert_maxtime(maxtime)
-
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p,
-                                          prob.ucons === nothing ? 0 : length(prob.ucons))
-
-    f.cons_j === nothing &&
-        error("This optimizer requires derivative definitions for nonlinear constraints. If the problem does not have nonlinear constraints, choose a different optimizer. Otherwise define the derivative for cons using OptimizationFunction either directly or automatically generate them with one of the autodiff backends")
-
     _loss = function (θ)
-        x = f.f(θ, prob.p, cur...)
+        x = cache.f.f(θ, cache.p, cur...)
         __x = first(x)
-        return prob.sense === Optimization.MaxSense ? -__x : __x
+        return cache.sense === Optimization.MaxSense ? -__x : __x
     end
     fg! = function (G, θ)
         if G !== nothing
-            f.grad(G, θ, cur...)
-            if prob.sense === Optimization.MaxSense
+            cache.f.grad(G, θ, cur...)
+            if cache.sense === Optimization.MaxSense
                 G .*= false
             end
         end
         return _loss(θ)
     end
     gg = function (G, θ)
-        f.grad(G, θ, cur...)
-        if prob.sense === Optimization.MaxSense
+        cache.f.grad(G, θ, cur...)
+        if cache.sense === Optimization.MaxSense
             G .*= false
         end
     end
 
     hh = function (H, θ)
-        f.hess(H, θ, cur...)
-        if prob.sense === Optimization.MaxSense
+        cache.f.hess(H, θ, cur...)
+        if cache.sense === Optimization.MaxSense
             H .*= false
         end
     end
-    F = eltype(prob.u0)
-    optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, prob.u0, real(zero(F)),
-                                        Optim.NLSolversBase.alloc_DF(prob.u0,
-                                                                     real(zero(F))),
-                                        isnothing(f.hess_prototype) ?
-                                        Optim.NLSolversBase.alloc_H(prob.u0,
-                                                                    real(zero(F))) :
-                                        convert.(F, f.hess_prototype))
+    u0_type = eltype(cache.u0)
+    optim_f = Optim.TwiceDifferentiable(_loss, gg, fg!, hh, cache.u0,
+                                        real(zero(u0_type)),
+                                        Optim.NLSolversBase.alloc_DF(cache.u0,
+                                                                     real(zero(u0_type))),
+                                        isnothing(cache.f.hess_prototype) ?
+                                        Optim.NLSolversBase.alloc_H(cache.u0,
+                                                                    real(zero(u0_type))) :
+                                        convert.(u0_type, cache.f.hess_prototype))
 
     cons_hl! = function (h, θ, λ)
         res = [similar(h) for i in 1:length(λ)]
-        f.cons_h(res, θ)
+        cache.f.cons_h(res, θ)
         for i in 1:length(λ)
             h .+= λ[i] * res[i]
         end
     end
 
-    lb = prob.lb === nothing ? [] : prob.lb
-    ub = prob.ub === nothing ? [] : prob.ub
-    optim_fc = Optim.TwiceDifferentiableConstraints(f.cons, f.cons_j, cons_hl!, lb, ub,
-                                                    prob.lcons, prob.ucons)
+    lb = cache.lb === nothing ? [] : cache.lb
+    ub = cache.ub === nothing ? [] : cache.ub
+    optim_fc = Optim.TwiceDifferentiableConstraints(cache.f.cons, cache.f.cons_j, cons_hl!,
+                                                    lb, ub,
+                                                    cache.lcons, cache.ucons)
 
-    opt_args = __map_optimizer_args(prob, opt, callback = _cb, maxiters = maxiters,
-                                    maxtime = maxtime, abstol = abstol, reltol = reltol;
-                                    kwargs...)
+    opt_args = __map_optimizer_args(cache, cache.opt, callback = _cb,
+                                    maxiters = cache.solver_args.maxiters,
+                                    maxtime = cache.solver_args.maxtime,
+                                    abstol = cache.solver_args.abstol,
+                                    reltol = cache.solver_args.reltol;
+                                    cache.solver_args...)
 
     t0 = time()
-    opt_res = Optim.optimize(optim_f, optim_fc, prob.u0, opt, opt_args)
+    opt_res = Optim.optimize(optim_f, optim_fc, cache.u0, cache.opt, opt_args)
     t1 = time()
     opt_ret = Symbol(Optim.converged(opt_res))
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt,
+    SciMLBase.build_solution(cache, cache.opt,
                              opt_res.minimizer, opt_res.minimum;
                              original = opt_res, retcode = opt_ret)
 end
diff --git a/lib/OptimizationOptimJL/test/runtests.jl b/lib/OptimizationOptimJL/test/runtests.jl
index 229dc17e1..270671a0c 100644
--- a/lib/OptimizationOptimJL/test/runtests.jl
+++ b/lib/OptimizationOptimJL/test/runtests.jl
@@ -40,11 +40,14 @@ using Test
     optprob = OptimizationFunction(rosenbrock, Optimization.AutoModelingToolkit();
                                    cons = cons)
 
-    prob = OptimizationProblem(optprob, x0, _p, lcons = [-Inf], ucons = [Inf])
+    prob = OptimizationProblem(optprob, x0, _p, lcons = [-5.0], ucons = [10.0])
     sol = solve(prob, IPNewton())
     @test 10 * sol.objective < l1
 
-    prob = OptimizationProblem(optprob, x0, _p, lcons = [-5.0], ucons = [10.0])
+    optprob = OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff();
+                                   cons = cons)
+
+    prob = OptimizationProblem(optprob, x0, _p, lcons = [-Inf], ucons = [Inf])
     sol = solve(prob, IPNewton())
     @test 10 * sol.objective < l1
 
@@ -115,4 +118,19 @@ using Test
 
     sol = solve(prob, Optim.KrylovTrustRegion())
     @test 10 * sol.objective < l1
+
+    @testset "cache" begin
+        objective(x, p) = (p[1] - x[1])^2
+        x0 = zeros(1)
+        p = [1.0]
+
+        prob = OptimizationProblem(objective, x0, p)
+        cache = Optimization.init(prob, Optim.NelderMead())
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[1.0] atol=1e-3
+
+        cache = Optimization.reinit!(cache; p = [2.0])
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[2.0] atol=1e-3
+    end
 end
diff --git a/lib/OptimizationOptimisers/src/OptimizationOptimisers.jl b/lib/OptimizationOptimisers/src/OptimizationOptimisers.jl
index 6c6f22906..f75aa353a 100644
--- a/lib/OptimizationOptimisers/src/OptimizationOptimisers.jl
+++ b/lib/OptimizationOptimisers/src/OptimizationOptimisers.jl
@@ -4,42 +4,107 @@ using Reexport, Printf, ProgressLogging
 @reexport using Optimisers, Optimization
 using Optimization.SciMLBase
 
-function SciMLBase.__solve(prob::OptimizationProblem, opt::AbstractRule,
-                           data = Optimization.DEFAULT_DATA;
-                           maxiters::Number = 0, callback = (args...) -> (false),
-                           progress = false, save_best = true, kwargs...)
-    if data != Optimization.DEFAULT_DATA
-        maxiters = length(data)
-    else
-        maxiters = Optimization._check_and_convert_maxiters(maxiters)
-        data = Optimization.take(data, maxiters)
+const OptimisersOptimizers = Union{Descent, Adam, Momentum, Nesterov, RMSProp,
+                                   AdaGrad, AdaMax, AdaDelta, AMSGrad, NAdam, RAdam, OAdam,
+                                   AdaBelief,
+                                   WeightDecay, ClipGrad, ClipNorm, OptimiserChain}
+
+struct OptimisersOptimizationCache{F <: OptimizationFunction, RC, LB, UB, S, D, O} <:
+       SciMLBase.AbstractOptimizationCache
+    f::F
+    reinit_cache::RC
+    lb::LB
+    ub::UB
+    sense::S
+    opt::O
+    data::D
+    solver_args::NamedTuple
+end
+
+function Base.getproperty(cache::OptimisersOptimizationCache, x::Symbol)
+    if x in fieldnames(Optimization.ReInitCache)
+        return getfield(cache.reinit_cache, x)
     end
+    return getfield(cache, x)
+end
 
-    θ = copy(prob.u0)
+function OptimisersOptimizationCache(prob::OptimizationProblem, opt, data; kwargs...)
+    reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit`
+    f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype)
+    return OptimisersOptimizationCache(f, reinit_cache, prob.lb, prob.ub, prob.sense, opt,
+                                       data, NamedTuple(kwargs))
+end
+
+SciMLBase.supports_opt_cache_interface(opt::OptimisersOptimizers) = true
+SciMLBase.has_reinit(cache::OptimisersOptimizationCache) = true
+function SciMLBase.reinit!(cache::OptimisersOptimizationCache; p = missing, u0 = missing)
+    if p === missing && u0 === missing
+        p, u0 = cache.p, cache.u0
+    else # at least one of them has a value
+        if p === missing
+            p = cache.p
+        end
+        if u0 === missing
+            u0 = cache.u0
+        end
+        if (eltype(p) <: Pair && !isempty(p)) || (eltype(u0) <: Pair && !isempty(u0)) # one is a non-empty symbolic map
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :ps) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `p` keyword argument as a vector of values, paying attention to parameter order."))
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :states) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
+            p, u0 = SciMLBase.process_p_u0_symbolic(cache, p, u0)
+        end
+    end
+
+    cache.reinit_cache.p = p
+    cache.reinit_cache.u0 = u0
+
+    return cache
+end
+
+function SciMLBase.__init(prob::OptimizationProblem, opt::OptimisersOptimizers,
+                          data = Optimization.DEFAULT_DATA;
+                          maxiters::Number = 0, callback = (args...) -> (false),
+                          progress = false, save_best = true, kwargs...)
+    return OptimisersOptimizationCache(prob, opt, data; maxiters, callback, progress,
+                                       save_best, kwargs...)
+end
+
+function SciMLBase.__solve(cache::OptimisersOptimizationCache)
+    if cache.data != Optimization.DEFAULT_DATA
+        maxiters = length(cache.data)
+        data = cache.data
+    else
+        maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters)
+        data = Optimization.take(cache.data, maxiters)
+    end
+    opt = cache.opt
+    θ = copy(cache.u0)
     G = copy(θ)
 
     local x, min_err, min_θ
-    min_err = typemax(eltype(real(prob.u0))) #dummy variables
+    min_err = typemax(eltype(real(cache.u0))) #dummy variables
     min_opt = 1
-    min_θ = prob.u0
+    min_θ = cache.u0
 
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
     state = Optimisers.setup(opt, θ)
 
     t0 = time()
-    Optimization.@withprogress progress name="Training" begin for (i, d) in enumerate(data)
-        f.grad(G, θ, d...)
-        x = f.f(θ, prob.p, d...)
-        cb_call = callback(θ, x...)
+    Optimization.@withprogress cache.solver_args.progress name="Training" begin for (i, d) in enumerate(data)
+        cache.f.grad(G, θ, d...)
+        x = cache.f.f(θ, cache.p, d...)
+        cb_call = cache.solver_args.callback(θ, x...)
         if !(typeof(cb_call) <: Bool)
             error("The callback should return a boolean `halt` for whether to stop the optimization process. Please see the sciml_train documentation for information.")
         elseif cb_call
             break
         end
         msg = @sprintf("loss: %.3g", x[1])
-        progress && ProgressLogging.@logprogress msg i/maxiters
+        cache.solver_args.progress && ProgressLogging.@logprogress msg i/maxiters
 
-        if save_best
+        if cache.solver_args.save_best
             if first(x) < first(min_err)  #found a better solution
                 min_opt = opt
                 min_err = x
@@ -49,7 +114,7 @@ function SciMLBase.__solve(prob::OptimizationProblem, opt::AbstractRule,
                 opt = min_opt
                 x = min_err
                 θ = min_θ
-                callback(θ, x...)
+                cache.solver_args.callback(θ, x...)
                 break
             end
         end
@@ -58,8 +123,7 @@ function SciMLBase.__solve(prob::OptimizationProblem, opt::AbstractRule,
 
     t1 = time()
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt, θ,
-                             x[1], solve_time = t1 - t0)
+    SciMLBase.build_solution(cache, opt, θ, x[1], solve_time = t1 - t0)
     # here should be build_solution to create the output message
 end
 
diff --git a/lib/OptimizationOptimisers/test/runtests.jl b/lib/OptimizationOptimisers/test/runtests.jl
index 30245a0b0..8850ae568 100644
--- a/lib/OptimizationOptimisers/test/runtests.jl
+++ b/lib/OptimizationOptimisers/test/runtests.jl
@@ -16,7 +16,7 @@ using Zygote
     @test 10 * sol.objective < l1
 
     prob = OptimizationProblem(optprob, x0, _p)
-    sol = solve(prob, Optimisers.ADAM(), maxiters = 1000, progress = false)
+    sol = solve(prob, Optimisers.Adam(), maxiters = 1000, progress = false)
     @test 10 * sol.objective < l1
 
     x0 = 2 * ones(ComplexF64, 2)
@@ -28,6 +28,23 @@ using Zygote
 
     prob = OptimizationProblem(optprob, x0, _p)
 
-    sol = solve(prob, Optimisers.ADAM(), maxiters = 1000)
+    sol = solve(prob, Optimisers.Adam(), maxiters = 1000)
     @test 10 * sol.objective < l1
+
+    @testset "cache" begin
+        objective(x, p) = (p[1] - x[1])^2
+        x0 = zeros(1)
+        p = [1.0]
+
+        prob = OptimizationProblem(OptimizationFunction(objective,
+                                                        Optimization.AutoForwardDiff()), x0,
+                                   p)
+        cache = Optimization.init(prob, Optimisers.Adam(0.1), maxiters = 1000)
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[1.0] atol=1e-3
+
+        cache = Optimization.reinit!(cache; p = [2.0])
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[2.0] atol=1e-3
+    end
 end
diff --git a/lib/OptimizationSpeedMapping/src/OptimizationSpeedMapping.jl b/lib/OptimizationSpeedMapping/src/OptimizationSpeedMapping.jl
index 8f5542d8c..d8ec43053 100644
--- a/lib/OptimizationSpeedMapping/src/OptimizationSpeedMapping.jl
+++ b/lib/OptimizationSpeedMapping/src/OptimizationSpeedMapping.jl
@@ -11,7 +11,61 @@ struct SpeedMappingOpt end
 SciMLBase.allowsbounds(::SpeedMappingOpt) = true
 SciMLBase.allowscallback(::SpeedMappingOpt) = false
 
-function __map_optimizer_args(prob::OptimizationProblem, opt::SpeedMappingOpt;
+struct SpeedMappingOptimizationCache{F <: OptimizationFunction, RC, LB, UB, O, P} <:
+       SciMLBase.AbstractOptimizationCache
+    f::F
+    reinit_cache::RC
+    lb::LB
+    ub::UB
+    opt::O
+    progress::P
+    solver_args::NamedTuple
+end
+
+function Base.getproperty(cache::SpeedMappingOptimizationCache, x::Symbol)
+    if x in fieldnames(Optimization.ReInitCache)
+        return getfield(cache.reinit_cache, x)
+    end
+    return getfield(cache, x)
+end
+
+function SpeedMappingOptimizationCache(prob::OptimizationProblem, opt; progress, kwargs...)
+    reinit_cache = Optimization.ReInitCache(prob.u0, prob.p) # everything that can be changed via `reinit`
+    f = Optimization.instantiate_function(prob.f, reinit_cache, prob.f.adtype)
+    return SpeedMappingOptimizationCache(f, reinit_cache, prob.lb, prob.ub, opt, progress,
+                                         NamedTuple(kwargs))
+end
+
+SciMLBase.supports_opt_cache_interface(opt::SpeedMappingOpt) = true
+SciMLBase.has_reinit(cache::SpeedMappingOptimizationCache) = true
+function SciMLBase.reinit!(cache::SpeedMappingOptimizationCache; p = missing, u0 = missing)
+    if p === missing && u0 === missing
+        p, u0 = cache.p, cache.u0
+    else # at least one of them has a value
+        if p === missing
+            p = cache.p
+        end
+        if u0 === missing
+            u0 = cache.u0
+        end
+        if (eltype(p) <: Pair && !isempty(p)) || (eltype(u0) <: Pair && !isempty(u0)) # one is a non-empty symbolic map
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :ps) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `p` keyword argument as a vector of values, paying attention to parameter order."))
+            hasproperty(cache.f, :sys) && hasfield(typeof(cache.f.sys), :states) ||
+                throw(ArgumentError("This cache does not support symbolic maps with `remake`, i.e. it does not have a symbolic origin." *
+                                    " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
+            p, u0 = SciMLBase.process_p_u0_symbolic(cache, p, u0)
+        end
+    end
+
+    cache.reinit_cache.p = p
+    cache.reinit_cache.u0 = u0
+
+    return cache
+end
+
+function __map_optimizer_args(cache::SpeedMappingOptimizationCache, opt::SpeedMappingOpt;
                               callback = nothing,
                               maxiters::Union{Number, Nothing} = nothing,
                               maxtime::Union{Number, Nothing} = nothing,
@@ -41,39 +95,44 @@ function __map_optimizer_args(prob::OptimizationProblem, opt::SpeedMappingOpt;
     return mapped_args
 end
 
-function SciMLBase.__solve(prob::OptimizationProblem, opt::SpeedMappingOpt;
-                           maxiters::Union{Number, Nothing} = nothing,
-                           maxtime::Union{Number, Nothing} = nothing,
-                           abstol::Union{Number, Nothing} = nothing,
-                           reltol::Union{Number, Nothing} = nothing,
-                           progress = false,
-                           kwargs...)
-    local x
-
-    maxiters = Optimization._check_and_convert_maxiters(maxiters)
-    maxtime = Optimization._check_and_convert_maxtime(maxtime)
+function SciMLBase.__init(prob::OptimizationProblem, opt::SpeedMappingOpt;
+                          maxiters::Union{Number, Nothing} = nothing,
+                          maxtime::Union{Number, Nothing} = nothing,
+                          abstol::Union{Number, Nothing} = nothing,
+                          reltol::Union{Number, Nothing} = nothing,
+                          progress = false,
+                          kwargs...)
+    return SpeedMappingOptimizationCache(prob, opt; maxiters, maxtime, abstol, reltol,
+                                         progress, kwargs...)
+end
 
-    f = Optimization.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)
+function SciMLBase.__solve(cache::SpeedMappingOptimizationCache)
+    local x
 
     _loss = function (θ)
-        x = f.f(θ, prob.p)
+        x = cache.f.f(θ, cache.p)
         return first(x)
     end
 
-    if isnothing(f.grad)
+    if isnothing(cache.f.grad)
         @info "SpeedMapping's ForwardDiff AD backend is used to calculate the gradient information."
     end
 
-    opt_args = __map_optimizer_args(prob, opt, maxiters = maxiters, maxtime = maxtime,
-                                    abstol = abstol, reltol = reltol; kwargs...)
+    maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters)
+    maxtime = Optimization._check_and_convert_maxtime(cache.solver_args.maxtime)
+    opt_args = __map_optimizer_args(cache, cache.opt, maxiters = maxiters,
+                                    maxtime = maxtime,
+                                    abstol = cache.solver_args.abstol,
+                                    reltol = cache.solver_args.reltol; cache.solver_args...)
 
     t0 = time()
-    opt_res = SpeedMapping.speedmapping(prob.u0; f = _loss, (g!) = f.grad, lower = prob.lb,
-                                        upper = prob.ub, opt_args...)
+    opt_res = SpeedMapping.speedmapping(cache.u0; f = _loss, (g!) = cache.f.grad,
+                                        lower = cache.lb,
+                                        upper = cache.ub, opt_args...)
     t1 = time()
     opt_ret = Symbol(opt_res.converged)
 
-    SciMLBase.build_solution(SciMLBase.DefaultOptimizationCache(prob.f, prob.p), opt,
+    SciMLBase.build_solution(cache, cache.opt,
                              opt_res.minimizer, _loss(opt_res.minimizer);
                              original = opt_res, retcode = opt_ret, solve_time = t1 - t0)
 end
diff --git a/lib/OptimizationSpeedMapping/test/runtests.jl b/lib/OptimizationSpeedMapping/test/runtests.jl
index de501b578..33cb82561 100644
--- a/lib/OptimizationSpeedMapping/test/runtests.jl
+++ b/lib/OptimizationSpeedMapping/test/runtests.jl
@@ -23,4 +23,19 @@ using Test
     prob = OptimizationProblem(f, x0, _p; lb = [-1.0, -1.0], ub = [1.5, 1.5])
     sol = solve(prob, SpeedMappingOpt())
     @test 10 * sol.objective < l1
+
+    @testset "cache" begin
+        objective(x, p) = (p[1] - x[1])^2
+        x0 = zeros(1)
+        p = [1.0]
+
+        prob = OptimizationProblem(objective, x0, p)
+        cache = Optimization.init(prob, SpeedMappingOpt())
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[1.0] atol=1e-3
+
+        cache = Optimization.reinit!(cache; p = [2.0])
+        sol = Optimization.solve!(cache)
+        @test sol.u≈[2.0] atol=1e-3
+    end
 end
diff --git a/src/function/finitediff.jl b/src/function/finitediff.jl
index a7b97a280..3903181ff 100644
--- a/src/function/finitediff.jl
+++ b/src/function/finitediff.jl
@@ -161,7 +161,7 @@ function instantiate_function(f, cache::ReInitCache,
 
     if f.hv === nothing
         hv = function (H, θ, v, args...)
-            res = ArrayInterfaceCore.zeromatrix(θ)
+            res = ArrayInterface.zeromatrix(θ)
             hess(res, θ, args...)
             H .= res * v
         end
diff --git a/src/function/forwarddiff.jl b/src/function/forwarddiff.jl
index fb36743ea..2df0769d3 100644
--- a/src/function/forwarddiff.jl
+++ b/src/function/forwarddiff.jl
@@ -134,7 +134,7 @@ function instantiate_function(f::OptimizationFunction{true}, cache::ReInitCache,
 
     if f.hv === nothing
         hv = function (H, θ, v, args...)
-            res = ArrayInterfaceCore.zeromatrix(θ)
+            res = ArrayInterface.zeromatrix(θ)
             hess(res, θ, args...)
             H .= res * v
         end
diff --git a/src/function/function.jl b/src/function/function.jl
index bd413992b..2d6810dbf 100644
--- a/src/function/function.jl
+++ b/src/function/function.jl
@@ -43,7 +43,8 @@ function that is not defined, an error is thrown.
 For more information on the use of automatic differentiation, see the
 documentation of the `AbstractADType` types.
 """
-function instantiate_function(f, x, ::AbstractADType, p, num_cons = 0)
+function instantiate_function(f, x, ::AbstractADType,
+                              p, num_cons = 0)
     grad = f.grad === nothing ? nothing : (G, x, args...) -> f.grad(G, x, p, args...)
     hess = f.hess === nothing ? nothing : (H, x, args...) -> f.hess(H, x, p, args...)
     hv = f.hv === nothing ? nothing : (H, x, v, args...) -> f.hv(H, x, v, p, args...)
diff --git a/src/function/mtk.jl b/src/function/mtk.jl
index 24dbff233..2a67324f8 100644
--- a/src/function/mtk.jl
+++ b/src/function/mtk.jl
@@ -109,7 +109,7 @@ end
 
 function instantiate_function(f, cache::ReInitCache,
                               adtype::AutoModelingToolkit, num_cons = 0)
-    p = isnothing(p) ? SciMLBase.NullParameters() : cache.p
+    p = isnothing(cache.p) ? SciMLBase.NullParameters() : cache.p
 
     sys = ModelingToolkit.modelingtoolkitize(OptimizationProblem(f, cache.u0, cache.p;
                                                                  lcons = fill(0.0,
@@ -132,33 +132,7 @@ function instantiate_function(f, cache::ReInitCache,
         H .= res * v
     end
 
-    expr = symbolify(ModelingToolkit.Symbolics.toexpr(ModelingToolkit.equations(sys)))
-    pairs_arr = if p isa SciMLBase.NullParameters
-        [Symbol(_s) => Expr(:ref, :x, i)
-         for (i, _s) in enumerate(ModelingToolkit.states(sys))]
-    else
-        vcat([Symbol(_s) => Expr(:ref, :x, i)
-              for (i, _s) in enumerate(ModelingToolkit.states(sys))],
-             [Symbol(_p) => Expr(:ref, :p, i)
-              for (i, _p) in enumerate(ModelingToolkit.parameters(sys))])
-    end
-    rep_pars_vals!(expr, pairs_arr)
-
-    if f.cons === nothing
-        cons = nothing
-        cons_exprs = nothing
-    else
-        cons = (res, θ) -> f.cons(res, θ, cache.p)
-        cons_oop = (x, p) -> (_res = zeros(eltype(x), num_cons); f.cons(_res, x, p); _res)
-
-        cons_sys = ModelingToolkit.modelingtoolkitize(NonlinearProblem(cons_oop, cache.u0, cache.p)) # 0 == f(x)
-        cons_eqs_lhss = ModelingToolkit.lhss(ModelingToolkit.Symbolics.canonical_form.(ModelingToolkit.equations(cons_sys))) # -f(x) == 0
-        cons_exprs = map(cons_eqs_lhss) do lhs
-            e = symbolify(ModelingToolkit.Symbolics.toexpr(-lhs)) # 0 == f(x)
-            rep_pars_vals!(e, pairs_arr)
-            return Expr(:call, :(==), e, :0)
-        end
-    end
+    cons = (res, θ) -> f.cons(res, θ, cache.p)
 
     cons_j = (J, θ) -> f.cons_j(J, θ, cache.p)
 
diff --git a/test/minibatch.jl b/test/minibatch.jl
index 07090857b..b0937dace 100644
--- a/test/minibatch.jl
+++ b/test/minibatch.jl
@@ -59,7 +59,7 @@ optfun = OptimizationFunction((θ, p, batch, time_batch) -> loss_adjoint(θ, bat
                               Optimization.AutoZygote())
 optprob = OptimizationProblem(optfun, pp)
 using IterTools: ncycle
-res1 = Optimization.solve(optprob, Optimisers.ADAM(0.05), ncycle(train_loader, numEpochs),
+res1 = Optimization.solve(optprob, Optimisers.Adam(0.05), ncycle(train_loader, numEpochs),
                           callback = callback, maxiters = numEpochs)
 @test 10res1.objective < l1
 
@@ -68,7 +68,7 @@ optfun = OptimizationFunction((θ, p, batch, time_batch) -> loss_adjoint(θ, bat
                               Optimization.AutoForwardDiff())
 optprob = OptimizationProblem(optfun, pp)
 using IterTools: ncycle
-res1 = Optimization.solve(optprob, Optimisers.ADAM(0.05), ncycle(train_loader, numEpochs),
+res1 = Optimization.solve(optprob, Optimisers.Adam(0.05), ncycle(train_loader, numEpochs),
                           callback = callback, maxiters = numEpochs)
 @test 10res1.objective < l1
 
@@ -77,7 +77,7 @@ optfun = OptimizationFunction((θ, p, batch, time_batch) -> loss_adjoint(θ, bat
                               Optimization.AutoModelingToolkit())
 optprob = OptimizationProblem(optfun, pp)
 using IterTools: ncycle
-@test_broken res1 = Optimization.solve(optprob, Optimisers.ADAM(0.05),
+@test_broken res1 = Optimization.solve(optprob, Optimisers.Adam(0.05),
                                        ncycle(train_loader, numEpochs),
                                        callback = callback, maxiters = numEpochs)
 # @test 10res1.objective < l1
@@ -96,6 +96,6 @@ optfun = OptimizationFunction((θ, p, batch, time_batch) -> loss_adjoint(θ, bat
                               grad = loss_grad)
 optprob = OptimizationProblem(optfun, pp)
 using IterTools: ncycle
-res1 = Optimization.solve(optprob, Optimisers.ADAM(0.05), ncycle(train_loader, numEpochs),
+res1 = Optimization.solve(optprob, Optimisers.Adam(0.05), ncycle(train_loader, numEpochs),
                           callback = callback, maxiters = numEpochs)
 @test 10res1.objective < l1