Hi All,
I want to optimize the Model created based on ModelingToolkitStandardLibrary.
So I have a small RC-Model copied from the official documentation to ModelingToolkitStandardLibrary.
module Grid
export get_model
using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
function get_model(; resistance::Float64=1.0, capacitance::Float64=1.0, voltage::Float64=1.0)
@variables t
@named R = Resistor(R=resistance)
@named C = Capacitor(C=capacitance)
@named V = Voltage()
@named ground = Ground()
@named source = Constant(k=voltage)
rc_eqs = [
connect(V.p, R.p)
connect(R.n, C.p)
connect(C.n, V.n, ground.g)
connect(V.V, source.output)
]
@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model, [R, C, V, source, ground])
sys = structural_simplify(rc_model)
return sys, states(rc_model)
end
end
The full code to run optimization is:
using Revise
includet("./grid.jl")
using .Grid
using ModelingToolkit, DifferentialEquations
using Optimization, OptimizationOptimisers, DiffEqSensitivity, Zygote, OptimizationOptimJL
using CSV, Plots, DataFrames, Statistics
function get_obs(fname::String, obs_vars::Vector{String})
df = DataFrame(CSV.File(fname))
return df[:, "time"], df[:, obs_vars]
end
function get_param_vars(model::ModelingToolkit.AbstractSystem)
params = parameters(model)
return params
end
function loss(p::Vector{Float64})
pred_sol = solve(prob, Tsit5(), p=p, saveat=t_save)
pred = Array(pred_sol')
mse = mean((pred - obs_vals) .^ 2)
return mse
end
# observation variable names
obs_vars = ["C₊v(t)"]
obs_data_file = "data/res.csv"
t_save, obs = get_obs(obs_data_file, obs_vars)
model, model_state = get_model()
param_vars = get_param_vars(model)
u0 = [s => 0.0 for s in states(model)]
tspan = (t_save[1], last(t_save))
lb = [0.0, 0.0, 0.99]
ub = [3.0, 3.0, 1.01]
prob = ODAEProblem(model, u0, tspan)
obs_vals = Array(obs)
p_init = [1.0, 1.0, 1.0]
adtype = Optimization.AutoZygote()
optf = Optimization.OptimizationFunction((x, p) -> loss(x), adtype)
optprob = Optimization.OptimizationProblem(optf, p_init, lb=lb, ub=ub)
result = Optimization.solve(optprob, LBFGS(), maxiters=300)
That code works well, and I can get the correct solution.
My problem is in the ordering of the parameters. The official doc says “symbol ordering is not guaranteed” but if I do something like this:
function loss(p::Vector{Float64})
p1 = [param_vars[1] => p[1], param_vars[2] => p[2], param_vars[3] => p[3]]
p2 = ModelingToolkit.varmap_to_vars(p1, param_vars)
pred_sol = solve(prob, Tsit5(), p=p2, saveat=t_save)
pred = Array(pred_sol')
mse = mean((pred - obs_vals) .^ 2)
return mse
end
I get the next error:
ERROR: Compiling Tuple{Type{Dict}, Vector{Pair{Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}}: try/catch is not supported.
Refer to the Zygote documentation for fixes.
So how can I do optimization and guarantee the correct order of the parameters?