# PDE Parameter Estimation

**URL:** <https://discourse.julialang.org/t/pde-parameter-estimation/120724>\
**Category:** Optimization (Mathematical)\
**Tags:** question, diffeq, optimization, modelingtoolkit\
**Created:** [September 30, 2024, 2:00pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724 "2024-09-30T14:00:33Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![cartiery](https://avatars.discourse-cdn.com/v4/letter/c/e47c2d/32.png) [@cartiery](https://discourse.julialang.org/u/cartiery)\
**Post date:** [September 30, 2024, 2:00pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724/1 "2024-09-30T14:00:33Z")

</div>

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](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/16). 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.

```julia
##### 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](https://docs.sciml.ai/DiffEqParamEstim/stable/getting_started/):

```julia
##### 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:

```julia
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](https://docs.sciml.ai/SciMLSensitivity/stable/examples/pde/pde_constrained/) for PDE. The code:

```julia
##### 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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![cartiery](https://avatars.discourse-cdn.com/v4/letter/c/e47c2d/32.png) [@cartiery](https://discourse.julialang.org/u/cartiery)\
**Post date:** [November 8, 2024, 2:59pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724/2 "2024-11-08T14:59:19Z")

</div>

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:

```julia
##### 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** :

```julia
##### 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:

```julia
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:

```julia
##### 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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 8, 2024, 4:32pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724/3 "2024-11-08T16:32:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![cartiery](https://avatars.discourse-cdn.com/v4/letter/c/e47c2d/32.png) [@cartiery](https://discourse.julialang.org/u/cartiery)\
**Post date:** [November 8, 2024, 5:20pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724/4 "2024-11-08T17:20:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 8, 2024, 5:33pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724/5 "2024-11-08T17:33:39Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dmetivie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmetivie/32/6926_2.png) [@dmetivie](https://discourse.julialang.org/u/dmetivie)\
**Post date:** [November 12, 2024, 2:37pm UTC](https://discourse.julialang.org/t/pde-parameter-estimation/120724/6 "2024-11-12T14:37:15Z")

</div>

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

```julia
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

```julia
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

```julia
##### 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](https://github.com/SciML/MethodOfLines.jl/issues/427)