ModelingToolkit connector error with Arrays

General Usage
1

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

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:

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:

@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!

1 Like
2

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:

@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.

@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, …).

@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.

@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 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.

1 Like
3

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

"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:

"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:

"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.

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.

4

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…

5

Are you aware of

?

6

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.

7

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 :slight_smile:

8

@fracpl, depending on what you’re trying to do, there is also RigidBodyDynamics, part of the 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.

1 Like
9

I use 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.

1 Like
10

(post deleted by author)

11

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)

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:

	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")