SDE Tutorial

An IJulia notebook corresponding to this tutorial will be available on GitHub soon.

Topics covered in this tutorial include:

the swing equation
fixpoint search of ODE systems
constructing an SDE problem with NetworkDynamics.jl

Example: Fluctuations in Power Grids

This tutorial explains the use of stochastic differential equations (SDE's) in NetworkDynamics.jl. As an example system we will simulate fluctuations in a simple power grid model.

The phase and frequency dynamics in power grids can be modeled by the swing equation. In the complex systems community this is also known as the Kuramoto model with inertia.

\[\begin{aligned} \dot \theta_i &= \omega_i \\ M_i \dot \omega_i &= P_i(t) - D_i \omega_i - \sum_{i=1}^N K_{ij} \sin(\theta_i - \theta_j) \end{aligned}\]

Here, $M_i$ and $D_i$ are inertia and damping parameters, $K_{ij}$ is the power capacitiy of the line and $P_i(t)$ is the power infeed at the node. We assume that this power can be separated into a constant power setpoint $P^\circ_i$ and a stochastic fluctuation.

\[\begin{aligned} P_i(t) = P^\circ_i + \sigma_i \cdot \xi_i(t) \end{aligned}\]

In this tutorial we assume the fluctuations $\xi_i(t)$ to be white Gaussian noise. Separating the deterministic and stochastic part, we can write this problem as a stochastic differential equation.

\[\begin{aligned} d\omega_i = \left[P_i - D_i \omega_i - \sum_{i=1}^N K_{ij} \sin(\theta_i - \theta_j)\right] dt + \sigma_i dW_i \end{aligned}\]

Here, $dW_i = \xi_i dt$ is the infinitesimal increment of the Wiener process. Problems of this type can be numerically solved in Julia as described in the SDE Tutorial of DifferentialEquations.

Implementing the Swing Equation

First we will implement the node and edge functions for the deterministic case without the fluctuations. We assume the inertia and damping parameters to be $M_i = 1.0$ and $D_i = 0.1$.

function swing_equation!(dv, v, edges, P, t)
    dv[1] = v[2]
    dv[2] = P - 0.1 * v[2]
    for e in edges
        dv[2] += e[1]
    end
end

function powerflow!(e, v_s, v_d, K, t)
    e[1] = K * sin(v_s[1] - v_d[1])
end

Contructing the Deterministic Dynamics

For the graph structure we will use a simple 4 node ring network.

using Graphs
g = watts_strogatz(4, 2, 0.0)

Then we can construct the ODEFunction of the deterministic system by using network_dynamics().

using NetworkDynamics

swing_vertex = ODEVertex(; f=swing_equation!, dim=2, sym=[:θ, :ω])
powerflow_edge = StaticEdge(; f=powerflow!, dim=1)

nd = network_dynamics(swing_vertex, powerflow_edge, g)

[ Info: Reconstructing EdgeFunction with :undefined coupling type..

Fixpoint Search

Now we need to define the dynamic parameters of vertices and edges. For simplicity we assume homogeneous capacities on the lines. For the nodes we assume half of them to be net producers ($P = 1.0$) and half of them to be net consumers ($P = -1.0$) of power.

K = 6.0
P = [1.0, -1.0, 1.0, -1.0]
p = (P, K)

We want to simulate fluctuations around an equilibrium state of our model system. Therefore, we need to find a fixpoint of the determinitic system which can be done by using the utility function find_fixpoint(). As an initial guess we take all variables equal to zero.

u0 = find_fixpoint(nd, p, zeros(8))

using StochasticDiffEq, OrdinaryDiffEq
ode_prob = ODEProblem(nd, u0, (0.0, 500.0), p)
ode_sol = solve(ode_prob, Tsit5())

using Plots, LaTeXStrings
plot(ode_sol; vars=syms_containing(nd, "ω"), ylims=(-1.0, 1.0), ylabel=L"\omega", legend=false, fmt=:png)

We see that this is in fact a fixpoint solution. We will later use this as an initial condition for the numerical integration of the SDE system.

Adding a Stochastic Layer

For adding the stochastic part of the dynamics we have to define a second graph layer. In our example, the fluctuations at different nodes are independent of each other. Therefore, we define a second graph with the same number of vertices but without any edges.

h = SimpleGraph(4, 0)

The dynamics at the nodes has to have the same dimension as in the deterministic case. In our example we only have fluctuations in the second variable.

function fluctuation!(dx, x, edges, p, t)
    dx[1] = 0.0
    dx[2] = 0.05
end

Now we can construct the dynamics of the second layer by using network_dynamics(). Since the graph structure of the stochastic layer has no edges we can take the edge function of the deterministic case as a placeholder.

fluctuation_vertex = ODEVertex(; f=fluctuation!, dim=2)
nd_noise = network_dynamics(fluctuation_vertex, powerflow_edge, h)

[ Info: Reconstructing EdgeFunction with :undefined coupling type..

Simulating the SDE

Finally, we can create an SDEProblem and solve it with DifferentialEquations.

sde_prob = SDEProblem(nd, nd_noise, u0, (0.0, 500.0), p)
sde_sol = solve(sde_prob, SOSRA())
plot(sde_sol; vars=syms_containing(nd, "ω"), ylims=(-1.0, 1.0), ylabel=L"\omega", legend=false, fmt=:png)

More details on SDE problems, e.g. how to include correlations or how to define an EnsembleProblem, can be found in the documentation of DifferentialEquations.