Error with ForwardDiff no method matching Float64

I am trying to use forwardDiff to compute the gradient of a function in which I am solving an ODE, but I am getting an error that there is no method matching Float64. I have been stumped with this error for awile and would greatly appreciate any assistance that could be provided.

using DifferentialEquations
using ForwardDiff


function CostFun(x)
    
    function SpringEqu!(du, u, x, t)
        du[1] = u[2]
        du[2] = -(x[1] / x[3]) * u[2] - (x[2] / x[3]) * u[1] + 50 / x[3]
    end
    
    u0 = [2.0, 0.0]
    tspan = (0.0, 1.0)
    prob = ODEProblem(SpringEqu!, u0, tspan, x)
    sol = solve(prob)


    Simpos = zeros(length(sol.t))
    Simvel = zeros(length(sol.t))
    tout = zeros(length(sol.t))
    for i = 1:length(sol.t)
        tout[i] = sol.t[i]
        Simpos[i] = sol[1, i]
        Simvel[i] = sol[2, i]
    end

    totalCost = sum(Simpos)
    return totalCost
end


xin = [2000.0, 20000.0, 80.0]

g = ForwardDiff.gradient(CostFun, xin)

The error from the function is the following:

LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3})
Closest candidates are:
  Float64(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  Float64(::T) where T<:Number at boot.jl:715
  Float64(::Int8) at float.jl:60
  ...
Stacktrace:
 [1] convert(::Type{Float64}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3}) 
at .\number.jl:7
 [2] setindex!(::Array{Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3}, ::Int64) at .\array.jl:826
 [3] CostFun(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3},1}) at C:\Users\spenc\github\MetaOpt\autodiffTest.jl:23
 [4] vector_mode_dual_eval at C:\Users\spenc\.juliapro\JuliaPro_v1.4.1-1\packages\ForwardDiff\cXTw0\src\apiutils.jl:37 [inlined]
 [5] vector_mode_gradient(::typeof(CostFun), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3},1}}) at C:\Users\spenc\.juliapro\JuliaPro_v1.4.1-1\packages\ForwardDiff\cXTw0\src\gradient.jl:97
 [6] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3},1}}, ::Val{true}) at C:\Users\spenc\.juliapro\JuliaPro_v1.4.1-1\packages\ForwardDiff\cXTw0\src\gradient.jl:17
 [7] gradient(::Function, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(CostFun),Float64},Float64,3},1}}) at C:\Users\spenc\.juliapro\JuliaPro_v1.4.1-1\packages\ForwardDiff\cXTw0\src\gradient.jl:15 (repeats 2 times)
 [8] top-level scope at C:\Users\spenc\github\MetaOpt\autodiffTest.jl:34
in expression starting at C:\Users\spenc\github\MetaOpt\autodiffTest.jl:34

You need to allow numbers of type ForwardDiff.Dual to propagate through CostFun, something like this:

julia> function CostFun(x::AbstractVector{T}) where T
           
           function SpringEqu!(du, u, x, t)
               du[1] = u[2]
               du[2] = -(x[1] / x[3]) * u[2] - (x[2] / x[3]) * u[1] + 50 / x[3]
           end
           
           u0 = T[2.0, 0.0]
           tspan = (0.0, 1.0)
           prob = ODEProblem(SpringEqu!, u0, tspan, x)
           sol = solve(prob)


           Simpos = zeros(T, length(sol.t))
           Simvel = zeros(T, length(sol.t))
           tout = zeros(T, length(sol.t))
           for i = 1:length(sol.t)
               tout[i] = sol.t[i]
               Simpos[i] = sol[1, i]
               Simvel[i] = sol[2, i]
           end

           totalCost = sum(Simpos)
           return totalCost
       end
CostFun (generic function with 2 methods)

julia> g = ForwardDiff.gradient(CostFun, xin)
3-element Array{Float64,1}:
  0.0013186230010380866
 -0.00014123587229879185
  0.001994165580567926
6 Likes

Thanks for your help.

Why did you make this change Simpos = zeros(T, length(sol.t))? I am having trouble figuring out why this was done.

It needs to be be able to fit the dual numbers, otherwise you’ll get the same error as before.

Perhaps clearer without gradient:

using ForwardDiff: Dual
xin = [2000.0 + Dual(0, (1,0,0)), 20000.0 + Dual(0, (0,1,0)), 80.0]
T = eltype(xin) # Dual{Nothing,Float64,3}
u0 = T[2.0, 0.0]
tspan = (0.0, 1.0)
prob = ODEProblem(SpringEqu!, u0, tspan, xin)
sol = solve(prob)

gives

julia> sol[1, :]
19-element Array{Dual{Nothing,Float64,3},1}:
  Dual{Nothing}(2.0,0.0,0.0,0.0)
  Dual{Nothing}(1.9999999989997657,8.341478310022593e-18,-5.0074311889727496e-14,0.0)
  Dual{Nothing}(1.9999999506887869,2.887412874324848e-15,-2.4686464570508713e-12,0.0)
...