Is it possible to simulate a resistive network using ModelingToolkit?

Specific Domains Modelling & Simulations
1

Hi,

I’m trying to simulate a simple circuit of a constant voltage source connected to a resistor, using MTK. But I’m getting empty results when I solve the problem. Bear in mind this is a simple, toy problem, which I could solve using Symbolics or SymPy, but Iḿ looking ahead where I need to simulate the system with varying voltages, etc.

Below is the toy problem:

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
using ModelingToolkit: t_nounits as t

@mtkmodel simple_resistor begin
    @parameters begin
        V = 800.0
        R1 = 4.0*499.0e3
    end
    @components begin
        r1    = Resistor(R = R1)
        source = Voltage()
        constant = Constant(k = V)
        ground = Ground()
    end
    @equations begin
        connect(constant.output, source.V)
        connect(source.p, r1.p)
        connect(source.n, r1.n, ground.g)
    end
end

@mtkbuild sys = simple_resistor()

prob=ODEProblem(sys, [], (0, 10.0))

sol = solve(prob)

The result of solve is:

retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
  0.0
 10.0
u: 2-element Vector{Vector{Float64}}:
 []
 []

The return of using Pkg, Pkg.status() is

Status `~/.julia/environments/v1.11/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [a93c6f00] DataFrames v1.7.0
⌃ [82cc6244] DataInterpolations v6.5.2
⌃ [0c46a032] DifferentialEquations v7.14.0
  [0b43b601] Groebner v0.8.2
⌃ [7073ff75] IJulia v1.25.0
  [a98d9a8b] Interpolations v0.15.1
⌃ [961ee093] ModelingToolkit v9.49.0
  [16a59e39] ModelingToolkitStandardLibrary v2.17.0
⌅ [8913a72c] NonlinearSolve v3.15.1
⌃ [1dea7af3] OrdinaryDiffEq v6.89.0
⌃ [f0f68f2c] PlotlyJS v0.18.14
⌃ [91a5bcdd] Plots v1.40.8
  [24249f21] SymPy v2.2.0
  [2efcf032] SymbolicIndexingInterface v0.3.34

The return of versioninfo() is:

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 24 × AMD Ryzen 9 5900X 12-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 24 virtual cores)
2

This system is linear with no differential equations so it analytically solves it. You can use symbolic indexing to recover any of the values.

julia> sol[sys.r1.v]
2-element Vector{Float64}:
 800.0
 800.0

julia> sol(5.0; idxs = sys.r1.v)
800.0
1 Like
3

Many thanks! That solved the issue!

4

I’m finding another issue when I place a second resistor in series with r1. The updated script is below:

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
using ModelingToolkit: t_nounits as t

@mtkmodel simple_resistor begin
    @parameters begin
        V = 800.0
        R1 = 4.0*499.0e3
        R2 = 4.0*499.0e3
    end
    @components begin
        r1    = Resistor(R = R1)
        r2    = Resistor(R = R2)
        source = Voltage()
        constant = Constant(k = V)
        ground = Ground()
    end
    @equations begin
        connect(constant.output, source.V)
        connect(source.p, r1.p)
        connect(r1.n, r2.p)
        connect(source.n, r2.n, ground.g)
    end
end

@mtkbuild sys = simple_resistor()

prob=ODEProblem(sys, [], (0, 10.0))

sol = solve(prob)

The “prob=ODEProblem(sys, , (0, 10.0))” fails with:

┌ Warning: Did not converge after `maxiters = 0` substitutions. Either there is a cycle in the rules or `maxiters` needs to be higher.
└ @ Symbolics /home/ricardo/.julia/packages/Symbolics/6CYZh/src/variable.jl:528

{
	"name": "ArgumentError",
	"message": "ArgumentError: SymbolicUtils.BasicSymbolic{Real}[1.996e6r1₊i(t)] are either missing from the variable map or missing from the system's unknowns/parameters list.",
	"stack": "ArgumentError: SymbolicUtils.BasicSymbolic{Real}[1.996e6r1₊i(t)] are either missing from the variable map or missing from the system's unknowns/parameters list.

Stacktrace:
  [1] throw_missingvars_in_sys(vars::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/utils.jl:771
  [2] promote_to_concrete(vs::Vector{SymbolicUtils.BasicSymbolic{Real}}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/utils.jl:790
  [3] better_varmap_to_vars(varmap::Dict{Any, Any}, vars::Vector{SymbolicUtils.BasicSymbolic{Real}}; tofloat::Bool, use_union::Bool, container_type::Type, toterm::Function, promotetoconcrete::Nothing, check::Bool, allow_symbolic::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/problem_utils.jl:288
  [4] process_SciMLProblem(constructor::Type, sys::NonlinearSystem, u0map::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}, pmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}; build_initializeprob::Bool, implicit_dae::Bool, t::Nothing, guesses::Dict{Any, Any}, warn_initialize_determined::Bool, initialization_eqs::Vector{Any}, eval_expression::Bool, eval_module::Module, fully_determined::Bool, check_initialization_units::Bool, tofloat::Bool, use_union::Bool, u0_constructor::typeof(identity), du0map::Nothing, check_length::Bool, symbolic_u0::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/problem_utils.jl:508
  [5] process_SciMLProblem
    @ ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/problem_utils.jl:413 [inlined]
  [6] (NonlinearProblem{true})(sys::NonlinearSystem, u0map::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}, parammap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}; check_length::Bool, kwargs::@Kwargs{eval_expression::Bool, eval_module::Module})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/nonlinear/nonlinearsystem.jl:499
  [7] NonlinearProblem(::NonlinearSystem, ::Dict{SymbolicUtils.BasicSymbolic{Real}, Float64}, ::Vararg{Any}; kwargs::@Kwargs{eval_expression::Bool, eval_module::Module})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/nonlinear/nonlinearsystem.jl:490
  [8] ModelingToolkit.InitializationProblem{true, SciMLBase.AutoSpecialize}(sys::ODESystem, t::Int64, u0map::Dict{Any, Any}, parammap::Dict{Any, Any}; guesses::Dict{Any, Any}, check_length::Bool, warn_initialize_determined::Bool, initialization_eqs::Vector{Any}, fully_determined::Bool, check_units::Bool, kwargs::@Kwargs{eval_expression::Bool, eval_module::Module})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:1330
  [9] (ModelingToolkit.InitializationProblem{true})(::ODESystem, ::Int64, ::Vararg{Any}; kwargs::@Kwargs{guesses::Dict{Any, Any}, warn_initialize_determined::Bool, initialization_eqs::Vector{Any}, eval_expression::Bool, eval_module::Module, fully_determined::Bool, check_units::Bool})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:1232
 [10] ModelingToolkit.InitializationProblem(::ODESystem, ::Int64, ::Vararg{Any}; kwargs::@Kwargs{guesses::Dict{Any, Any}, warn_initialize_determined::Bool, initialization_eqs::Vector{Any}, eval_expression::Bool, eval_module::Module, fully_determined::Bool, check_units::Bool})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:1220
 [11] process_SciMLProblem(constructor::Type, sys::ODESystem, u0map::Vector{Any}, pmap::SciMLBase.NullParameters; build_initializeprob::Bool, implicit_dae::Bool, t::Int64, guesses::Dict{Any, Any}, warn_initialize_determined::Bool, initialization_eqs::Vector{Any}, eval_expression::Bool, eval_module::Module, fully_determined::Bool, check_initialization_units::Bool, tofloat::Bool, use_union::Bool, u0_constructor::typeof(identity), du0map::Nothing, check_length::Bool, symbolic_u0::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/problem_utils.jl:466
 [12] process_SciMLProblem
    @ ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/problem_utils.jl:413 [inlined]
 [13] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Any}, tspan::Tuple{Int64, Float64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, eval_expression::Bool, eval_module::Module, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:811
 [14] ODEProblem
    @ ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:799 [inlined]
 [15] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Any}, tspan::Tuple{Int64, Float64})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:799
 [16] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:786
 [17] (ODEProblem{true})(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:785
 [18] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:775
 [19] ODEProblem(::ODESystem, ::Vector{Any}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/eiNg3/src/systems/diffeqs/abstractodesystem.jl:774
 [20] top-level scope
    @ ~/jupyter/isolation_sensing/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W3sZmlsZQ==.jl:1"
}
5

In your second case, there is one unknown:

julia> @mtkbuild sys = simple_resistor()
Model sys:
Equations (1):
  1 standard: see equations(sys)
Unknowns (1): see unknowns(sys)
  r2₊i(t)
Parameters (6): see parameters(sys)
  R1 [defaults to 1.996e6]
  constant₊k [defaults to V]: Constant output value of block
  r2₊R [defaults to R2]: Resistance
  R2 [defaults to 1.996e6]
  r1₊R [defaults to R1]: Resistance
  V [defaults to 800.0]
Observed (21): see observed(sys)

So you need to supply an initial condition for it. In this case, you can supply anything you want and it will be treated as an initial guess for the solver.

julia> prob=ODEProblem(sys, [0.0], (0, 10.0))
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 0.0

julia> sol = solve(prob);

julia> sol(0, idxs=[sys.r1.v, sys.r2.v])
2-element Vector{Float64}:
 399.99999994373604
 400.00000005626396
6

And you might as well use NonlinearProblem instead of ODEProblem in situations like this with no ODEs.