Merge pull request #222 from JuliaDynamics/hw/datainterp

use DataInterpolations in gas network example

Author: hexaeder

docs/Project.toml

@@ -3,6 +3,7 @@ name = "ND.jl docs"
 [deps]
 Bonito = "824d6782-a2ef-11e9-3a09-e5662e0c26f8"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
@@ -12,7 +13,6 @@ Electron = "a1bb12fb-d4d1-54b4-b10a-ee7951ef7ad3"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-LinearInterpolations = "b20c7882-2490-4592-a606-fbbfe9e745e8"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NetworkDynamics = "22e9dc34-2a0d-11e9-0de0-8588d035468b"
@@ -39,6 +39,7 @@ NetworkDynamicsInspector = {path = "../NetworkDynamicsInspector"}
 [compat]
 Bonito = "≥0.0.1"
 CairoMakie = "0.13.1"
+DataInterpolations = "7"
 DiffEqCallbacks = "4.2.2"
 Distributions = "0.25.117"
 Documenter = "1.8.0"
@@ -48,7 +49,6 @@ Electron = "≥0.0.1"
 GraphMakie = "0.5.13"
 Graphs = "≥0.0.1"
 LaTeXStrings = "1.4.0"
-LinearInterpolations = "0.1.4"
 Literate = "2.20.1"
 ModelingToolkit = "≥0.0.1"
 NetworkDynamics = "≥0.0.1"

docs/examples/gas_network.jl

@@ -1,13 +1,13 @@
 #=
-# Dynamic Flow in simple Gas Network
+# Dynamic Flow in Simple Gas Network
 
-This Example is based on the paper
+This example is based on the paper
 
 > Albertus J. Malan, Lukas Rausche, Felix Strehle, Sören Hohmann,
 > Port-Hamiltonian Modelling for Analysis and Control of Gas Networks,
 > IFAC-PapersOnLine, Volume 56, Issue 2, 2023, https://doi.org/10.1016/j.ifacol.2023.10.193.
 
-and tries replicate a simple simulation of flow in a 3-node gas network.
+and tries to replicate a simple simulation of flow in a 3-node gas network.
 
 This example can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md
 
@@ -19,21 +19,21 @@ using DynamicQuantities
 using ModelingToolkit: D as Dt, t as t
 using Test
 using StaticArrays
-using LinearInterpolations
+using DataInterpolations
 using OrdinaryDiffEqTsit5
 using CairoMakie
 CairoMakie.activate!(type="svg") #hide
 
 #=
 ## Node Models
 
-In this example, we use the equation based modeling using `ModelingToolkit.jl`.
-To verify the equations on a basic level we also provide units to eveything to
+In this example, we use equation-based modeling using `ModelingToolkit.jl`.
+To verify the equations on a basic level, we also provide units to everything to
 perform dimensionality checks.
 
 There are 2 node models used in the paper. The first node type has a constant
 pressure.
-Additionally, we ad some "internal" state `q̃_inj` which we want to plot later (see also [Observables](@ref)).
+Additionally, we add some "internal" state `q̃_inj` which we want to plot later (see also [Observables](@ref)).
 =#
 @mtkmodel ConstantPressureNode begin
     @parameters begin
@@ -55,12 +55,15 @@ nothing #hide
 The second node model is a variable pressure node. It has one output state (the pressure) and one input state,
 the aggregated flows from the connected pipes.
 As an internal state we have the injected flow from our source/load.
-The source/load behaviour itself is provided via a time dependent function.
+The source/load behavior itself is provided via a time-dependent function.
 =#
 @mtkmodel VariablePressureNode begin
     @structural_parameters begin
         load_profile # time dependent load profile
     end
+    @constants begin
+        load_unit = 1, [description="unit of the load profile", unit=u"m^3/s"]
+    end
     @parameters begin
         C, [description="Lumped capacitance of connected pipes", unit=u"m^4 * s^2 / kg"]
     end
@@ -70,7 +73,7 @@ The source/load behaviour itself is provided via a time dependent function.
         q̃_nw(t), [description="aggregated flow from pipes into node", unit=u"m^3/s", input=true]
     end
     @equations begin
-        q̃_inj ~ load_profile(t)
+        q̃_inj ~ load_profile(t) * load_unit
         C * Dt(p) ~ q̃_inj + q̃_nw # (30)
     end
 end
@@ -79,10 +82,10 @@ nothing #hide
 #=
 ## Pipe Model
 
-The pipe is modeld as a first order ODE for the volumetric flow at the `dst` end. It has two inputs:
-the pressure at the source and and the pressure at the destination end.
+The pipe is modeled as a first-order ODE for the volumetric flow at the `dst` end. It has two inputs:
+the pressure at the source and the pressure at the destination end.
 Later on, we'll specify the model to be antisymmetric, thus the flow is calculated explicitly for the
-destination end, but the source end will just recive just that times (-1).
+destination end, but the source end will just receive that times (-1).
 =#
 @mtkmodel Pipe begin
     @parameters begin
@@ -157,7 +160,6 @@ T  = 278u"K"           # simulation temperature
 r  = 0.012u"mm"        # pipe roughness
 D  = 0.6u"m"           # pipe diameter
 
-## TODO: here is switched the lenths l12 and l13. The results are better. Is this a mistake in the paper?
 L₁₂ = 90u"km"
 L₁₃ = 80u"km"
 L₂₃ = 100u"km"
@@ -175,7 +177,7 @@ sinθ₂₃ = 0.0
 nothing # hide
 
 #=
-Lastly, we need to calculate the compressibility factor, the speed of sound and
+Lastly, we need to calculate the compressibility factor, the speed of sound, and
 the density at standard conditions:
 =#
 Z̃ = 1 - 3.52 * p̃/pc * exp(-2.26*(T̃/Tc)) + 0.274 * (p̃/pc)^2 * exp(-1.878*(T̃/Tc)) # (5)
@@ -188,56 +190,54 @@ nothing # hide
 The equivalent "pressure capacity" at the nodes is calculated as a sum of the connected
 pipe parameters according to (28).
 
-Here use defintions based on the speed and "standard" conditions.
+Here we use definitions based on the speed and "standard" conditions.
 =#
 C₂ = L₁₂*A/(2*ρ̃*c̃^2) + L₂₃*A/(2*ρ̃*c̃^2) # (28)
 C₃ = L₁₃*A/(2*ρ̃*c̃^2) + L₂₃*A/(2*ρ̃*c̃^2) # (28)
 nothing #hide
 
 #=
-Alternatively, we could calculate `Z2` and `Z3` based on the actuel pressure and simulation temperature.
-Then we could calculated the speed of sound for the "correct" conditions at the node.
-It seems to have very little effect on the actual results so I kept it simple.
+Alternatively, we could calculate `Z2` and `Z3` based on the actual pressure and simulation temperature.
+Then we could calculate the speed of sound for the "correct" conditions at the node.
+It seems to have very little effect on the actual results, so I kept it simple.
 =#
 nothing #hide
 
 #=
 ## Load Profile
 
-The paper specifies the load profile at two nodes. We use the package [`LinearInterpolations.jl`](https://github.com/jw3126/LinearInterpolations.jl)
-to get a callable object which represents this picewise linear interpolation.
+The paper specifies the load profile at two nodes. We use the package [`DataInterpolations.jl`](https://github.com/SciML/DataInterpolations.jl)
+to get a callable object which represents this piecewise linear interpolation.
 
-However, this function is not Symbolics.jl compatible, so we need to stop Symbolics.jl/ModelingToolkit.jl
-from tracing it. To do so, we use `@register_symbolic` to declare it as a symbolic function which is treated
-as a blackbox.
+Currently, the linear interpolation does not support any units yet. To satisfy the static unit check,
+we multipy the interpolation output by a constant 1 of that unit.
 
-Additionally, we need to tell ModelingToolkit about the units of this object. This is just used
-for the static unit check during construction of the model. Later one, when we generate the
-Julia code from the symbolic reepresentation all units will be stripped.
+Note however, that the unit check is only performed at the construction of the
+model. Later on, when the nummeric code will be generated from the symbolic
+representation, all units will be stripped.
 
 !!! note "Discontinuities in RHS"
-    The picewise linear interpolated function creates discontinuities in the RHS of the system.
-    However since we know the times exactly, we can handle this by simply giving a list of explicit
+    The piecewise linear interpolated function creates discontinuities in the RHS of the system.
+    However, since we know the times exactly, we can handle this by simply giving a list of explicit
     tstops to the solve command, to make sure those are hit exactly.
 =#
-load2(t) = -Interpolate(SA[0, 4, 12, 20, 24]*3600, SA[20, 30, 10, 30, 20], extrapolate=LinearInterpolations.Constant(20))(t)
-load3(t) = -Interpolate(SA[0, 4, 12, 20, 24]*3600, SA[40, 50, 30, 50, 40], extrapolate=LinearInterpolations.Constant(40))(t)
-@register_symbolic load2(t)
-@register_symbolic load3(t)
-ModelingToolkit.get_unit(op::typeof(load2), _) = u"m^3/s"
-ModelingToolkit.get_unit(op::typeof(load3), _) = u"m^3/s"
+load2 = LinearInterpolation(-1*Float64[20, 30, 10, 30, 20], [0, 4, 12, 20, 24]*3600.0; extrapolation=ExtrapolationType.Constant)
+load3 = LinearInterpolation(-1*Float64[40, 50, 30, 50, 40], [0, 4, 12, 20, 24]*3600.0; extrapolation=ExtrapolationType.Constant)
+ModelingToolkit.get_unit(::LinearInterpolation, _ ) = 1.0 # type piracy!
 nothing #hide
 
 #=
+As a workaround we had to explicitly define `LinearInterpolations` as unitless, which is type piracy! Don't to this in any package code!
+
 ## Building the Network
 
-To bild the Network we first need to define the components. This is a two step process:
+To build the network, we first need to define the components. This is a two-step process:
 
 - first create the symbolic `ODESystem` using ModelingToolkit
 - secondly build a NetworkDynamics component model ([`VertexModel`](@ref)/[`EdgeModel`](@ref)) based on the symbolic system.
 
 In the first step we can use the keyword arguments to pass "default" values for our parameters and states.
-Those values will be automaticially transfered to the metadata of the component model the second step.
+Those values will be automatically transferred to the metadata of the component model in the second step.
 
 The second step requires to define the interface variables, i.e. what are the "input" states of your
 component model and what are the "output" states.
@@ -255,10 +255,10 @@ v2 = VertexModel(v2_mtk, [:q̃_nw], [:p]; name=:v2, vidx=2)
 v3 = VertexModel(v3_mtk, [:q̃_nw], [:p]; name=:v3, vidx=3)
 
 #=
-For the edge Model we have two inputs: the pressure on both source and destination end.
-There is a single output state: the volumetric flow. However we also need to tell
+For the edge model we have two inputs: the pressure on both source and destination end.
+There is a single output state: the volumetric flow. However, we also need to tell
 NetworkDynamics about the coupling type. In this case we use `AntiSymmetric`, which
-meas that the source end will recieve the same flow, just inverted sign.
+means that the source end will receive the same flow, just with inverted sign.
 =#
 @named e12_mtk = Pipe(; L=L₁₂, sinθ=sinθ₁₂, D, A, γ, η, r, g, T, Tc, pc, Rs, c̃, ρ̃, p̃)
 @named e13_mtk = Pipe(; L=L₁₃, sinθ=sinθ₁₃, D, A, γ, η, r, g, T, Tc, pc, Rs, c̃, ρ̃, p̃)
@@ -271,10 +271,10 @@ e23 = EdgeModel(e23_mtk, [:p_src], [:p_dst], AntiSymmetric([:q̃]); name=:e23, s
 #=
 To build the network object we just need to pass the vertices and edges to the constructor.
 
-Note that we've used the `vidx` and `src`/`dst` keywords in the constructors to define for
-each component to which "part" of the network it belongs.
+Note that we've used the `vidx` and `src`/`dst` keywords in the constructors to define
+for each component to which "part" of the network it belongs.
 
-This means, the constructor can automaticially construct a graph based on those informations and
+This means the constructor can automatically construct a graph based on that information and
 we don't need to pass it explicitly.
 =#
 
@@ -286,38 +286,28 @@ structurally different, because both functions capure a unique loadprofile funct
 
 ## Finding a Steady State
 
-To simulate the systme, we first need to find a steadystate. As a "guess" for that
+To simulate the system, we first need to find a steady state. As a "guess" for that,
 we create a `NWState` object from the network.
 This will allocate flat arrays for states `u` and parameters `p` and fill them with the
 default values.
 =#
 uguess = NWState(nw)
 
 #=
-This is not a steadystate of the system however. To find a true steadystate we want to ensure
-that the lhs of the system is zero.
-We can solve for a steady state numerically by defining a Nonlinear Rootfind problem.
-
-To do so, we need to wrap the Network object in a closure.
-=#
-nwwrap = (du, u, p) -> begin
-    nw(du, u, p, 0)
-    nothing
-end
-initprob = NonlinearProblem(nwwrap, uflat(uguess), pflat(uguess))
-initsol = solve(initprob)
+This is not a steady state of the system however. To find a true steady state, we want to ensure
+that the LHS of the system is zero.
 
-#=
-We can create a new `NWState` object by wrapping the solution from the nonlinear problem and the
-original prameters in a new `NWState` object.
+We can use the [`find_fixpoint`](@ref) from `NetworkDynamics.jl` to initialize the system.
+Internally, this uses a nummerical solve for the rootfind problem `0 = rhs`.
+The result is automaticially wrapped as a [`NWState`](@ref) object.
 =#
-u0 = NWState(nw, initsol.u, uguess.p)
+u0 = find_fixpoint(nw, uguess)
 
 #=
 ## Solving the ODE
 
 Using this as our initial state we can create the actual `ODEProblem`.
-Since the ode allways operates on flat state and aprameter arrays we use `uflat` and `pflat` to extract them.
+Since the ODE always operates on flat state and parameter arrays, we use `uflat` and `pflat` to extract them.
 =#
 prob = ODEProblem(nw, uflat(u0), (0.0,24*3600), copy(pflat(u0)))
 sol = solve(prob, Tsit5(), tstops=[0,4,12,20,24]*3600)
@@ -358,8 +348,8 @@ Notably, the "internal" states defined in the symbolic models are not "states" i
 For example, we captured the load profile in the `q̃_inj` state of the `VariablePressureNode`.
 The only dynamic state of the model however is `p`.
 Using the "observables" mechanism from SciML, which is implemented by NetworkDynamics, we can reconstruct
-those "optimized" states which have been removed symbolicially.
-Here we plot the reconstructed load profile of nodes 2 and 3. Also, we know that node 1 is infinetly stiff,
+those "optimized" states which have been removed symbolically.
+Here we plot the reconstructed load profile of nodes 2 and 3. Also, we know that node 1 is infinitely stiff,
 acting as an infinite source of volumetric flow. We can reconstruct this flow too.
 =#
 fig = begin

docs/src/mtk_integration.md

@@ -36,7 +36,7 @@ as tearing of thousands of equations.
 
 ## RC-Circuit Example
 In good [MTK tradition](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/acausal_components/), this feature will be explained along a simple RC circuit example.
-The [Dynamic Flow in simple Gas Network](@ref) example is another showcase of the MTK constructors.
+The [Dynamic Flow in Simple Gas Network](@ref) example is another showcase of the MTK constructors.
 
 The system to model is 2 node, 1 edge network. The node output states are the voltage (to ground), the edge output sates are the currents at both ends.
 ```

ext/MTKExt.jl

@@ -298,7 +298,7 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     params = setdiff(allparams, Set(allinputs))
 
     # extract the main equations and observed equations
-    eqs::Vector{Equation} = full_equations(sys)
+    eqs::Vector{Equation} = ModelingToolkit.subs_constants(full_equations(sys))
     fix_metadata!(eqs, sys);
 
     # assert the ordering of states and equations
@@ -325,7 +325,8 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     # extract observed equations. They might depend on eachother so resolve them
     obs_subs = Dict(eq.lhs => eq.rhs for eq in observed(sys))
     obseqs = map(observed(sys)) do eq
-        eq.lhs ~ fixpoint_sub(eq.rhs, obs_subs)
+        expanded_rhs = fixpoint_sub(eq.rhs, obs_subs)
+        eq.lhs ~ ModelingToolkit.subs_constants(expanded_rhs)
     end
     fix_metadata!(obseqs, sys);
     # obs can only depend on parameters (including allinputs) or states

src/utils.jl

@@ -118,20 +118,28 @@ end
 
 function rand_inputs_fg(rng, cf)
     @argcheck hasindim(cf) "ComponentModel has no specified input dimensions/syms"
-    du = rand(rng, dim(cf))
+    du = [NaN for _ in 1:dim(cf)]
     u = rand(rng, dim(cf))
     p = rand(rng, pdim(cf))
     ins = Tuple(rand(rng, l) for l in values(indim(cf)))
     if has_external_input(cf)
         ext = rand(rng, extdim(cf))
         ins = (ins..., ext)
     end
-    outs = Tuple(rand(rng, l) for l in values(outdim(cf)))
+    outs = Tuple([NaN for _ in 1:l] for l in values(outdim(cf)))
     t = NaN
-    (outs, du, u, ins, p, t)
+    (outs, du, u, ins, p, t) # fg extrancs outs and ints internally!
 end
 rand_inputs_fg(cf) = rand_inputs_fg(Random.default_rng(), cf)
 
+function rand_inputs_obsf(rng, cf)
+    (_, u, p, ins, p, t) = rand_inputs_fg(rng, cf)
+    N = length(cf.obssym)
+    outs = [NaN for _ in 1:N]
+    (outs, u, ins..., p, t)
+end
+rand_inputs_obsf(cf) = rand_inputs_obsf(Random.default_rng(), cf)
+
 
 # abstract symbolic index types
 abstract type SymbolicIndex{C,S} end

test/MTK_test.jl

@@ -257,3 +257,38 @@ v = VertexModel(fullyimplicit, [:u], [:z])
 data = NetworkDynamics.rand_inputs_fg(v)
 b = @b $(NetworkDynamics.compfg(v))($data...)
 @test b.allocs == 0
+
+@testset "Test constants in MTK models" begin
+    @mtkmodel DQSwing_Constants begin
+        @extend DQBus()
+        @variables begin
+            ω(t) = 0.0, [description = "Rotor frequency"]
+            θ(t) = 0.0, [description = "Rotor angle"]
+            Pel(t), [description = "Electrical Power"]
+        end
+        @constants begin
+            V = 1, [description = "Voltage magnitude"]
+            useless = 0
+        end
+        @parameters begin
+            M = 1, [description = "Inertia"]
+            D = 0.1, [description = "Damping"]
+            Pmech, [description = "Mechanical Power"]
+        end
+        @equations begin
+            Dt(θ) ~ ω
+            Dt(ω) ~ 1/M * (Pmech - D*ω + Pel)
+            Pel ~ real((u_r + im*u_i) * (i_r - im*i_i)) + useless
+            u_r ~ V*cos(θ)
+            u_i ~ V*sin(θ)
+        end
+    end
+    @named withconst = DQSwing_Constants()
+    v = VertexModel(withconst, [:i_r, :i_i], [:u_r, :u_i])
+
+    data = NetworkDynamics.rand_inputs_fg(v)
+    NetworkDynamics.compfg(v)(data...) # no error
+
+    data = NetworkDynamics.rand_inputs_obsf(v)
+    v.obsf(data...) # no error
+end