Network Dynamics without NetworkDynamics.jl

In this example we model a Kuramoto system on a complex network.

using Graphs
using BlockSystems
using ModelingToolkit
using OrdinaryDiffEq
using Plots

Dynamical systems on complex networks are best modelled on directed graphs, because in general an interaction of node i and j may be asymmetric, e.g. when node i influences node j but not vice versa.

Symmetries of the system such as (anti-)symmetric coupling functions may justify other approaches, but those are of limited scope, for a detailed discussion see our paper, Section II.D.

In the following each edge is represented by a function:

     e₁₂ = f(1,2)
     .--->---.
   (1)       (2)
     ˙---<---˙
     e₂₁ = f(2,1)

Each edge function sees the values of the connected nodes, its source and its destination. Nodes are modeled as functions as well and see all their incoming edges.

The goal is to generate IOBlocks for edges and vertices based on a given graph.

While defining the edge block is straightforward, we have to make the vertex blocks a bit special: since MTK does not support vector inputs yet we need a special IOBlock which depends on the number of connected edges.

function gen_edge_block(name)
    @parameters t src(t) dst(t) K
    @variables out(t)
    IOBlock([out ~ K*sin(src-dst)], [src, dst], [out], name=Symbol(name))
end

# append subscript, `sybscript(:foo, 2) => :foo₂`
subscript(s, i) = Symbol(s, Char(0x02080 + i))

function gen_vertex_block(n_edges, name)
    @parameters t ω
    @variables ϕ(t)
    D = Differential(t)
    # the way array variables work changed. This is a hack to retrieve the old behavior of
    # this closely mimics the old @parameters edge[1:n_edges](t)
    edge = Num[]
    for i in 1:n_edges
        symname = subscript(:edge, i)
        append!(edge, @parameters $symname(t))
    end

    IOBlock([D(ϕ) ~ ω + (+)(edge...)],
            [edge...],
            [ϕ],
            name=Symbol(name))
end

For simplicity our graph will be a simple ring network specified with Graphs.jl

N = 8
g = SimpleDiGraph(watts_strogatz(N,2,0)) # ring network

First we generate a list of all edge-blocks because they don't depend on the vertices.

edgelist = [(i=i, src=e.src, dst=e.dst, block=gen_edge_block("e_$(e.src)_$(e.dst)"))
            for (i, e) in enumerate(edges(g))]
edge_blocks = [e.block for e in edgelist]

Now we can generate vertex blocks based on their number of incoming edges. We will also create the connections

e_i_k.out => node.edge₁
e_j_k.out => node.edge₂
...

for all edges that go from vertices i or j to vertex k.

vert_blocks = IOBlock[]
connections = Pair[]

for i in vertices(g)
    # collect the incoming edges for each node
    edges = filter(e -> e.dst == i, edgelist)

    node = gen_vertex_block(length(edges), "v$i")
    push!(vert_blocks, node)

    # each node has the open inputs edge₁, edge₂, ...
    # we need to connect the outputs of the edge-blocks to the
    # inputs of the node like edge_j_to_1.out => node.edge₁
    for (i, edge) in enumerate(edges)
        node_input_i = getproperty(node, subscript(:edge, i))
        push!(connections, edge.block.out => node_input_i)
    end
end

Once the vertices are generated we can plug the edges' src and dst to the output of the # corresponding vertex, in this case its the oscillators angle ϕ

for edge in edgelist
    push!(connections, vert_blocks[edge.src].ϕ => edge.block.src)
    push!(connections, vert_blocks[edge.dst].ϕ => edge.block.dst)
end

We want the connect_system to get rid of the algebraic states for the edges. Therefore we have to provide a list of outputs which only contains the outputs of the vertices. By doing so the edge outputs will become internal istates of the IOSystem and upon connection may be reduced.

outputs = [block.ϕ for block in vert_blocks]

network = IOSystem(connections, vcat(vert_blocks, edge_blocks), outputs=outputs)

networkblock = connect_system(network, verbose=false)

As the output shows the system has be reduced to just N equations. Well now we can generate the functions...

gen = generate_io_function(networkblock,
                           f_states=[v.ϕ for v in vert_blocks],
                           f_params=vcat([v.ω for v in vert_blocks],
                                         [e.K for e in edge_blocks]),
                           warn=false);

... enclose the f_ip to get rid of the empty inputs field...

odefun(du, u, p, t) = gen.f_ip(du, u, (), p, t)

... set the starting conditions ...

ω = collect(1:N)./N
ω .-= sum(ω)/N
K = [3.0 for i in edge_blocks]
p = (ω..., K...)

x0 = collect(1:N)./N
x0 .-= sum(x0)./N

... and solve the system!

tspan = (0., 10.)
prob = ODEProblem(odefun, x0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol, ylabel="ϕ")

This page was generated using Literate.jl.