function f(u, p, t)
    return -u
end

prob = ODEProblem(f, ones(5), (0.0, 100.0))
prob = SteadyStateProblem(prob)
 solve(prob, DynamicSS(Tsit5()))

type Tsit5ConstantCache has no field tmp

Stacktrace:
 [1] allDerivPass at /home/henrik/.julia/v0.6/DiffEqCallbacks/src/terminatesteadystate.jl:11 [inlined]
 [2] #48 at /home/henrik/.julia/v0.6/DiffEqCallbacks/src/terminatesteadystate.jl:25 [inlined]
 [3] apply_discrete_callback! at /home/henrik/.julia/v0.6/OrdinaryDiffEq/src/callbacks.jl:189 [inlined]
 [4] handle_callbacks!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5,Array{Float64,1},Float64,Void,Float64,Float64,Float64,Array{Array{Float64,1},1},DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,false,Void,#f,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},#f,Void,OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqCallbacks.##48#50{Float64,Float64,DiffEqCallbacks.#allDerivPass},DiffEqCallbacks.##49#51,DiffEqBase.#INITIALIZE_DEFAULT}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1}}) at /home/henrik/.julia/v0.6/OrdinaryDiffEq/src/integrators/integrator_utils.jl:306
 [5] loopfooter!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5,Array{Float64,1},Float64,Void,Float64,Float64,Float64,Array{Array{Float64,1},1},DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,false,Void,#f,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},#f,Void,OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqCallbacks.##48#50{Float64,Float64,DiffEqCallbacks.#allDerivPass},DiffEqCallbacks.##49#51,DiffEqBase.#INITIALIZE_DEFAULT}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1}}) at /home/henrik/.julia/v0.6/OrdinaryDiffEq/src/integrators/integrator_utils.jl:264
 [6] solve!(::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5,Array{Float64,1},Float64,Void,Float64,Float64,Float64,Array{Array{Float64,1},1},DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,false,Void,#f,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#f,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},#f,Void,OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,DiffEqBase.#ODE_DEFAULT_NORM,DiffEqBase.CallbackSet{Tuple{},Tuple{DiffEqBase.DiscreteCallback{DiffEqCallbacks.##48#50{Float64,Float64,DiffEqCallbacks.#allDerivPass},DiffEqCallbacks.##49#51,DiffEqBase.#INITIALIZE_DEFAULT}}},DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN,DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE,DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK,DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Void,Void,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,1}}) at /home/henrik/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:342
 [7] #solve#1263(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Array{Float64,1},Float64,false,Void,#f,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/henrik/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:7
 [8] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{Array{Float64,1},Float64,false,Void,#f,Void,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)
 [9] #solve#16(::Array{Any,1}, ::Function, ::DiffEqBase.SteadyStateProblem{Array{Float64,1},false,Void,Void,#f,UniformScaling{Int64}}, ::SteadyStateDiffEq.DynamicSS{OrdinaryDiffEq.Tsit5,Float64,Float64}) at /home/henrik/.julia/v0.6/SteadyStateDiffEq/src/solve.jl:45
 [10] solve(::DiffEqBase.SteadyStateProblem{Array{Float64,1},false,Void,Void,#f,UniformScaling{Int64}}, ::SteadyStateDiffEq.DynamicSS{OrdinaryDiffEq.Tsit5,Float64,Float64}) at /home/henrik/.julia/v0.6/SteadyStateDiffEq/src/solve.jl:42`

using FiniteDiff
using ForwardDiff
using Optimization
using OptimizationOptimJL
using NonlinearSolve
using DifferentialEquations
using SteadyStateDiffEq
using Plots

function sys!(du, u, (p,p_fix,kf,plateau,kd,start_time,end_time), t)
    #@show t,forcing_function(t,plateau,kd,start_time,end_time),p,p_fix,kf
    du[1] = kf[1]*(1.0 - forcing_function(t,plateau,kd,start_time,end_time)) - p[1]*u[1]
    du[2] = p[1]*u[1] - p_fix[1]*u[2]
    du[3] = p_fix[1]*u[2] - p[2]*u[3]
    if t > 50
     #   @show u
      #  @show du
      @show t
    end
    return nothing
end

function forcing_function(t,plateau,kd,start_time,end_time)
    if t < start_time
        return plateau*t
    else
        if t <= end_time
            return plateau
        else
            if t > end_time && t < Inf
                return plateau*(exp(-kd*(t-end_time)))
            else
                return 0.0
            end
        end
    end
end



function set_kf_outer(probss,u03)
    function set_kf_inner(kf_guess,(u0,p_guess,p_fix,plateau,kd,start_time,end_time))
        # Solve for steady state with current p_guess and an initial guess of kf and u0
        plateau = 0.0
        sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol = 1e-14,reltol = 1e-14))
        # Compare kf at steady state to measured u0[3]
        # @show sol.u
        # @show sol.resid
        # @show sol.prob
        @show (sol[3]-u03)^2
        return (sol[3]-u03)^2
    end
end

function main()
    #### Parameters and state variables
    kf0= [10.0] #true value of zero order production rate
    p0 = [1.0,0.25,0.5] #true values of parameters
    p_fix = [p0[2]] # a parameter that is known, note that the system is otherwise unidentifiable
    p_var = [p0[1],p0[3]] # Unknown parameters
    u0_guess = [1000.0,750.0,20.0] #initial guess at steady state of state variables, except u0[3] is known to be 20.0
    u03 = 20.0 # measured value
    plateau = 50.0 # Zero order production rate changes at time = 0, reaches a fractional change given by plateau at start_time
    kd = 1.0 # Rate constant for exponential decay back to original value starting at end_time
    start_time = 1.0
    end_time = 10.0

    #### Set up ODE problem
    t_end = 36.0 # Length of simulation
    tspan = (0.0,t_end) # timespan of simulation
    probODE = ODEProblem(sys!,u0_guess,tspan,(p_var,p_fix,kf0,plateau,kd,start_time,end_time))
    
    #### Set up steady state problem
    probss = SteadyStateProblem(probODE)
    sol_ss = solve(probss)
    u0_ss = sol_ss.u
    @show "u0_ss true",sol_ss.u
    
    #### Initial guesses
    p_guess = [0.1,3.0] #inital guess at values in p_val
    kf_guess = [100.0] #initial guess at kf0
    

    #### Set up optimization problem to determine kf for current value of p_guess
    set_kf = set_kf_outer(probss,u03)
    optkf = OptimizationFunction(set_kf,Optimization.AutoForwardDiff())
    prob_set_kf = OptimizationProblem(optkf,kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    
    @show FiniteDiff.finite_difference_gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess)
    @show ForwardDiff.gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess )

    return nothing
end

main() ```

using FiniteDiff
using ForwardDiff
using Optimization
using OptimizationOptimJL
using NonlinearSolve
using DifferentialEquations
using SteadyStateDiffEq
using Plots

function sys!(du, u, (p,p_fix,kf,plateau,kd,start_time,end_time), t)
    #@show t,forcing_function(t,plateau,kd,start_time,end_time),p,p_fix,kf
    du[1] = kf[1]*(1.0 - forcing_function(t,plateau,kd,start_time,end_time)) - p[1]*u[1]
    du[2] = p[1]*u[1] - p_fix[1]*u[2]
    du[3] = p_fix[1]*u[2] - p[2]*u[3]
    
     #   @show u
      #  @show du
      @show t,dt
    
    return nothing
end

function forcing_function(t,plateau,kd,start_time,end_time)
    if t < start_time
        return plateau*t
    else
        if t <= end_time
            return plateau
        else
            if t > end_time && t < Inf
                return plateau*(exp(-kd*(t-end_time)))
            else
                return 0.0
            end
        end
    end
end



function set_kf_outer(probss,u03)
    function set_kf_inner(kf_guess,(u0,p_guess,p_fix,plateau,kd,start_time,end_time))
        # Solve for steady state with current p_guess and an initial guess of kf and u0
        plateau = 0.0
        sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol = 1e-16,reltol = 1e-16))
        # Compare kf at steady state to measured u0[3]
        # @show sol.u
        # @show sol.resid
        # @show sol.prob
        @show (sol[3]-u03)^2
        return (sol[3]-u03)^2
    end
end

function main()
    #### Parameters and state variables
    kf0= [10.0] #true value of zero order production rate
    p0 = [1.0,0.25,0.5] #true values of parameters
    p_fix = [p0[2]] # a parameter that is known, note that the system is otherwise unidentifiable
    p_var = [p0[1],p0[3]] # Unknown parameters
    u0_guess = [1000.0,750.0,20.0] #initial guess at steady state of state variables, except u0[3] is known to be 20.0
    u03 = 20.0 # measured value
    plateau = 0.5 # Zero order production rate changes at time = 0, reaches a fractional change given by plateau at start_time
    kd = 1.0 # Rate constant for exponential decay back to original value starting at end_time
    start_time = 1.0
    end_time = 10.0

    #### Set up ODE problem
    t_end = 36.0 # Length of simulation
    tspan = (0.0,t_end) # timespan of simulation
    probODE = ODEProblem(sys!,u0_guess,tspan,(p_var,p_fix,kf0,plateau,kd,start_time,end_time))
    
    #### Set up steady state problem
    probss = SteadyStateProblem(probODE)
    sol_ss = solve(probss)
    u0_ss = sol_ss.u
    @show "u0_ss true",sol_ss.u
    
    #### Initial guesses
    p_guess = [0.1,3.0] #inital guess at values in p_val
    kf_guess = [100.0] #initial guess at kf0
    

    #### Set up optimization problem to determine kf for current value of p_guess
    set_kf = set_kf_outer(probss,u03)
    optkf = OptimizationFunction(set_kf,Optimization.AutoForwardDiff())
    prob_set_kf = OptimizationProblem(optkf,kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    
    
    @show FiniteDiff.finite_difference_gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess)
    @show ForwardDiff.gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess )
    # @show set_kf(kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
    # @show set_kf(kf_guess,(u0_guess,[p_guess[1]*1.01,p_guess[2]],p_fix,plateau,kd,start_time,end_time))
    # @show set_kf(kf_guess,(u0_guess,[p_guess[1],p_guess[2]*1.01],p_fix,plateau,kd,start_time,end_time))
    
   

    return nothing
end

main()

using DifferentialEquations, Optimization, OptimizationPolyalgorithms, OptimizationOptimJL,
      SciMLSensitivity, Zygote, Plots, BenchmarkTools,ForwardDiff,NLsolve, Statistics
      

function SimpleKinetics!(du, u, p, t)
    x, y,x2,y2 = u
    β, δ = p

    du[1] =   (1.0-0.1*t) - β*x
    du[2] =  -δ*y + β*x
    du[3] =   0.1*t - β*x2
    du[4] =  -δ*y2 + β*x2
end

function SimpleKineticsSS(du, u, p, t)
    t = 0.0
    SimpleKinetics!(du, u, p, t)
end

function loss_wrapper(prob,tsteps,probSS)
    function loss_closure(p)
        u0=solve(remake(probSS,p=p),SSRootfind())

        sol = solve(remake(prob,u0=u0,p=p), Tsit5(), saveat = tsteps)

        loss = sum(abs2, sol.-0.5) 
        @show loss
        return loss, sol
    end
end

callback = function (p, l, pred)
  
  plt = plot(pred, ylim = (0, 6))
  display(plt)
  
  return false
end


function main()
    tspan = (0.0, 10.0)
    tsteps = 0.0:0.1:10.0
    p = [1.0, 3.0]
    u0 = [2.0, 1.0,0.0,0.0]
   
    prob = ODEProblem(SimpleKinetics!, u0, tspan, p)
   
    probSS = SteadyStateProblem{true}(SimpleKineticsSS, u0, p)
 
    solODE = solve(remake(prob,u0=u0,p=p), Tsit5(), saveat = tsteps)
    display(plot(solODE))
    
    u0=solve(probSS,SSRootfind())
    solODE = solve(remake(prob,u0=u0,p=p), Tsit5(), saveat = tsteps)
    display(plot(solODE))

    loss = loss_wrapper(prob,tsteps,probSS)

    adtype = Optimization.AutoForwardDiff()
    optf = Optimization.OptimizationFunction((x,p)->loss(x), adtype)
    optprob = Optimization.OptimizationProblem(optf, p)

    sol = Optimization.solve(optprob, PolyOpt(), callback = callback, maxiters = 100)
                                   
    p =sol.u
    @show p
    u0=solve(remake(probSS,p=p),SSRootfind())
    solODE = solve(remake(prob,u0=u0,p=p), Tsit5(), saveat = tsteps)
    display(plot(solODE))

    @show sum(abs2, solODE.-0.5) 
end

main()