ModelingToolkit, SciMLSensitivity: Computed parameter?

Specific Domains Modelling & Simulations
, ,
1

I am working on a model of a leaky bucket using ModelingToolkit.jl. I was experimenting with fitting parameters of the model based on randomized data using SciMLSensitivity.jl. In the attached example, I define my constants and parameters to be optimized as well as the diff equation of my simple model.

I derived the model using bond graphs. The generalized compliance C of a regularly shaped bucket is C = (A_b/(ρ*g). My question is, how can I correctly use this computation inside of the @mtkmodel macro. Before I started to experiment with fitting parameters, it caused no problems in @parameters, but in the example below, I am not sure where to calculate it from constants and optimized parameters.

Here is the example:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

using RecursiveArrayTools, Plots
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO

using OrdinaryDiffEq, SciMLSensitivity, Zygote

@mtkmodel FOL begin
    @constants begin
        ρ = 1000  # kg/m³
        g = 9.81 # m/s²
    end
    @parameters begin
        A_b = 0.2 # m
        R = 10^5 
    end
    @variables begin
        V(t) = 1.0 # m³
    end
    @equations begin
        D(V) ~ -V/( (A_b/(ρ*g)) *R) # C = (A_b/(ρ*g)
    end
end

@mtkbuild fol = FOL()

# Simulation interval and intermediary points
tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

prob = ODEProblem(fol, [], (0.0, 10.0), [])
sol = solve(prob, Tsit5())

# generate synthetic data
tt = collect(range(0, stop = 10, length = 101))
randomized = VectorOfArray([vec(sol(tt[i])) + 0.1randn(length(sol(tt[i]))) for i in 1:length(tt)])
data = convert(Array, randomized)

# define loss function
function loss(p)
    sol = solve(prob, Tsit5(), p = p, saveat = tsteps)
    loss = sum(abs2, sol .- data)
    return loss
end

# fit parameters
adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)

p = [0.25, 10^5 + 5*10^4] # initial conditions
optprob = Optimization.OptimizationProblem(optf, p, lb = [0.1, 10^4], ub = [0.3, 10^6])
result_ode = Optimization.solve(optprob, BFGS(), maxiters = 100)

remade_solution = solve(remake(prob, p = result_ode.u), Tsit5(),
saveat = tsteps)

# plot result
plot(remade_solution, label="Remade Solution")
plot!(sol, label="Original solution")
scatter!(tt, data', label="Randomized Data", marker=:dot)
2

It looks like you’re fitting all parameters. You probably want to follow the docs on just fitting the tunables.

1 Like