Array-valued parameters are not interpreted correctly when creating ODESystem

[johannesnauta]

Describe the bug 🐞

When defining array-valued parameters using the @parameters macro, I noticed odd behavior that the fact that the parameter is array-valued is lost when creating an ODESystem. More interestingly, it worked just fine until today. While I ran Pkg.update(), I actually pinned the version of ModelingToolkit, so I am unsure what changed (recently) in regards to this, as I am using v9.23.0, which is also the stable documentation.

Minimal Reproducible Example 👇

julia> using ModelingToolkit
julia> S = 2
julia> @variables t x(t)[1:S]
2-element Vector{Any}:
 t
  (x(t))[1:2]

julia> @parameters a[1:S]
1-element Vector{Symbolics.Arr{Num, 1}}:
 a[1:2]

julia> D = Differential(t)
Differential(t)

julia> eqns = [D(x[i]) ~ a[i] * x[i] for i in 1:S]
2-element Vector{Equation}:
 Differential(t)((x(t))[1]) ~ a[1]*(x(t))[1]
 Differential(t)((x(t))[2]) ~ a[2]*(x(t))[2]

julia> @named sys = ModelingToolkit.ODESystem(eqns, t)
Model sys with 2 equations
Unknowns (2):
  (x(t))[1]
  (x(t))[2]
Parameters (1):
  a

Expected behavior
I expect this:

Parameters (2):
  a[1]
  a[2]

Interestingly, it seems to work when the parameters are matrix-valued, i.e.,

julia> @parameters b[1:S,1:S]
1-element Vector{Symbolics.Arr{Num, 2}}:
 b[1:2,1:2]

julia> eqns = [D(x[i]) ~ b[i,i] * x[i] for i in 1:S]
2-element Vector{Equation}:
 Differential(t)((x(t))[1]) ~ b[1, 1]*(x(t))[1]
 Differential(t)((x(t))[2]) ~ b[2, 2]*(x(t))[2]

julia> @named sys = ModelingToolkit.ODESystem(eqns, t)
Model sys with 2 equations
Unknowns (2):
  (x(t))[1]
  (x(t))[2]
Parameters (2):
  b[1, 1]
  b[2, 2]

Error & Stacktrace ⚠️
No relevant error, simply missing parameters.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [6e4b80f9] BenchmarkTools v1.5.0
  [336ed68f] CSV v0.10.14
  [13f3f980] CairoMakie v0.12.5
  [479239e8] Catalyst v14.0.1
  [5ae59095] Colors v0.12.11
  [861a8166] Combinatorics v1.0.2
  [a93c6f00] DataFrames v1.6.1
  [0c46a032] DifferentialEquations v7.13.0
  [31c24e10] Distributions v0.25.109
  [ffbed154] DocStringExtensions v0.9.3
  [b7d42ee7] Einsum v0.4.1
  [d872a56f] ElectronDisplay v1.0.1
  [d4d017d3] ExponentialUtilities v1.26.1
  [68837c9b] FHist v0.11.5
  [86223c79] Graphs v1.11.2
  [033835bb] JLD2 v0.4.50
  [ccbc3e58] JumpProcesses v9.11.1
  [b964fa9f] LaTeXStrings v1.3.1
  [23fbe1c1] Latexify v0.16.4
⌃ [961ee093] ModelingToolkit v9.23.0 ⚲
  [2774e3e8] NLsolve v4.5.1
  [77ba4419] NaNMath v1.0.2
  [67456a42] OhMyThreads v0.5.3
  [d96e819e] Parameters v0.12.3
  [91a5bcdd] Plots v1.40.5
  [1fd47b50] QuadGK v2.9.4
  [0bca4576] SciMLBase v2.42.0
⌃ [a9a3c162] SparseArrayKit v0.2.1
  [276daf66] SpecialFunctions v2.4.0
  [0c5d862f] Symbolics v5.33.0
  [bd369af6] Tables v1.12.0
  [6aa20fa7] TensorOperations v5.0.0
  [3a884ed6] UnPack v1.0.2
  [8ba89e20] Distributed
  [37e2e46d] LinearAlgebra
  [9a3f8284] Random
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0

• Output of versioninfo()

Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × 12th Gen Intel(R) Core(TM) i5-1235U
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, alderlake)
Threads: 8 default, 0 interactive, 4 GC (on 12 virtual cores)
Environment:
  JULIA_NUM_THREADS = 4

[johannesnauta]

The problem seems to persist even when going back to month-old version (e.g., ModelingToolkit@9.19.0), even though I had successfully ran the above example a week ago. It seems that only reverting back to Julia v1.9.4 results in the example above to run correctly, but in this case I cannot use some other, newer, ModelingToolkit functionality that my code relies on (specifically, the setu and setp functions).

Edit: the above code does seem to work on ModelingToolkit version 8.76.0

[SebastianM-C]

The bug might be with the matrix parameters, the vector parameters work correctly. The fact that you get

Parameters (1):
  a

is not a bug. The a parameter is 1 vector parameter, as MTK doesn't scalarize vectors anymore on v9 (it did that on v8).

[johannesnauta]

So how should one then define array-valued parameters? When I use @parameters a[1:2] and use the ODESystem to create and ODEProblem I need to be able to supply a[1] and a[2] (e.g., by using params = [a[1] => 1.0, a[2] => 1.0]), but currently (in v9+) I get an error that a is Missing from the variable map.

[ChrisRackauckas]

params = [a => [1.0,1.0]]

[johannesnauta]

Ok that makes a lot of sense. Must have missed this in the docs at some point, and I guess I was using v8 for some time. Thanks for the swift replies.