I have the simple system from the documentation.
module Grid1
export get_model
using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
function get_model(; resistance::Float64=0.4, capacitance::Float64=0.6, 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, parameters(sys), states(sys), C
end
end
I want to run optimization of that system wrt parameters [R and C]. The following code works well.
using Revise
includet("./grid1.jl")
using .Grid1
using ModelingToolkit, DifferentialEquations
using Optimization, OptimizationOptimisers, DiffEqSensitivity, ForwardDiff, Zygote, OptimizationOptimJL
using CSV, DataFrames, Statistics
function get_obs(fname::String, obs_vars::Vector{String})
df = DataFrame(CSV.File(fname))
return df[:, "time"], df[:, obs_vars]
end
function loss(p::Vector{Float64})
pred_sol = solve(prob, Tsit5(), p=p, saveat=t_save)
pred = Array(pred_sol[model_state[1]])
# pred = Array(pred_sol')
mse = mean((pred - obs_vals) .^ 2)
return mse
end
function callback(p, l)
println(l, " ", p)
return false
end
obs_var = ["Câ‚Šv(t)"]
obs_data_file = "working_example/res1.csv"
t_save, obs = get_obs(obs_data_file, obs_var)
model, param_vars, model_state, target_var = get_model()
u0 = [s => 0.0 for s in model_state]
tspan = (t_save[1], last(t_save))
lb = zeros(length(param_vars))
ub = zeros(length(param_vars)) .+ 10
prob = ODAEProblem(model, u0, tspan)
obs_vals = Array(obs)
p_init = rand(length(param_vars))
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(), callback=callback, maxiters=300)
Sometimes I need to point out the name of the output variable explicitly. I mean in the cost function I want to do something like that:
function loss(p::Vector{Float64})
pred_sol = solve(prob, Tsit5(), p=p, saveat=t_save)
pred = Array(pred_sol[model_state[1]]) # this doesn't work
# pred = Array(pred_sol') # this works
mse = mean((pred - obs_vals) .^ 2)
return mse
end
the model_state is:
julia> model_state
1-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
Câ‚Šv(t)
I have variable C but if I do pred = Array(pred_sol[C.v]) I get the same error.
The error:
ERROR: ArgumentError: invalid index: Câ‚Šv(t) of type Term{Real, Base.ImmutableDict{DataType, Any}}
It’s important to me to have that option to point name explicitly because sometimes it can be variable not from the state. For e.g. it can be R.v etc.
If I solve this outside of loss-function
sol = solve(prob, Tsit5(), p=p, saveat=t_save)
I can get sol[C.v] without any errors.
Any idea how I can use the names of variables in the loss function?
Thanks!