Merge branch 'SciML:main' into bgc/revolute

Author: bradcarman

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkitStandardLibrary"
 uuid = "16a59e39-deab-5bd0-87e4-056b12336739"
 authors = ["Chris Rackauckas and Julia Computing"]
-version = "1.6.0"
+version = "1.7.0"
 
 [deps]
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"

docs/make.jl

@@ -7,6 +7,8 @@ using ModelingToolkitStandardLibrary.Magnetic.FluxTubes
 using ModelingToolkitStandardLibrary.Electrical
 using ModelingToolkitStandardLibrary.Thermal
 
+ENV["GKSwstype"] = "100"
+
 include("pages.jl")
 
 makedocs(sitename = "ModelingToolkitStandardLibrary.jl",

src/Blocks/continuous.jl

@@ -414,16 +414,29 @@ function LimPID(; name, k = 1, Ti = false, Td = false, wp = 1, wd = 1,
 end
 
 """
-    StateSpace(A, B, C, D=0; x_start=zeros(size(A,1)), name)
+    StateSpace(A, B, C, D=0; x_start=zeros(size(A,1)), u0=zeros(size(B,2)), y0=zeros(size(C,1)), name)
 
 A linear, time-invariant state-space system on the form.
-```
+```math
 ẋ = Ax + Bu
 y = Cx + Du
 ```
 Transfer functions can also be simulated by converting them to a StateSpace form.
+
+`y0` and `u0` can be used to set an operating point, providing them changes the dynamics from an LTI system to the affine system 
+```math
+ẋ = Ax + B(u - u0)
+y = Cx + D(u - u0) + y0
+```
+For a nonlinear system
+```math
+ẋ = f(x, u)
+y = h(x, u)
+```
+linearized around the operating point `x₀, u₀`, we have `y0, u0 = h(x₀, u₀), u₀`.
 """
-function StateSpace(; A, B, C, D = nothing, x_start = zeros(size(A, 1)), name)
+function StateSpace(; A, B, C, D = nothing, x_start = zeros(size(A, 1)), name,
+                    u0 = zeros(size(B, 2)), y0 = zeros(size(C, 1)))
     nx, nu, ny = size(A, 1), size(B, 2), size(C, 1)
     size(A, 2) == nx || error("`A` has to be a square matrix.")
     size(B, 1) == nx || error("`B` has to be of dimension ($nx x $nu).")
@@ -442,9 +455,11 @@ function StateSpace(; A, B, C, D = nothing, x_start = zeros(size(A, 1)), name)
     # pars = @parameters A=A B=B C=C D=D # This is buggy
     eqs = [ # FIXME: if array equations work
         [Differential(t)(x[i]) ~ sum(A[i, k] * x[k] for k in 1:nx) +
-                                 sum(B[i, j] * input.u[j] for j in 1:nu) for i in 1:nx]..., # cannot use D here
+                                 sum(B[i, j] * (input.u[j] - u0[j]) for j in 1:nu)
+         for i in 1:nx]..., # cannot use D here
         [output.u[j] ~ sum(C[j, i] * x[i] for i in 1:nx) +
-                       sum(D[j, k] * input.u[k] for k in 1:nu) for j in 1:ny]...,
+                       sum(D[j, k] * (input.u[k] - u0[k]) for k in 1:nu) + y0[j]
+         for j in 1:ny]...,
     ]
     compose(ODESystem(eqs, t, vcat(x...), [], name = name), [input, output])
 end

src/Mechanical/Translational/Translational.jl

@@ -17,7 +17,7 @@ include("components.jl")
 export PositionSensor
 include("sensors.jl")
 
-#export SineForce
+export Force
 include("sources.jl")
 
 end

src/Mechanical/Translational/sources.jl

@@ -1,15 +1,13 @@
-#=
-function SineForce(; name, amp, freq)
+using ModelingToolkitStandardLibrary.Blocks
+
+function Force(; name)
     @named flange = MechanicalPort()
+    @named f = RealInput()
 
-    pars = @parameters begin
-        amp = amp
-        freq = freq
-    end
-    vars = []
+    vars = pars = []
     eqs = [
-        flange.f ~ amp * sin(2 * π * t * freq),
+        flange.f ~ -f.u,
     ]
-    compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]), flange)
+    compose(ODESystem(eqs, t, vars, pars; name = name, defaults = [flange.v => 0]),
+            [flange, f])
 end
-=#

test/Blocks/continuous.jl

@@ -108,6 +108,25 @@ end
     # equilibrium point is at [1, 0]
     @test sol[ss.x[1]][end]≈1 atol=1e-3
     @test sol[ss.x[2]][end]≈0 atol=1e-3
+
+    # non-zero operating point
+    u0 = [1] # This causes no effective input to the system since c.k = 1
+    y0 = [2]
+    @named ss = StateSpace(; A, B, C, D, x_start = zeros(2), u0, y0)
+    @named model = ODESystem([
+                                 connect(c.output, ss.input),
+                             ],
+                             t,
+                             systems = [ss, c])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, Pair[], (0.0, 100.0))
+    sol = solve(prob, Rodas4())
+    @test sol.retcode == :Success
+
+    @test sol[ss.x[1]][end ÷ 2]≈0 atol=1e-3 # Test that x did not move
+    @test sol[ss.x[1]][end]≈0 atol=1e-3 # Test that x did not move
+    @test sol[ss.x[2]][end]≈0 atol=1e-3
+    @test sol[ss.output.u[1]][end]≈y0[] atol=1e-3 # Test that the output equals the operating point
 end
 
 """Second order demo plant"""

test/Mechanical/translational.jl

@@ -1,5 +1,6 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
+using ModelingToolkitStandardLibrary.Blocks
 import ModelingToolkitStandardLibrary.Mechanical.Translational as TV
 import ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition as TP
 
@@ -44,3 +45,50 @@ D = Differential(t)
     @test solp[bp.v][1] == 1.0
     @test solp[bp.v][end]≈0.0 atol=1e-4
 end
+
+@testset "driven spring damper mass" begin
+    @named dv = TV.Damper(d = 1, v_a_0 = 1)
+    @named dp = TP.Damper(d = 1, v_a_0 = 1, s_a_0 = 3, s_b_0 = 1)
+
+    @named sv = TV.Spring(k = 1, v_a_0 = 1, delta_s_0 = 1)
+    @named sp = TP.Spring(k = 1, s_a_0 = 3, s_b_0 = 1, l = 1)
+
+    @named bv = TV.Mass(m = 1, v_0 = 1)
+    @named bp = TP.Mass(m = 1, v_0 = 1, s_0 = 3)
+
+    @named gv = TV.Fixed()
+    @named gp = TP.Fixed(s_0 = 1)
+
+    @named fv = TV.Force()
+    @named fp = TP.Force()
+
+    @named source = Sine(frequency = 3, amplitude = 2)
+
+    function simplify_and_solve(damping, spring, body, ground, f, source)
+        eqs = [connect(f.f, source.output)
+               connect(f.flange, body.flange)
+               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, f, source])
+
+        sys = structural_simplify(model)
+
+        println.(full_equations(sys))
+
+        prob = ODEProblem(sys, [], (0, 20.0), [])
+        sol = solve(prob, Rodas4())
+
+        return sol
+    end
+
+    solv = simplify_and_solve(dv, sv, bv, gv, fv, source)
+    solp = simplify_and_solve(dp, sp, bp, gp, fp, source)
+
+    for sol in (solv, solp)
+        lb, ub = extrema(solv(15:0.05:20, idxs = bv.v).u)
+        @test -lb≈ub atol=1e-2
+        @test -0.11 < lb < -0.1
+    end
+end