Problems with stiff solvers for MTK model

ohmsweetohm1 · April 20, 2021, 12:44pm

I have a MTK model (copy paste from docs):

using ModelingToolkit, Plots, DifferentialEquations

@parameters t

function Pin(;name)
    @variables v(t) i(t)
    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
end

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

function Resistor(;name, R = 1.0)
    val = R
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters R
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           v ~ p.i * R
          ]
    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
end

function Capacitor(; name, C = 1.0)
    val = C
    @named p = Pin()
    @named n = Pin()
    @variables v(t)
    @parameters C
    D = Differential(t)
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           D(v) ~ p.i / C
          ]
    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
end

function ConstantVoltage(;name, V = 1.0)
    val = V
    @named p = Pin()
    @named n = Pin()
    @parameters V
    eqs = [
           V ~ p.v - n.v
           0 ~ p.i + n.i
          ]
    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
end

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

function connect_pins(ps...)
    eqs = [
        0 ~ sum(p->p.i, ps) # KCL
    ]
    # KVL
    for i in 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end

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

@named rc_model = ODESystem(rc_eqs, t,
    systems=[resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
    capacitor.v => 0.0
    capacitor.p.i => 0.0
]
prob = ODAEProblem(sys, u0, (0, 10.0))

Which runs fine with Tsit5

sol = solve(prob, Tsit5())

But errors with Rodas5, or Rosenbrock23 (I did not test others)

sol = solve(prob, Rodas5())

ERROR: LoadError: DimensionMismatch("row lengths must match")
Stacktrace:
  [1] hcat(::BitMatrix, ::BitMatrix)
    @ Base .\bitarray.jl:1792
  [2] generate_chunked_partials(x::Vector{Float64}, colorvec::Vector{Int64}, #unused#::Val{1})
    @ SparseDiffTools ~\.julia\packages\SparseDiffTools\HWfTE\src\differentiation\compute_jacobian_ad.jl:57
  [3] SparseDiffTools.ForwardColorJacCache(f::Function, x::Vector{Float64}, _chunksize::Nothing; dx::Nothing, colorvec::Vector{Int64}, sparsity::SparseArrays.SparseMatrixCSC{Bool, Int64})
    @ SparseDiffTools ~\.julia\packages\SparseDiffTools\HWfTE\src\differentiation\compute_jacobian_ad.jl:34
  [4] build_jac_config
    @ ~\.julia\packages\OrdinaryDiffEq\la32G\src\derivative_wrappers.jl:122 [inlined]
  [5] build_jac_config
    @ ~\.julia\packages\OrdinaryDiffEq\la32G\src\derivative_wrappers.jl:117 [inlined]
  [6] alg_cache(alg::Rodas5{0, true, DefaultLinSolve, DataType}, u::Vector{Float64}, rate_prototype::Vector{Float64}, uEltypeNoUnits::Type, uBottomEltypeNoUnits::Type, tTypeNoUnits::Type, uprev::Vector{Float64}, uprev2::Vector{Float64}, f::ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#265"), Symbol("##arg#266"), Symbol("##arg#267"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0415a075, 0xa1d706f6, 0x283e10f1, 0xf60795b9, 0xb10e7014)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#53"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, t::Float64, dt::Float64, reltol::Float64, p::Vector{Float64}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~\.julia\packages\OrdinaryDiffEq\la32G\src\caches\rosenbrock_caches.jl:586
  [7] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#265"), Symbol("##arg#266"), Symbol("##arg#267"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0415a075, 0xa1d706f6, 0x283e10f1, 0xf60795b9, 0xb10e7014)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#53"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Rodas5{0, true, DefaultLinSolve, DataType}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Rational{Int64}, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, 
Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~\.julia\packages\OrdinaryDiffEq\la32G\src\solve.jl:293
  [8] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#265"), Symbol("##arg#266"), Symbol("##arg#267"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0415a075, 0xa1d706f6, 0x283e10f1, 0xf60795b9, 0xb10e7014)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#53"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Rodas5{0, true, DefaultLinSolve, DataType}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq ~\.julia\packages\OrdinaryDiffEq\la32G\src\solve.jl:66
  [9] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#265"), Symbol("##arg#266"), Symbol("##arg#267"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0415a075, 0xa1d706f6, 0x283e10f1, 0xf60795b9, 0xb10e7014)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#53"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Rodas5{0, true, DefaultLinSolve, DataType}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~\.julia\packages\OrdinaryDiffEq\la32G\src\solve.jl:4
 [10] __solve
    @ ~\.julia\packages\OrdinaryDiffEq\la32G\src\solve.jl:4 [inlined]
 [11] #solve_call#56
    @ ~\.julia\packages\DiffEqBase\qntkj\src\solve.jl:61 [inlined]
 [12] solve_call
    @ ~\.julia\packages\DiffEqBase\qntkj\src\solve.jl:48 [inlined]
 [13] #solve_up#58
    @ ~\.julia\packages\DiffEqBase\qntkj\src\solve.jl:82 [inlined]
 [14] solve_up
    @ ~\.julia\packages\DiffEqBase\qntkj\src\solve.jl:75 [inlined]
 [15] #solve#57
    @ ~\.julia\packages\DiffEqBase\qntkj\src\solve.jl:70 [inlined]
 [16] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#265"), Symbol("##arg#266"), Symbol("##arg#267"), :t), ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", ModelingToolkit.StructuralTransformations.var"#_RGF_ModTag", (0x0415a075, 0xa1d706f6, 0x283e10f1, 0xf60795b9, 0xb10e7014)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Bool, Int64}, Nothing, Nothing, Nothing, Vector{Symbol}, Nothing, ModelingToolkit.StructuralTransformations.var"#generated_observed#53"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Rodas5{0, true, DefaultLinSolve, DataType})
    @ DiffEqBase ~\.julia\packages\DiffEqBase\qntkj\src\solve.jl:68
 [17] top-level scope

Is this a bug or am I doing something wrong?

Sciemon · April 20, 2021, 1:12pm

Just saying I also had this error message using a similar electrical model. I thought it was because of using complex numbers, but it isn’t apparently. I tested many other (stiff) solvers but most didn’t work.

ChrisRackauckas · April 20, 2021, 4:56pm

]st and versioninfo()?

ohmsweetohm1 · April 21, 2021, 6:53am

pkg> st
      Status `.\Project.toml`
  [0c46a032] DifferentialEquations v6.16.0
  [961ee093] ModelingToolkit v5.15.0
  [91a5bcdd] Plots v1.12.0

julia> versioninfo()
Julia Version 1.6.0
Commit f9720dc2eb (2021-03-24 12:55 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 4

ohmsweetohm1 · April 25, 2021, 8:14am

Works with

(modeling_toolkit) pkg> st
      Status `.\Project.toml`
  [0c46a032] DifferentialEquations v6.16.0
  [961ee093] ModelingToolkit v5.16.0
  [91a5bcdd] Plots v1.12.0