|
| 1 | +# ModelingToolkit Integration |
| 2 | + |
| 3 | +NetworkDynamics.jl is compatible with [`ModelingTookit.jl`](https://github.com/SciML/ModelingToolkit.jl) (MTK). |
| 4 | +The general idea is to use MTK to define *component models* (i.e. edge and vertex dynamics) |
| 5 | +which are then connected on network scale using NetworkDynamics. |
| 6 | + |
| 7 | +The main entry point for this interop are the constructors |
| 8 | +```julia |
| 9 | +VertexFunction(::ODESystem, inputs, outputs) |
| 10 | +EdgeFunction(::ODESystem, srcin, dstin, [srscout], dstout) |
| 11 | +``` |
| 12 | +whose docstrings can be found in the [Component Functions with MTK](@ref) section in the API. |
| 13 | + |
| 14 | +These constructors will: |
| 15 | +- transforming the states marked as input to parameters and `structural_simplify`ing the system, |
| 16 | +- generating the `f` and `g` functions, |
| 17 | +- generate code for observables, |
| 18 | +- port all supported [Metadata](@ref) from MTK symbols to component symbols and |
| 19 | +- output a `Vertex-`/`EdgeFunction` function compatible with NetworkDynamics.jl. |
| 20 | + |
| 21 | +The main usecase for this feature is when you want to build relatively complex |
| 22 | +component models but interconnect them in a very homogeneous way (i.e. having the |
| 23 | +same output/input pairings in the whole system). |
| 24 | + |
| 25 | +In theory, you can achieve everything you want to do with plain MTK. The idea of |
| 26 | +combining the two is, that NetworkDynamics offers *far less flexibility* when in |
| 27 | +comes to interconnection of subsystems on the network level. This might allow ND |
| 28 | +to exploit more knowledge of the structure without very expensive operations such |
| 29 | +as tearing of thousands of equations. |
| 30 | + |
| 31 | +!!! warning |
| 32 | + ModelingToolkit is a fast paced library with lots of functionality and ever |
| 33 | + growing complexity. As such the provided interface is kinda experimental. |
| 34 | + Some features of MTK are straight up unsupported, for example events within |
| 35 | + models or delay differential equations. |
| 36 | + |
| 37 | +## RC-Circuit Example |
| 38 | +In good [MTK tradition](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/acausal_components/), this feature will be explained along a simple RC circuit example. |
| 39 | +The [Dynamic Flow in simple Gas Network](@ref) example is another showcase of the MTK constructors. |
| 40 | + |
| 41 | +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. |
| 42 | +``` |
| 43 | +
|
| 44 | +ideal v source Resistor Capacitor |
| 45 | + v1 o─←────MMM────→─o v2 |
| 46 | + │ ┴ |
| 47 | + (↗) ┬ |
| 48 | + │ │ |
| 49 | + ⏚ ⏚ |
| 50 | +``` |
| 51 | +Obviously there is no need in modeling such a small system using NetworkDynamics, however |
| 52 | +the method extends quite easily to construct large electrical networks reusing the same |
| 53 | +fundamental building blocks. |
| 54 | + |
| 55 | +```@example mtk |
| 56 | +using NetworkDynamics |
| 57 | +using ModelingToolkit |
| 58 | +using ModelingToolkit: t_nounits as t, D_nounits as D |
| 59 | +using OrdinaryDiffEqTsit5 |
| 60 | +using CairoMakie |
| 61 | +``` |
| 62 | + |
| 63 | +All our components have "terminals", which have a voltage and current. We don't use the `@connector` from MTK here because our pins mark the interface towards the network and do not follow the MTK connector semantics. |
| 64 | + |
| 65 | +```@example mtk |
| 66 | +@mtkmodel NWTerminal begin |
| 67 | + @variables begin |
| 68 | + v(t), [description="Voltage at node"] |
| 69 | + i(t), [description="Current flowing into node"] |
| 70 | + end |
| 71 | +end |
| 72 | +nothing #hide |
| 73 | +``` |
| 74 | + |
| 75 | +An ideal voltage source is just a model which pins its output voltage to a fixed parameter. |
| 76 | +The source ejects whatever current is necessary. We introduce another variable `i(t)` |
| 77 | +to "capture" this current. This variable will be removed during structural simplify, but will |
| 78 | +be available for plotting through the [Observables](@ref) mechanism. |
| 79 | +The `VertexFunction` can be generated from an `ODESystem` by providing names of the input and output states: |
| 80 | + |
| 81 | +```@example mtk |
| 82 | +@mtkmodel VoltageSource begin |
| 83 | + @components begin |
| 84 | + p = NWTerminal() |
| 85 | + end |
| 86 | + @parameters begin |
| 87 | + V = 1.0 |
| 88 | + end |
| 89 | + @variables begin |
| 90 | + i(t), [description="produced current by ideal voltage source (observable)"] |
| 91 | + end |
| 92 | + @equations begin |
| 93 | + i ~ -p.i |
| 94 | + p.v ~ V |
| 95 | + end |
| 96 | +end |
| 97 | +@named vs = VoltageSource() |
| 98 | +vs_vertex = VertexFunction(vs, [:p₊i], [:p₊v]; vidx=1) |
| 99 | +``` |
| 100 | + |
| 101 | +A capacitor is a slightly more complicated model. Its voltage is defined as an differential |
| 102 | +equation based on the inflowing current. |
| 103 | + |
| 104 | +```@example mtk |
| 105 | +@mtkmodel Capacitor begin |
| 106 | + @components begin |
| 107 | + p = NWTerminal(;v=0) |
| 108 | + end |
| 109 | + @parameters begin |
| 110 | + C = 1.0 |
| 111 | + end |
| 112 | + @equations begin |
| 113 | + D(p.v) ~ p.i / C |
| 114 | + end |
| 115 | +end |
| 116 | +@named cap = Capacitor() |
| 117 | +cap_vertex = VertexFunction(cap, [:p₊i], [:p₊v], vidx=2) |
| 118 | +``` |
| 119 | + |
| 120 | +For the resistor we need two pins, one for the `src` and one for the `dst` side. |
| 121 | +The equations are straight forward. |
| 122 | + |
| 123 | +```@example mtk |
| 124 | +@mtkmodel Resistor begin |
| 125 | + @components begin |
| 126 | + src = NWTerminal() |
| 127 | + dst = NWTerminal() |
| 128 | + end |
| 129 | + @parameters begin |
| 130 | + R = 1.0 |
| 131 | + end |
| 132 | + @equations begin |
| 133 | + dst.i ~ (src.v - dst.v)/R |
| 134 | + src.i ~ -dst.i |
| 135 | + end |
| 136 | +end |
| 137 | +@named resistor = Resistor() |
| 138 | +resistor_edge = EdgeFunction(resistor, [:src₊v], [:dst₊v], [:src₊i], [:dst₊i]; src=:vs, dst=:cap) |
| 139 | +``` |
| 140 | + |
| 141 | +Having all those components defined, we can build the network. We don't need to provide a graph |
| 142 | +object here because we specified the placement in the graph on a per component basis. |
| 143 | + |
| 144 | +```@example mtk |
| 145 | +nw = Network([vs_vertex, cap_vertex], [resistor_edge]) |
| 146 | +``` |
| 147 | + |
| 148 | +We can see, that NetworkDynamics internally is able to reduce all of the "output" states. We end up with a plain ODE of a single state. |
| 149 | + |
| 150 | +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` |
| 151 | +and solve the system. |
| 152 | + |
| 153 | +```@example mtk |
| 154 | +u0 = NWState(nw) # generate state based on default values |
| 155 | +prob = ODEProblem(nw, uflat(u0), (0, 10.0), pflat(u0)) |
| 156 | +sol = solve(prob, Tsit5()) |
| 157 | +
|
| 158 | +# plot the solution |
| 159 | +fig, ax1, p = plot(sol; idxs=VIndex(:cap, :p₊v), label="Capacitor Voltage"); |
| 160 | +axislegend(ax1) |
| 161 | +ax2 = Axis(fig[2,1]) |
| 162 | +plot!(ax2, sol; idxs=VIndex(:vs, :i), label="Current produced by ideal v source") |
| 163 | +axislegend(ax2) |
| 164 | +fig # hide |
| 165 | +``` |
| 166 | + |
0 commit comments