PDE Parameter Estimation

Specific Domains Optimization (Mathematical)
, , ,
1

Hello, I’m new a DifferentialEquations/ModelingToolkit user that would like to perform parameter estimation of source terms in advection-dispersion equation.

I’ve defined the PDE using ModelingToolkit, based on this post. I turned the PDE system into an ODE problem using MethodOfLines, and created synthetic data solving the problem for a given source amplitude. The code is below.

##### Define PDE problem #####
using ModelingToolkit 
using DifferentialEquations
using MethodOfLines
# Variables
@variables x 
@independent_variables t 
@variables c(..) 
∂t = Differential(t)
∂x = Differential(x)
∂²x = Differential(x)^2
# Parameters
@parameters λ₁ 
# Constants
u = 0.5
K = 0.01
D = 100
# Sources terms
function δ(x,width) #gaussian
    return 1/(width*sqrt(π/18))*exp(-18*x^2/(width^2))
end
source=λ₁*δ(t-1,1)*δ(x-1,1)
# PDE
eqn = ∂t(c(x, t)) ~ D * ∂²x(c(x, t)) - u * ∂x(c(x, t)) - K * c(x, t) + source
# Space and time domains
L=3.0 
T=3.0 
dom = [x ∈ (0, L), t ∈ (0, T)]
# Initial and boundary conditions
bcs = [
       c(0,t) ~ 0., 
       c(x,0) ~ 0., 
       ∂t(c(x,0)) ~ 0., 
       ∂x(c(0,t)) ~ 0.
       ] 
# Define PDE system
@named sys = PDESystem(eqn,
                       bcs,
                       dom,
                       [x, t],
                       [c(x, t)],
                       [λ₁=>1]
                       )
# Define ODE problem
prob = discretize(sys, 
                  MOLFiniteDifference([x => 0.05],
                                      t,
                                      advection_scheme=scheme=UpwindScheme(2)
                                      )             
                  );
##### Generate synthetic data #####
function solve_prob(prob, p::Vector{Pair{Num, Float64}}; kwargs...)
    newprob=remake(prob, p=p)
    sol = solve(newprob,
                   Tsit5(),
                  ;kwargs...
                   )
    return(sol)
end
@time sol=solve_prob(prob,
                [λ₁=>1.]
                ;
                saveat=0.1, 
                )
# Collect solution
key=collect(keys(sol.u))[1]
data=sol.u[key];

From this point, I would like to estimate λ₁ from synthetic data. I first tried using DiffEqParamEstim, adapting ODE parameter example:

##### Parameter estimation using DiffEqParamEstim ######
using Optimization
using DiffEqParamEstim
# Create cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(sol.t, data),
                                     Optimization.AutoForwardDiff(),
                                     maxiters = 1000000, verbose = true)
# Test cost function
cost_function(1.)

but I get the following error:

MethodError: getindex(::SciMLBase.PDETimeSeriesSolution{Float64, 1, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Float64, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing}, Nothing, Vector{Float64}, Tuple{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Vector{Num}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Float64, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Dict{Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}, SciMLBase.DEStats}, ::Int64) is ambiguous.

Candidates:
  getindex(A::AbstractDiffEqArray, i::Int64)
    @ RecursiveArrayTools deprecated.jl:103
  getindex(A::Union{SciMLBase.PDENoTimeSolution{T, N, S, D}, SciMLBase.PDETimeSeriesSolution{T, N, S, D}}, sym) where {T, N, S, D<:MethodOfLines.MOLMetadata}
    @ MethodOfLines C:\Users\yoann.cartier\.julia\packages\MethodOfLines\iIYN1\src\interface\solution\common.jl:36
  getindex(A::AbstractVectorOfArray, I::Int64)
    @ RecursiveArrayTools deprecated.jl:103
  getindex(A::Union{SciMLBase.PDENoTimeSolution{T, N, S, D}, SciMLBase.PDETimeSeriesSolution{T, N, S, D}}, sym, args...) where {T, N, S, D<:MethodOfLines.MOLMetadata}
    @ MethodOfLines C:\Users\yoann.cartier\.julia\packages\MethodOfLines\iIYN1\src\interface\solution\common.jl:64

Possible fix, define
  getindex(::SciMLBase.PDETimeSeriesSolution{T, N, S, D}, ::Int64) where {T, N, S, D<:MethodOfLines.MOLMetadata}

So I tried using SciMLSensitivity instead, following a SciML example for PDE. The code:

##### Parameter estimation using SciMLSensitivity ######
# Define optimization problem
using SciMLSensitivity
# Optimization methods
using Zygote
using OptimizationPolyalgorithms
using OptimizationOptimJL

function solve_prob(prob, p::Vector{<:Real}; kwargs...)
    newprob=remake(prob, p=p)
    sol = solve(newprob,
                   Tsit5(),
                  ;kwargs...
                   )
    return(sol)
end
# Define loss function, here based on RMSE
# This should return a scalar, the loss value, as the first return output and if any additional outputs are returned, they will be passed to the callback function
function loss(θ)
    pred_sol=solve_prob(prob,θ,saveat=0.1)
    return sqrt(sum(abs2.(pred_sol.u[key].-data))),pred_sol.u[key]
end

# Callback function to observe training
function CB(θ, l, pred) 
    display(l)
    append!(PRED, [pred])
    append!(LOSS, l)
    append!(PARS, [θ])
    false
end

# Define the optimization problem: 
adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x,_) -> loss(x),                                        adtype)
#First argument: f(u,p,args...): the function to optimize. u are the optimization variables and p are parameters used in definition of the objective, even if no such parameters are used in the objective it should be an argument in the function. This can also take any additional arguments that are relevant to the objective function                                  
optprob = Optimization.OptimizationProblem(optf,[1.1]);
# Perform optimization
res = Optimization.solve(optprob, PolyOpt(), callback = CB)
@show res.u

But I get the error:

MethodError: no method matching SciMLBase.PDETimeSeriesSolution(::SciMLBase.PDETimeSeriesSolution{Float64, 1, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, @NamedTuple{callback::Nothing}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing}, Nothing, Vector{Float64}, Tuple{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Vector{Num}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, @NamedTuple{callback::Nothing}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Dict{Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}, SciMLBase.DEStats}, ::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization})

Closest candidates are:
  SciMLBase.PDETimeSeriesSolution(!Matched::SciMLBase.AbstractODESolution{T}, ::MethodOfLines.MOLMetadata) where T
   @ MethodOfLines C:\Users\yoann.cartier\.julia\packages\MethodOfLines\iIYN1\src\interface\solution\timedep.jl:11

So I’m stuck for now. Does anyone know how to make parameter estimation work for PDEs using ModelingToolkit and Optimization ?

Many thanks in advance for your ideas!

2

Hello again, this self-answer is to provide a shorter MWE based on heat equation, trying to estimate heat transfert coefficient.

My main questions are: how to use DiffEqParamEstim and a PDE created with ModelingToolKit ? Is this optimal in terms of performance?

I first define PDE problem using ModelingToolKit:

##### Define PDE problem #####
using ModelingToolkit: Interval
using ModelingToolkit 
using DifferentialEquations
using MethodOfLines
# Variables
@variables x 
@independent_variables t 
@variables c(..) 
∂t = Differential(t)
∂²x = Differential(x)^2
# Parameters
@parameters D
# PDE
eqn = ∂t(c(x, t)) ~ D * ∂²x(c(x, t))
# Space and time domains
L=3.0 
T=3.0 
dom = [x ∈ Interval(0, L), t ∈ Interval(0, T)]
# Initial and boundary conditions
bcs = [
       c(0,t) ~ 1.0, 
       c(x,0) ~ 0.0,
       ] 
# Define PDE system
@named sys = PDESystem(eqn,
                       bcs,
                       dom,
                       [x, t],
                       [c(x, t)],
                       [D=>0.1]
                       )
# Define ODE problem
prob = discretize(sys, 
                  MOLFiniteDifference([x => 0.05],
                                      t,
                                      advection_scheme=scheme=UpwindScheme(2)
                                      )             
                  );
##### Generate synthetic data #####
function solve_prob(prob, p::Vector{Pair{Num, Float64}}; kwargs...)
    newprob=remake(prob, p=p)
    sol = solve(newprob,
                   Tsit5(),
                  ;kwargs...
                   )
    return(sol)
end
D₀=0.1
@time sol=solve_prob(prob,
                [D=>D₀]
                ;
                saveat=0.1, 
                )
# Collect solution
key=collect(keys(sol.u))[1]
data=sol.u[key];

I would like to estimate D from synthetic data. I first tried using DiffEqParamEstim:

##### Parameter estimation using DiffEqParamEstim ######
using Optimization
using DiffEqParamEstim
# Create cost function
cost_function = build_loss_objective(prob, Tsit5(), L2Loss(sol.t, data),
                                     Optimization.AutoForwardDiff(),
                                     maxiters = 1000000, verbose = true)
# Test cost function
cost_function(0.2)

But I get the error:

MethodError: getindex(::SciMLBase.PDETimeSeriesSolution{Float64, 1, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Float64, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing}, Nothing, Vector{Float64}, Tuple{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Vector{Num}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Float64, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Float64, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Dict{Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}, SciMLBase.DEStats}, ::Int64) is ambiguous.

Candidates:
  getindex(A::AbstractDiffEqArray, i::Int64)
    @ RecursiveArrayTools deprecated.jl:103
  getindex(A::Union{SciMLBase.PDENoTimeSolution{T, N, S, D}, SciMLBase.PDETimeSeriesSolution{T, N, S, D}}, sym) where {T, N, S, D<:MethodOfLines.MOLMetadata}
    @ MethodOfLines C:\Users\yoann.cartier\.julia\packages\MethodOfLines\iIYN1\src\interface\solution\common.jl:36
  getindex(A::AbstractVectorOfArray, I::Int64)
    @ RecursiveArrayTools deprecated.jl:103
  getindex(A::Union{SciMLBase.PDENoTimeSolution{T, N, S, D}, SciMLBase.PDETimeSeriesSolution{T, N, S, D}}, sym, args...) where {T, N, S, D<:MethodOfLines.MOLMetadata}
    @ MethodOfLines C:\Users\yoann.cartier\.julia\packages\MethodOfLines\iIYN1\src\interface\solution\common.jl:64

Possible fix, define
  getindex(::SciMLBase.PDETimeSeriesSolution{T, N, S, D}, ::Int64) where {T, N, S, D<:MethodOfLines.MOLMetadata}

I also tried using SciMLSensitivity instead:

##### Parameter estimation using SciMLSensitivity ######
# Define optimization problem
using SciMLSensitivity
# Optimization methods
using Zygote
using OptimizationPolyalgorithms
using OptimizationOptimJL

function solve_prob(prob, p::Vector{<:Real}; kwargs...)
    newprob=remake(prob, p=p)
    sol = solve(newprob,
                   Tsit5(),
                  ;kwargs...
                   )
    return(sol)
end
# Define loss function, here based on RMSE
# This should return a scalar, the loss value, as the first return output and if any additional outputs are returned, they will be passed to the callback function
function loss(θ)
    pred_sol=solve_prob(prob,θ,saveat=0.1)
    return sqrt(sum(abs2.(pred_sol.u[key].-data))),pred_sol.u[key]
end

# Define the optimization problem: 
adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x,_) -> loss(x), adtype)
#First argument: f(u,p,args...): the function to optimize. u are the optimization variables and p are parameters used in definition of the objective, even if no such parameters are used in the objective it should be an argument in the function. This can also take any additional arguments that are relevant to the objective function                                  
optprob = Optimization.OptimizationProblem(optf,[0.2]);
# Perform optimization
res = Optimization.solve(optprob, PolyOpt())
@show res.u

But I get the error:

MethodError: no method matching SciMLBase.PDETimeSeriesSolution(::SciMLBase.PDETimeSeriesSolution{Float64, 1, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, @NamedTuple{callback::Nothing}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing}, Nothing, Vector{Float64}, Tuple{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Vector{Num}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, @NamedTuple{callback::Nothing}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Dict{Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}, SciMLBase.DEStats}, ::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization})

Closest candidates are:
  SciMLBase.PDETimeSeriesSolution(!Matched::SciMLBase.AbstractODESolution{T}, ::MethodOfLines.MOLMetadata) where T
   @ MethodOfLines C:\Users\yoann.cartier\.julia\packages\MethodOfLines\iIYN1\src\interface\solution\timedep.jl:11

Many thanks in advancefor your help!

3

DiffEqParamEstim isn’t compatible with MethodOfLines at this moment. You’d need to more directly write a loss function for now.

1 Like
4

Ok, and this is also the case for SciMLSensitivity from what I understand?
I’ll write the loss function for now then, thanks!

5

SciMLSensitivity.jl is just the system for getting gradients. It would be used with any gradient-based local optimization.

6

To get an interface as close as possible to DiffEqParamEstim.jl I overloaded L2Loss for PDETimeSeriesSolution

import ModelingToolkit.SciMLBase: PDETimeSeriesSolution
function (f::L2Loss)(sol::PDETimeSeriesSolution)
    xidx, data = f.data
    weight = f.data_weight
    diff_weight = f.differ_weight
    colloc_grad = f.colloc_grad
    dudt = f.dudt

    if sol isa DiffEqBase.AbstractEnsembleSolution
        failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol)
    else
        failure = !SciMLBase.successful_retcode(sol.retcode)
    end
    failure && return Inf

    sum(abs2, data - sol[c(x, t)][x_idx, :])
end

It looks a bit hacky (maybe suboptimal somewhere?) as I have to provided observed points indexes. However, it works here

x_obs = 0:0.05:L
x_idx = findall(x -> x ∈ x_obs, sol[x])

cost_function = build_loss_objective(prob, Tsit5(), L2Loss(sol.t, (x_idx, data)),
                                     Optimization.AutoForwardDiff(),
                                     maxiters = 1000000, verbose = true)

cost_function(0.2)
cost_function(0.1) # = 0 as expected
optprob = Optimization.OptimizationProblem(cost_function, [0.3])
using OptimizationOptimJL
@time "optimization" optsol = solve(optprob, BFGS())

Here is complete code

##### Define PDE problem #####
using ModelingToolkit: Interval
using ModelingToolkit 
using DifferentialEquations
using MethodOfLines
# Variables
@variables x 
@independent_variables t 
@variables c(..) 
∂t = Differential(t)
∂²x = Differential(x)^2
# Parameters
@parameters D
# PDE
eqn = ∂t(c(x, t)) ~ D * ∂²x(c(x, t))
# Space and time domains
L=3.0 
T=3.0 
dom = [x ∈ Interval(0, L), t ∈ Interval(0, T)]
# Initial and boundary conditions
bcs = [
       c(0,t) ~ 1.0, 
       c(x,0) ~ 0.0,
       ] 
# Define PDE system
@named sys = PDESystem(eqn,
                       bcs,
                       dom,
                       [x, t],
                       [c(x, t)],
                       [D];
                       defaults=Dict(D => 0.1)
                       )

# Define ODE problem
prob = discretize(sys, 
                  MOLFiniteDifference([x => 0.05],
                                      t,
                                      advection_scheme=scheme=UpwindScheme(2)
                                      )             
                  );
##### Generate synthetic data #####
@time "PDE Solve" sol = solve(prob, Tsit5(), saveat=0.1)

key=collect(keys(sol.u))[1]
data=sol.u[key]

##### Parameter estimation using DiffEqParamEstim ######
using Optimization
using DiffEqParamEstim
# # Create cost function
# cost_function = build_loss_objective(prob, Tsit5(), L2Loss(sol.t, data),
#                                      Optimization.AutoForwardDiff(),
#                                      maxiters = 1000000, verbose = true)
# # Test cost function
# cost_function(0.2)

import ModelingToolkit.SciMLBase: PDETimeSeriesSolution
function (f::L2Loss)(sol::PDETimeSeriesSolution)
    xidx, data = f.data
    weight = f.data_weight
    diff_weight = f.differ_weight
    colloc_grad = f.colloc_grad
    dudt = f.dudt

    if sol isa DiffEqBase.AbstractEnsembleSolution
        failure = any(!SciMLBase.successful_retcode(s.retcode) for s in sol)
    else
        failure = !SciMLBase.successful_retcode(sol.retcode)
    end
    failure && return Inf

    sum(abs2, data - sol[c(x, t)][x_idx, :])
end

x_obs = 0:0.05:L
x_idx = findall(x -> x ∈ x_obs, sol[x])

cost_function = build_loss_objective(prob, Tsit5(), L2Loss(sol.t, (x_idx, data)),
                                     Optimization.AutoForwardDiff(),
                                     maxiters = 1000000, verbose = true)

cost_function(0.2)
optprob = Optimization.OptimizationProblem(cost_function, [0.3])
using OptimizationOptimJL
@time "optimization" optsol = solve(optprob, BFGS())

Notice I get a few warning (that I don’t understand all) during last solve.
Plus I changed compare to original @cartiery code the PDESystem definition according to latest notation, see Discretization failed (no method matching getmetadata(...)) · Issue #427 · SciML/MethodOfLines.jl · GitHub