Tool for composable circuit simulation?

Specific Domains Modelling & Simulations
,
1

I’m looking for a (free) tool to do some rather basic electrical circuit simulations describing a unit with a couple of current sources, a couple of resistors, switches and some external connections. Having that I would like to instantiate several of the units, controlling aspects of each one separately, add external loads and connections between the units and try out scenarios. Essentially a DC simulation.

I have done a little but of exploration and found SIMS.jl (now FunctionalModels.jl) and ACME.jl
My feeling so far is that FunctionalModels.jl seems more promising, but according to its documentation it is being re-aligned with ModelingToolkit.jl, documentation is not up to date and last commit was 2021…
So what I wonder is what Julians typically use in similar cases (I guess I have no unique needs)

2 Likes
2

Based on your description, it sounds like ModelingToolkit.jl and ModelingToolkitStandardLibrary.jl should have you covered. There are some basic Circuit examples in the stdlib docs, do they look like something you’re looking for?

3

It sure looks like it should do fine, but I have to admit that using and extending it is a bit more of a challeng than I initially expected. Is anyone aware of some examples (both defining eg an open/close switch and creating a small “instantiable” component would be helpful) beyond the ubiquitous RC-circuit?

4

How would you describe your switch? Does it switch at a particular, pre-defined time? Or does it trigger when the input reaches a certain value?

5

Defined as a bool at construction time, (isClosed=true or isClosed=false) comparable to R in the Resistor model. Then using ifelse to define the equations for closed and open states.
It seems I run into issues with Symbolics.jl and that @parameters wrap things as Num.
The ifelse work before I define isClosed as a @parameter, but not after.
I guess that either passing 0 or 1 as parameter and comparing that in the ifelse is a possible but rather ugly an “un-julian” workaround, as is using a Resistor with a very high or very low R - but I would prefer do do it “right”.

6

Could you post an example where you have tried it “the right way”? Maybe we can figure out what the problems with Symbolics is?

7

Here is a seemingly working example.
Moving the equations definition to before the parameter definition makes the ifelse evaluate OK.
I guess that there are still problems left here…
To go back to your question about how the switches are defined, I guess it would be better to be able to define the state (open/closed) with parameter when the ODAEProblem is defined? Not sure if that would work with the equations being defined as they are now though…
Anyway, code below, laced with quite a few debug printouts in the IdealSwitch
I realise the isequal(isClosed,True) is redundant at start, but it is needed (but failing)
after the @parameter step.

using Revise
using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Electrical: OnePort, Ground
using OrdinaryDiffEq
using Plots
using IfElse: ifelse

@parameters t

# No switch... need to define that.

@component function IdealSwitch(; name, isClosed=true)
    println("Creating switch, isClosed=$(isClosed)")
    println(" Before param assign: $(isClosed)  -> t? $(isequal(isClosed,true)) -> f? $(isequal(isClosed, false))")
    println("typeof(isClosed): $(typeof(isClosed))")
    @named oneport = OnePort()
    @unpack v, i = oneport
    @variables closed  # Dummy to include in eqs
    eqs = [ v ~ ifelse(isequal(isClosed, true), 0, v)
            i ~ ifelse(isequal(isClosed, true), i, 0)
            closed ~ isequal(isClosed, true)
    ]
    println(" eqs: $(eqs)")
    
    ps = @parameters isClosed = isClosed
    println("ps: $ps")

    # Original placement of eqs definition - ie _after_ params: here isequal always evaluate to false
    #eqs = [ v ~ ifelse(isequal(isClosed, true), 0, v)
    #        i ~ ifelse(isequal(isClosed, true), i, 0)
    #        closed ~ isequal(isClosed, true)
    #]
    #println(" eqs: $(eqs)")

    println("After param assign: $(isClosed)  -> t? $(isequal(isClosed,true)) -> f? $(isequal(isClosed, false))")
    println("  typeof(isClosed): $(typeof(isClosed))")
    println("  Switch is $(ifelse(isequal(isClosed, true), "CLOSED", "OPEN"))")
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

@named s_closed = IdealSwitch(;isClosed=true)
# equations(s_closed)
# @named s_open = IdealSwitch(;isClosed=false)
# equations(s_open)

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

@named vcc = ConstantVoltage()
@named resistor = Resistor(R=2)
@named switch = IdealSwitch(isClosed=true)
@named ground = Ground()

eqs = [
    connect(vcc.p, resistor.p)
    connect(resistor.n, switch.p)
    connect(vcc.n, switch.n, ground.g)]

@named switch_sys = ODESystem(eqs, t;
    systems = [vcc, resistor, switch, ground])
sys = structural_simplify(switch_sys)
prob = ODAEProblem(sys, [], (0, 10.0))
sol = solve(prob, Tsit5())

plot(sol, idxs = [resistor.v, resistor.i],
    title = "Switch Demonstration",
    labels = ["Resistor Voltage" "Resistor Current"])
8

I’d suggest dropping isClosed parameter and closed variable and simplifying IdealSwitch to

@component function IdealSwitch(; name, isClosed=true)
    @named oneport = OnePort()
    @unpack v, i = oneport
    eqs = [ v ~ ifelse(isClosed, 0, v)
            i ~ ifelse(isClosed, i, 0)
    ]
    extend(ODESystem(eqs, t, [], []; name = name), oneport)
end

However I can see why isClosed can be a parameter. In that case, define the @parameter isClosed = isClosed before equations and use ifelse(ModelingToolkit.getdefault(isClosed), 0, v)

In commented equations in the implementaion here, a boolean is compared with symbolic var, which will always return false.

@parameters a = true b = true

isequal(a, true) # false. As `a` is a symbolic var while true is a bool
isequal(a, b) # false; as a and b are different symbolic vars
isequal(ModelingToolkit.getdefault(a), true) # true; as default value is a boolean
isequal(a, a) # true
9

Here is a bit of a hack that I think works. You have to turn the Bool into a term (Num in Symbolics) to get the correct ifelse dispatch. Th rest is carefully arranging the equations to make simplification happy.

@component function IdealSwitch(; name, isClosed=true)
    @named n = Pin()
    @named p = Pin()
    ps = @parameters isClosed=isClosed
    eqs = [ 0 ~ ifelse(isClosed == true, n.v - p.v, n.i)
            0 ~ ifelse(isClosed == true, n.i + p.i, p.i)
    ]
    ODESystem(eqs, t, [], ps; name, systems=[p, n])
end

FYI, IfElse is no longer needed since #27343 was merged.

10

And a more complete example just for fun.

using Revise
using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks: RealInput, Step
using ModelingToolkitStandardLibrary.Electrical
using OrdinaryDiffEq
using Plots

@parameters t

@component function IdealSwitch(; name, isClosed=true)
    @named n = Pin()
    @named p = Pin()
    ps = @parameters isClosed=isClosed
    eqs = [ 0 ~ ifelse(isClosed == true, n.v - p.v, n.i)
            0 ~ ifelse(isClosed == true, n.i + p.i, p.i)
    ]
    ODESystem(eqs, t, [], ps; name, systems=[p, n])
end

"""
A switch whose state is controlled by a RealInput pin.

The switching threshold is 0.5, open when input is less closed when input is greater.
"""
@component function ControlledSwitch(; name)
    @named n = Pin()
    @named p = Pin()
    @named u = RealInput()
    eqs = [ 0 ~ ifelse(u.u > 0.5, n.v - p.v, n.i)
            0 ~ ifelse(u.u > 0.5, n.i + p.i, p.i)
    ]
    ODESystem(eqs, t, [], []; name, systems=[p, n, u])
end

function plot_ideal_switch(;isClosed=true)
    @named vcc = Voltage()
    @named resistor = Resistor(R=2)
    @named switch = IdealSwitch(;isClosed)
    @named ground = Ground()

    ps = @parameters supply_voltage=12.0

    eqs = [
        connect(vcc.p, resistor.p)
        connect(resistor.n, switch.p)
        connect(vcc.n, switch.n, ground.g)
        vcc.V.u ~ supply_voltage
       ]

    @named switch_sys = ODESystem(eqs, t, [], ps;
        systems = [vcc, resistor, switch, ground])
    sys = structural_simplify(switch_sys)
    prob = ODAEProblem(sys, [], (0, 10.0))
    sol = solve(prob, Tsit5())

    plot(sol, idxs = [resistor.v, resistor.i],
        title = "Switch Demonstration",
        labels = ["Resistor Voltage" "Resistor Current"])
end

"""
With no dynamics (no time-derivative terms), no time stepping is performed.
This example adds a parallel RC circuit some stepping happens and the switch
time can be plotted.
'"""
function plot_controlled_switch()
    @named vcc = Voltage()
    @named resistor = Resistor(R=2)
    @named r2 = Resistor(R=75e3)
    @named cap = Capacitor(C=10e-6)
    @named switch = ControlledSwitch()
    @named ground = Ground()
    @named control = Step(;start_time=5.0)

    ps = @parameters supply_voltage=12.0

    eqs = [
        connect(vcc.p, resistor.p, r2.p)
        connect(resistor.n, switch.p)
        connect(r2.n, cap.p)
        connect(vcc.n, switch.n, cap.n, ground.g)
        vcc.V.u ~ supply_voltage
        connect(switch.u, control.output)
    ]

    @named switch_sys = ODESystem(eqs, t, [], ps;
        systems = [vcc, resistor, switch, r2, cap, ground, control])
    sys = structural_simplify(switch_sys)
    prob = ODAEProblem(sys, [], (0, 10.0))
    sol = solve(prob, Tsit5())

    plot(sol, idxs = [resistor.v, resistor.i, control.output.u],
        title = "Switch Demonstration",
        labels = ["Resistor Voltage" "Resistor Current" "Switch"])
end
11

A big thank you to all who have helped so far!
ModelingToolkit sure is a powerful thing, but to be honest with the mixing in of Symbolics I feel a little bit as falling fown the rabbit hole. Things are not as clear as they usually are in Julia, at least not to me.
With the help provided so far its fine to model an IdealSwitch.
But as in the seminal Julia blogpost “I want more”. Specifically I would like to set up a nested structure of components - including switches - defining the structure of the problem, do the structural_simplify step and then solve for different parameter (eg switch position) combinations. Essentially “wire it all up, then play with the switches and observe the results”. I would also like to define the state ot the switches (or parameters in general) on the parent, ie composed level, ie the SwR component below and as descibed Here (Inheritance and Combine)
The code below fails, not finding swr.i in the solution when I’m trying to plot it. I had kind of expected to be able to plot anything that is listed as an observable in the structurally simplified system but that is apparently not the case. I guess there is a way to list the “indexes” of a solution, but I have failed to find how.

Am I doing something fundamentally wrong here or is this just a case of the devil in the details? I did not expect struggling so much with composition in ModelingToolkit.

using Revise
using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Electrical: OnePort, Ground
using OrdinaryDiffEq
using Plots

@parameters t

@component function IdealSwitch(; name, isClosed=true)
    @named n = Pin()
    @named p = Pin()
    ps = @parameters isClosed=isClosed
    eqs = [ 0 ~ ifelse(isClosed == true, n.v - p.v, n.i)
            0 ~ ifelse(isClosed == true, n.i + p.i, p.i)
    ]
    ODESystem(eqs, t, [], ps; name, systems=[p, n])
end


# Simple aggregated component - Switch + resistor in series
@component function SwR(; name, sw_state=true, R=1.0)
    ps = @parameters (sw_state=sw_state), (R=R) 
    @variables i(t), v(t)

    @named p = Pin()
    @named sw = IdealSwitch()
    @named resistor = Resistor()
    @named n = Pin()

    eqs = [
        connect(p, sw.p)
        connect(sw.n, resistor.p)
        connect(resistor.n, n)
        i ~ p.i
        v ~ p.v - n.v
    ]  
    # Trying to couple the swr level sw_state to the sw (IdealSwicth)'s isClosed parameter
    par_defaults = [resistor.R => R, sw.isClosed => sw_state]  
    ODESystem(eqs, t, [], [R, sw_state];
        systems = [p, sw, resistor, n], name=name, defaults=par_defaults)
end

@named swr = SwR(;sw_state=true, R=2.0)
@named vcc = Voltage()
@named gnd = Ground()
eqs = [
    vcc.V.u ~ 5
    connect(vcc.p, swr.p)
    connect(vcc.n, swr.n, gnd.g)
]
@named swr_sys = ODESystem(eqs, t; systems = [vcc, swr, gnd])

sys = structural_simplify(swr_sys)
prob = ODAEProblem(sys, [], (0, 10.0), [swr.sw_state => true])
sol = solve(prob, Tsit5())
plot(sol, idxs = [swr.i],
    title = "Compound Switch Demonstration",
    labels = ["Switch branch Current"])
12

i is simply missing from the list of variables of the ODESystem, I think.

13 
# Simple aggregated component - Switch + resistor in series
@component function SwR(; name, sw_state=true, R=1.0)
    ## you assign, ps here, might as well use it when creating the ODESystem
    ps = @parameters (sw_state=sw_state), (R=R) 
    ## Also need to keep track of variables
    vs = @variables i(t), v(t)

    @named p = Pin()
    @named sw = IdealSwitch()
    @named resistor = Resistor()
    @named n = Pin()

    eqs = [
        connect(p, sw.p)
        connect(sw.n, resistor.p)
        connect(resistor.n, n)
        i ~ p.i
        v ~ p.v - n.v
    ]  
    ## I would leave this out as it is redundant with the defaults assigned above when creating the parameters.
    # Trying to couple the swr level sw_state to the sw (IdealSwicth)'s isClosed parameter
    ## par_defaults = [resistor.R => R, sw.isClosed => sw_state]
    # Use the assigned arrays of variables and parameters when creating the system.
    ODESystem(eqs, t, vs, ps; systems = [p, sw, resistor, n], name=name)
end

When nesting components like this it is often helpful to move parameters up and down the hierarchy using the scope tools.

14

After looking at JuliaCon 2023 presentation including @baggepinnen and Venkateshprasad Bhat here, I decided to take a step back and try to use the @mtkmodel macro and friends, as this is described as the “new way”, looks cleaner and potentially also enables the use of GUI builder for the resulting component(s).
I also wanted to try to model the IdealSwitch directly as an extension of OnePort as that seems to be more in line with the standard lib.
I further wanted to extend the example a bit to look more like what I’m actually trying to build, which is a circuit with a Voltage source, a current-limiting resistor and two switches in series. I finally want to define parameters at the Gadget level, ie I want to have two booleans at the Gadget level to set, rather than one boolean on each switch, one level further down in the hierarchy.
Ultimately I’d like to be able to instantiate more than one Gadget, connect them together and to external loads and then set up problems parameterised something like
[gadget1.b1 => true, gadget1.b2 => false, gadget2.b1 => true, gadget2.b2 => true] etc.

My code is below, but it has at least two problems:
When I try to simplify a created Gadget, I get an unbalanced system error:

ERROR: ExtraEquationsSystemException: The system is unbalanced. There are 27 highest order derivative variables and 29 equations.
More equations than variables, here are the potential extra equation(s):
 0 ~ r₊i(t) - s1₊i(t)
 0 ~ s2₊i(t) - s1₊i(t)

Since both the indicated equations include the IdealSwitch that is the likely culprit but I fail to see the error. Defined the way it is now it’s equations should end up as either as redundant i ~ i or trivial i ~ 0 depending on the state of the switch - and similar for v. Essentially just a connection or a break. I’d though the structural simplification would deal with that.

I also expect that there are issues with the way I try to apply the scope tools as suggested by @contradict in the previous post. I looked at the tools page but don’t quite see how the example there should be applied to composing a system like I do here. I hope the code is clear bout how I expect it to work and hope somene can explain how it actually works :wink: Since this thread seems to get a bit of views, I guess there are more people than me interested in this.

using ModelingToolkit
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Electrical: Pin, OnePort, Ground
using OrdinaryDiffEq
using Plots

@parameters t

@mtkmodel IdealSwitch begin
    @parameters begin
        isClosed = true, [description = "State of switch"]
    end
    @extend v, i = one_port = OnePort()
    @equations begin
        i ~ ifelse(isClosed == true, i, 0)
        v ~ ifelse(isClosed == true, 0, v)
    end
end

# Simple test circuit _and_ dumbed-down version of actual application
# Real application will eg define 3 connectors located between r-s1, s1-s2 and s2-src
#  for connection of externl loads.
@mtkmodel Gadget begin
    @parameters begin
        b1 = true, [description = "State of switch s1"]
        b2 = true, [description = "State of switch s2"]
    end
    @components begin
        src = Voltage()
        r = Resistor(R=2)
        s1 = IdealSwitch()
        s2 = IdealSwitch()
        gnd = Ground()
    end
    @equations begin
        src.V.u ~ 10
        connect(src.p, r.p)
        connect(r.n, s1.p)
        connect(s1.n, s2.p)
        connect(s2.n, src.n, gnd.g)
    end
    # Trying to "promote" switch states to be parameters on Gadget level
    begin
        b1 = ParentScope(s1.isClosed)
        b2 = ParentScope(s2.isClosed)
    end
end

@named gadget = Gadget()
equations(gadget)
sys = structural_simplify(gadget)
prob = ODAEProblem(sys, [], (0, 10.0), [gadget.b1 => true, gadget.b2 => true])
sol = solve(prob, Tsit5())
plot(sol, idxs = [gadget.r.i])
15

I think the problem here is that structural_simplify cannot look inside ifelse. Looking at the equations defined by gadget, for example s1₊i(t) is defined to be some opaque function of itself (the ifelse from IdealSwitch), that would be fine if that were the only equation determining the value of that variable, the DAE solvers are good at this kind of thing. The problem is there is also another simple linear expression of other variables that defines that variable. It is trying to tell you that with the extra equations error, but those are often difficult to figure out why it picked the particular once it did to be extra. In this case, it is a bit easier to see when you print out the full system of equations and see that some values in these equations are also defined as fixed-point functions of themselves. Of course, that can’t be the whole story, but I still don’t really understand how simplify works well enough to explain exactly why this fails. In my IdealSwitch, I tried to avoid these cases so that there is an obvious transform to a DAE. I couldn’t figure out how to do that and extend an OnePort, however there is no reason old and new style components shouldn’t integrate nicely.

ParentScope didn’t do what you expect because you have to use the result of it in the subcomponent. The variable has already been namespaced when you extract it from the subcomponent. You probably don’t want to move the isClosed values up the hierarchy since they need separate names anyway. The block at the end re-assignes the names b1 and b2 to the values of the parent-scoped sn₊isClosed parameters and then they are discarded as you can see in the output of parameters(gadget). If you want to create new parameters that provide the default values for isClosed, you can use them when defining the switch

s1 = IdealSwitch(;isClosed=b1)

No namespace trickery needed.

16

Using the following defnition for IdealSwitch seems to work:

@mtkmodel IdealSwitch begin
    @parameters begin
        isClosed = true, [description = "State of switch"]
    end
    @components begin
        p = Pin()
        n = Pin()
    end
    @equations begin
        0 ~ ifelse(isClosed == true, n.v - p.v, n.i)
        0 ~ ifelse(isClosed == true, n.i + p.i, p.i)
    end
end

However, I still run into problems when i try to pass parameters while setting up the problem.
I have updated to pass on the b1 and b2 parametrs as suggested: s1 = IdealSwitch(;isClosed=b1) etc but it seems that the default value set in the @parameters section in always used, regardless what i try to pass as parameters to the ODEAProblem. Not giving a default value at parameter definition gives a complaint about missing variables. So I assume that I am not passing the right keys for parameters to the problem creation?

@mtkmodel Gadget begin
    @parameters begin
        b1 = true, [description = "State of switch s1"]
        b2 = true, [description = "State of switch s2"]
    end
    @components begin
        src = Voltage()
        r = Resistor(R=2)
        s1 = IdealSwitch(; isClosed=b1)
        s2 = IdealSwitch(; isClosed=b2)
        gnd = Ground()
    end
    @equations begin
        src.V.u ~ 10
        connect(src.p, r.p)
        connect(r.n, s1.p)
        connect(s1.n, s2.p)
        connect(s2.n, src.n, gnd.g)
    end
end

@named gadget = Gadget()
sys = structural_simplify(gadget)

# Run with both switches closed, expecting current 10/2 = 5 Amps. Works.
prob = ODAEProblem(sys, [], (0, 10.0), [gadget.b1 => true, gadget.b2 => true])
sol = solve(prob, Tsit5())
sol[gadget.r.i, 1]

# Run with one switch open, one closed, expecting current 0 (open circuit) Fails (still get 5...)
prob = ODAEProblem(sys, [], (0, 10.0), [gadget.b1 => true, gadget.b2 => false])
sol = solve(prob, Tsit5())
sol[gadget.r.i, 1]

Oh, I’d like to add a `modelingtoolkit* tag to this post but I don’t seem to be able to do that post–creation. If anyone can, please do.

17

This is a namespacing problem that I don’t understand.

It is easy to see that all the parameters are there:

julia> parameters(sys)
5-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 b1
 b2
 r₊R
 s1₊isClosed
 s2₊isClosed

And that the defaults are what you wanted:

julia> ModelingToolkit.get_defaults(sys)
Dict{Any, Any} with 10 entries:
  b1          => true
  s2₊isClosed => b2
  r₊i(t)      => 0.0
  b2          => true
  src₊v(t)    => 0.0
  s1₊isClosed => b1
  src₊V₊u(t)  => 0.0
  r₊R         => 2
  r₊v(t)      => 0.0
  src₊i(t)    => 0.0

But setting parameters the intuitive way you did results in no change to the defaults (compare to parameters(sys) to see the names):

julia> prob.p
5-element Vector{Float64}:
 1.0
 1.0
 2.0
 1.0
 1.0

Using p=[:b1=>true, :b2=>false] also does not work. I think the problem is that

julia> gadget.b1
gadget₊b1

instead of just b1.

I have found two workarounds.
The first is to find the indices into parameters(sys) with something like

"""
getparam(sys, var)

Retrieve a parameter object by symbol from a system.
"""
function getparam(sys::ODESystem, name)
    ps = parameters(sys)
    idx = findfirst(ps) do p
        string(p) == string(name)
    end
    ps[idx]
end

Which can be called like p=[getparam(sys, :b1)=>true, getparam(sys, :b2)=>false].

The other is to create new parameters with the correct names.

@paramters b1 b2
p = [b1 => true, b2 => false]
prob = ODAEProblem(sys, [], (0, 10.0), p)

After which you can see the values have taken effect:

julia> prob.p
5-element Vector{Float64}:
 1.0
 0.0
 2.0
 1.0
 0.0
18

Good job finding not just one, but two workarounds!
Still, this means defining the system, running struct_simplify and then finding the namea of the parameter(s) to set. That feels as a somewhat broken workflow.
In the tutorial here it is shown that variables in the original system can still be accessed in the solution, even if they were optimised away by struct_simplify. Can this be seen as the same thing (but missing) for parameters?

@baggepinnen At the end of the presentation linked above you mention a couple of things in Future plans - including “if-else conditions” and “Symbolic defaults for subpomponents”. This sounds a lot like what we are struggling with here - are we (I?) just trying this too early? What I try to do is very much like the DC motor in your last example (the example shows multi-domain elecrical-mechanical modeling of a DC motor with a mechanical load) but if one builds a car with two DC motors I guess one would not like to need to run structural-simplify to find out what the parameters for each DC motor are called and then edit the source to assign using the right parameter names? (Sorry to keep tagging you, but you seem to be the one closest to ModelingToolkit development in this thread)

19

You are experiencing an unfortunate namespacing issue. When you refer to gadget.b1, this is adding the namespace gadget to the parameter, but the parameters are expected without the top-level namespace attached to them. There are two ways to work around this, one is to

@unpack b1, b2 = gadget

which grabs the parameters b1, b2 without namespace from gadget, the other one is to complete the model, after which you can refer to parameters within it without appending the top-level namespace. The solution using this method is shown below

using Test
# Run with one switch open, one closed, expecting current 0 (open circuit) 
cgadget = complete(gadget)
prob = ODAEProblem(sys, [], (0, 10.0), [cgadget.b1 => true, cgadget.b2 => false])
sol = solve(prob, Tsit5())
@test sol[gadget.r.i, 1] ≈ 0 atol=1e-6

This is a very common mistake to make, and we have an issue open about it

1 Like
20

Thanks, that makes some sense.
However, using the complete way and upping the bets with a system including two gadgets (just placed in the same system, not even connected trips me up again:

@mtkmodel G2 begin
    @components begin
        g1 = Gadget()
        g2 = Gadget()
    end
end

@named gsys = G2()
cgsys = complete(gsys)
sys = structural_simplify(gsys)

# Run g1 with one switch open, one closed, expecting current 0 (open circuit) 
prob = ODAEProblem(sys, [], (0, 10.0), [cgsys.g1.b1 => true, cgsys.g1.b2 => false])
collect(zip(parameters(sys), prob.p))
sol = solve(prob, Tsit5())
sol[gsys.g1.r.i, 1]

This works at least as far as parameter setting, as " (g1₊b2, 0.0)" can be seen in the output from the parameter dumping. In fact this may well have worked perfectly, but I am so far unable to access the interesting variable in the solution. I have tried all of sys, gsys, cgsys as prefixes (and no prefix at all) when indexing into the solution - all tries failing. Intuitively I’d expect the notation in the code above to be the correct one?
General question is actually how to get all the usable indexes for an ODESolution?