# ModelingToolkit connector error with Arrays

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565>\
**Category:** General Usage\
**Tags:** modelingtoolkit\
**Created:** [December 7, 2024, 12:27pm UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565 "2024-12-07T12:27:18Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![fracpl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fracpl/32/214015_2.png) [@fracpl](https://discourse.julialang.org/u/fracpl)\
**Post date:** [December 7, 2024, 12:27pm UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/1 "2024-12-07T12:27:18Z")

</div>

Hi, I am trying to implement from scratch some 3D multibody dynamics with ModelingToolkit.jl. When I define my connector (Frame) using arrays:

```julia
using ModelingToolkit
using LinearAlgebra
using ModelingToolkit: t_nounits as t, D_nounits as D

@connector Frame begin
    r(t)[1:3]
    q(t)[1:4]
    f(t)[1:3], [connect=Flow]
    τ(t)[1:3], [connect=Flow]
end

@named f = Frame()

```

I get the following error:

```julia
ERROR: MethodError: no method matching CartesianIndices(::Tuple{SubArray{Any, 1, Vector{Any}, Tuple{UnitRange{Int64}}, true}})
The type `CartesianIndices` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  CartesianIndices(::Tuple{})
   @ Base multidimensional.jl:271
  CartesianIndices(::RecursiveArrayTools.AbstractVectorOfArray)
   @ RecursiveArrayTools C:\Users\francesco capolupo\.julia\packages\RecursiveArrayTools\DX3HL\src\vector_of_array.jl:684       
  CartesianIndices(::AbstractArray)
   @ Base multidimensional.jl:281
  ...

Stacktrace:
 [1] eachindex(A::Symbolics.Arr{Num, 1})
   @ Symbolics C:\Users\francesco capolupo\.julia\packages\Symbolics\YbNrd\src\arrays.jl:589
 [2] iterate(A::Symbolics.Arr{Num, 1})
   @ Base .\abstractarray.jl:1207
 [3]
   @ Main .\none:0
 [4] __Frame__
   @ .\none:0 [inlined]
 [5] #_#441
   @ C:\Users\francesco capolupo\.julia\packages\ModelingToolkit\0O7FS\src\systems\model_parsing.jl:25 [inlined]
 [6] top-level scope
   @ C:\Users\francesco capolupo\.julia\packages\ModelingToolkit\0O7FS\src\systems\abstractsystem.jl:2113
Some type information was truncated. Use `show(err)` to see complete types.

```

Everything works fine if, instead of using arrays, I use scalar components:

```julia
@connector Frame2 begin
    x(t)
    y(t)
    z(t)
    qs(t)
    qx(t)
    qy(t)
    qz(t)
    fx(t), [connect=Flow]
    fy(t), [connect=Flow]
    fz(t), [connect=Flow]
    τx(t), [connect=Flow]
    τy(t), [connect=Flow]
    τz(t), [connect=Flow]
end
@named f2 = Frame2()

┌ Warning: f2 contains 6 flow variables, yet 7 regular (non-flow, non-stream, non-input, non-output) variables. This could lead to imbalanced model that are difficult to debug. Consider marking some of the regular variables as input/output variables.      
└ @ ModelingToolkit C:\Users\francesco capolupo\.julia\packages\ModelingToolkit\0O7FS\src\systems\connectors.jl:54
Model f2:
Equations (13):
  13 connecting: see equations(expand_connections(f2))
Unknowns (13): see unknowns(f2)
  x(t)
  y(t)
  z(t)
  qs(t)
  qx(t)
  qy(t)
  ⋮

```

Am I doing something wrong with the arrays? Thanks!

---

<div class="post-metadata">

**Author:** ![Walter\_Dal\_Maz\_Silva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/walter_dal_maz_silva/32/202490_2.png) [@Walter\_Dal\_Maz\_Silva](https://discourse.julialang.org/u/Walter_Dal_Maz_Silva)\
**Post date:** [December 19, 2024, 7:28am UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/2 "2024-12-19T07:28:47Z")

</div>

I am facing the same problem, which I think is related to the way `@variables` is parsing the brackets in the expression; to illustrate this, let me start by defining a simple `@connector` with symbolic metadata:

```julia
@connector FlowPort begin
    @variables begin
        ṁ(t) , [unit = u"kg/s", connect = Flow]
        Y(t) , [unit = u"kg/kg"]
    end
end

@named f = FlowPort()

```

The above works simply fine. Now let’s add a `@structural_parameters` for array sizing (in my case I intend to work with generic chemical kinetics, thus parametrizing is important). It fails exactly the way you are describing.

```julia
@connector FlowPort begin
    @structural_parameters begin
        Nc = 2
    end	
    @variables begin
        ṁ(t) , [unit = u"kg/s", connect = Flow]
        Y(t)[1:Nc] , [unit = u"kg/kg"]
    end
end

@named f = FlowPort(Nc = 3)

```

Yesterday I had the same error in another context (building a simpler prototype without reusable components) and someway (I forgot why, still trying to figure out) I realized the trouble was related to the square brackets, so I tried the following which works but does not meet my requirements (having units, description, …).

```julia
@connector FlowPort begin
    @structural_parameters begin
        Nc = 2
    end
    @variables begin
        ṁ(t) , [unit = u"kg/s", connect = Flow]
        Y(t)[1:Nc]
    end
end
	
@named f = FlowPort(Nc=3)	

```

Next I tried a few things regarding array initialization and in _some cases_ it works; I still get the warning but (I think) that is normal because a connector is not expected to be instantiated as a standalone object.

```julia
@connector FlowPort begin
    @structural_parameters begin
        Nc
    end
    @parameters begin
        # FAILS: Yc[1:Nc] = zeros(Nc), [description = "test", unit = u"kg/kg"]
        # WORKS: Yc[1:Nc] = zeros(Nc), [description = "test"]
        # WORKS: Yc[1:Nc] = zeros(Nc)
    end
    @variables begin
        ṁ(t) , [unit = u"kg/s", connect = Flow]
        Y(t)[1:Nc] = 0u"kg/kg", [unit = u"kg/kg"]
    end
end

@named f = FlowPort(Nc=3)

```

What I find funny here is that [DynamicQuantity Literals](https://docs.sciml.ai/ModelingToolkit/stable/basics/Validation/#DynamicQuantities-Literals) are supposed not to work in this context, but well, I just entered a literal there. But that is another bug unrelated to our root problem and if you are not using units that is normal syntax.

I will keep on developing to see if it works when integrated and I will update a comment here (wrote at this stage because I will forget all the steps to get to this point).

PS: please notice that I used different combinations of standard and constructor values of `Nc` and a case without standard value.

UPDATE 1: The warnings are not related to incomplete model. They persist when assembling the system of equations. Investigation ongoing.

UPDATE 2: Using `0u"kg/kg"` or a parameter `Yc` to initialize the variables does not produce the same outcomes: if using the parameter (even though it has a default value) a free parameter is identified by the system.

---

<div class="post-metadata">

**Author:** ![Walter\_Dal\_Maz\_Silva](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/walter_dal_maz_silva/32/202490_2.png) [@Walter\_Dal\_Maz\_Silva](https://discourse.julialang.org/u/Walter_Dal_Maz_Silva)\
**Post date:** [December 20, 2024, 8:53am UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/3 "2024-12-20T08:53:50Z")

</div>

After trying for the entire day, I realized the `@connector` API is not fully ready to support array variables yet. While waiting for its full development, if any, I came out with a solution that would allow me later to easily convert my model into a full-MTK one.

(1) Created a standardized parameter/variable generator function

```julia
"Shared state of any mixture element."
function mixture_state(kind, Nc)
	state = if kind == :variables
		@variables begin
			ṁ(t) , [unit = u"kg/s"]
			T(t) , [unit = u"K"]
			p(t) , [unit = u"Pa"]
			Y(t)[1:Nc] , [unit = u"kg/kg"]
		end
	else
		@parameters begin
			ṁ , [unit = u"kg/s"]
			T , [unit = u"K"]
			p , [unit = u"Pa"]
			Y[1:Nc] , [unit = u"kg/kg", tunable = false]
		end
	end

	return state
end

```

(2) Created a replacement to `connect` using the arguments generated by the above function:

```julia
"Replace `connect` for elements using `MixtureFlowPorts`."
function plug(d, u)
	# Notice that the sign convention here is intended to produce models
	# that remain compatible with a @connector (0 ~ ustream.ṁ + dstream.ṁ)
	# because there is nothing handling the I/O directionality here as in
	# fields with metadata `connect = Flow`.
	return [
		d.ṁ ~ -u.ṁ
		d.T ~ u.T
		d.p ~ u.p
		scalarize(d.Y ~ u.Y)...
	]
end

```

(3) Used old-styled functions to parametrically build what previously was a connector:

```julia
"A multicomponent mass flow port with state variables."
function MixtureFlowPort(; name, Nc)
    sts = mixture_state(:variables, Nc)
    return ODESystem(Equation[], t, sts, []; name)
end

```

The difficulty that you need to create a `plug` function for each family of variables you may need to connect (if several streams exist in the same model). I am also thinking about making `mixture_state` a structure and then create typed versions of `plug`.

From that you can get an idea of what the specific models (blocks) look like. You can check the full code [here](https://github.com/wallytutor/WallyToolbox.jl/blob/main/docs/src/Pluto/acausal-modeling.jl).

PS: Please notice that I use my own toolbox in extremely specific ways and if you want to run the code it is simpler to clone my repository and use script `wjulia` at its root (tweaking to match your system’s paths) and launch Pluto from the opened session. That is required so that Julia find all my library files. I created this monster upon frustration for the lack of the equivalent of `pip install -e .` in Julia.

---

<div class="post-metadata">

**Author:** ![fracpl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fracpl/32/214015_2.png) [@fracpl](https://discourse.julialang.org/u/fracpl)\
**Post date:** [December 24, 2024, 8:28am UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/4 "2024-12-24T08:28:39Z")

</div>

Thanks for your help, I made a few additional tests  
on my side too and got to the same conclusion: ModelingToolkit is not fully ready to handle arrays yet. I hope they’ll fix this soon…

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 24, 2024, 8:40am UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/5 "2024-12-24T08:40:54Z")

</div>

Are you aware of

> **[GitHub - JuliaComputing/Multibody.jl: Model and simulate multibody systems in Julia](https://github.com/JuliaComputing/Multibody.jl)**
>
> Model and simulate multibody systems in Julia

?

---

<div class="post-metadata">

**Author:** ![fracpl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fracpl/32/214015_2.png) [@fracpl](https://discourse.julialang.org/u/fracpl)\
**Post date:** [January 6, 2025, 4:27pm UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/6 "2025-01-06T16:27:28Z")

</div>

Yes, but I understand that Multibody.jl requires a (quite expensive!) JuliaSim licence. I was looking for an open source and free of charge multibody library.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [January 6, 2025, 4:32pm UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/7 "2025-01-06T16:32:25Z")

</div>

> [@fracpl](#):
>
> but I understand that [Multibody.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Multibody) requires a (quite expensive!) JuliaSim licence

Yes, a license is required for professional, non-academic use.

I wish you the best of luck if you want to implement this yourself, in this case, you can view the existence of Multibody.jl as simply a proof that it is possible 🙂

---

<div class="post-metadata">

**Author:** ![RobbesU](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robbesu/32/212459_2.png) [@RobbesU](https://discourse.julialang.org/u/RobbesU)\
**Post date:** [January 31, 2025, 12:33pm UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/8 "2025-01-31T12:33:36Z")

</div>

@fracpl, depending on what you’re trying to do, there is also [RigidBodyDynamics](https://juliarobotics.org/RigidBodyDynamics.jl/stable/), part of the [JuliaRobotics.org](https://juliarobotics.org/) ecosystem.  
Very much robotics oriented, causal and limited to rigid bodies (as the name suggests) although can define springs, dampers etc also. Was dormant for a few years but appears that @ferrolho has picked things up again.

---

<div class="post-metadata">

**Author:** ![ferrolho](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ferrolho/32/213665_2.png) [@ferrolho](https://discourse.julialang.org/u/ferrolho)\
**Post date:** [January 31, 2025, 1:02pm UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/9 "2025-01-31T13:02:32Z")

</div>

I use [RigidBodyDynamics.jl](https://github.com/JuliaRobotics/RigidBodyDynamics.jl) on a regular basis for my robotics projects. Although the package has been dormant for a while, it is still very much in working condition.

---

<div class="post-metadata">

**Author:** ![RobbesU](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/robbesu/32/212459_2.png) [@RobbesU](https://discourse.julialang.org/u/RobbesU)\
**Post date:** [February 3, 2025, 11:22am UTC](https://discourse.julialang.org/t/modelingtoolkit-connector-error-with-arrays/123565/11 "2025-02-03T11:22:08Z")

</div>

This works (`Julia 1.11.3`, `MTK 9.61.0` and `9.62.0`)

This little test system comprises a mass-flange-inputforce system with gravity, a forced sinusoidal input function via the flange, and a constant force direct on the mass.  
(edit: added the packages and constants used)

```julia
using ModelingToolkit
using DifferentialEquations
using LaTeXStrings
using Plots
import LinearAlgebra as LA
import ModelingToolkitStandardLibrary.Mechanical.Translational as T
import ModelingToolkitStandardLibrary.Blocks as B

using ModelingToolkit: t_nounits as t, D_nounits as D

# a few constants for easy test switching
const G_EARTH = [0.0, 0.0, -9.81] # m/s²
const G_MOON = [0.0, 0.0, -1.625]
const G_MARS = [0.0, 0.0, -3.72076]
const Earth="Earth"
const Moon="Moon"
const Mars="Mars"
const gravities=Dict(Earth => G_EARTH, Moon=>G_MOON, Mars=>G_MARS)

@connector function Flange3D(;name)

	vars = @variables begin
		f(t)[1:3], [connect=Flow, description = "Cut force resolved in connector"]
		v(t)[1:3], [description = "Cut velocities in connector"]
	end

	pars = []
	
	v,f = collect.((v,f))
	vars = [v;f]
	eqs = [1~1] # This avoids "BoundsError: attempt to access 0-element Vector{Vector{Any}} at index [0]"
	eq = eqs==[] ? Equation[] : eqs
	
	ODESystem(eq, t, vars, pars; name)
end

@component function Mass3D(; name, m, g=G_EARTH) 

	pars = @parameters begin # this allows the params to be set in mtkbuild call OR as system parameters
		m=m
		g[1:3]=g  
		f0[1:3]=[0.,0.,0.]
	end
	
    vars = @variables begin
        s(t)[1:3] # =[0.,0.,0.], [description = "Positions in [m] along X, Y and Z axes of the mass"]
        v(t)[1:3] # =[0.,0.,0.], [description = "Velocities in [m/s] along X, Y and Z axes of the mass"]
		a(t)[1:3], [description = "Acceleration in [m/s²] along the X,Y,Z axes of the mass"]
		f(t)[1:3], [description = "Forces in [N] along X,Y,Z axes of the mass"]
    end

    systems = @named begin
		flange = Flange3D() # connects v & f
	end

    eqs = [
        collect(D.(s) .~ v)...
        collect(D.(v) .~ a)...
        collect(v .~ flange.v)...
        collect(f .~ flange.f)...
        collect(a .~ (f .+f0)./m .+ g)...
	]

	initialization_eqs=[]
	
	return ODESystem(eqs, t, vars, pars; systems, name, initialization_eqs)
end

@component function ForceFunction3D(; name, a0=[0.,0.,0.])
	pars = @parameters a0[1:3]=a0
	vars = @variables f(t) 
    systems = @named begin
        output = B.RealOutputArray(; nout=3)
    end
	@constants π=pi # this makes for 'pretty' equations in display
	
	f = a0 .* sin(π*t) # this gets passed through!
    eqs = [
        collect(output.u .~ -f)...
	]
	
	return ODESystem(eqs, t, vars, pars; systems, name)
end

@component function Force3D(; name)
	pars = []
	vars = []
    systems = @named begin
        input = B.RealInputArray(; nin=3, guess=[0.,0.,0.])
		flange = Flange3D()
    end
    eqs = [
		collect(flange.f .~ input.u)...
	]
	
	return ODESystem(eqs, t, vars, pars; name, systems) 
end

@component function ForcedMass3D(;name, m, g=G_EARTH, a0=[0.,0.,0.] )
	pars = []
	vars = []
	systems = @named begin
		force = Force3D()
		force_func = ForceFunction3D(a0=a0)
		mass = Mass3D(m=m,g=g)
	end
	eqs = [
		ModelingToolkit.connect(force_func.output, force.input )
		ModelingToolkit.connect(force.flange, mass.flange)
	]

	return ODESystem(eqs, t, vars, pars; systems, name)
end

@mtkbuild sys1 = ForcedMass3D(m=1.0); # = @named + structural_simplity()

parameters(sys1) # mass.f0, mass.g, mass.m, force_func.a0

```

Simulate it with the variable parameters available without recreating the ODESystem:

```julia
	dt = 0.02
	duration = 5.0
	tspan, ts = time_span(duration, dt) # duration, dt
	
	# MTK unknows are pos s & vel v: unknowns(simple_sys) 
	s00 = [0.,0.,0.]
	v00 = [0.,0.,0.]
	s0 = [10., 20., 10.]
	v0 = [5., 2., 3.]
	a0 = [0.,10.,+20] # amplitude of input force function
	f0 = [1.,0.,10.] # constant forces on mass over and above gravity
	z0 = [0.,0.,0.]
	
	p = Earth
	g = gravities[p]
	m = 1.
	
	u0 = []  
	p0 = Dict( :mass₊m=>m, :mass₊g=>g, :mass₊f0=>f0, :force_func₊a0=>a0) # ... this works :-)
	
	forcedmass_init_eqs = [
		collect(sys1.mass.s .~ s0 )...
		collect(sys1.mass.v .~ v0 )...
	]
	
	prob = ODEProblem(sys1, u0, tspan, p0; initialization_eqs=forcedmass_init_eqs)

	sol = solve(prob, Rodas5(), dt=dt, saveat=ts)
	plot(sol, title="$p? trajectory")

```