Is it possible to simulate a resistive network using ModelingToolkit?

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)

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

Many thanks! That solved the issue!

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"
}

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

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

Many thanks for the tip!

By the way, is there a method to print the nice model information you quoted there?

Model sys:
Equations (1):
  1 standard: see equations(sys)
Unknowns (1): see unknowns(sys)
...

cool, thanks, good to know!

I checked in this link that the @mtkmodel macro does not support NonlinearSystem yet, so I’ll stick with the ODEProblem for the time being.

If you mean exactly what I quoted there, that is the show method as used by the REPL display. I evaluated the expression from your code in the REPL and quoted what was printed.

If you mean digging in to more details, you can use unknowns(sys), equations(sys), or observed(sys) to show the full list of each.

1 Like

You just do NonlinearProblem on an ODE and it’ll work.

By that, I guess you mean the following, which works fine:

using NonlinearSolve
prob_NL = NonlinearProblem(prob)

NonlinearProblem with uType Nothing. In-place: true
u0: nothing

but when I solve, it returns empty:

sol_NL = solve(prob_NL)
retcode: Success
u: Float64[]

show(sol_NL)
Float64[]

Yes because the solution is trivial (analytically solved, no numeric), but again symbolic indexing works fine.

by that, I guess you mean:

sol_NL[sys.r1.v]

but then I get:

BoundsError: attempt to access Tuple{Vector{Float64}, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}} at index [3]

Stacktrace:
 [1] getindex
   @ ./tuple.jl:31 [inlined]
 [2] macro expansion
   @ ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:162 [inlined]
 [3] macro expansion
   @ ./none:0 [inlined]
 [4] generated_callfunc
   @ ./none:0 [inlined]
 [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#5964805160111424296"), :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x280425a4, 0x4a32d630, 0x943b2a62, 0x38dff8e0, 0x05cd040d), Nothing})(::Vector{Float64}, ::MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}})
   @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/M9ZX8/src/RuntimeGeneratedFunctions.jl:150
 [6] (::SymbolicIndexingInterface.TimeIndependentObservedFunction{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#5964805160111424296"), :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x280425a4, 0x4a32d630, 0x943b2a62, 0x38dff8e0, 0x05cd040d), Nothing}})(::SymbolicIndexingInterface.NotTimeseries, prob::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Nothing, NonlinearProblem{Nothing, true, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#37#46"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#911"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x47755f63, 0x54879625, 0xa6ce3993, 0x677dc31f, 0x0ea26d1b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xe4768855, 0xc75a2b9e, 0xbdbf2985, 0x3b9edfca, 0xfaebd888), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing, Nothing, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, Nothing, Nothing, Nothing, Nothing, Nothing})
   @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/sGNlE/src/state_indexing.jl:142
 [7] (::SymbolicIndexingInterface.TimeIndependentObservedFunction{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#5964805160111424296"), :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x280425a4, 0x4a32d630, 0x943b2a62, 0x38dff8e0, 0x05cd040d), Nothing}})(prob::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Nothing, NonlinearProblem{Nothing, true, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#37#46"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#911"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x47755f63, 0x54879625, 0xa6ce3993, 0x677dc31f, 0x0ea26d1b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xe4768855, 0xc75a2b9e, 0xbdbf2985, 0x3b9edfca, 0xfaebd888), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing, Nothing, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, Nothing, Nothing, Nothing, Nothing, Nothing})
   @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/sGNlE/src/value_provider_interface.jl:166
 [8] getindex(A::SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Nothing, NonlinearProblem{Nothing, true, MTKParameters{Vector{Float64}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#37#46"{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#911"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x47755f63, 0x54879625, 0xa6ce3993, 0x677dc31f, 0x0ea26d1b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xe4768855, 0xc75a2b9e, 0xbdbf2985, 0x3b9edfca, 0xfaebd888), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing, Nothing, Nothing, Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing}, @Kwargs{}, SciMLBase.StandardNonlinearProblem}, Nothing, Nothing, Nothing, Nothing, Nothing}, sym::Num)
   @ SciMLBase ~/.julia/packages/SciMLBase/hJh6T/src/solutions/solution_interface.jl:117
 [9] top-level scope
   @ ~/jupyter/isolation_sensing/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_W4sZmlsZQ==.jl:1

Oh, when we nest we problems we might need to pass on SII @cryptic.ax .

Using SteadyStateProblem here is fine:

using NonlinearSolve
prob_NL = SteadyStateProblem(sys, [sys.r2.i => 0.0])
sol_NL = solve(prob_NL, TrustRegion())
julia> sol_NL[sys.r2.i]
0.0002004008016032064

Oh, when we nest we problems we might need to pass on SII @cryptic.ax .

That’s not the issue here. The system inside sol_NL is an ODESystem. It generates observed functions of the form f(u, p, t), but since sol_NL is not time dependent, SII calls the function has f(u, p) which leads to the above issue. In SII terminology, the index provider is time-dependent but the value provider isn’t.

Ahh we neeed to make that conversion better

Cool, thanks for the reply!

I opened a bug issue here, even though this seems to me more like an improvement (there was no “improvement” issue type):