Skip to content

Fix initialization for TranslationalPosition and add more tests #392

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
May 13, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 7 additions & 7 deletions docs/src/connectors/connections.md
Original file line number Diff line number Diff line change
Expand Up @@ -172,7 +172,7 @@ Now, let's consider the position-based approach. We can build the same model wi
const TP = ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition

systems = @named begin
damping = TP.Damper(d = 1, va = 1, vb = 0.0)
damping = TP.Damper(d = 1)
body = TP.Mass(m = 1, v = 1)
ground = TP.Fixed(s_0 = 0)
end
Expand All @@ -191,7 +191,7 @@ nothing # hide
As can be seen, we get exactly the same result. The only difference here is that we are solving an extra equation, which allows us to plot the body position as well.

```@example connections
prob = ODEProblem(sys, [], (0, 10.0), [])
prob = ODEProblem(sys, [], (0, 10.0), [], fully_determined=true)
sol_p = solve(prob)

p1 = plot(sol_p, idxs = [body.v])
Expand Down Expand Up @@ -226,21 +226,21 @@ In this problem, we have a mass, spring, and damper which are connected to a fix

#### Damper

The damper will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point. For both position- and velocity-based domains, we set the damping constant `d=1` and `va=1` and leave the default for `v_b_0` at 0. For the position domain, we also need to set the initial positions for `flange_a` and `flange_b`.
The damper will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point. For both position- and velocity-based domains, we set the damping constant `d=1` and `va=1` and leave the default for `v_b_0` at 0.

```@example connections
@named dv = TV.Damper(d = 1)
@named dp = TP.Damper(d = 1, va = 1, vb = 0.0, flange_a__s = 3, flange_b__s = 1)
@named dp = TP.Damper(d = 1)
nothing # hide
```

#### Spring

The spring will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point. For both position- and velocity-based domains, we set the spring constant `k=1`. The velocity domain then requires the initial velocity `va` and initial spring stretch `delta_s`. The position domain instead needs the initial positions for `flange_a` and `flange_b` and the natural spring length `l`.
The spring will connect the flange/flange 1 (`flange_a`) to the mass, and flange/flange 2 (`flange_b`) to the fixed point. For both position- and velocity-based domains, we set the spring constant `k=1`. The velocity domain then requires the initial velocity `va` and initial spring stretch `delta_s`. The position domain instead needs the natural spring length `l`.

```@example connections
@named sv = TV.Spring(k = 1)
@named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
@named sp = TP.Spring(k = 1, l = 1)
nothing # hide
```

Expand Down Expand Up @@ -281,7 +281,7 @@ function simplify_and_solve(damping, spring, body, ground; initialization_eqs =

println.(full_equations(sys))

prob = ODEProblem(sys, [], (0, 10.0), []; initialization_eqs)
prob = ODEProblem(sys, [], (0, 10.0), []; initialization_eqs, fully_determined=true)
sol = solve(prob; abstol = 1e-9, reltol = 1e-9)

return sol
Expand Down
39 changes: 15 additions & 24 deletions src/Mechanical/TranslationalPosition/components.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ Flange fixed in housing at a given position.
s_0 = 0
end
@components begin
flange = Flange(; s = s_0)
flange = Flange()
end
@equations begin
flange.s ~ s_0
Expand Down Expand Up @@ -53,7 +53,7 @@ Sliding mass with inertia
f(t)
end
@components begin
flange = Flange(; s = s)
flange = Flange()
end
@equations begin
flange.s ~ s
Expand All @@ -65,7 +65,7 @@ end

const REL = Val(:relative)
@component function Spring(::Val{:relative}; name, k, va = 0.0, vb = 0.0,
delta_s = 0, flange_a__s = 0, flange_b__s = 0)
delta_s = 0)
pars = @parameters begin
k = k
end
Expand All @@ -87,56 +87,47 @@ const REL = Val(:relative)
flange_b.f ~ -f]

return compose(
ODESystem(eqs, t, vars, pars; name = name,
defaults = [
flange_a.s => flange_a__s,
flange_b.s => flange_b__s,
flange_a.f => +delta_s * k,
flange_b.f => -delta_s * k
]),
ODESystem(eqs, t, vars, pars; name = name),
flange_a,
flange_b)
end

const ABS = Val(:absolute)

"""
Spring(; name, k, flange_a__s = 0, flange_b__s = 0, l=0)
Spring(; name, k, l=0)

Linear 1D translational spring

# Parameters:

- `k`: [N/m] Spring constant
- `l`: Unstretched spring length
- `flange_a__s`: [m] Initial value of absolute position of flange_a
- `flange_b__s`: [m] Initial value of absolute position of flange_b

# Connectors:

- `flange_a: 1-dim. translational flange on one side of spring`
- `flange_b: 1-dim. translational flange on opposite side of spring` #default function
"""
function Spring(; name, k, flange_a__s = 0, flange_b__s = 0, l = 0)
Spring(ABS; name, k, flange_a__s, flange_b__s, l)
function Spring(; name, k, l = 0)
Spring(ABS; name, k, l)
end #default function

@component function Spring(::Val{:absolute};
name, k, flange_a__s = 0,
flange_b__s = 0, l = 0)
name, k, l = 0)
pars = @parameters begin
k = k
l = l
end
vars = @variables begin
f(t) = k * (flange_a__s - flange_b__s - l)
f(t)
end

@named flange_a = Flange(; s = flange_a__s, f = k * (flange_a__s - flange_b__s - l))
@named flange_b = Flange(; s = flange_b__s, f = -k * (flange_a__s - flange_b__s - l))
@named flange_a = Flange()
@named flange_b = Flange()

eqs = [
# delta_s ~ flange_a.s - flange_b.s
# delta_s ~ flange_a.s - flange_b.s
f ~ k * (flange_a.s - flange_b.s - l) #delta_s
flange_a.f ~ +f
flange_b.f ~ -f]
Expand Down Expand Up @@ -166,12 +157,12 @@ Linear 1D translational damper
@variables begin
va(t)
vb(t)
f(t) = +(va - vb) * d
f(t)
end

@components begin
flange_a = Flange(; s = 0.0, f = (va - vb) * d)
flange_b = Flange(; s = 0.0, f = -(va - vb) * d)
flange_a = Flange()
flange_b = Flange()
end

@equations begin
Expand Down
23 changes: 12 additions & 11 deletions test/Mechanical/translational.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,10 @@ end

@testset "Spring, Damper, Mass, Fixed" begin
@named dv = TV.Damper(d = 1)
@named dp = TP.Damper(d = 1, va = 1, vb = 0.0, flange_a.s = 3, flange_b.s = 1)
@named dp = TP.Damper(d = 1)

@named sv = TV.Spring(k = 1)
@named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
@named sp = TP.Spring(k = 1, l = 1)

@named bv = TV.Mass(m = 1)
@named bp = TP.Mass(m = 1, v = 1, s = 3)
Expand All @@ -52,7 +52,8 @@ end

sys = structural_simplify(model)

prob = ODEProblem(sys, [], (0, 20.0), []; initialization_eqs)
prob = ODEProblem(
sys, [], (0, 20.0), []; initialization_eqs, fully_determined = true)
sol = solve(prob; abstol = 1e-9, reltol = 1e-9)

return sol
Expand All @@ -71,10 +72,10 @@ end

@testset "driven spring damper mass" begin
@named dv = TV.Damper(d = 1)
@named dp = TP.Damper(d = 1, va = 1.0, vb = 0.0, flange_a.s = 3, flange_b.s = 1)
@named dp = TP.Damper(d = 1)

@named sv = TV.Spring(k = 1)
@named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
@named sp = TP.Spring(k = 1, l = 1)

@named bv = TV.Mass(m = 1)
@named bp = TP.Mass(m = 1, v = 1, s = 3)
Expand All @@ -101,12 +102,13 @@ end

model = System(dv, sv, bv, gv, fv, source)
sys = structural_simplify(model)
prob = ODEProblem(sys, [bv.s => 0], (0, 20.0), [])
prob = ODEProblem(
sys, [bv.s => 0, sv.delta_s => 1], (0, 20.0), [], fully_determined = true)
solv = solve(prob, Rodas4())

model = System(dp, sp, bp, gp, fp, source)
sys = structural_simplify(model)
prob = ODEProblem(sys, [], (0, 20.0), [])
prob = ODEProblem(sys, [], (0, 20.0), [], fully_determined = true)
solp = solve(prob, Rodas4())

for sol in (solv, solp)
Expand Down Expand Up @@ -186,8 +188,7 @@ end
function mass_spring(; name)
systems = @named begin
fixed = TP.Fixed()
spring = TP.Spring(;
k = 10.0, l = 1.0, flange_a__s = 0.0, flange_b__s = 2.0)
spring = TP.Spring(; k = 10.0, l = 1.0)
mass = TP.Mass(; m = 100.0, s = 2.0, v = 0.0)
pos_sensor = TP.PositionSensor()
force_sensor = TP.ForceSensor()
Expand All @@ -207,7 +208,7 @@ end
@named model = mass_spring()
sys = structural_simplify(model)

prob = ODEProblem(sys, [], (0.0, 1.0))
prob = ODEProblem(sys, [], (0.0, 1.0), fully_determined = true)
sol = solve(prob, Tsit5())

@test all(sol[sys.spring.flange_a.f] .== sol[sys.force_value.u])
Expand All @@ -230,7 +231,7 @@ end
@named sys = ODESystem(
eqs, t, [], []; systems = [force, source, mass, acc, acc_output])
s = complete(structural_simplify(sys))
prob = ODEProblem(s, [mass.s => 0], (0.0, pi))
prob = ODEProblem(s, [mass.s => 0], (0.0, pi), fully_determined = true)
sol = solve(prob, Tsit5())
@test sol[sys.acc_output.u] ≈ (sol[sys.mass.f] ./ m)
end
Expand Down
Loading