Reproduce Acausal RC Example with Mass Spring Model

[traversaro]

I am trying to reproduce the following Modelica Model with ModellingToolkit.jl:

model OscillatorExample
  Modelica.Mechanics.Translational.Components.Fixed fixed annotation(
    Placement(visible = true, transformation(origin = {30, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Mechanics.Translational.Components.Mass mass(v(fixed = true, start = 1))  annotation(
    Placement(visible = true, transformation(origin = {-50, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  Modelica.Mechanics.Translational.Components.Spring spring(c = 100)  annotation(
    Placement(visible = true, transformation(origin = {-10, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
equation
  connect(spring.flange_b, fixed.flange) annotation(
    Line(points = {{0, 10}, {30, 10}}, color = {0, 127, 0}));
  connect(spring.flange_a, mass.flange_b) annotation(
    Line(points = {{-20, 10}, {-40, 10}}, color = {0, 127, 0}));

annotation(
    uses(Modelica(version = "4.0.0")));
end OscillatorExample;

That gives this simple results using OpenModelica 1.18 :

[traversaro]

Inspired by https://mtk.sciml.ai/stable/tutorials/acausal_components/ and the Mechanics.Translational compontent of the Modelica Standard Library, I tried this example:

using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Translational_Flange(;name)
    sts = @variables s(t)=0.0 f(t)=0.0
    ODESystem(Equation[], t, sts, []; name=name)
end

function Translational_Fixed(;name, s0 = 0.0)
    @named flange = Translational_Flange()
    ps = @parameters s0=s0
    eqs = [flange.s ~ s0]
    compose(ODESystem(eqs, t, [], ps; name=name), flange)
end

function Translational_TwoFlanges(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    compose(ODESystem([], t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialRigid
function Translational_TwoFlangesRigid(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        flange_a.s ~ flange_b.s
       ]
    compose(ODESystem(eqs, t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialCompliant
function Translational_TwoFlangesCompliant(;name)
    sts = @variables s_rel(t)=0.0 f(t)=0.0 
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        s_rel ~ flange_b.s - flange_a.s
        flange_b.f ~ f
        flange_a.f ~ -f
       ]
    compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Components.Mass
function Translational_Mass(;name, m = 1.0)
    @named twoflanges = Translational_TwoFlangesRigid()
    @unpack flange_a, flange_b = twoflanges
    sts = @variables s(t)=0.0 v(t) a(t)
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
           s ~ flange_a.s
           v ~ D(s)
           a ~ D(v)
           m*a ~ flange_a.f + flange_b.f
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), twoflanges)
end

# Inspired from Modelica.Mechanics.Translational.Components.Spring
function Translational_Spring(;name, c = 100.0, s_rel0 = 0.0)
    @named twoflanges = Translational_TwoFlangesCompliant()
    @unpack s_rel, f = twoflanges
    ps = @parameters c=c s_rel0=s_rel0
    eqs = [
           f ~ c * (s_rel - s_rel0)
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), twoflanges)
end

#=
m = 1.0
@named mass = Translational_Mass(m=m)

ms_eqs = [
         mass.flange_a.f ~ 0
         mass.flange_b.f ~ 0
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass])
sys = structural_simplify(ms_model)
# Mass starts with velocity = 1
u0 = [
    mass.v => 1.0
     ]
prob = ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob, Tsit5())
plot(sol, vars=[mass.v,mass.s,mass.a,mass.flange_b.f])
=#

m = 1.0
s0 = 0.0
c = 100.0
@named mass = Translational_Mass(m=m)
@named spring = Translational_Spring(c=c)
@named fixed = Translational_Fixed(s0=s0)

ms_eqs = [
          connect(mass.flange_b, spring.flange_a)
          connect(spring.flange_b, fixed.flange)
          mass.flange_a.f ~ 0
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass, spring, fixed])
sys = structural_simplify(ms_model)
# Mass starts with velocity = 1
u0 = [
    mass.v => 1.0
     ]
prob = ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob, Tsit5())
plot(sol, vars=[mass.v,mass.s])

The problem is that for some reason the position does not change at all with this model, i.e. the result is:

State of environment:

(@v1.7) pkg> st
      Status `~/.julia/environments/v1.7/Project.toml`
  [0c46a032] DifferentialEquations v7.0.0
  [961ee093] ModelingToolkit v8.0.0
  [91a5bcdd] Plots v1.25.4

[traversaro]

I updated, and I found a MWE for my problem. In this model, I had the same problem:

using ModelingToolkit, Plots, DifferentialEquations

@variables t

function Translational_Mass(;name, m = 1.0)
    sts = @variables s(t) v(t) a(t)
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
           v ~ D(s)
           a ~ D(v)
          ]
    ODESystem(eqs, t, sts, ps; name=name)
end

m = 1.0
@named mass = Translational_Mass(m=m)

ms_eqs = [
         mass.a ~ 1
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass])
# Mass starts with velocity = 1
u0 = [
    mass.s => 0.0
    mass.v => 0.0
     ]

sys = structural_simplify(ms_model)
prob_complex =  ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob_complex, Tsit5())
plot(sol, vars=[mass.v,mass.s,mass.a])

while this one is working fine (where I just changed the order of the equations):

using ModelingToolkit, Plots, DifferentialEquations

@variables t

function Translational_Mass(;name, m = 1.0)
    sts = @variables s(t) v(t) a(t)
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
           D(s) ~ v
           D(v) ~ a
          ]
    ODESystem(eqs, t, sts, ps; name=name)
end

m = 1.0
@named mass = Translational_Mass(m=m)

ms_eqs = [
         mass.a ~ 1
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass])
# Mass starts with velocity = 1
u0 = [
    mass.s => 0.0
    mass.v => 0.0
     ]

sys = structural_simplify(ms_model)
prob_complex =  ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob_complex, Tsit5())
plot(sol, vars=[mass.v,mass.s,mass.a])

The underlying issue is https://github.com/SciML/ModelingToolkit.jl/issues/860 , and the PR that will fix the problem is https://github.com/SciML/ModelingToolkit.jl/pull/1164.

[traversaro]

Slight more complex model that works as well:

using ModelingToolkit, Plots, DifferentialEquations

@variables t

@connector function Translational_Flange(;name)
    sts = @variables s(t)=0.0 f(t)=0.0
    ODESystem(Equation[], t, sts, []; name=name)
end

function Translational_Fixed(;name, s0 = 0.0)
    @named flange = Translational_Flange()
    ps = @parameters s0=s0
    eqs = [flange.s ~ s0]
    compose(ODESystem(eqs, t, [], ps; name=name), flange)
end

@connector function Translational_TwoFlanges(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    compose(ODESystem([], t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialRigid
@connector function Translational_TwoFlangesRigid(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        flange_a.s ~ flange_b.s
       ]
    compose(ODESystem(eqs, t, [], []; name=name), flange_a, flange_b)
end

function Translational_Mass(;name, m = 1.0)
    @named twoflanges = Translational_TwoFlangesRigid()
    @unpack flange_a, flange_b = twoflanges
    sts = @variables s(t) v(t) a(t)
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
           s ~ flange_a.s
           D(s) ~ v 
           D(v) ~ a
           m*a ~ flange_a.f + flange_b.f
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name),twoflanges)
end

m = 1.0
@named mass = Translational_Mass(m=m)

ms_eqs = [
          mass.flange_a.f ~ 1
          mass.flange_b.f ~ 0
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass])
# Mass starts with velocity = 1
u0 = [
    mass.s => 0.0
    mass.v => 0.0
     ]

sys = structural_simplify(ms_model)
prob_complex =  ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob_complex, Tsit5())
plot(sol, vars=[mass.v,mass.s,mass.a])

[traversaro]

With a bit of manual modification of the equations to what I guess are bugs, I was able to get a solution similar to the one of OpenModelica. This is the model:

using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Translational_Flange(;name)
    sts = @variables s(t)=0.0 f(t)=0.0
    ODESystem(Equation[], t, sts, []; name=name)
end

function Translational_Fixed(;name, s0 = 0.0)
    @named flange = Translational_Flange()
    ps = @parameters s0=s0
    eqs = [flange.s ~ s0]
    compose(ODESystem(eqs, t, [], ps; name=name), flange)
end

@connector function Translational_TwoFlanges(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    compose(ODESystem([], t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialRigid
@connector function Translational_TwoFlangesRigid(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        flange_a.s ~ flange_b.s
       ]
    compose(ODESystem(eqs, t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialCompliant
@connector function Translational_TwoFlangesCompliant(;name)
    sts = @variables s_rel(t)=0.0 f(t)=0.0 
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        s_rel ~ flange_b.s - flange_a.s
        flange_b.f ~ f
        flange_a.f ~ -f
       ]
    compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Components.Mass
function Translational_Mass(;name, m = 1.0)
    @named twoflanges = Translational_TwoFlangesRigid()
    @unpack flange_a, flange_b = twoflanges
    sts = @variables s(t)=0.0 v(t)=0.0 a(t)=0.0
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
           s ~ flange_a.s
           D(s) ~ v
           D(v) ~ a
           a ~ (flange_a.f + flange_b.f) / m
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), twoflanges)
end

# Inspired from Modelica.Mechanics.Translational.Components.Spring
function Translational_Spring(;name, c = 100.0, s_rel0 = 0.0)
    @named twoflanges = Translational_TwoFlangesCompliant()
    @unpack s_rel, f = twoflanges
    ps = @parameters c=c s_rel0=s_rel0
    eqs = [
           f ~ c * (s_rel - s_rel0)
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), twoflanges)
end

m = 1.0
s0 = 0.0
c = -100.0
@named mass = Translational_Mass(m=m)
@named spring = Translational_Spring(c=c)
@named fixed = Translational_Fixed(s0=s0)

ms_eqs = [
          connect(mass.flange_b, spring.flange_a)
          connect(spring.flange_b, fixed.flange)
          mass.flange_a.f ~ 0
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass, spring, fixed])
sys = structural_simplify(ms_model)
# Mass starts with velocity = 1
u0 = [
    mass.s => 0.0
    mass.v => 1.0
     ]
prob = ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob, Tsit5())
plot(sol, vars=[mass.v,mass.s])

Env:

(@v1.7) pkg> st
      Status `~/.julia/environments/v1.7/Project.toml`
  [0c46a032] DifferentialEquations v7.1.0
  [961ee093] ModelingToolkit v8.3.0
  [91a5bcdd] Plots v1.25.6

Result:

Open points:

• W.r.t. to the equations in MSL, I had to move manually the derivative from the RHS to the LHS (this should be fixed by https://github.com/SciML/ModelingToolkit.jl/pull/1164
• W.r.t. to the equations in MSL, I had to write Newton Second Law as a ~ (flange_a.f + flange_b.f) / m instead of m*a ~ (flange_a.f + flange_b.f)
• W.r.t. to the equations in MSL and in the Modelica model, I had to set to the spring a value of -100.0 instead of 100.0, probably there is a some other - in the equations somewhere else

[traversaro]

W.r.t. to the equations in MSL, I had to write Newton Second Law as a ~ (flange_a.f + flange_b.f) / m instead of m*a ~ (flange_a.f + flange_b.f)

I opened the issue https://github.com/SciML/ModelingToolkit.jl/issues/1426 for this problem.

[traversaro]

W.r.t. to the equations in MSL and in the Modelica model, I had to set to the spring a value of -100.0 instead of 100.0, probably there is a some other - in the equations somewhere else

The problem here is for some reason in the Modelica model spring.f and mass.flange_b.f are equal, while in ModelingToolkit model they are opposite:

[traversaro]

W.r.t. to the equations in MSL and in the Modelica model, I had to set to the spring a value of -100.0 instead of 100.0, probably there is a some other - in the equations somewhere else

The problem here is for some reason in the Modelica model spring.f and mass.flange_b.f are equal, while in ModelingToolkit model they are opposite:

The problem is that in the connector definition, f in Modelica was marked as flow, while on ModelingToolkit not. By setting it to flow also in the ModellingToolkit model, now everything works as expected. The latest version of the model is:

using ModelingToolkit, Plots, DifferentialEquations

@variables t
@connector function Translational_Flange(;name)
    sts = @variables s(t)=0.0 f(t)=0.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Translational_Fixed(;name, s0 = 0.0)
    @named flange = Translational_Flange()
    ps = @parameters s0=s0
    eqs = [flange.s ~ s0]
    compose(ODESystem(eqs, t, [], ps; name=name), flange)
end

@connector function Translational_TwoFlanges(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    compose(ODESystem([], t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialRigid
@connector function Translational_TwoFlangesRigid(;name)
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        flange_a.s ~ flange_b.s
       ]
    compose(ODESystem(eqs, t, [], []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Interfaces.PartialCompliant
@connector function Translational_TwoFlangesCompliant(;name)
    sts = @variables s_rel(t)=0.0 f(t)=0.0 
    @named flange_a = Translational_Flange()
    @named flange_b = Translational_Flange()
    eqs = [
        s_rel ~ flange_b.s - flange_a.s
        flange_b.f ~ f
        flange_a.f ~ -f
       ]
    compose(ODESystem(eqs, t, sts, []; name=name), flange_a, flange_b)
end

# Inspired from Modelica.Mechanics.Translational.Components.Mass
function Translational_Mass(;name, m = 1.0)
    @named twoflanges = Translational_TwoFlangesRigid()
    @unpack flange_a, flange_b = twoflanges
    sts = @variables s(t) v(t) a(t)
    ps = @parameters m=m
    D = Differential(t)
    eqs = [
           s ~ flange_a.s
           D(s) ~ v
           D(v) ~ a
           a ~ (flange_a.f + flange_b.f) / m
          ]
    extend(ODESystem(eqs, t, sts, ps; name=name), twoflanges)
end

# Inspired from Modelica.Mechanics.Translational.Components.Spring
function Translational_Spring(;name, c = 100.0, s_rel0 = 0.0)
    @named twoflanges = Translational_TwoFlangesCompliant()
    @unpack s_rel, f = twoflanges
    ps = @parameters c=c s_rel0=s_rel0
    eqs = [
           f ~ c * (s_rel - s_rel0)
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), twoflanges)
end

m = 1.0
s0 = 0.0
c = 100.0
@named mass = Translational_Mass(m=m)
@named spring = Translational_Spring(c=c)
@named fixed = Translational_Fixed(s0=s0)

ms_eqs = [
          connect(mass.flange_b, spring.flange_a)
          connect(spring.flange_b, fixed.flange)
          mass.flange_a.f ~ 0
         ]

@named _ms_model = ODESystem(ms_eqs, t)
@named ms_model = compose(_ms_model,
                          [mass, spring, fixed])
sys = structural_simplify(ms_model)
# Mass starts with velocity = 1
u0 = [
    mass.s => 0.0
    mass.v => 1.0
     ]
prob = ODAEProblem(sys, u0, (0, 1.0))
sol = solve(prob, Tsit5())
plot(sol, vars=[mass.v,mass.s,spring.f,mass.flange_b.f])

[traversaro]

As the example was reproduced and the remaining problem are tracked in different issues, I think I can close this issue.