Matrix inversion in ModellingToolkit

Specific Domains Modelling & Simulations
1

Hi all, I am trying to perform matrix inversion on a 2-dimensional variable in ModellingToolkit, but I get a MethodError that there is no matching method. Is there some implementation of matrix inversion that I could use in ModellingToolkit?

Here is the MWE:

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

myRange = 1:4
rangeLength = length(myRange)

@parameters A[myRange,myRange] = rand(rangeLength, rangeLength)
@variables X(t)[myRange,myRange] = rand(rangeLength,rangeLength)
@variables Y(t)[myRange,myRange] = rand(rangeLength,rangeLength)

eqs = [
[D(X)[i,j] ~ X[i,j] * A[i,j] for i=myRange for j=myRange]...
[Y[i,j] ~ inv(X)[i,j] for i=myRange for j=myRange]...
]

@named sys = ODESystem(eqs, t)
simple = structural_simplify(sys)
prob = ODEProblem(simple, [], (0,50))
sol = solve(prob, Euler(), dt=1, dtmax = 1)

Below is the stacktrace:

ERROR: MethodError: no method matching inv(::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true})

Closest candidates are:
  inv(::Crayons.ANSIStyle)
   @ Crayons C:\Users\robva\.julia\packages\Crayons\u3AH8\src\crayon.jl:76
  inv(::Complex{Num})
   @ Symbolics C:\Users\robva\.julia\packages\Symbolics\Eas9m\src\num.jl:54
  inv(::ComplexF64)
   @ Base complex.jl:476
  ...

Stacktrace:
  [1] macro expansion
    @ C:\Users\robva\.julia\packages\SymbolicUtils\c0xQb\src\code.jl:375 [inlined]
  [2] macro expansion
    @ C:\Users\robva\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:163 [inlined]
  [3] macro expansion
    @ .\none:0 [inlined]
  [4] generated_callfunc
    @ .\none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…}, ::Vector{…}, ::Int64)
    @ RuntimeGeneratedFunctions C:\Users\robva\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:150
  [6] (::ModelingToolkit.var"#f#711"{…})(du::Vector{…}, u::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Int64)
    @ ModelingToolkit C:\Users\robva\.julia\packages\ModelingToolkit\PtnSY\src\systems\diffeqs\abstractodesystem.jl:344
  [7] (::SciMLBase.Void{ModelingToolkit.var"#f#711"{…}})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\robva\.julia\packages\SciMLBase\hSv8d\src\utils.jl:482
  [8] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Vector{…}, arg2::Vector{…}, arg3::ModelingToolkit.MTKParameters{…}, arg4::Int64)
    @ FunctionWrappers C:\Users\robva\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
  [9] macro expansion
    @ C:\Users\robva\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined]
 [10] do_ccall
    @ C:\Users\robva\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined]
 [11] FunctionWrapper
    @ C:\Users\robva\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined]
 [12] _call
    @ C:\Users\robva\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:12 [inlined]
 [13] FunctionWrappersWrapper
    @ C:\Users\robva\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:10 [inlined]
 [14] ODEFunction
    @ C:\Users\robva\.julia\packages\SciMLBase\hSv8d\src\scimlfunctions.jl:2296 [inlined]
 [15] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.EulerCache{…})
    @ OrdinaryDiffEq C:\Users\robva\.julia\packages\OrdinaryDiffEq\0tf1M\src\perform_step\fixed_timestep_perform_step.jl:76
 [16] __init(prob::ODEProblem{…}, alg::Euler, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; 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::Int64, dtmin::Int64, dtmax::Int64, force_dtmin::Bool, adaptive::Bool, gamma::Int64, abstol::Nothing, reltol::Nothing, qmin::Int64, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Int64, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(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), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEq C:\Users\robva\.julia\packages\OrdinaryDiffEq\0tf1M\src\solve.jl:518
 [17] __init (repeats 5 times)
    @ C:\Users\robva\.julia\packages\OrdinaryDiffEq\0tf1M\src\solve.jl:11 [inlined]
 [18] #__solve#766
    @ C:\Users\robva\.julia\packages\OrdinaryDiffEq\0tf1M\src\solve.jl:6 [inlined]
 [19] __solve
    @ C:\Users\robva\.julia\packages\OrdinaryDiffEq\0tf1M\src\solve.jl:1 [inlined]
 [20] solve_call(_prob::ODEProblem{…}, args::Euler; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\robva\.julia\packages\DiffEqBase\NaUtB\src\solve.jl:612
 [21] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::ModelingToolkit.MTKParameters{…}, args::Euler; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\robva\.julia\packages\DiffEqBase\NaUtB\src\solve.jl:1080
 [22] solve_up
    @ C:\Users\robva\.julia\packages\DiffEqBase\NaUtB\src\solve.jl:1066 [inlined]
 [23] solve(prob::ODEProblem{…}, args::Euler; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\robva\.julia\packages\DiffEqBase\NaUtB\src\solve.jl:1003
 [24] top-level scope
    @ c:\Users\robva\ToBe\ToBe-Model\sandbox\test_matrix_inverse.jl:24

Thanks in advance!

2

Do you need the actual inverse or can you write it as a linear solve? Normally it can be written as a linear solve () and that would fix expression growth some. Use solve_for if you can

3

I need the actual inverse of the matrix. I implemented the inverse using solve_for (see code below), and that implementation works. However, I think it won’t scale well as the equations become very long. It works relatively fast on a 4x4 matrix, but when I tried a 10x10 matrix it took a very long time.

Is there a more efficient way of implementing the inverse of a matrix in ModelingToolkit?

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

myRange = 1:4
rangeLength = length(myRange)

@parameters A[myRange,myRange] = rand(rangeLength, rangeLength)
@parameters Z[myRange,myRange] = Matrix{Float64}(I, rangeLength, rangeLength)
@variables X(t)[myRange,myRange] = rand(rangeLength,rangeLength)
@variables Y(t)[myRange,myRange] = rand(rangeLength,rangeLength)

invertEquation = [Z[i,j] ~ sum([X[i,k] * Y[k,j] for k=myRange]) for i=myRange for j=myRange]
inverseX = ModelingToolkit.solve_for(invertEquation, [Y[i,j] for i=myRange for j=myRange])
inverseX = reshape(inverseX, rangeLength, rangeLength)

eqs = [
    [D(X)[i,j] ~ X[i,j] * A[i,j] for i=myRange for j=myRange]...
    [Y[i,j] ~ inverseX[i,j] for i=myRange for j=myRange]...
]

@named sys = ODESystem(eqs, t)

simple = structural_simplify(sys)
prob = ODEProblem(simple, [], (0,50))
sol = solve(prob, Euler(), dt=1, dtmax = 1)
4

No that expression right there does not need the matrix inverse. Keep it lazy if you want to scale. There is no matrix inverse expression without expression growth