PDE Parameter Estimation

Specific Domains Optimization (Mathematical)
, , ,
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