Optimization and Modelingtoolkit

Specific Domains Modelling & Simulations
1

Note: I am brand new to Julia and system modeling in general so I dont really have any idea of what I am doing haha.

Ive got a simple voltage divider system (my intended model is more complex but I figured Id start small) modeled using modelingtoolkit. What I am trying to do next is the inverse problem of taking a given set of outputs and finding the model parameters that best match those outputs (given voltages and currents, find the R’s) Is there any magic I can do with modelingtoolkit to accomplish this?

I tried this using SciPy and OMPython (python bindings for open modelica) and used one of the SciPy minimization functions to find the parameters. This eventually worked but took forever (OMPython, or maybe even modelica in general generates and compiles a bunch of intermediate C files prior to solving the model which slowed things down a lot).

Side question: Whats the best way to solve non-dynamic models such as this one? Since the outputs are all the same at every time step I can just use a timespan of (0.0, 0.0) and be fine, but is there a better way?

For reference, My Julia code for this model (based on the RC circuit example from modelingtoolkit) is below:

using ModelingToolkit, DifferentialEquations, Plots

@parameters t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0
    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 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

function ModelingToolkit.connect(::Type{Pin}, pins...)
    eqs = [
           0 ~ sum(p->p.i, pins) # KCL
          ]
    # KVL
    for i in 1:length(pins)-1
        push!(eqs, pins[i].v ~ pins[i+1].v)
    end

    return eqs
end

@named resistor1 = Resistor(R=200.0)
@named resistor2 = Resistor(R=100.0)
@named source = ConstantVoltage(V=10.0)
@named ground = Ground()

connections = [
    connect(source.p, resistor1.p)
    connect(resistor1.n, resistor2.p)
    connect(resistor2.n, source.n, ground.g)
]


@named voltage_divider_model = ODESystem(connections, t, systems=[resistor1, resistor2, source, ground])
sys = structural_simplify(voltage_divider_model)
tspan = (0.0, 2.0)
prob = ODEProblem(sys, [], tspan)
sol = solve(prob, Rodas4())
plot(sol, vars=[source.p.v, resistor1.v, resistor2.v])
2

Regarding the “took forever”: it sounds like your modelica model recompiled at every optimization step. Are you aware that if you make the parameters you vary a modelica parameter then you don’t have to recompile, just run the Model again with the updated parameter value?

3

Use a NonlinearSystem type instead. Note that in the near future we’re going to make the lowering process automatically detect this, but for now you have to specify it.

Yes, use https://diffeqflux.sciml.ai/dev/ for solving inverse problems with the generated ODEProblem.

1 Like
4

Initially, I was under the impression that yes I was only generating the model once, but you were right! I was rebuilding each time :joy:! I fixed that and its much faster now! Thanks for the insight.

1 Like
5

Hey Chris thanks for replying. I will try using DiffEqFlux and seeing if I can get something working.

Regarding the Nonlinear model, I tried making all my component models NonlinearSystems instead of ODESystems but ran into an issue when using the extend function with a nonlinear model with no independent variable (none of my states are a function of t). It encounters a problem when making the convert_system call in /src/systems/abstractsystem.jl:

function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys))
    T = SciMLBase.parameterless_type(basesys)
    iv = independent_variable(basesys)
    # iv is ::Nothing here
    sys = convert_system(T, sys, iv)
   ...

I was able to make it work by adding an if statement and calling it without the iv variable:

function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys))
    T = SciMLBase.parameterless_type(basesys)
    iv = independent_variable(basesys)
    if iv === nothing
        sys = convert_system(T, sys)
   else
       sys = convert_system(T, sys, iv)
   ...

Is that the correct fix or was the issue that I did something wrong on my end? Should I make an issue/submit a PR or something?

6

That looks like the correct fix and worth a PR.

7

awesome! I will submit it tonight when I get the chance

8

I was able to incorporate the DiffEqFlux package and solve the inverse problem : ). Very speedy compared to my original open modelica based solution

using ModelingToolkit, DifferentialEquations, DiffEqFlux, Plots

@parameters t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0
    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 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

function ModelingToolkit.connect(::Type{Pin}, pins...)
    eqs = [
           0 ~ sum(p->p.i, pins) # KCL
          ]
    # KVL
    for i in 1:length(pins)-1
        push!(eqs, pins[i].v ~ pins[i+1].v)
    end

    return eqs
end

@named resistor1 = Resistor(R=200.0)
@named resistor2 = Resistor(R=100.0)
@named source = ConstantVoltage(V=10.0)
@named ground = Ground()

connections = [
    connect(source.p, resistor1.p)
    connect(resistor1.n, resistor2.p)
    connect(resistor2.n, source.n, ground.g)
]


@named voltage_divider_model = ODESystem(connections, t, systems=[resistor1, resistor2, source, ground])
sys = structural_simplify(voltage_divider_model)
tspan = (0.0, 2.0)
prob = ODEProblem(sys, [], tspan)
true_sol = solve(prob, Rodas4())
plot(true_sol, vars=[source.p.v, resistor1.v, resistor2.v])

function loss(p)
    sol = solve(prob, Rodas4(), p=p)
    trues = [true_sol[resistor1.v]; true_sol[resistor2.v]; true_sol[resistor1.i]; ]
    sols = [sol[resistor1.v]; sol[resistor2.v]; sol[resistor1.i]; ]
    loss = sum(abs2, sols .- trues)
    return loss, sol
end

callback = function (p, l, pred)
    display(l)
    display(p)
    # Tell sciml_train to not halt the optimization. If return true, then
    # optimization stops.
    return false
end

ps = [5.0, 5.0, 1.0]
result_ode = DiffEqFlux.sciml_train(loss, ps,
                                    cb = callback)

sol = solve(prob, Rodas4(), p=result_ode.u)
plot(sol, vars=[source.p.v, resistor1.v, resistor2.v])
2 Likes