Maximum of matrix vector product for "large" matrices

General Usage
1

I am using MTK and noticed that when I construct an algebraic variable defined as the max of the matrix-vector product of the state vector, the results are accurate for small matrices, but go inaccurate rapidly as the matrix size increases.
Specifically, I get accurate results up 100x10. For 500x10, it’s always inaccurate.
Copy-paste example (you can simply uncomment the definition of C):

using LinearAlgebra
using ModelingToolkit, Plots, DifferentialEquations
using Symbolics: scalarize

@variables t
D = Differential(t)

function F(;name, α0 = fill(0.0, 10), C = fill(1.0, 10, 10))

    sts = @variables T(t)=0.0 α(t)[1:length(α0)]=α0

    eqs = [ scalarize(D.(α) .~ 0.0)
            T ~ maximum(C*α)
            ]

    ODESystem(eqs, t, [α..., T], []; name=name)
end

α0 = rand(10)

#C = rand(100, 10)
C = rand(500, 10)

@named mySys = F(α0=α0, C=C)

sys = structural_simplify(mySys)

prob = ODEProblem(sys, [], (0.0, 1.0))
sol = solve(prob)

println(maximum(C*α0))
println(sol[mySys.T][end])

Any clues how to get around that using MTK? If i build the same system directly in DifferentialEquations with an algebraic state defined the same manner, I don’t run into the same problem.

2

@YingboMa @shashi could this be fixed by using the MuladdMacro expression changing on the generated code?

3

I see very similar values for the last two lines

julia> println(maximum(C*α0))
4.421761822317552

julia> println(sol[mySys.T][end])
4.421761822317553

I started a bigger system (1000x100) running but it hasn’t finished yet…

(MTK) pkg> st
Project MTK v0.1.0
Status `~/Desktop/Julia/MTK/Project.toml`
  [0c46a032] DifferentialEquations v7.6.0
  [961ee093] ModelingToolkit v8.30.0
  [91a5bcdd] Plots v1.35.5
  [731186ca] RecursiveArrayTools v2.32.1
  [0c5d862f] Symbolics v4.13.0

Edit
1000x100 finished!

julia> println(maximum(C*α0))
33.386913594597175

julia> println(sol[mySys.T][end])
33.38691359459718
4

Ah, good find! I checked and made sure the package versions were the same, but it didn’t make a difference in Julia 1.7.2. When I tried in Julia 1.8, I got nearly the same result (not equal but good enough)! Thank you for pointing that out!

5

Ok, now for part II now that the max issue is resolved by switching to Julia 1.8.2 :slight_smile: :

My original code looks a lot more like:

using LinearAlgebra
using ModelingToolkit, Plots, DifferentialEquations
using Symbolics: scalarize

@variables t
D = Differential(t)


function myComponent(;name, Nm = 10, α0=zeros(Nm))

    sts = @variables α(t)[1:Nm]=α0
    eqs = scalarize(D.(α) .~ 0.0)

    ODESystem(eqs, t, [α...], []; name)
end

function F(;name, α0s = [fill(0.0, 10) for i in 1:20], C = fill(1.0, 10, 200))

    components = [myComponent(name = Symbol(string(i)), α0 = α0s[i] ) for i in 1:length(α0s)]

    α = reduce(vcat, [component.α for component in components])

    @variables T(t)=0.0

    eqs = [  T ~ maximum(C*α)
            ]

    compose(ODESystem(eqs, t, [T], [],; name), components...)

end

α0 = fill(rand(10), 20)

#C = rand(100, 10)
C = rand(500, 200)

@named mySys = F(α0s=α0, C=C)

sys = structural_simplify(mySys)

prob = ODEProblem(sys, [], (0.0, 1.0))
@time solve(prob)

When I run that on Julia 1.7.2. I get:

5.666656 seconds (10.51 M allocations: 634.450 MiB, 1.81% gc time, 99.94% compilation time)
retcode: Success

When i try to run that on Julia 1.8.2, it takes > 5min (I’m still waiting for the answer and will edit).

Now I’m sure I’m doing something dumb with this concatenation (if this wasn’t symbolic, i’d at least pre-allocate). Any tips for a newbie at symbolic math/MTK?

Thanks for the help!

6

It’s been 15min now, so I give up.