Merge branch 'SciML:main' into heating

Author: langestefan

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkitStandardLibrary"
 uuid = "16a59e39-deab-5bd0-87e4-056b12336739"
 authors = ["Chris Rackauckas and Julia Computing"]
-version = "2.18.0"
+version = "2.19.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

README.md

@@ -50,30 +50,34 @@ using ModelingToolkitStandardLibrary.Electrical
 using ModelingToolkitStandardLibrary.Blocks: Constant
 using ModelingToolkit: t_nounits as t
 
-R = 1.0
-C = 1.0
-V = 1.0
-systems = @named begin
-    resistor = Resistor(R = R)
-    capacitor = Capacitor(C = C, v = 0.0)
-    source = Voltage()
-    constant = Constant(k = V)
-    ground = Ground()
+@mtkmodel RC begin
+    @parameters begin
+        R = 1.0
+        C = 1.0
+        V = 1.0
+    end
+    @components begin
+        resistor = Resistor(R = R)
+        capacitor = Capacitor(C = C, v = 0.0)
+        source = Voltage()
+        constant = Constant(k = V)
+        ground = Ground()
+    end
+    @equations begin
+        connect(constant.output, source.V)
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n, ground.g)
+    end
 end
 
-rc_eqs = [connect(constant.output, source.V)
-          connect(source.p, resistor.p)
-          connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)]
-
-@named rc_model = ODESystem(rc_eqs, t; systems)
-sys = structural_simplify(rc_model)
+@mtkbuild sys = RC()
 prob = ODEProblem(sys, Pair[], (0, 10.0))
-sol = solve(prob, Tsit5())
-plot(sol, idxs = [capacitor.v, resistor.i],
+sol = solve(prob)
+
+plot(sol, idxs = [sys.capacitor.v, sys.resistor.i],
     title = "RC Circuit Demonstration",
     labels = ["Capacitor Voltage" "Resistor Current"])
-savefig("plot.png")
 ```
 
 ![](https://user-images.githubusercontent.com/1814174/164912983-c3f73628-0e19-4e42-b085-4f62ba6f23d1.png)

docs/Project.toml

@@ -13,7 +13,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 [compat]
 ControlSystemsBase = "1.1"
 DataFrames = "1.7"
-DataInterpolations = "6.4"
+DataInterpolations = "6.4, 7"
 DifferentialEquations = "7.6"
 Documenter = "1"
 IfElse = "0.1"

docs/make.jl

@@ -31,6 +31,7 @@ makedocs(sitename = "ModelingToolkitStandardLibrary.jl",
         ModelingToolkitStandardLibrary.Hydraulic,
         ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible],
     clean = true, doctest = false, linkcheck = true,
+    linkcheck_ignore = ["https://www.mathworks.com/help/simscape/ug/basic-principles-of-modeling-physical-networks.html#bq89sba-6"],
     warnonly = [:docs_block, :missing_docs, :cross_references],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/"),

docs/src/API/blocks.md

@@ -36,6 +36,11 @@ Add
 Add3
 Product
 Division
+UnaryMinus
+Power
+Modulo
+Floor
+Ceil
 StaticNonLinearity
 Abs
 Sign

src/Blocks/Blocks.jl

@@ -10,7 +10,7 @@ using ModelingToolkit: getdefault, t_nounits as t, D_nounits as D
 export RealInput, RealInputArray, RealOutput, RealOutputArray, SISO
 include("utils.jl")
 
-export Gain, Sum, MatrixGain, Feedback, Add, Add3, Product, Division
+export Gain, Sum, MatrixGain, Feedback, Add, Add3, Product, Division, Power, Modulo, UnaryMinus, Floor, Ceil
 export Abs, Sign, Sqrt, Sin, Cos, Tan, Asin, Acos, Atan, Atan2, Sinh, Cosh, Tanh, Exp
 export Log, Log10
 include("math.jl")
@@ -29,6 +29,8 @@ include("continuous.jl")
 
 export AnalysisPoint, get_sensitivity, get_comp_sensitivity,
        get_looptransfer, open_loop
-include("analysis_points.jl")
+@static if !isdefined(ModelingToolkit, :AnalysisPoint)
+    include("analysis_points.jl")
+end
 
 end

src/Blocks/math.jl

@@ -208,6 +208,110 @@ Output first input divided by second input.
     end
 end
 
+"""
+    Power(; name)
+
+Output the exponential with base as the first input and exponent as second input i.e u1^u2
+
+# Connectors:
+
+  - `base`
+  - `exponent`
+  - `output`
+"""
+@mtkmodel Power begin
+    @components begin
+        base = RealInput()
+        exponent = RealInput()
+        output = RealOutput()
+    end
+    @equations begin
+        output.u ~ base.u ^ exponent.u
+    end
+end
+
+"""
+    Modulo(; name)
+
+Output the remainder when the first input is divided by second input.
+
+# Connectors:
+
+  - `dividend`
+  - `divisor`
+  - `remainder`
+"""
+@mtkmodel Modulo begin
+    @components begin
+        dividend = RealInput()
+        divisor = RealInput(guess = 1.0) # denominator can not be zero
+        remainder = RealOutput()
+    end
+    @equations begin
+        remainder.u ~ mod(dividend.u, divisor.u)
+    end
+end
+
+"""
+    UnaryMinus(; name)
+
+Output the product of -1 and the input.
+
+# Connectors:
+
+  - `input`
+  - `output`
+"""
+@mtkmodel UnaryMinus begin
+    @components begin
+        input = RealInput()
+        output = RealOutput()
+    end
+    @equations begin
+        output.u ~ -(input.u)
+    end
+end
+
+"""
+    Floor(; name)
+
+Output the floor rounding of the input.
+
+# Connectors:
+
+  - `input`
+  - `output`
+"""
+@mtkmodel Floor begin
+    @components begin
+        input = RealInput()
+        output = RealOutput()
+    end
+    @equations begin
+        output.u ~ floor(input.u)
+    end
+end
+
+"""
+    Ceil(; name)
+
+Output the ceiling rounding of the input.
+
+# Connectors:
+
+  - `input`
+  - `output`
+"""
+@mtkmodel Ceil begin
+    @components begin
+        input = RealInput()
+        output = RealOutput()
+    end
+    @equations begin
+        output.u ~ ceil(input.u)
+    end
+end
+
 """
     StaticNonLinearity(func; name)
 

src/Blocks/sources.jl

@@ -672,10 +672,7 @@ data input component.
     eqs = [
         output.u ~ get_sampled_data(t, p)
     ]
-    return ODESystem(eqs, t, vars, [pars; p]; name, systems,
-        defaults = [
-            output.u => get_sampled_data(0.0, p)
-        ])
+    return ODESystem(eqs, t, vars, [pars; p]; name, systems)
 end
 
 """
@@ -701,8 +698,7 @@ data input component.
     eqs = [
         output.u ~ get_sampled_data(t, buffer)
     ]
-    return ODESystem(eqs, t, vars, pars; name, systems,
-        defaults = [output.u => get_sampled_data(0.0, buffer)])
+    return ODESystem(eqs, t, vars, pars; name, systems)
 end
 
 SampledData(x::SampledDataType.Option; kwargs...) = SampledData(Val(x); kwargs...)

src/Electrical/Analog/ideal_components.jl

@@ -357,10 +357,9 @@ R = R_const + pos * R_ref * (1 + alpha * (port.T - T_ref))
 
     @parameters begin
         R_ref = 1.0,
-        [description = "Resistance at temperature T_ref when fully closed (pos=1.0)",
-            unit = "Ω"]
-        T_ref = 300.15, [description = "Reference temperature", unit = "K"]
-        R_const = 1e-3, [description = "Constant resistance between p and n", unit = "Ω"]
+        [description = "Resistance at temperature T_ref when fully closed (pos=1.0) (Ω)"]
+        T_ref = 300.15, [description = "Reference temperature (K)"]
+        R_const = 1e-3, [description = "Constant resistance between p and n (Ω)"]
     end
 
     @components begin
@@ -369,13 +368,12 @@ R = R_const + pos * R_ref * (1 + alpha * (port.T - T_ref))
 
     @variables begin
         pos(t), [description = "Position of the wiper (normally 0-1)"]
-        R(t), [description = "Resistance", unit = "Ω"]
+        R(t), [description = "Resistance (Ω)"]
     end
 
     if T_dep
         @parameters begin
-            alpha = 1e-3,
-            [description = "Temperature coefficient of resistance", unit = "K^-1"]
+            alpha = 1e-3, [description = "Temperature coefficient of resistance (K^-1)"]
         end
         @components begin
             port = HeatPort()

test/Blocks/math.jl

@@ -157,6 +157,103 @@ end
     @test sol[prod.output.u] ≈ 2 * sin.(2 * pi * sol.t)
 end
 
+@testset "Power" begin
+    @named c1 = Sine(; frequency = 1)
+    @named c2 = Constant(; k = 2)
+    @named pow = Power(;)
+    @named int = Integrator(; k = 1)
+    @named model = ODESystem(
+        [
+            connect(c1.output, pow.base),
+            connect(c2.output, pow.exponent),
+            connect(pow.output, int.input)
+        ],
+        t,
+        systems = [int, pow, c1, c2])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 1.0))
+    sol = solve(prob, Rodas4())
+    @test isequal(unbound_inputs(sys), [])
+    @test sol.retcode == Success
+    @test sol[pow.output.u] ≈ sin.(2 * pi * sol.t) .^ 2
+end
+
+@testset "Modulo" begin
+    @named c1 = Ramp(height = 2, duration = 1, offset = 1, start_time = 0, smooth = false)
+    @named c2 = Constant(; k = 1) 
+    @named modl = Modulo(;)
+    @named model = ODESystem(
+        [
+            connect(c1.output, modl.dividend),
+            connect(c2.output, modl.divisor)
+        ],
+        t,
+        systems = [modl, c1, c2])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, [], (0.0, 1.0))
+    sol = solve(prob, Rodas4())
+    @test isequal(unbound_inputs(sys), [])
+    @test sol.retcode == Success
+    @test sol[modl.remainder.u] ≈ mod.(2 * sol.t,1)
+end
+
+@testset "UnaryMinus" begin
+    @named c1 = Sine(; frequency = 1)
+    @named minu = UnaryMinus(;)
+    @named int = Integrator(; k = 1)
+    @named model = ODESystem(
+        [
+            connect(c1.output, minu.input),
+            connect(minu.output, int.input)
+        ],
+        t,
+        systems = [int, minu, c1])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 1.0))
+    sol = solve(prob, Rodas4())
+    @test isequal(unbound_inputs(sys), [])
+    @test sol.retcode == Success
+    @test sol[minu.output.u] ≈ - sin.(2 * pi * sol.t)
+end
+
+@testset "Floor" begin
+    @named c1 = Sine(; frequency = 1)
+    @named flr = Floor(;)
+    @named int = Integrator(; k = 1)
+    @named model = ODESystem(
+        [
+            connect(c1.output, flr.input),
+            connect(flr.output, int.input)
+        ],
+        t,
+        systems = [int, flr, c1])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 1.0))
+    sol = solve(prob, Rodas4())
+    @test isequal(unbound_inputs(sys), [])
+    @test sol.retcode == Success
+    @test sol[flr.output.u] ≈ floor.(sin.(2 * pi * sol.t))
+end
+
+@testset "Ceil" begin
+    @named c1 = Sine(; frequency = 1)
+    @named cel = Ceil(;)
+    @named int = Integrator(; k = 1)
+    @named model = ODESystem(
+        [
+            connect(c1.output, cel.input),
+            connect(cel.output, int.input)
+        ],
+        t,
+        systems = [int, cel, c1])
+    sys = structural_simplify(model)
+    prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 1.0))
+    sol = solve(prob, Rodas4())
+    @test isequal(unbound_inputs(sys), [])
+    @test sol.retcode == Success
+    @test sol[cel.output.u] ≈ ceil.(sin.(2 * pi * sol.t))
+end
+
 @testset "Division" begin
     @named c1 = Sine(; frequency = 1)
     @named c2 = Constant(; k = 2)