Merge pull request #303 from JuliaDynamics/hw/mtk10take2

bump MTK to version 10, take 2

Author: hexaeder

.github/workflows/tests.yml

@@ -161,7 +161,6 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/julia-buildpkg@v1
       - name: Load Julia packages from cache
         id: julia-cache
         uses: julia-actions/cache@v2
@@ -202,6 +201,9 @@ jobs:
             using Pkg
             using Test
             pkg"dev ."
+            # pkg"update"
+            pkg"instantiate"
+            pkg"build"
             @testset "PowerDynamics Downstream Tests" begin
                 @testset "Normal Package tests" begin
                     Pkg.test("PowerDynamics"; coverage=false)

NEWS.md

@@ -1,8 +1,16 @@
 # NetworkDynamics Release Notes
 
+## v0.10.4 Changelog
+- [#303](https://github.com/JuliaDynamics/NetworkDynamics.jl/pull/303) update for ModelingToolkit.jl v10 compatibility:
+  - rename all `ODESystem` -> `System` (follows MTK v10 API)
+  - MTK extension now uses `mtkcompile` instead of `structural_simplify` internally
+  - Add new `implicit_output` function to handle fully implicit output variables in MTK models
+  - Add documentation for handling fully implicit outputs in MTK integration
+  - Update minimum ModelingToolkit.jl requirement from v9.67 to v10
+
 ## v0.10.3 Changelog
 - [#301](https://github.com/JuliaDynamics/NetworkDynamics.jl/pull/301) improve callback system performance and flexibility:
-  - Add callback batching for better DiscreteComponentCallback performance 
+  - Add callback batching for better DiscreteComponentCallback performance
   - Allow `EIndex(1=>2)` as standalone edge index with relaxed type constraints
   - Optimize CallbackSet construction to prevent performance bottlenecks
   - Add important documentation warning about parameter array copying in callbacks

Project.toml

@@ -1,7 +1,7 @@
 name = "NetworkDynamics"
 uuid = "22e9dc34-2a0d-11e9-0de0-8588d035468b"
 authors = ["Frank Hellmann <[email protected]>, Michael Lindner <[email protected]>, Hans Würfel <[email protected]"]
-version = "0.10.3"
+version = "0.10.4"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
@@ -67,7 +67,7 @@ KernelAbstractions = "0.9.18"
 LinearAlgebra = "1"
 MacroTools = "0.5.15"
 Mixers = "0.1.2"
-ModelingToolkit = "9.82"
+ModelingToolkit = "10"
 NNlib = "0.9.13"
 NonlinearSolve = "4"
 OrderedCollections = "1.8.0"

docs/examples/gas_network.jl

@@ -233,7 +233,7 @@ As a workaround we had to explicitly define `LinearInterpolations` as unitless,
 
 To build the network, we first need to define the components. This is a two-step process:
 
-- first create the symbolic `ODESystem` using ModelingToolkit
+- first create the symbolic `System` 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.

docs/examples/init_tutorial.jl

@@ -61,15 +61,19 @@ end
 nothing #hide
 
 #=
-**B) A static prosumer node which forces a certain flow (pressure is fully implicit)**
+**B) A static prosumer node which forces a certain flow**
+
+!!! note "Fully Implicit Output"
+    We need to use `implicit_output(p)` to handle the fully implicit pressure
+    output. See [fully implicit outputs](@ref Fully-Implicit-Outputs) for details.
 =#
 @mtkmodel StaticProsumerNode begin
     @extend GasNode()
     @parameters begin
         q̃_prosumer, [description="flow injected by prosumer"]
     end
     @equations begin
-        -q̃_nw ~ q̃_prosumer
+        -q̃_nw ~ q̃_prosumer + implicit_output(p)
     end
 end
 nothing #hide

docs/src/API.md

@@ -18,9 +18,9 @@ EdgeModel()
 
 ## Component Models with MTK
 ```@docs
-VertexModel(::ModelingToolkit.ODESystem, ::Any, ::Any)
-EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::Any, ::Any)
-EdgeModel(::ModelingToolkit.ODESystem, ::Any, ::Any, ::Any)
+VertexModel(::ModelingToolkit.System, ::Any, ::Any)
+EdgeModel(::ModelingToolkit.System, ::Any, ::Any, ::Any, ::Any)
+EdgeModel(::ModelingToolkit.System, ::Any, ::Any, ::Any)
 ```
 
 ### Output Function Helpers/Wrappers
@@ -245,6 +245,7 @@ save_parameters!
 ff_to_constraint
 Base.copy(::NetworkDynamics.ComponentModel)
 extract_nw
+implicit_output
 ```
 
 ## NetworkDynamicsInspector API

docs/src/metadata.md

@@ -33,7 +33,7 @@ Special cases for symbol metadata are:
 
 For those, there are special functions `has_*`, `get_*`, `set_*!`, `delete_*!` and `strip_*!`. The `strip_*!` functions remove all metadata of a specific type from all symbols in a component. See [Per Symbol Metadata API](@ref).
 
-These are closely aligned with the [metadata use in ModelingToolkit](@extref ModelingToolkit symbolic_metadata). They are automatically copied from the `ODESystem` if you use MTK models to create NetworkDynamics models.
+These are closely aligned with the [metadata use in ModelingToolkit](@extref ModelingToolkit symbolic_metadata). They are automatically copied from the `System` if you use MTK models to create NetworkDynamics models.
 
 ## Metadata Utils
 Accessing metadata (especially defaults) of states and parameters is a very

docs/src/mtk_integration.md

@@ -6,13 +6,13 @@ which are then connected on network scale using NetworkDynamics.
 
 The main entry point for this interop are the constructors
 ```julia
-VertexModel(::ODESystem, inputs, outputs)
-EdgeModel(::ODESystem, srcin, dstin, [srscout], dstout)
+VertexModel(::System, inputs, outputs)
+EdgeModel(::System, srcin, dstin, [srscout], dstout)
 ```
 whose docstrings can be found in the [Component Models with MTK](@ref) section in the API.
 
 These constructors will:
-- transform the states marked as input to parameters and `structural_simplify`ing the system,
+- transform the states marked as input to parameters and `mtkcompile`ing the system,
 - generate the `f` and `g` functions,
 - generate code for observables,
 - port all supported [Metadata](@ref) from MTK symbols to component symbols and
@@ -43,7 +43,7 @@ The system to model is 2 node, 1 edge network. The node output states are the vo
 
 ideal v source      Resistor     Capacitor
             v1 o─←────MMM────→─o v2
-               │               ┴ 
+               │               ┴
               (↗)              ┬
                │               │
                ⏚               ⏚
@@ -76,7 +76,7 @@ An ideal voltage source is just a model which pins its output voltage to a fixed
 The source ejects whatever current is necessary. We introduce another variable `i(t)`
 to "capture" this current. This variable will be removed during structural simplify, but will
 be available for plotting through the [Observables](@ref) mechanism.
-The `VertexModel` can be generated from an `ODESystem` by providing names of the input and output states:
+The `VertexModel` can be generated from an `System` by providing names of the input and output states:
 
 ```@example mtk
 @mtkmodel VoltageSource begin
@@ -140,7 +140,7 @@ resistor_edge = EdgeModel(resistor, [:src₊v], [:dst₊v], [:src₊i], [:dst₊
 
 Having all those components defined, we can build the network. We don't need to provide a graph
 object here because we specified the placement in the graph on a per component basis.
- 
+
 ```@example mtk
 nw = Network([vs_vertex, cap_vertex], [resistor_edge])
 ```
@@ -149,7 +149,7 @@ We can see, that NetworkDynamics internally is able to reduce all of the "output
 
 Now we can simulate the system. For that we generate the `u0` object. Since the metadata (such as default values) was automatically transferred, we can straight away construct the `ODEProblem`
 and solve the system.
- 
+
 ```@example mtk
 u0 = NWState(nw) # generate state based on default values
 prob = ODEProblem(nw, uflat(u0), (0, 10.0), pflat(u0))
@@ -164,3 +164,67 @@ axislegend(ax2)
 fig # hide
 ```
 
+## Fully Implicit Outputs
+When working with MTK systems in NetworkDynamics, you may encounter situations where
+your desired output variables don't explicitly appear in the equations. This creates **fully
+implicit outputs** - variables that are determined by the system's constraints but aren't
+directly computed.
+
+!!! tip "tl;dr"
+    Introduce "fake" dependencies to your input-forcing equations `0 ~ in + implicit_output(y)`.
+    Which is mathematically equivalent to `0 ~ in` but helps MTK to reason about dependencies.
+
+Consider a system with a fully implicit output:
+```
+   u┌───────┐y
+  ─→┤ 0 ~ u ├→─
+    └───────┘
+```
+Here, $y$ does not appear in the equations at all. In general, that doesn't make too much sense.
+During simplification, MTK will potentially get rid of the equation as it does not contribute to the system's state.
+
+However, in NetworkDynamics, we're always dealing with **open loop models** on the equation level, which is not exactly what MTK was made for.
+If you build a closed loop between a subsystem A which **has input forcing** and a subsystem
+B which has **input feed forward**, the resulting system can be solved:
+```
+    (system with input forcing)
+          ua┌─────────┐ya
+        ┌──→┤  0 ~ ua ├→──┐
+        │   └─────────┘   │
+        │ yb┌─────────┐ub │
+        └──←┤ yb ~ ub ├←──┘
+            └─────────┘
+(system with input feed forward)
+```
+
+Since MTK does not know about the closed loop (which is only introduced on the NetworkDynamics level once we leave the equation based domain) we need to help MTK to figure out those dependencies.
+We can do so by introducing "fake" dependencies using [`implicit_output`](@ref).
+This function is defined as
+```julia
+implicit_output(x) = 0
+ModelingToolkit.@register_symbolic implicit_output(x)
+```
+which makes it numerically equivalent to zero (no effect on the simulation) but is
+opaque to the Symbolic Simplification.
+
+### Example
+
+Consider a "Kirchhoff Node" between multiple resistors:
+- the currents through the resistors directly depend on the voltage output of the node (input feed forward) and
+- the Kirchhoff node requires the sum of all inflowing currents to be zero (input forcing).
+
+We can model this type of node like this:
+```@example mtk
+@mtkmodel KirchhoffNode begin
+    @variables begin
+        v(t), [description="Node voltage", output=true]
+        i_sum(t), [description="Sum of incoming currents", input=true]
+    end
+    @equations begin
+        0 ~ i_sum + implicit_output(v)  # Kirchhoff's current law
+    end
+end
+@named kirchhoff = KirchhoffNode()
+VertexModel(kirchhoff, [:i_sum], [:v])
+```
+where we "trick" MTK into believing that the input forcing equation depends on the output too.

ext/MTKExt_utils.jl

@@ -152,7 +152,7 @@ function warn_missing_features(sys)
     cev = ModelingToolkit.get_continuous_events(sys)
     dev = ModelingToolkit.get_discrete_events(sys)
     if !isempty(cev) || !isempty(dev)
-        @warn "Model has attached events, which is not supportet."
+        @warn "Model has attached events, which is not supported."
     end
 
     if !isempty(ModelingToolkit.initialization_equations(sys))
@@ -215,3 +215,16 @@ function fix_metadata!(invalid_eqs, sys)
     end
     invalid_eqs .= fixedeqs
 end
+
+function remove_implicit_output_fn!(eqs)
+    r = SymbolicUtils.@rule implicit_output(~~x) => 0
+    chain = SymbolicUtils.Chain([r])
+    rewriter = SymbolicUtils.Prewalk(chain)
+
+    for i in eachindex(eqs)
+        eq = eqs[i]
+        eqs[i] = rewriter(eq.lhs) ~ rewriter(eq.rhs)
+    end
+
+    eqs
+end