`dae_index_reduction` not working for symbolic arrays

[AayushSabharwal]

Reproducible example:

julia> @parameters t m=1.0

julia> @variables x[1:2](t)=[1.,-1.] v[1:2](t)=[0., 0.]  λ(t)

julia> eqs = [
 Differential(t)(x[1](t)) ~ v[1](t)
 Differential(t)(x[2](t)) ~ v[2](t)
 Differential(t)(v[1](t)) ~ 2λ(t)*x[1](t)*(m^-1)
 Differential(t)(v[2](t)) ~ (m^-1)*(2λ(t)*x[2](t) - (9.81m))
 x[1](t)^2 + x[2](t)^2 - 2 ~ 0
]

julia> @named sys = ODESystem(eqs, t)
Model sys with 5 equations
States (5):
  x[1](t) [defaults to 1.0]
  x[2](t) [defaults to -1.0]
  v[1](t) [defaults to 0.0]
  v[2](t) [defaults to 0.0]
  λ(t)
Parameters (1):
  m [defaults to 1.0]


julia> structural_simplify(sys)
ERROR: ArgumentError: A differentiated state's operation must be a `Sym`, so states like `D(u + u)` are disallowed. Got `x[1]`.
Stacktrace:
 [1] diff2term(O::Term{Real, Nothing})
   @ Symbolics ~/.julia/packages/Symbolics/fd3w9/src/utils.jl:114
 [2] pantelides_reassemble(sys::ODESystem, eqassoc::Vector{Int64}, assign::Vector{Int64})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/1IpUm/src/structural_transformation/pantelides.jl:28
 [3] dae_index_lowering(sys::ODESystem; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/1IpUm/src/structural_transformation/pantelides.jl:153
 [4] dae_index_lowering
   @ ~/.julia/packages/ModelingToolkit/1IpUm/src/structural_transformation/pantelides.jl:150 [inlined]
 [5] structural_simplify(sys::ODESystem; simplify::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/1IpUm/src/systems/abstractsystem.jl:811
 [6] structural_simplify(sys::ODESystem)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/1IpUm/src/systems/abstractsystem.jl:808
 [7] top-level scope
   @ REPL[61]:1

It appears as if dae_index_reduction specifically doesn't like symbolic arrays. This also occurs on the following (equivalent) system of equations:

3-element Vector{Equation}:
 Differential(t)(Differential(t)(x[1](t))) ~ 2λ(t)*x[1](t)*(m^-1)
 Differential(t)(Differential(t)(x[2](t))) ~ (m^-1)*(2λ(t)*x[2](t) - (9.81m))
 x[1](t)^2 + x[2](t)^2 - 2 ~ 0

My environment:

      Status `~/d/Julia/MTK_Circuits/Project.toml`
  [0c46a032] DifferentialEquations v6.19.0
  [961ee093] ModelingToolkit v6.4.8
  [91a5bcdd] Plots v1.21.1