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.
