Write acausal components for scalar, vector and rotations integrator

[traversaro]

While working on https://github.com/traversaro/MechanicsMTKExperiments/issues/2, I had several issue, that I tought it was easier to debug by building some simpler components: 1) A scalar integrator 2) A vector integrator 3) A rotation integrator (a component/system in which the input is an angular velocity $\omega$ and the output is $R$ such that $\dot{R}R^T = skew(\omega)$

[traversaro]

Scalar Integrator

This was quite easy:

using ModelingToolkit, Plots, DifferentialEquations

@variables t

@connector function ScalarConnector(;name)
    sts = @variables u(t)=0.0
    ODESystem(Equation[], t, sts, []; name=name)
end

function ScalarIntegrator(;name)
    @named input = ScalarConnector()
    @named output = ScalarConnector()
    D = Differential(t)
    eqs = [
           D(output.u) ~ input.u
          ]
    compose(ODESystem(eqs, t, [], []; name=name), input, output)
end

function ScalarConstant(;name)
    @named output = ScalarConnector()
    eqs = [
            output.u ~ 1.0
          ]
    compose(ODESystem(eqs, t, [], []; name=name), output)
end

@named constant = ScalarConstant()
@named integrator = ScalarIntegrator()

rc_eqs = [
          connect(constant.output, integrator.input)
         ]
@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
                          [constant, integrator])
sys = structural_simplify(rc_model)
u0 = [
        integrator.output.u => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
plot(sol)

[traversaro]

Vector Integrator

I tried to port the scalar integrato to a vector one, but I obtained some errors even in a much simpler problem:

using ModelingToolkit, Plots, DifferentialEquations

@variables t

@connector function VectorConnector(;name, N = 10)
    sts = @variables u[1:N](t)
    ODESystem(Equation[], t, sts, []; name=name)
end

function VectorConstant(;name, N = 10)
    @named output = VectorConnector(N=N)
    compose(ODESystem(Equation[], t, [], []; name=name), output)
end

N = 10
@named constant = VectorConstant(N = N)

This fails with error:

ERROR: LoadError: MethodError: no method matching hasmetadata(::Vector{Num}, ::Type{Symbolics.VariableDefaultValue})
Closest candidates are:
  hasmetadata(::SymbolicUtils.Symbolic, ::Any) at ~/.julia/packages/SymbolicUtils/FvmiV/src/types.jl:15
  hasmetadata(::Complex{Num}, ::Any) at ~/.julia/packages/Symbolics/HDE84/src/Symbolics.jl:147
  hasmetadata(::Num, ::Any) at ~/.julia/packages/Symbolics/HDE84/src/Symbolics.jl:147
Stacktrace:
  [1] hasdefault(v::Vector{Num})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/KItSa/src/utils.jl:190
  [2] collect_defaults!
    @ ~/.julia/packages/ModelingToolkit/KItSa/src/utils.jl:201 [inlined]
  [3] process_variables!(var_to_name::Dict{Any, Any}, defs::Dict{Any, Any}, vars::Vector{Vector{Num}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/KItSa/src/utils.jl:195
  [4] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Symbolics.Arr{Num, 1}}, ps::Vector{Any}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Any, Any}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, checks::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/KItSa/src/systems/diffeqs/odesystem.jl:146
  [5] macro expansion
    @ ~/jltest/vector_connector2.jl:7 [inlined]
  [6] (::var"#f#5"{Symbol, Int64})()
    @ Main ~/.julia/packages/ModelingToolkit/KItSa/src/systems/connectors.jl:17
  [7] VectorConnector(; name::Symbol, N::Int64)
    @ Main ~/.julia/packages/ModelingToolkit/KItSa/src/systems/connectors.jl:19
  [8] VectorConstant(; name::Symbol, N::Int64)
    @ Main ~/jltest/vector_connector2.jl:11
  [9] top-level scope
    @ ~/jltest/vector_connector2.jl:16
 [10] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [11] top-level scope
    @ REPL[2]:1
in expression starting at /home/traversaro/jltest/vector_connector2.jl:16

I am not an expert, but perhaps there is some problem with vectors in connectors (I did not find any example of that kind).

[traversaro]

I am not an expert, but perhaps there is some problem with vectors in connectors (I did not find any example of that kind).

I opened an issue for this in https://github.com/SciML/ModelingToolkit.jl/issues/1432 .

[traversaro]

By the way, a related library of ModelingToolkit components can be found in https://github.com/SciML/ModelingToolkitStandardLibrary.jl, in particular the "Block" module.

[traversaro]

By the way, a related library of ModelingToolkit components can be found in https://github.com/SciML/ModelingToolkitStandardLibrary.jl, in particular the "Block" component.

This was indeed interesting. Something learnt (from code like https://github.com/SciML/ModelingToolkitStandardLibrary.jl/blob/2d7a4e891abaa8182e8c31b252f2053e303bb718/src/Blocks/continuous.jl#L188):

• Modelica-style input and output annotations also exist for variables in ModelingTookit, even outside connectors
• You can define components that have input and outputs (even of vector variables!) without using connectors, but just by defining variables in the system. This was added in https://github.com/SciML/ModelingToolkit.jl/pull/1132 .

[traversaro]

I was able to write a vector integrator:

using ModelingToolkit, Plots, DifferentialEquations

import Symbolics.scalarize

@variables t

function VectorConstant(;N=1, y0=zeros(N), name)
    @variables y[1:N](t)=y0 [output=true]
    y = collect(y)
    eqs = y .~ ones(N)
    ODESystem(eqs, t, name=name)
end

function VectorIntegrator(;N=1, y0=zeros(N), name)
    @variables u[1:N](t)=0 [input=true] y[1:N](t)=y0 [output=true]
    u = collect(u) # https://github.com/JuliaSymbolics/Symbolics.jl/issues/379
    y = collect(y)
    D = Differential(t)
    eqs = D.(y) .~ (1:N).*u
    ODESystem(eqs, t, name=name)
end

N = 10
@named constant = VectorConstant(N = N)
@named integrator = VectorIntegrator(N = N)

rc_eqs = constant.y ~ integrator.u
@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
                          [constant, integrator])
sys = structural_simplify(rc_model)
u0 = [
        integrator.y => zeros(N,1)
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())
plot(sol)

Output:

However, dealing with vector values it is still a bit cumbersome and requires some understanding of ModelingToolkit and Symbolics internals.

[traversaro]

I am not an expert, but perhaps there is some problem with vectors in connectors (I did not find any example of that kind).

I opened an issue for this in SciML/ModelingToolkit.jl#1432 .

This was fixed.