DiffEqParamEstim with Modellingtoolkit

Specific Domains Modelling & Simulations
7

Of course below should demonstrate my problem just taken the RC-Circuit from the MTK tutorial

using ModelingToolkit, Plots, DifferentialEquations, DiffEqParamEstim

@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function ConstantVoltage(;name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           V ~ v
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = ConstantVoltage(V=V)
@named ground = Ground()

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
          connect(capacitor.n, ground.g)
         ]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
                          [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
      capacitor.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
#sol = solve(prob, Tsit5(), save_idxs = [capacitor.v]) can not index a solution using symbolics, can also not index at observed values 
t = collect(range(0,stop=10,length=200))
sol = solve(prob, Tsit5(), saveat = t)

plot(sol, idxs = [resistor.p.i,capacitor.v ]) #here resisitor is observed variable 

V = [Vector{Float64}(undef,2) for _ in 1:200]
for i in 1:200
    V[i][1] = sol[resistor.p.i][i]
end
for i in 1:200
    V[i][2] = sol[capacitor.v ][i]
end
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(V[:][i] + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
cost_func = build_loss_objective(prob,Tsit5(),L2Loss(t,data),save_idxs = [resistor.p.i,capacitor.v ])
prange = 0.8:0.01:1.2
plot(prange, [cost_func([i, 1.0, 1.0]) for i in prange], lw = 3, title = "R Liklihood (Rest at base state)", xlabel = "Emax", ylabel = "Objective Function Value - LS")
## Working cost function ##
randomized = VectorOfArray([(sol(t[i]) + .01randn(1)) for i in 1:length(t)])
data = convert(Array,randomized)
cost_func = build_loss_objective(prob,Tsit5(),L2Loss(t,data))
prange = 0.8:0.01:1.2
plot(prange, [cost_func([i, 1.0, 1.0]) for i in prange], lw = 3, title = "R Liklihood (Rest at base state)", xlabel = "Emax", ylabel = "Objective Function Value - LS")

## Appears the symbolic indexing in the cost function causes problems