Merge pull request #334 from SciML/bgc/siso_defaults

Remove defaults from Blocks.SISO

Author: ChrisRackauckas

docs/src/connectors/connections.md

@@ -136,8 +136,8 @@ using ModelingToolkitStandardLibrary
 const TV = ModelingToolkitStandardLibrary.Mechanical.Translational
 
 systems = @named begin
-    damping = TV.Damper(d = 1, flange_a.v = 1)
-    body = TV.Mass(m = 1, v = 1)
+    damping = TV.Damper(d = 1)
+    body = TV.Mass(m = 1)
     ground = TV.Fixed()
 end
 
@@ -155,7 +155,8 @@ nothing # hide
 As expected, we have a similar solution…
 
 ```@example connections
-prob = ODEProblem(sys, [], (0, 10.0), [])
+prob = ODEProblem(
+    sys, [], (0, 10.0), []; initialization_eqs = [sys.body.s ~ 0, sys.body.v ~ 1])
 sol_v = solve(prob)
 
 p1 = plot(sol_v, idxs = [body.v])
@@ -228,7 +229,7 @@ In this problem, we have a mass, spring, and damper which are connected to a fix
 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`.
 
 ```@example connections
-@named dv = TV.Damper(d = 1, flange_a.v = 1)
+@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)
 nothing # hide
 ```
@@ -238,7 +239,7 @@ nothing # hide
 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`.
 
 ```@example connections
-@named sv = TV.Spring(k = 1, flange_a__v = 1, delta_s = 1)
+@named sv = TV.Spring(k = 1)
 @named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
 nothing # hide
 ```
@@ -248,7 +249,7 @@ nothing # hide
 For both position- and velocity-based domains, we set the mass `m=1` and initial velocity `v=1`. Like the damper, the position domain requires the position initial conditions set as well.
 
 ```@example connections
-@named bv = TV.Mass(m = 1, v = 1)
+@named bv = TV.Mass(m = 1)
 @named bp = TP.Mass(m = 1, v = 1, s = 3)
 nothing # hide
 ```
@@ -270,7 +271,7 @@ As can be seen, the position-based domain requires more initial condition inform
 Let's define a quick function to simplify and solve the 2 different systems. Note, we will solve with a fixed time step and a set tolerance to compare the numerical differences.
 
 ```@example connections
-function simplify_and_solve(damping, spring, body, ground)
+function simplify_and_solve(damping, spring, body, ground; initialization_eqs = Equation[])
     eqs = [connect(spring.flange_a, body.flange, damping.flange_a)
            connect(spring.flange_b, damping.flange_b, ground.flange)]
 
@@ -280,8 +281,8 @@ function simplify_and_solve(damping, spring, body, ground)
 
     println.(full_equations(sys))
 
-    prob = ODEProblem(sys, [], (0, 10.0), [])
-    sol = solve(prob; dt = 0.1, adaptive = false, reltol = 1e-9, abstol = 1e-9)
+    prob = ODEProblem(sys, [], (0, 10.0), []; initialization_eqs)
+    sol = solve(prob; abstol = 1e-9, reltol = 1e-9)
 
     return sol
 end
@@ -291,7 +292,10 @@ nothing # hide
 Now let's solve the velocity domain model
 
 ```@example connections
-solv = simplify_and_solve(dv, sv, bv, gv);
+initialization_eqs = [bv.s ~ 3
+                      bv.v ~ 1
+                      sv.delta_s ~ 1]
+solv = simplify_and_solve(dv, sv, bv, gv; initialization_eqs);
 nothing # hide
 ```
 
@@ -346,12 +350,11 @@ By definition, the spring stretch is
 \Delta s = s - s_{b_0} - l
 ```
 
-Which means both systems are actually solving the same exact system.  We can plot the numerical difference between the 2 systems and see the result is negligible.
+Which means both systems are actually solving the same exact system.  We can plot the numerical difference between the 2 systems and see the result is negligible (much less than the tolerance of 1e-9).
 
 ```@example connections
 plot(title = "numerical difference: vel. vs. pos. domain", xlabel = "time [s]",
     ylabel = "solv[bv.v] .- solp[bp.v]")
 time = 0:0.1:10
 plot!(time, (solv(time)[bv.v] .- solp(time)[bp.v]), label = "")
-Plots.ylims!(-1e-15, 1e-15)
 ```

src/Blocks/nonlinear.jl

@@ -93,20 +93,30 @@ Initial value of state `Y` can be set with `int.y`
   - `input`
   - `output`
 """
-@mtkmodel SlewRateLimiter begin
-    @parameters begin
-        rising = 1.0, [description = "Maximum rising slew rate of SlewRateLimiter"]
-        falling = -rising, [description = "Derivative time constant of SlewRateLimiter"]
-        Td = 0.001, [description = "Derivative time constant"]
-    end
-    begin
-        getdefault(rising) ≥ getdefault(falling) ||
-            throw(ArgumentError("`rising` must be smaller than `falling`"))
-        getdefault(Td) > 0 ||
-            throw(ArgumentError("Time constant `Td` must be strictly positive"))
+@component function SlewRateLimiter(;
+        name, y_start = 0.0, rising = 1.0, falling = -rising, Td = 0.001)
+    pars = @parameters begin
+        rising = rising, [description = "Maximum rising slew rate of SlewRateLimiter"]
+        falling = falling, [description = "Derivative time constant of SlewRateLimiter"]
+        Td = Td, [description = "Derivative time constant"]
+        y_start = y_start
     end
-    @extend u, y = siso = SISO(; y_start)
-    @equations begin
+
+    getdefault(rising) ≥ getdefault(falling) ||
+        throw(ArgumentError("`rising` must be smaller than `falling`"))
+    getdefault(Td) > 0 ||
+        throw(ArgumentError("Time constant `Td` must be strictly positive"))
+
+    @named siso = SISO(; y_start)
+    @unpack y, u = siso
+
+    eqs = [
         D(y) ~ max(min((u - y) / Td, rising), falling)
-    end
+    ]
+
+    initialization_eqs = [
+        y ~ y_start
+    ]
+
+    return extend(ODESystem(eqs, t, [], pars; name, initialization_eqs), siso)
 end

src/Blocks/utils.jl

@@ -130,8 +130,8 @@ Single input single output (SISO) continuous system block.
         y_start = 0.0
     end
     @variables begin
-        u(t) = u_start, [description = "Input of SISO system"]
-        y(t) = y_start, [description = "Output of SISO system"]
+        u(t), [guess = u_start, description = "Input of SISO system"]
+        y(t), [guess = y_start, description = "Output of SISO system"]
     end
     @components begin
         input = RealInput(guess = u_start)

src/Mechanical/Translational/components.jl

@@ -65,9 +65,9 @@ Sliding mass with inertia
     @named flange = MechanicalPort()
 
     vars = @variables begin
-        s(t)
-        v(t)
-        f(t)
+        s(t), [guess = 0]
+        v(t), [guess = 0]
+        f(t), [guess = 0]
     end
 
     eqs = [flange.v ~ v
@@ -98,22 +98,21 @@ Linear 1D translational spring
   - `flange_a`: 1-dim. translational flange on one side of spring
   - `flange_b`: 1-dim. translational flange on opposite side of spring
 """
-@component function Spring(; name, k, delta_s = 0.0, flange_a__v = 0.0, flange_b__v = 0.0)
-    Spring(REL; name, k, delta_s, flange_a__v, flange_b__v)
+@component function Spring(; name, k)
+    Spring(REL; name, k)
 end # default
 
-@component function Spring(::Val{:relative}; name, k, delta_s = 0.0, flange_a__v = 0.0,
-        flange_b__v = 0.0)
+@component function Spring(::Val{:relative}; name, k)
     pars = @parameters begin
         k = k
     end
     vars = @variables begin
-        delta_s(t) = delta_s
-        f(t)
+        delta_s(t), [guess = 0]
+        f(t), [guess = 0]
     end
 
-    @named flange_a = MechanicalPort(; v = flange_a__v)
-    @named flange_b = MechanicalPort(; v = flange_b__v)
+    @named flange_a = MechanicalPort()
+    @named flange_b = MechanicalPort()
 
     eqs = [D(delta_s) ~ flange_a.v - flange_b.v
            f ~ k * delta_s
@@ -125,20 +124,19 @@ end # default
 end
 
 const ABS = Val(:absolute)
-@component function Spring(::Val{:absolute}; name, k, sa = 0, sb = 0, flange_a__v = 0.0,
-        flange_b__v = 0.0, l = 0)
+@component function Spring(::Val{:absolute}; name, k, l = 0)
     pars = @parameters begin
         k = k
         l = l
     end
     vars = @variables begin
-        sa(t) = sa
-        sb(t) = sb
-        f(t)
+        sa(t), [guess = 0]
+        sb(t), [guess = 0]
+        f(t), [guess = 0]
     end
 
-    @named flange_a = MechanicalPort(; v = flange_a__v)
-    @named flange_b = MechanicalPort(; v = flange_b__v)
+    @named flange_a = MechanicalPort()
+    @named flange_b = MechanicalPort()
 
     eqs = [D(sa) ~ flange_a.v
            D(sb) ~ flange_b.v
@@ -169,8 +167,8 @@ Linear 1D translational damper
         d
     end
     @variables begin
-        v(t)
-        f(t)
+        v(t), [guess = 0]
+        f(t), [guess = 0]
     end
 
     @components begin

test/Blocks/utils.jl

@@ -1,5 +1,7 @@
 using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkit
+using OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @testset "Array Guesses" begin
     for (block, guess) in [
@@ -23,6 +25,22 @@ end
     end
 end
 
+@testset "SISO Check" begin
+    k, w, d = 1.0, 1.0, 0.5
+    @named c = Constant(; k = 1)
+    @named so = SecondOrder(; k = k, w = w, d = d, xd = 1)
+    @named iosys = ODESystem(connect(c.output, so.input), t, systems = [so, c])
+    sys = structural_simplify(iosys)
+
+    initsys = ModelingToolkit.generate_initializesystem(sys)
+    initsys = structural_simplify(initsys)
+    initprob = NonlinearProblem(initsys, [t => 0])
+    initsol = solve(initprob)
+
+    @test initsol[sys.so.xd] == 1.0
+    @test initsol[sys.so.u] == 1.0
+end
+
 @test_deprecated RealInput(; name = :a, u_start = 1.0)
 @test_deprecated RealInput(; name = :a, nin = 2, u_start = ones(2))
 @test_deprecated RealOutput(; name = :a, u_start = 1.0)

test/Mechanical/translational.jl

@@ -34,7 +34,7 @@ end
     @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 sv = TV.Spring(k = 1, flange_a__v = 1, delta_s = 1)
+    @named sv = TV.Spring(k = 1)
     @named sp = TP.Spring(k = 1, flange_a__s = 3, flange_b__s = 1, l = 1)
 
     @named bv = TV.Mass(m = 1)
@@ -43,21 +43,23 @@ end
     @named gv = TV.Fixed()
     @named gp = TP.Fixed(s_0 = 1)
 
-    function simplify_and_solve(damping, spring, body, ground)
+    function simplify_and_solve(
+            damping, spring, body, ground; initialization_eqs = Equation[])
         eqs = [connect(spring.flange_a, body.flange, damping.flange_a)
                connect(spring.flange_b, damping.flange_b, ground.flange)]
 
         @named model = ODESystem(eqs, t; systems = [ground, body, spring, damping])
 
         sys = structural_simplify(model)
 
-        prob = ODEProblem(sys, [body.s => 0], (0, 20.0), [])
-        sol = solve(prob, ImplicitMidpoint(), dt = 0.01)
+        prob = ODEProblem(sys, [], (0, 20.0), []; initialization_eqs)
+        sol = solve(prob; abstol = 1e-9, reltol = 1e-9)
 
         return sol
     end
 
-    solv = simplify_and_solve(dv, sv, bv, gv)
+    solv = simplify_and_solve(
+        dv, sv, bv, gv; initialization_eqs = [bv.s ~ 3, bv.v ~ 1, sv.delta_s ~ 1])
     solp = simplify_and_solve(dp, sp, bp, gp)
 
     @test solv[bv.v][1] == 1.0