Add translational library

[kianmolani]

This PR is linked to issue #84

[ValentinKaisermayer]

Doc strings are missing and actual tests

[baggepinnen]

Have a look at https://github.com/baggepinnen/AutomaticDocstrings.jl to quickly generate docstring stubs, helps a lot when you have a lot of docstrings to generate

[kianmolani]

Can you add the connectors in the doc strings as well? And name the states s(t) to be consistent with the rest of the library.

If you want to add images, simply drag them into a comment or an issue.

Thank you for the suggestions, Valentin. The connectors in utils.js already have doc strings added (I mirrored the documentation from the same file in the Rotational lib). Are you referring to something else? Also, what do you mean about naming the states s(t)? The position state is already denoted by s.

[ValentinKaisermayer]

I meant a # Connectors: section in the doc strings, listing the connectors of the components. You can have a look at the doc strings in the other modules.

[kianmolani]

I meant a # Connectors: section in the doc strings, listing the connectors of the components. You can have a look at the doc strings in the other modules.

Okay, that's clear. Will do, and thanks!

[kianmolani]

There appear to be issues with our simulated results for a mass-spring-damper model. Julia and Modelica simulation results do not match, and the difference is very apparent. I've attached below screenshots of Julia and Modelica simulation results of our mass position vs. time for underdamped, overdamped, and critically damped systems. For the underdamped system, the difference appears to be accounted for with an offset of about 0.5, as well as a slight phase shift. For the overdamped and critically damped systems, the solutions possess different equilibrium states.

I've ensured that both Modelica and Julia simulations are initialized in the same way. Note that simulation results between Modelica and Julia for a mass-damper model are very close (tested across a range of parameter values).

This last screenshot shows how I've defined my system connections for a mass-spring-damper model:

[baggepinnen]

Could the difference between julia and modelica lie in the initialization? Is there a difference already at the first point saved in the solution?

[kianmolani]

Could the difference between julia and modelica lie in the initialization? Is there a difference already at the first point saved in the solution?

I believe that they've been initialized in the same way. Though I will triple-check this (the solutions match for the mass-damper system, tested across a range of param values, so perhaps the culprit here is the spring). I've attached the Modelica simulation settings below, though I don't know if they'd be of use to you (the simulation interval settings are the same, and so is the integration tolerance level).

[baggepinnen]

Let me see if I understand this correctly, the fixed point is at s = 4.5 and the spring has resting length s_rel0 = 1, then I'd expect the steady-state position to be either 5.5 or 3.5 depending on the direction the spring is pointing in. It thus looks like the Julia solution is correct since it comes to rest at 3.5? I assume there's no gravity or other external forces acting anywhere? Is the length of the spring and all other params the same in modelica?

[kianmolani]

Let me see if I understand this correctly, the fixed point is at s = 4.5 and the spring has resting length s_rel0 = 1, then I'd expect the steady-state position to be either 5.5 or 3.5 depending on the direction the spring is pointing in. It thus looks like the Julia solution is correct since it comes to rest at 3.5? I assume there's no gravity or other external forces acting anywhere? Is the length of the spring and all other params the same in modelica?

No other forces and all other params confirmed to be the same.

The issue might be that Modelica takes in an additional parameter that we do not: mass length (set to 1 in this simulation). Because Modelica also has the parameter s_start (the initial value of the absolute position of the mass), and because this parameter is set to 3, it means that Modelica is tracking the center of this mass. Therefore, the equilibrium state at x = 3 makes sense for the Modelica simulation.

However, again, we do not take in a parameter for mass length, and so I guess that we're treating the mass as a point mass? In any case, initializing s_start to also equal 3 in our simulation (given that we do not have a mass length) might be problematic in that there's a discontinuity in the connections. I've drawn two schematics to illustrate this point:

For Modelica:

For Julia:

[kianmolani]

Indeed, I've just confirmed that setting the spring length to 1.5m in our Julia simulation in order to account for this "gap" yields matching simulation results.

[kianmolani]

Perhaps mass length is a property that we should include as well.

[codecov[bot]]

Codecov Report

Merging #85 (f7f188e) into main (db37027) will decrease coverage by 5.56%. The diff coverage is 0.00%.

@@            Coverage Diff             @@
##             main      #85      +/-   ##
==========================================
- Coverage   68.59%   63.03%   -5.57%     
==========================================
  Files          23       26       +3     
  Lines         952     1036      +84     
==========================================
  Hits          653      653              
- Misses        299      383      +84     

Impacted Files	Coverage Δ
src/Mechanical/Translational/components.jl	0.00% <0.00%> (ø)
src/Mechanical/Translational/sources.jl	0.00% <0.00%> (ø)
src/Mechanical/Translational/utils.jl	0.00% <0.00%> (ø)

:mega: We’re building smart automated test selection to slash your CI/CD build times. Learn more

[YingboMa]

You should remove .vscode and .DS_Store, then the PR is good to go after formatting.

[kianmolani]

You should remove .vscode and .DS_Store, then the PR is good to go after formatting.

All .vscode and .DS_Store files have now been deleted.

[YingboMa]

Format failure.

[bradcarman]

Hi, I'm reviewing the code and I noticed that the Flange connector uses absolute position for the state variable. Might I suggest changing this to velocity? The reason for this is so we have less matching parameters to set. Let me show an example...

The below code requires that each of the 3 components have a position parameter set, and s_rel0 is really just s0-s_start, so it can already be seen as redundant.

s0 = 4.5 # fixed offset position of flange housing
s_start = 3 # initial value of absolute position of sliding mass
s_rel0 = 1.5 # unstretched spring length
v_start = 10 # initial value of absolute linear velocity of sliding mass

@named fixed = Fixed(; s0)
@named mass = Mass(; m=1, s_start, v_start)
@named spring = Spring(; c=10, s_rel0)

eqs = [
    connect(mass.flange_a, spring.flange_a)
    connect(spring.flange_b, fixed.flange)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01)

Plotting sol[mass.s]...

If I remove the position parameters from Fixed and Spring, you can see we get the incorrect result...

@named fixed = Fixed()
@named mass = Mass(; m, s_start, v_start)
@named spring = Spring(; c)

Now, for comparison (using #107) we can see that only the Body position is needed, which to me makes more sense because we really only need to track absolute positions of actual moving rigid bodies. Additionally the Fixed component is defined by simply setting the port velocity to 0.

x0 = 3 # initial value of absolute position of sliding mass
v0 = 10 # initial value of absolute linear velocity of sliding mass

function Fixed(; name)
    @named port = Port()
    eqs = [port.v ~ 0]
    return compose(ODESystem(eqs, t, [], []; name = name), port)
end

@named fixed = Fixed()
@named mass = Body(; m=1, x0, v0)
@named spring = Spring(; k=10)

eqs = [
    connect(mass.port, spring.port_a)
    connect(spring.port_b, fixed.port)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0, 10.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01)

This is also how other acausal modeling tools I know of define their thru and across (flow) variables for translational components.

[baggepinnen]

Might I suggest changing this to velocity?

The intention is to follow modelica stdlib, if we're deviating from them, it's probably an oversight and should be fixed :)

[ValentinKaisermayer]

Although, Simscape uses velocity and angular velocity for its mechanical domains.

In the example you show, I think you do not need to specify s_rel0 for the spring. It should be fully determined by the mass and fixed position.

[bradcarman]

In this example the spring starts at it's natural length (i.e. uncompressed), therefore one should not need to specify s_rel0. By using position rather than velocity, this overcomplicates the model and now one must specify both the position of the Fixed component as well as s_rel0 for the spring.

Additionally, because the Flange is defined with s defaulting to 0, we now introduce a requirement that the model must solve for it's initial states. For a small model, this is not such a big deal, but for the large models I'm dealing with, initialization is something I'd like to avoid for 1) I need my model to start immediately, waiting for initialization is not an option for real time control, 2) I'd like to guarantee my states.

If we run this simple model without initialization and simplification, it can be seen that the result is not good.

Full code below...

using ModelingToolkitStandardLibrary.Mechanical.Translational, ModelingToolkit,
      OrdinaryDiffEq, Test
using Setfield

@parameters t
D = Differential(t)

s0 = 4.5 # fixed offset position of flange housing
s_start = 3 # initial value of absolute position of sliding mass
s_rel0 = 1.5 # unstretched spring length
v_start = 10 # initial value of absolute linear velocity of sliding mass

@named fixed = Fixed(; s0)
@named mass = Mass(; m=1, s_start, v_start)
@named spring = Spring(; c=10, s_rel0)

eqs = [
    connect(mass.flange_a, spring.flange_a)
    connect(spring.flange_b, fixed.flange)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = expand_connections(model)
@set! sys.eqs = [(typeof(eq.lhs) != ModelingToolkit.Connection) & !ModelingToolkit._iszero(eq.lhs) & !ModelingToolkit.isdifferential(eq.lhs) ? 0 ~ eq.rhs - eq.lhs : eq for eq in sys.eqs]

prob = ODEProblem(sys, [], (0, 2.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01, initializealg=NoInit())

using CairoMakie
begin
    f,a,=lines(sol.t, sol[mass.s], label="sol[mass.s]");    
    lines!(a, sol.t, sol[mass.flange_a.s], label="sol[mass.flange_a.s]")
    lines!(a, sol.t, sol[spring.flange_a.s], label="sol[spring.flange_a.s]", linestyle=:dash)
    lines!(a, sol.t, sol[spring.flange_b.s], label="sol[spring.flange_b.s]")
    lines!(a, sol.t, sol[fixed.flange.s], label="sol[fixed.flange.s]", linestyle=:dash)
    a.xlabel="time [s]"; a.ylabel="position [m]"; 
    Legend(f[1,2], a)
    f
end

Additionally we can see that spring.s_rel is not automatically resolved

In contrast, by using velocity and force for the translational Connection, we only need to specify the position of the Body/Mass, and no initialization is required (using #107)...

using ModelingToolkitStandardLibrary.Mechanical.Translational, ModelingToolkit,
      OrdinaryDiffEq, Test
using Setfield

@parameters t
D = Differential(t)

x0 = 3 # initial value of absolute position of sliding mass
v0 = 10 # initial value of absolute linear velocity of sliding mass

function Fixed(; name)
    @named port = Port()
    eqs = [port.v ~ 0]
    return compose(ODESystem(eqs, t, [], []; name = name), port)
end

@named fixed = Fixed()
@named mass = Body(; m=1, x0, v0)
@named spring = Spring(; k=10)

eqs = [
    connect(mass.port, spring.port_a)
    connect(spring.port_b, fixed.port)
]

@named model = ODESystem(eqs, t, [], []; systems=[fixed, mass, spring])

sys = expand_connections(model)
@set! sys.eqs = [(typeof(eq.lhs) != ModelingToolkit.Connection) & !ModelingToolkit._iszero(eq.lhs) & !ModelingToolkit.isdifferential(eq.lhs) ? 0 ~ eq.rhs - eq.lhs : eq for eq in sys.eqs]
prob = ODEProblem(sys, [], (0, 2.0))
sol = solve(prob, ImplicitMidpoint(), dt=0.01, initializealg=NoInit())

using CairoMakie
begin
    f,a,=lines(sol.t, sol[mass.x], label="sol[mass.x]");    
    a.xlabel="time [s]"; a.ylabel="position [m]"; 
    Legend(f[1,2], a)
    f
end

[bradcarman]

Making my point in another way, one should be able to define the Fixed component such that position is not needed, we simply specify that the change in position is zero...

function Fixed(; name)
    @named flange = Flange()
    eqs = [D(flange.s) ~ 0]
    return compose(ODESystem(eqs, t, [], []; name = name), flange)
end

But, unfortunately this will not work in the end because connect( fixed.flange, component.flange ) will give the equation...

fixed.flange.s ~ component.flange.s

Therefore, we are forced to specify the position for every single Flange in the model. But in theory, it should not be necessary to specify the position for a Fixed component, this is really overly defining the model.

[Firionus]

@bradcarman How would you express in your approach a spring that has non-zero length in its relaxed state? That would certainly be a important for someone that wants to integrate with visualization or if someone wants to model something precisely the way they imagine their coordinate system.

[bradcarman]

One final point (then I'll shut up, ha ha), the necessity that one must specify a spring length for the Spring component is silly when one is considering models that really just deal with stiffness. For example, defining a Hard Stop component that defines the stiffness of the boundary. Should the length/thickness of the boundary really need to be defined? Mathematically this information is irrelevant. My recommend spring definition then is...

function Spring(;name, k = 1e3, delta0 = 0.0)
    pars = @parameters k=k delta0=delta0
    vars = @variables begin
        Δx(t) = delta0
        f(t) = k*delta0
    end
    
    @named port_a = Port()
    @named port_b = Port()
    
    eqs = [
        D(Δx) ~ port_a.v - port_b.v
        f ~ k*Δx
        port_a.f ~ +f
        port_b.f ~ -f
    ]
    return compose(ODESystem(eqs, t, vars, pars; name = name), port_a, port_b)
end

Then, we only need to define 2 parameters, the stiffness, and the initial compression state (delta0).

And to answer @Firionus, for visualization, if we are plotting the position states, the only states tracked by the model are actual moving bodies/masses. So if we have 2 bodies connected by a spring, then we can plot both body positions.

[Firionus]

And to answer @Firionus, for visualization, if we are plotting the position states, the only states tracked by the model are actual moving bodies/masses. So if we have 2 bodies connected by a spring, then we can plot both body positions.

This doesn't hold true if the spring is connected to a fixed element. Also, visualizing springs with a finite relaxed length is much closer to a real mechanical setup than springs that are infinitely small sometimes, don't you think?

One final point (then I'll shut up, ha ha), the necessity that one must specify a spring length for the Spring component is silly when one is considering models that really just deal with stiffness.

If you find it silly you can just leave it at the default of s_rel0 = 0.0 and achieve that. I realize you'll still have to define the position of fixed points but I think this is a small price to pay for gaining better modelling of real life with the absolute position approach.

Ultimately, it comes down to the choice of programming interface. Should a "mechanical component" only keep track of velocities and forces, i.e. its energy exchange, or should it also track its absolute position. The Modelica Standard Library (MSL 3.2.3) goes with the second approach:

connector Flange "One-dimensional translational flange"

  SI.Position s "Absolute position of flange";
  flow SI.Force f "Cut force directed into flange";
  annotation (...);
end Flange;

And since @baggepinnen said:

The intention is to follow modelica stdlib

This PR gets that right, since it follows the MSL:

@connector function Flange(; name)
    sts = @variables begin
        s(t)
        f(t), [connect = Flow]
    end
    ODESystem(Equation[], t, sts, [], name = name, defaults = Dict(s => 0.0, f => 0.0))
end

source