Network Construction

Building a Network

The main type of NetworkDynamics.jl is a Network. A network bundles various component models (edge and vertex models) together with a graph to form a callable object which represents the right hand side (RHS) of the overall dynamical system, see Mathematical Model.

A Network is build by passing a graph g, vertex models vertexm and edge models edgem to the Network constructor:.

nw = Network(g, vertexm, edgem; kwargs...)

Two important keywords for the Network constructor are:

execution: Defines the ExecutionStyle of the coreloop, e.g. SequentialExecution{true}(). A execution style is a special Julia object, which tells the backend how to parallelize (e.g. ThreadedExecution{true}() will use native Julia threads to parallelize the RHS call). A list of available executions styles can be found under Execution Types in the API.
aggregator: Instructs the backend how to perform the aggregation and which aggregation function to use. Aggregation is the process of creating a single vertex input by reducing over the outputs of adjecent edges of said vertex. The aggregator contains both the function and the algorithm. E.g. SequentialAggregator(+) is a sequential aggregation by summation. A list of availabe Aggregators can be found under Aggregators in the API.

Graphless Constructor

If each of the network components has a "graphelement" metadata, we may omit the explicit graph.

nw = Network(vertexm, edgem)

The graphelement metadata can be set using the following syntax:

VertexModel(; ..., vidx=1)         # places vertex at position 1
EdgeModel(; ..., src=1, dst=2)     # places edge between 1 and 2
EdgeModel(; ..., src=:v1, dst=:v2) # places edge between vertices with names `:v1` and `:v2`

Parallel Edges

Several edge models may connect the same pair of vertices, for example parallel lines between two buses. Each of them is an independent edge with its own states, parameters and outputs.

nw = Network(vertexm, [EdgeModel(; ..., src=1, dst=2, name=:line_a),
                       EdgeModel(; ..., src=1, dst=2, name=:line_b)])

In that case the graphless constructor builds a NetworkDynamics.ComponentGraph instead of a SimpleGraph or SimpleDiGraph, and the edge index is simply the position in the provided list. Without parallel edges, the edges are sorted into the order of edges(g) of a SimpleGraph/SimpleDiGraph as before. Pass legacy_graph=false to always get the input order, which keeps indices stable if parallel edges are added later.

Since EIndex(1 => 2) would be ambiguous for parallel edges it throws an error there; address such edges by their integer index or by a unique name instead.

Building `VertexModel`s

This chapter will walk you through the most important aspects of defining a custom vertex model. For a list of all keyword arguments please check out the docstring of VertexModel.

As an example, we'll construct an second order kuramoto model.

function kuramoto_f!(dv, v, esum, p, t)
    M, P, D = p
    dv[1] = v[2]
    dv[2] = (P - D*v[2] + esum[1])/M
    nothing
end
function kuramoto_g!(y, v, esum, p, t)
    y[1] = v[1]
    nothing
end
VertexModel(; f=kuramoto_f!, g=kuramoto_g!, dim=2, pdim=3, outdim=1)

VertexModel :VertexM FeedForward()
 ├─ 2 states: [v₁, v₂]
 ├─ 1 output: [o]
 └─ 3 params: [p₁, p₂, p₃]

Those keywords are the minimum metadata we need to provide.

However there is a problem: the vertex is classified as a FeedForward vertex, which is unnecessary. We can improve the implementation of g according to the Feed Forward Behavior section.

function kuramoto_g_noff!(y, v, p, t)
    y[1] = v[1]
    nothing
end
VertexModel(; f=kuramoto_f!, g=kuramoto_g_noff!, dim=2, pdim=3, outdim=1)

VertexModel :VertexM NoFeedForward()
 ├─ 2 states: [v₁, v₂]
 ├─ 1 output: [o]
 └─ 3 params: [p₁, p₂, p₃]

To simplify your programming and avoid explicitly writing the above trivial output function you can use StateMask. By writing

VertexModel(; f=kuramoto_f!, g=StateMask(1:1), dim=2, pdim=3)

VertexModel :VertexM PureStateMap()
 ├─ 2 states: [v₁, v₂]
 ├─ 1 output: [v₁]
 └─ 3 params: [p₁, p₂, p₃]

we are instructing the vertex model, that the output is part of the states x[1:1]. This results in the following changes:

outdim is removed because it can be inferred from StateMask
outsym is not a generic :o any more but inferred from the state symbols.

We can be even less verbose by writing g=1:1 or just g=1.

Lastly, we define improved names for our states and parameters as well as assigning a position in the graph to enable the graphless network construction. Whenever you provide a sym keyword the corresponding dim keyword stops being neccessary. So, we end up with a relatively short definition

VertexModel(; f=kuramoto_f!, g=1,
              sym=[:θ, :ω], psym=[:M=>1, :P=>0.1, :D=>0],
              insym=[:P_nw], name=:swing, vidx=1)

VertexModel :swing PureStateMap() @ Vertex 1
 ├─ 1 input:  [P_nw]
 ├─ 2 states: [θ, ω]
 ├─ 1 output: [θ]
 └─ 3 params: [M=1, P=0.1, D=0]

Building `EdgeModel`s

This chapter walks you through the most important aspects when defining custom edge models. For a list of all keyword arguments please check the docstring of EdgeModel.

As an example edge model we define a standard sinusoidal coupling between the vertices in our network. The full definition is:

function edge_f!(de, e, vsrc, vdst, p, t)
    nothing
end
function edge_g!(ysrc, ydst, e, vsrc, vdst, p, t)
    ydst[1] = p[1] * sin(vsrc[1] - vdst[1])
    ysrc[1] = -ydst[1]
end
EdgeModel(; f=edge_f!, g=edge_g!, dim=0, pdim=1, outdim=1)

EdgeModel :StaticEdgeM FeedForward()
 ├─   0 states:  []  
 ├─ 1/1 outputs: src=[src₊o] dst=[dst₊o]
 └─   1 param:   [p]

This is a purely "static" edge without internal states. This means we can omit f and dim entirely. Also, we can define a variant of g without the e input

function edge_g_ff!(ysrc, ydst, vsrc, vdst, p, t)
    ydst[1] = p[1] * sin(vsrc[1] - vdst[1])
    ysrc[1] = -ydst[1]
end
EdgeModel(;g=edge_g_ff!, pdim=1, outdim=1)

EdgeModel :StaticEdgeM PureFeedForward()
 ├─   0 states:  []  
 ├─ 1/1 outputs: src=[src₊o] dst=[dst₊o]
 └─   1 param:   [p]

which classifies as a PureFeedForward edge. In cases like this, where the edge is actually anti-symmetrical we can define a single sided output function and wrap it in an AntiSymmetric object:

function edge_g_s!(ydst, vsrc, vdst, p, t)
    ydst[1] = p[1] * sin(vsrc[1] - vdst[1])
end
EdgeModel(;g=AntiSymmetric(edge_g_ff!), pdim=1, outdim=1)

EdgeModel :StaticEdgeM FeedForward()
 ├─   0 states:  []  
 ├─ 1/1 outputs: src=[₋o] dst=[o]
 └─   1 param:   [p]

This can also lead to briefer output naming. Available single sided wrappers are:

Directed (no coupling at src),
AntiSymmetric (same coupling at src and dst),
Symmetric (inverse coupling at dst) and
Fiducial (define separate g for both ends).

Once again we can add additonal data like defining a src and dst index

function edge_g_s!(ydst, vsrc, vdst, p, t)
    ydst[1] = p[1] * sin(vsrc[1] - vdst[1])
end
EdgeModel(;g=AntiSymmetric(edge_g_ff!), psym=:K=>1, outsym=:P, insym=:θ, src=1, dst=4)

EdgeModel :StaticEdgeM FeedForward() @ Edge 1=>4
 ├─ 1/1 inputs:  src=[src₊θ] dst=[dst₊θ]
 ├─   0 states:  []  
 ├─ 1/1 outputs: src=[₋P] dst=[P]
 └─   1 param:   [K=1]