Parameter estimation for Fractional differential equations

Hi, I am a new to Julia. I have been trying to fit my FDE model to a data set. I was trying to do it in MATLAB but there are no stiff solvers for fractional differential equations (not in my knowledge). Everything seems to run smoothly, i have have tried to see if the solver is able to sol or not and it provides the appropriate solution for some parameter set. Now when i try to build a cost function using build_loss_objective(), there seems to be a problem there. The code doesn’t throw an error at this stage though. I supply this cost function to Optimization.OptimizationProblem() and try to solve the problem of estimation using BFGS(). At this point, it throws an error which i cannot seem to understand how to fix. Please help me decipher that error.
Thank you.
Here is my MWE:

using DifferentialEquations, RecursiveArrayTools, Plots, DiffEqParamEstim, FractionalDiffEq
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO

p=[2.1528105468750001, 501317.0357499123201706, 0.0062344791199942,0.0443565286589824,0.0679359657330613,0.0294691281604861,0.0003138206637042,0.0002216919579355,0.9]

function SIQR(du, u, p, t)
    r, K, b, a1, b1, g, mu, d = p
    du[1] = r*u[1]*(1 - u[1]/K) -(b*u[2]*u[1])/(1 + a1*u[1] +b1*u[2])
    du[2] = (b*u[2]*u[1])/(1 + a1*u[1] +b1*u[2]) - (g+mu)*u[2]
    du[3] = g*u[2] -(d+mu)*u[3] 
    du[4] = d*u[3] -mu*u[4]
end
α = p[9] .* [1.0,1.0,1.0,1.0]
u1_0 = [1407600000.0, 81441.0, 43921.0, 11474683.0]
t1span = (0, 100)
prob1 = FODEProblem(SIQR, α, u1_0, t1span, p)

dataAux1 = [81441 89019 92998 103793 96557 115269 126315 131893 144829 152682 169914 160694 185248 199569 216850 233943 260778 275306 256947 294290 315802 332503 345147 349313 354531 319435 362902 379459 386888 402110]
data1 = collect(Float64,dataAux1)
t1 = collect(range(0, stop = 29, length = 30))

function loss(p)
    sol = solve(prob1, PIRect(),dt=0.01,p=p,saveat = t1)
    loss = sum(abs2, sol(t1)[3,:] .- data1')
    return loss
end

lowerbounds = collect(Float64,[-10,100000,0,0,0,0,0,0,0])
upperbounds = collect(Float64,[10,1407600000,1,1,1,1,1,1,1])

cost_function = build_loss_objective(prob1, PIRect(),dt=0.01, loss,
                                     Optimization.AutoForwardDiff(),
                                     maxiters = 10000, verbose = false)
optprob = Optimization.OptimizationProblem(cost_function, [2.1528105468750001, 501317.0357499123201706, 0.0062344791199942,0.0443565286589824,0.0679359657330613,0.0294691281604861,0.0003138206637042,0.0002216919579355,0.8],lb = lowerbounds,ub=upperbounds)
result_bfgs = solve(optprob, BFGS())

My initial guess was maybe i am trying to estimate the fractional order \alpha so i remade the code with a fixed \alpha, got the same error, as follows:

ERROR: Non-standard number type (i.e. not Float32, Float64,
ComplexF32, or ComplexF64) detected as the element type
for the initial condition or time span. These generic
number types are only compatible with the pure Julia
solvers which support generic programming, such as
OrdinaryDiffEq.jl. The chosen solver does not support
this functionality. Please double check that the initial
condition and time span types are correct, and check that
the chosen solver was correct.

Solver: PIRect()
u0 type: Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqParamEstim.var"#29#30"{Nothing, typeof(DiffEqParamEstim.STANDARD_PROB_GENERATOR), Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, @NamedTuple{dt::Float64, maxiters::Int64, verbose::Bool}}, FODEProblem{Vector{Float64}, Tuple{Int64, Int64}, Vector{Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(SIQR), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Float64}, FractionalDiffEq.StandardFODEProblem, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, PIRect, L2Loss{Vector{Float64}, Matrix{Float64}, Nothing, Nothing, Nothing}, Nothing, Tuple{}}, Float64}, Float64, 9}}
Timespan type: Tuple{Float64, Float64}

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Stacktrace:
  [1] check_prob_alg_pairing(prob::FODEProblem{…}, alg::PIRect)
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1598
  [2] solve_up(prob::FODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::Vector{…}, args::PIRect; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1195
  [3] solve(prob::FODEProblem{…}, args::PIRect; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1089
  [4] (::DiffEqParamEstim.var"#29#30"{…})(p::Vector{…}, ::SciMLBase.NullParameters)
    @ DiffEqParamEstim C:\Users\sanyam\.julia\packages\DiffEqParamEstim\6z5Ou\src\build_loss_objective.jl:12
  [5] (::DifferentiationInterface.FixTail{DiffEqParamEstim.var"#29#30"{…}, Tuple{…}})(args::Vector{ForwardDiff.Dual{…}})
    @ DifferentiationInterface C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\utils\context.jl:169
  [6] vector_mode_dual_eval!(f::DifferentiationInterface.FixTail{…}, cfg::ForwardDiff.GradientConfig{…}, x::Vector{…})
    @ ForwardDiff C:\Users\sanyam\.julia\packages\ForwardDiff\UBbGT\src\apiutils.jl:24
  [7] vector_mode_gradient!(result::Vector{…}, f::DifferentiationInterface.FixTail{…}, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff C:\Users\sanyam\.julia\packages\ForwardDiff\UBbGT\src\gradient.jl:98
  [8] gradient!(result::Vector{…}, f::DifferentiationInterface.FixTail{…}, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…}, ::Val{…})
    @ ForwardDiff C:\Users\sanyam\.julia\packages\ForwardDiff\UBbGT\src\gradient.jl:39
  [9] gradient!(f::DiffEqParamEstim.var"#29#30"{…}, grad::Vector{…}, prep::DifferentiationInterfaceForwardDiffExt.ForwardDiffGradientPrep{…}, backend::AutoForwardDiff{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceForwardDiffExt C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\ext\DifferentiationInterfaceForwardDiffExt\onearg.jl:436
 [10] (::OptimizationBase.var"#grad#16"{…})(res::Vector{…}, θ::Vector{…})
    @ OptimizationBase C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:28
 [11] (::OptimizationOptimJL.var"#19#23"{…})(G::Vector{…}, θ::Vector{…})
    @ OptimizationOptimJL C:\Users\sanyam\.julia\packages\OptimizationOptimJL\VaURt\src\OptimizationOptimJL.jl:285
 [12] value_gradient!!(obj::OnceDifferentiable{Float64, Vector{Float64}, Vector{Float64}}, x::Vector{Float64})
    @ NLSolversBase C:\Users\sanyam\.julia\packages\NLSolversBase\totsP\src\interface.jl:82
 [13] value_gradient!!(bw::Optim.BarrierWrapper{…}, x::Vector{…})
    @ Optim C:\Users\sanyam\.julia\packages\Optim\8dE7C\src\multivariate\solvers\constrained\fminbox.jl:90
 [14] initial_state(method::LBFGS{…}, options::Optim.Options{…}, d::Optim.BarrierWrapper{…}, initial_x::Vector{…})
    @ Optim C:\Users\sanyam\.julia\packages\Optim\8dE7C\src\multivariate\solvers\first_order\l_bfgs.jl:168
 [15] optimize(df::OnceDifferentiable{…}, l::Vector{…}, u::Vector{…}, initial_x::Vector{…}, F::Fminbox{…}, options::Optim.Options{…})
    @ Optim C:\Users\sanyam\.julia\packages\Optim\8dE7C\src\multivariate\solvers\constrained\fminbox.jl:451
 [16] __solve(cache::OptimizationCache{…})
    @ OptimizationOptimJL C:\Users\sanyam\.julia\packages\OptimizationOptimJL\VaURt\src\OptimizationOptimJL.jl:312
 [17] solve!(cache::OptimizationCache{…})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:227
 [18] solve(::OptimizationProblem{…}, ::LBFGS{…}; kwargs::@Kwargs{})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:129
 [19] solve(::OptimizationProblem{…}, ::LBFGS{…})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:126
 [20] top-level scope
    @ c:\Users\sanyam\Downloads\JuliaOnVSCode101\FractionPaperStuff\SIQRFracFittingTry1.jl:61
Some type information was truncated. Use `show(err)` to see complete types.

Hello and welcome to the community :waving_hand:

The error message seems pretty clear to me, the non-standard number type ForwardDiff.Dual appears because you are using AutoForwardDiff. You could try to make use of finite differences for computing the gradient instead of ForwardDiff, i.e., change

into Optimization.AutoFiniteDiff().

1 Like

Thank you, this has fixed the issue of number type.
I am now getting a different error. I have tried to fix it/ find any relevant information about it. I dont know if it is at all related to the original issue or is it a problem with something else?
Here is my error:

ERROR: MethodError: anyeltypedual(::ODESolution{…}, ::Type{…}) is ambiguous.

Candidates:
  anyeltypedual(sol::AbstractDiffEqArray, counter)
    @ DiffEqBaseForwardDiffExt C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\ext\DiffEqBaseForwardDiffExt.jl:360
  anyeltypedual(x, ::Type{Val{counter}}) where counter
    @ DiffEqBaseForwardDiffExt C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\ext\DiffEqBaseForwardDiffExt.jl:204

Possible fix, define
  anyeltypedual(::AbstractDiffEqArray, ::Type{Val{counter}}) where counter

Stacktrace:
  [1] promote_u0(u0::Vector{…}, p::ODESolution{…}, t0::Float64)
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\utils.jl:151
  [2] get_concrete_problem(prob::FODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1310
  [3] get_concrete_problem
    @ C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1301 [inlined]
  [4] solve_up(prob::FODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ODESolution{…}, args::PIRect; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1193
  [5] solve_up
    @ C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1177 [inlined]
  [6] solve(prob::FODEProblem{…}, args::PIRect; sensealg::Nothing, u0::Nothing, p::ODESolution{…}, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1089
  [7] loss(p::ODESolution{…})
    @ Main c:\Users\sanyam\Downloads\JuliaOnVSCode101\FractionPaperStuff\SIQRFracFittingTry1.jl:23
  [8] (::DiffEqParamEstim.var"#29#30"{…})(p::Vector{…}, ::SciMLBase.NullParameters)
    @ DiffEqParamEstim C:\Users\sanyam\.julia\packages\DiffEqParamEstim\6z5Ou\src\build_loss_objective.jl:18
  [9] prepare_gradient_nokwarg(strict::Val{…}, f::DiffEqParamEstim.var"#29#30"{…}, backend::AutoEnzyme{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterface C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\first_order\gradient.jl:92
 [10] prepare_gradient(f::DiffEqParamEstim.var"#29#30"{…}, backend::AutoEnzyme{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…}; strict::Val{…})
    @ DifferentiationInterface C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\first_order\gradient.jl:11
 [11] prepare_gradient(f::DiffEqParamEstim.var"#29#30"{…}, backend::AutoEnzyme{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterface C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\first_order\gradient.jl:8
 [12] instantiate_function(f::OptimizationFunction{…}, x::Vector{…}, adtype::AutoEnzyme{…}, p::SciMLBase.NullParameters, num_cons::Int64; g::Bool, h::Bool, hv::Bool, fg::Bool, fgh::Bool, cons_j::Bool, cons_vjp::Bool, cons_jvp::Bool, cons_h::Bool, lag_h::Bool)
    @ OptimizationBase C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:26
 [13] instantiate_function
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:17 [inlined]
 [14] #instantiate_function#40
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:274 [inlined]
 [15] instantiate_function
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:267 [inlined]
 [16] OptimizationCache(prob::OptimizationProblem{…}, opt::Fminbox{…}; callback::Function, maxiters::Nothing, maxtime::Nothing, abstol::Nothing, reltol::Nothing, progress::Bool, structural_analysis::Bool, manifold::Nothing, kwargs::@Kwargs{})
    @ OptimizationBase C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\cache.jl:60
 [17] OptimizationCache
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\cache.jl:25 [inlined]
 [18] #__init#2
    @ C:\Users\sanyam\.julia\packages\OptimizationOptimJL\VaURt\src\OptimizationOptimJL.jl:110 [inlined]
 [19] __init
    @ C:\Users\sanyam\.julia\packages\OptimizationOptimJL\VaURt\src\OptimizationOptimJL.jl:81 [inlined]
 [20] #init#744
    @ C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:213 [inlined]
 [21] init
    @ C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:208 [inlined]
 [22] solve(::OptimizationProblem{…}, ::BFGS{…}; kwargs::@Kwargs{})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:129
 [23] solve(::OptimizationProblem{…}, ::BFGS{…})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:126
 [24] top-level scope
    @ c:\Users\sanyam\Downloads\JuliaOnVSCode101\FractionPaperStuff\SIQRFracFittingTry1.jl:35
Some type information was truncated. Use `show(err)` to see complete types.

ERROR: MethodError: anyeltypedual(::ODESolution{…}, ::Type{…}) is ambiguous.

Candidates:
  anyeltypedual(sol::AbstractDiffEqArray, counter)
    @ DiffEqBaseForwardDiffExt C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\ext\DiffEqBaseForwardDiffExt.jl:360
  anyeltypedual(x, ::Type{Val{counter}}) where counter
    @ DiffEqBaseForwardDiffExt C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\ext\DiffEqBaseForwardDiffExt.jl:204

Possible fix, define
  anyeltypedual(::AbstractDiffEqArray, ::Type{Val{counter}}) where counter

Stacktrace:
  [1] promote_u0(u0::Vector{…}, p::ODESolution{…}, t0::Float64)
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\utils.jl:151
  [2] get_concrete_problem(prob::FODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1310
  [3] get_concrete_problem
    @ C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1301 [inlined]
  [4] solve_up(prob::FODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ODESolution{…}, args::PIRect; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1193
  [5] solve_up
    @ C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1177 [inlined]
  [6] solve(prob::FODEProblem{…}, args::PIRect; sensealg::Nothing, u0::Nothing, p::ODESolution{…}, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\sanyam\.julia\packages\DiffEqBase\5hvMq\src\solve.jl:1089
  [7] loss(p::ODESolution{…})
    @ Main c:\Users\sanyam\Downloads\JuliaOnVSCode101\FractionPaperStuff\SIQRFracFittingTry1.jl:23
  [8] (::DiffEqParamEstim.var"#29#30"{…})(p::Vector{…}, ::SciMLBase.NullParameters)
    @ DiffEqParamEstim C:\Users\sanyam\.julia\packages\DiffEqParamEstim\6z5Ou\src\build_loss_objective.jl:18
  [9] (::DifferentiationInterface.FixTail{DiffEqParamEstim.var"#29#30"{…}, Tuple{…}})(args::Vector{Float64})
    @ DifferentiationInterface C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\utils\context.jl:169
 [10] prepare_gradient_nokwarg(strict::Val{…}, f::Function, backend::AutoFiniteDiff{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceFiniteDiffExt C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\ext\DifferentiationInterfaceFiniteDiffExt\onearg.jl:261
 [11] #prepare_gradient#47
    @ C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\first_order\gradient.jl:11 [inlined]
 [12] prepare_gradient
    @ C:\Users\sanyam\.julia\packages\DifferentiationInterface\zJHX8\src\first_order\gradient.jl:8 [inlined]
 [13] instantiate_function(f::OptimizationFunction{…}, x::Vector{…}, adtype::AutoFiniteDiff{…}, p::SciMLBase.NullParameters, num_cons::Int64; g::Bool, h::Bool, hv::Bool, fg::Bool, fgh::Bool, cons_j::Bool, cons_vjp::Bool, cons_jvp::Bool, cons_h::Bool, lag_h::Bool)
    @ OptimizationBase C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:26
 [14] instantiate_function
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:17 [inlined]
 [15] #instantiate_function#40
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:274 [inlined]
 [16] instantiate_function
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\OptimizationDIExt.jl:267 [inlined]
 [17] OptimizationCache(prob::OptimizationProblem{…}, opt::Fminbox{…}; callback::Function, maxiters::Nothing, maxtime::Nothing, abstol::Nothing, reltol::Nothing, progress::Bool, structural_analysis::Bool, manifold::Nothing, kwargs::@Kwargs{})
    @ OptimizationBase C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\cache.jl:60
 [18] OptimizationCache
    @ C:\Users\sanyam\.julia\packages\OptimizationBase\yD5XZ\src\cache.jl:25 [inlined]
 [19] #__init#2
    @ C:\Users\sanyam\.julia\packages\OptimizationOptimJL\VaURt\src\OptimizationOptimJL.jl:110 [inlined]
 [20] __init
    @ C:\Users\sanyam\.julia\packages\OptimizationOptimJL\VaURt\src\OptimizationOptimJL.jl:81 [inlined]
 [21] #init#744
    @ C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:213 [inlined]
 [22] init
    @ C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:208 [inlined]
 [23] solve(::OptimizationProblem{…}, ::BFGS{…}; kwargs::@Kwargs{})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:129
 [24] solve(::OptimizationProblem{…}, ::BFGS{…})
    @ SciMLBase C:\Users\sanyam\.julia\packages\SciMLBase\Ha7rZ\src\solve.jl:126
 [25] top-level scope
    @ c:\Users\sanyam\Downloads\JuliaOnVSCode101\FractionPaperStuff\SIQRFracFittingTry1.jl:35
Some type information was truncated. Use `show(err)` to see complete types.