Matrix inversion in ODESystem in ModelingToolkit

Hi all, I’m trying to build an ODESystem in ModelingToolkit, with an equation that uses a matrix inversion. Basically, I initialize a random matrix A of size (industries,industries), and the matrix L of size (industries,industries) that represents the inverse. Both are of the same dimensions, but I get a DimensionMismatch error. Below you find my code, and my stacktrace.

using DifferentialEquations
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

industries = 4

function test_factory(;name)
    @parameters A[1:industries, 1:industries] = rand(industries, industries)
    @variables L(t)[1:industries, 1:industries]

    eq = [
        L ~ inv(A)
    ]

    return ODESystem(eq, t; name)
end

testFunction = true

if testFunction

    @named mysys = test_factory()

    eqsystem = compose(mysys)
    print(equations(eqsystem))
    eqsystem = complete(eqsystem)

    prob = ODEProblem(eqsystem, (0.0, 100.0))
    sol = solve(prob, Tsit5())

    Plots.display(plot(sol))
else

    @variables id[1:industries, 1:industries]
    @variables A[1:industries, 1:industries]
    @variables L[1:industries, 1:industries]

    id = I(industries)
    A = rand(industries, industries)
    L = inv(A)


    println(id)
    println(A)
    println(L)
end

Stacktrace:

Equation[L(t) ~ inv(A(t))]ERROR: DimensionMismatch: number of columns of each array must match (got (1, 4))
Stacktrace:
  [1] _typed_vcat(::Type{Equation}, A::Tuple{Vector{Equation}, Matrix{Equation}})
    @ Base .\abstractarray.jl:1702
  [2] typed_vcat
    @ .\abstractarray.jl:1716 [inlined]
  [3] vcat(::Vector{Equation}, ::Matrix{Equation})
    @ Base .\abstractarray.jl:1694
  [4] BottomRF
    @ .\reduce.jl:86 [inlined]
  [5] MappingRF
    @ .\reduce.jl:100 [inlined]
  [6] _foldl_impl(op::Base.MappingRF{…}, init::Vector{…}, itr::Vector{…})
    @ Base .\reduce.jl:58
  [7] foldl_impl
    @ .\reduce.jl:48 [inlined]
  [8] mapfoldl_impl
    @ .\reduce.jl:44 [inlined]
  [9] _mapreduce_dim
    @ .\reducedim.jl:362 [inlined]
 [10] mapreduce
    @ .\reducedim.jl:357 [inlined]
 [11] flatten_equations(eqs::Vector{Equation})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1526
 [12] TearingState(sys::ODESystem; quick_cancel::Bool, check::Bool)
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\systemstructure.jl:255
 [13] TearingState
    @ C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\systemstructure.jl:250 [inlined]
 [14] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Tuple{…}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), guesses::Dict{…}, t::Nothing, warn_initialize_determined::Bool, build_initializeprob::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:869
 [15] process_DEProblem
    @ C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:834 [inlined]
 [16] (ODEProblem{…})(sys::ODESystem, u0map::Tuple{…}, tspan::Nothing, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1085
 [17] ODEProblem
    @ C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1075 [inlined]
 [18] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Tuple{Float64, Float64})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1075
 [19] (ODEProblem{true})(sys::ODESystem, args::Tuple{Float64, Float64}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1062
 [20] (ODEProblem{true})(sys::ODESystem, args::Tuple{Float64, Float64})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1061
 [21] #ODEProblem#755
    @ C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1051 [inlined]
 [22] ODEProblem(sys::ODESystem, args::Tuple{Float64, Float64})
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\kByuD\src\systems\diffeqs\abstractodesystem.jl:1050
 [23] top-level scope
    @ c:\Users\robva\ToBe\Trials\Julia\Array tests.jl:28
Some type information was truncated. Use `show(err)` to see complete types.

What am I dong wrong? Or is this functionality not supported by ModelingToolkit?

Thanks!

This looks like a bug. Open an issue.