Open jleugeri opened 4 years ago
Thanks for your request! We have not forgotten about it. :wink:
One of our students has re-implemented a simple example for overload line failures in a five node power network as described in Schäfer et al. (2018) using NetworkDynamics
. The basic idea of this example is to set down the capacity parameter of an edge to zero if the flow on it exceeds a critical value. This is achieved by using a VectorContinuousCallback
. It is possible to interact with the state of the system by calling the network rhs with the signature (u, p, t, GetGD) as shown in the example. This will return a GraphData object, that allows for reading and manipulating the underlying data. It should be also straight forward to affect other nodes/edges than the one that triggered the event.
I'm not entirely sure if this approach solves your problem. We are currently working on a more simple example and a tutorial for event handling with NetworkDynamics
and eventually will also build an interface to make it more convenient for the user.
using NetworkDynamics
using LightGraphs
using DifferentialEquations
using Plots
struct SwingVertex
P::Float64 # power of node
I::Float64 # inertia of node
γ::Float64 # damping of node
end
function (sv::SwingVertex)(dv, v, e_s, e_d, p, t)
# v[1] -> δ, dv[1] = dγ = ω = v[2]
# v[2] -> ω, dv[2] -> dω
dv[1] = v[2]
dv[2] = sv.P - sv.γ*v[2] + flow_sum(e_s, e_d)
dv[2] = dv[2]/sv.I
nothing
end
function flow_sum(e_s, e_d)
sum = 0.0
for e in e_d
sum -= e[1]
end
for e in e_s
sum += e[1]
end
return sum
end
mutable struct SwingEdge
from::Int8
to::Int8
K::Float64
max_flow::Float64
end
function (se::SwingEdge)(e, v_s, v_d, p, t)
e[1] = se.K * sin(v_d[1] - v_s[1])
end
"""
Function returns callable object, which can be used as 'affect!' of a
callback. The affect! will set the coupling of edge idx to zero.
"""
function kill_line_affect(idx)
function affect!(integrator)
@info "Line $idx killed at t=$(integrator.t)."
integrator.f.f.edges![idx].f!.K = 0
u_modified!(integrator, false)
end
return affect!
end
"""
Function returns callable object, which can be used as 'affect!' of a
callback. The affect! will set the coupling of edge (source, destination)
in network to zero.
"""
function kill_line_affect(network, source, destination)
gs = network.f.graph_structure
s_ind = findall(x->x==source, gs.s_e)
d_ind = findall(x->x==destination, gs.d_e)
ind = intersect(s_ind, d_ind)
@assert length(ind) == 1
return kill_line_affect(ind[1])
end
"""
Function returns a VectorContinoursCallback which compares the max_flow of all
edges in the network an shuts them down once the flow is to high.
"""
function watch_line_limit_callback(network)
num_e = network.f.graph_structure.num_e
edges = network.f.edges!
function condition(out, u, t, integrator)
# get current edge values from integrator
graph_data = integrator.f.f(u, integrator.p, integrator.t, GetGD)
edge_values = graph_data.e_array
for i in 1:num_e
out[i] = edges[i].f!.max_flow - abs(edge_values[i])
end
end
function affect!(integrator, idx)
kill_line_affect(idx)(integrator)
end
return VectorContinuousCallback(condition, affect!, num_e)
end
function plot_flow(saved_edge_data, network)
t = saved_edge_data.t
vals = saved_edge_data.saveval
data = zeros(length(t), length(vals[1]))
for i in 1:length(vals)
# devide by 1.63 to normalize
data[i, :] = abs.(vals[i])/1.63
end
gs = network.f.graph_structure
sym = Array{String}(undef, (1,gs.num_e))
for i in 1:gs.num_e
sym[1,i] = string(gs.s_e[i]) * "->" * string(gs.d_e[i])
end
plot(t, data, label=sym, size=(800,600))
end
# Definition of nodes and edges according to schäfer18
I = 1.0
γ = 0.1
verticies = [
SwingVertex(-1.0, I, γ), #1
SwingVertex( 1.5, I, γ), #2
SwingVertex(-1.0, I, γ), #3
SwingVertex(-1.0, I, γ), #4
SwingVertex( 1.5, I, γ), #5
]
swingedges = [
SwingEdge(1, 2, 1.63, 0.6*1.63), #1
SwingEdge(1, 3, 1.63, 0.6*1.63), #2
SwingEdge(1, 5, 1.63, 0.6*1.63), #3
SwingEdge(2, 3, 1.63, 0.6*1.63), #4
SwingEdge(2, 4, 1.63, 0.6*1.63), #5
SwingEdge(3, 4, 1.63, 0.6*1.63), #6
SwingEdge(4, 5, 1.63, 0.6*1.63) #7
]
# BEWARE: ordering of edges ≠ ordering of add_edge calls!
# In this case I've ordered them manualy for this to match.
g = SimpleGraph(length(verticies))
for e in swingedges
add_edge!(g, e.from, e.to)
end
# Define nodes/edges and network
odeverts = [ODEVertex(f! = vert, dim=2, sym=[:d, :w]) for vert in verticies]
staticedges = [StaticEdge(f! = edge, dim=1, sym=[:F]) for edge in swingedges]
swing_network! = network_dynamics(odeverts, staticedges, g)
# x0 determined from static solution
# x0 = [0.0 for i in 1:2*nv(g)]
x0 = [-0.12692637482862684, -1.3649456633810975e-6, 0.14641121510104085, 4.2191082676726005e-7, -0.24376507587890778, 1.567589744768255e-6, -0.12692637482862684, -1.3649456633810975e-6, 0.35120661043511864, 7.403907552948938e-7]
# define callback to save edge values (for plotting)
function save_edges(u, t, integrator)
graph_data = integrator.f.f(u, integrator.p, integrator.t, GetGD)
return copy(graph_data.e_array)
end
saved_edgevalues = SavedValues(Float64, Array{Float64, 1})
save_callback = SavingCallback(save_edges, saved_edgevalues)
# define problem
prob = ODEProblem(swing_network!, x0, (0., 6))
# define kill_first callback, at t=1.0s remove edge(2,4)
kill_first = PresetTimeCallback( 1.0, kill_line_affect(swing_network!, 2, 4))
# define load_watch_callback
load_watch_callback = watch_line_limit_callback(swing_network!)
@info "Start solver"
sol = solve(prob, Tsit5(), callback=CallbackSet(save_callback, kill_first, load_watch_callback), dtmax=0.01)
# plot(sol, vars = idx_containing(swing_network!, :w))
# plot(sol, vars = idx_containing(swing_network!, :d))
plot_flow(saved_edgevalues, swing_network!)
Thanks for the reply and code! I'll try to get something basic implemented this way - if I run into major roadblocks or have some useful input, I'll let you know. In either case, I'd personally still leave the issue open until there's an actual "user-friendly" interface for this sort of stuff, but I also won't mind if you prefer to close it instead.
Hey, has there been an update on this issue? I have a similar use case for synaptic coupling between neurons, in which I need to detect an event, save the time of that event, and use that time to send pulses to coupled units. Specifically, when a neuron j
spikes (one of its variables increases above a threshold) at time t_s
, it sends a current to its neighboring neurons that is a function of the time difference t - t_s
. So the coupling term also needs to access the spike times of each neuron.
I would like to implement this with NetworkDynamics (nice package btw!), but I would first like to have an intuition if using NetworkDynamics for this will be efficient, or whether I should write a more specific code for this.
Thanks a lot!
Hey @KalelR, thank you for your interest in this package :)
I am not an expert on spiking neural networks but in general it sounds possible. In general this would need to be implemented with the DifferentialEquations.jl callback interface. One could add a further internal array of spiketimes for each neuron for a subtype of NetworkDiffEq
and an interface for conveniently building callbacks with access to these arrays.
There was a BoF session at this years JuliaCon 2022 on spiking neural nets. There seems to be a fairly large community of researchers who want to use this, here is a thread on discourse (slightly outdated): https://discourse.julialang.org/t/spiking-neural-networks/63270/17
One of the packages mentioned in this thread which is still actively developed would be: https://github.com/wsphillips/Conductor.jl
That said, if you want to build this into ND.jl i would gladly merge PRs and help with integrating it into the core of the package but i don't have the resources (or background on Spiking NNs) to do it myself at the moment
Hi @lindnemi, thanks for the answer! Indeed, I had a similar idea to your approach and it seems to work well. There's currently a small bug with the Callback interface when neurons are synchronized. I'll first work on fixing it and then at some point it would be nice to work in ND.jl. I quite like the flexibility of the package! Once I have time I can do this but no hurry for now :)
Hi!
During JuliaCon I mentioned this already in Slack, so I'm now just posting an Issue to make sure you don't forget about my feature request :smile:
It would be awesome to have an interface in
NetworkDynamics.jl
for triggering events based on the state of individual components (nodes or edges), but right now that seems to be a little bit complicated and un-documentedFor example for my work, I need to be able to trigger an event when the state of one node crosses above a critical value, which should then affect the state of a different node or edge. A clean way to achieve that might be to only allow interactions via edges, but then there should be a way to specify a callback that is triggered based on the state of (one of) the edge's connected nodes and can affect the states of (one of) the connected nodes. Or maybe there should be an
EventEdge
type?Since
DifferentialEquations.jl
already offers a specializedVectorContinuousCallback
, it'd probably be best ifNetworkDynamics.jl
provided some convenient way to define individual call-backs on a per-node / per-edge basis, but then aggregated all of the individual continuous callbacks behind the scenes into one vector-valued callback. (Probably similar for other callback types, but I haven't given that much thought.)For starters, I'd also be happy with a dirty hack, as this feature is the one thing currently keeping me from using
NetworkDynamics.jl
. My use-case is in neuro-science, but one of you mentioned a use-case for modeling failures in a power-grid, and I imagine there are many more, so I'm sure this would be a useful feature to have!