In [1]:
import networkx
import condynet.graph
import condynet.plot
import condynet.drawing
import h5obj
import glob
import bundle

Collective Nonlinear Dynamics of
Electricity Networks (CoNDyNet)


Determination of resonance frequencies of LC networks with binary link disorder


Daniel Jung

Jacobs University Bremen
School of Engineering and Science
Solid State Physics Group

December 10, 2014

Table of contents

  1. Motivation
    • LC circuit
    • Flow equations
  2. Model
    • Basic Graph Theory
    • The 2D grid model
    • The Small World Model
    • Adding Edge Attributes
  3. Numerical method
    • Resonance frequencies
  4. Resonance spectra
    • The 2D grid model
    • The Small World Model
  5. Summary

Motivation


\[Z = R + i X\]
\[\frac{1}{X} = \frac{1}{i\omega L} + i\omega C\]

A LC circuit has a resonance frequency of \(\omega_0 = \frac{1}{\sqrt{LC}}\)

Frequency \(\omega\) of AC current approaching the resonance frequency: \[\lim\limits_{\omega \rightarrow \omega_0} Z(\omega) = \infty\]

\(\Rightarrow\) No power transmission possible!

Normal operation: \(\omega\) should stay far below the smallest resonance.

Arbitrary LC network

  • Coupled LC oscillators
  • Set of resonance frequencies \(\omega_n\)

Flow equations I

An arbitrary one-phasic AC grid
with line impedances \(Z_{ij}\).

Generator: \(I_i > 0\)
Consumer: \(I_i < 0\)

\[\mathbf{Y} = \mathbf{Z}^{-1}\]

Combine Ohm's law

\[V_{ij} = Z_{ij} I_{ij} \quad\text{or}\quad I_{ij} = Y_{ij} V_{ij} \qquad\qquad\]

with Kirchhoff's laws for each node and mesh

\[I_i = \sum_j I_{ij} \qquad V_{ij} = V_i - V_j\]

to derive the current flow equations

\[I_i = \sum_j Y_{ij} (V_i - V_j) \quad\rm.\]

Name Symbol
Impedance matrix \(\mathbf{Z}\)
Admittance matrix \(\mathbf{Y}\)

R. Huang, G. Korniss, S. Nayak, 2009

Flow Equations II

An arbitrary one-phasic AC grid.

Generator: \(I_i > 0\)
Consumer: \(I_i < 0\)

\[I_i = \sum_j Y_{ij} (V_i - V_j)\]

To get a standard matrix-vector multiplication, reformulate to

\[I_i = \sum_j L_{ij} V_j \quad\mathrm{or}\quad \mathbf{I} = \mathbf{L} \mathbf{V}\]

with

\[L_{ij} = \delta_{ij} \sum_{k\neq i} Y_{ik} - (1-\delta_{ij}) Y_{ij}\]

Note:

  • \(\mathbf{L}\) is defined in analogy to the topological network Laplacian \(\mathbf{G}\).
  • \(\mathbf{L}\) is commonly referred to as admittance matrix as well.

R. Huang, G. Korniss, S. Nayak, 2009

In [2]:
figsize(4, 1.5)
g = networkx.Graph()
g.add_edges_from([(0, 1), (0, 2), (1, 2), (2, 3)])
networkx.draw(g, node_color='white', with_labels=True)
savefig('pic/small-graph.svg', bbox_inches='tight')
close(gcf())

Basic Graph Theory I

Example

\(N = 4\)

A graph is given by

  • a set of nodes \(i\)
  • a set of edges \((i, j)\)

Properties

Property Explanation
Order \(N\) Number of nodes
Size Number of edges
In [3]:
D = diag(g.degree().values())
print D
[[2 0 0 0]
 [0 2 0 0]
 [0 0 3 0]
 [0 0 0 1]]

In [4]:
import easytable
t = easytable.autotable(D, hc='&', hp=1)
tl = [r+'\\\\' for r in t.split('\n')]
print '\n'.join(tl)
 2 & 0 & 0 & 0 \\
 0 & 2 & 0 & 0 \\
 0 & 0 & 3 & 0 \\
 0 & 0 & 0 & 1 \\

In [5]:
E = networkx.adjacency_matrix(g).todense()
print E
[[0 1 1 0]
 [1 0 1 0]
 [1 1 0 1]
 [0 0 1 0]]

In [6]:
G = networkx.laplacian_matrix(g).todense()
print G
[[ 2 -1 -1  0]
 [-1  2 -1  0]
 [-1 -1  3 -1]
 [ 0  0 -1  1]]

Basic Graph Theory II

Example

\[D = \begin{pmatrix} 2 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 \\ 0 & 0 & 3 & 0 \\ 0 & 0 & 0 & 1 \\ \end{pmatrix}\]

\[E = \begin{pmatrix} 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ \end{pmatrix}\]

\[G = \begin{pmatrix} 2 & -1 & -1 & 0 \\ -1 & 2 & -1 & 0 \\ -1 & -1 & 3 & -1 \\ 0 & 0 & -1 & 1 \\ \end{pmatrix}\]

Node Properties

Property Explanation
Degree Number of incident edges

Matrices

Name Explanation
Degree matrix \(\mathbf{D}\) Diagonal matrix containing all node degrees
Adjacency matrix \(\mathbf{E}\) Nonzero only if nodes \(i\) and \(j\) are adjacent
Laplacian matrix \(\mathbf{G}\) \(\mathbf{G} = \mathbf{D} - \mathbf{E}\)

\[G_{ij} = \delta_{ij} \sum_{k\neq i} E_{ik} - (1-\delta_{ij}) E_{ij}\]

In [7]:
evs = networkx.laplacian_spectrum(g)
condynet.plot.pspectrum(evs, limits=[-.2, 4.2], figsize=[6, 1], width=1)
savefig('pic/laplacian-spectrum.svg', bbox_inches='tight')
print around(evs, 6)
[-0.  1.  3.  4.]

Basic Graph Theory III

Example

Laplacian spectrum:
[0.  1.  3.  4.]

\[G = \begin{pmatrix} 2 & -1 & -1 & 0 \\ -1 & 2 & -1 & 0 \\ -1 & -1 & 3 & -1 \\ 0 & 0 & -1 & 1 \\ \end{pmatrix}\]

Laplacian spectrum

Certain eigenvalues (EVs) of the Laplacian matrix \(\mathbf{G}\) have special properties:

  • All EVs are non-negative
    (as \(\mathbf{G}\) is positive semi-definite).
  • At least one eigenvalue is \(0\).
  • Number of eigenvalues equal
    to \(0\): Number of connected subgraphs.
  • Second-smallest EV:
    Algebraic connectivity.
In [8]:
g = networkx.grid_2d_graph(3, 3, periodic=True)
figsize(4, 4)
condynet.drawing.draw_grid2d_periodic(g, node_color='white', with_labels=True, node_size=2000)
savefig('pic/grid2d.svg')
close(gcf())

The regular 2D grid graph

\(N=9\), periodic boundary conditions
In [43]:
g = condynet.graph.small_world(12, p=1.5)
figsize(4, 4)
networkx.draw_circular(g, node_color='white', with_labels=True, node_size=500, width=3)
savefig('pic/small-world.svg', bbox_inches='tight')
close(gcf())

The Small World Model

Example

\(N=12 \qquad p=1.5\)

Closed ring with \(N\) nodes, plus \(S\) randomly chosen shortcuts.

Shortcut density \(\frac{p}{2} = \frac{S}{N}\)

Number of shortcuts \(S = \frac{Np}{2}\)

Two limiting cases:

Limiting case Resulting graph
\(p_{\rm min} = 0\) Closed ring
\(p_{\rm max} = N - 3\) Complete graph

Maximum number of shortcuts:

\[S_{\rm max} = \frac{N(N-3)}{2}\]


J. Travers, S. Milgram, 1969; M. Newman, D. Watts, 1999; R. Huang, G. Korniss, S. Nayak, 2009

In [39]:
q = .5
hvals = choose(rand(g.number_of_edges()) < q, [-1, 1])
networkx.set_edge_attributes(g, 'hval', dict(zip(g.edges(), hvals)))
In [41]:
edgelist, hvals = zip(*networkx.get_edge_attributes(g, 'hval').items())
edgecolors = choose(array(hvals) == 1, array(['r', 'g']))
networkx.draw_circular(g, edgelist=edgelist, edge_color=edgecolors, width=3,
                       node_color='white', node_size=500, with_labels=True)
savefig('pic/small-world-binary.svg', bbox_inches='tight')
close(gcf())

Adding Edge Attributes

Example

\(N=12\)
\(p=1.5 \qquad q=0.5\)

In order to describe LC networks, we attribute an impedance \(Z_{ij}\)
to each edge.

Here, we consider random impedances, using a binary distribution:

Edge type \(Z_{ij}\) Chance
Capacitance \((i\omega C)^{-1}\) \(q\)
Inductance \(i\omega L\) \(1-q\)

R. Huang, G. Korniss, S. Nayak, 2009

Edge type \(Z_{ij}\) Chance \(Y_{ij}\) \(h_{ij}\)
Capacitance \((i\omega C)^{-1}\) \(q\) \(y_1 = i\omega C\) \(-1\)
Inductance \(i\omega L\) \(1-q\) \(y_2 = (i\omega L)^{-1}\) \(1\)

Introduce matrix \(\mathbf{h}\):

\[h_{ij} = \begin{cases} -1 & \rm,\quad \text{edge $(i,j)$ carries capacitance $C$} \\ 1 & \rm,\quad \text{edge $(i,j)$ carries inductance $L$} \\ 0 & \rm,\quad \text{no edge between $i$ and $j$.} \end{cases}\]

\(\Rightarrow\) Similar to \(\mathbf{E}\), but also \(-1\) allowed.


R. Huang, G. Korniss, S. Nayak, 2009

Resonance frequencies I

Remember:

Flow equations:

\[I_i = \sum_j L_{ij} V_j \quad\mathrm{or}\quad \mathbf{I} = \mathbf{L} \mathbf{V}\]
with "Laplacian matrix" \(\mathbf{L}\).
  • Consider resonance case, \(I_i = 0\)
  • The system can be seen as a system of coupled LC oscillator circuits
  • The system has \(N\) resonance frequencies \(\omega_n\)

\[\sum_j L_{ij}(\omega)\, V_j = 0 \qquad\text{or}\qquad \mathbf{L} \mathbf{V} = \mathbf{0}\]


R. Huang, G. Korniss, S. Nayak, 2009

Resonance frequencies II

Remember: \[\begin{align} y_1 &= i\omega C \\ y_2 &= (i\omega L)^{-1} \end{align}\]
\[h_{ij} = \begin{cases} -1 &\rm,\quad Y_{ij} = y_1 \\ 1 &\rm,\quad Y_{ij} = y_2 \\ 0 &\rm,\quad Y_{ij} = 0 \end{cases}\]

Define \(\mathbf{H}\),

\[H_{ij} = \delta_{ij} \sum_{k\neq i} h_{ik} - (1-\delta_{ij}) h_{ij} \quad\rm,\]

and

\[\lambda = \frac{y_1 + y_2}{y_1 - y_2} \quad\rm,\]

so that the flow equations for the resonance case can be rewritten as

\[\mathbf{L} \mathbf{V} = \mathbf{0} \qquad\Rightarrow\qquad (\mathbf{H} - \lambda \mathbf{G}) \mathbf{V} = \mathbf{0} \quad\rm.\]

But this is not yet a regular eigenvalue problem...


R. Huang, G. Korniss, S. Nayak, 2009; Fyodorov 1999

Resonance frequencies III

Remember: \[\lambda = \frac{y_1 + y_2}{y_1 - y_2}\]
\[\begin{align} y_1 &= i\omega C \\ y_2 &= (i\omega L)^{-1} \end{align}\]

Define

\[\tilde{\mathbf{H}} = \mathbf{G}^{-1/2}\, \mathbf{H}\, \mathbf{G}^{-1/2}\]

(real symmetric) and

\[\tilde{\mathbf{V}} = \mathbf{G}^{1/2} \mathbf{V} \quad\rm.\]


Notes:

  1. \(\mathbf{G}\) is real and positive semi-definite, so its matrix square root \(\mathbf{G}^{1/2}\)
    is uniquely defined.
  2. \(\mathbf{G}\) is positive semi-definite, i.e. it has at least one eigenvalue \(0\) and hence is always singular. So its pseudo-inverse \(\mathbf{G}^{-1}\)
    has to be considered.

Fyodorov 1999

Resonance frequencies IV

Remember: \[\lambda = \frac{y_1 + y_2}{y_1 - y_2}\]
\[\begin{align} y_1 &= i\omega C \\ y_2 &= (i\omega L)^{-1} \end{align}\]

Define

\[\tilde{\mathbf{H}} = \mathbf{G}^{-1/2}\, \mathbf{H}\, \mathbf{G}^{-1/2}\]

(real symmetric) and

\[\tilde{\mathbf{V}} = \mathbf{G}^{1/2} \mathbf{V} \quad\rm.\]

Then, we can rewrite

\[(\mathbf{H} - \lambda \mathbf{G}) \mathbf{V} = \mathbf{0} \qquad\Rightarrow\qquad \mathbf{\tilde{H}}\, \mathbf{\tilde{V}}_n = \lambda_n \mathbf{\tilde{V}}_n\]

So we are facing a regular eigenvalue problem, with a known relationship between the eivenvalues \(\lambda_n\)
and the resonance frequencies \(\omega_n\):

\[\omega_n = \frac{1}{\sqrt{LC}} \sqrt{\frac{1 + \lambda_n}{1 - \lambda_n}}\]


Fyodorov 1999

In [12]:
# load data
data2d = {}
filenames = glob.glob('/home/djung/work/condynet/dorg2d/l????/q????.h5')
for filename in filenames:
    f = h5obj.File(filename)
    try:
        reson = f['reson']
        dens = f['dens']
        stderr = f['stderr']
        param = f['param']
    except:
        continue
    shape = param.shape
    shape0, shape1 = shape
    if shape0 != shape1:
        raise ValueError
    q = param.q
    if shape0 not in data2d:
        data2d[shape0] = {}
    if q in data2d[shape0]:
        raise ValueError('data for L=%i and q=%03f already exists' % (shape0, q))
    data2d[shape0][q] = bundle.Bundle(reson=reson, dens=dens, stderr=stderr, param=param)
In [13]:
figsize(8, 5)
L = 30
for q in [.2, .3, .4, .5]:
    d = data2d[L][q]
    errorbar(d.reson, d.dens, yerr=d.stderr, label='q=%.1f' % q)
axis([-1, 1, 0, None])
legend(loc=1)
title('Resonance spectrum of 2D grid with binary link disorder,\n'
      + r'$N\, =\, 30 \times\, 30\, =\, 900$, for different composite ratios $q$', fontsize=18)
xlabel(r'$\lambda$', fontsize=24)
ylabel(r'$\rho(\lambda)$', fontsize=24)
savefig('pic/ador-q05.svg', bbox_inches='tight')
close(gcf())

Resonance spectra I: 2D grid

Define the density of resonances (DOR):

\[\rho(\lambda) = \frac{1}{N} \sum_{n=1}^{N_{\mathrm{R}}} \delta(\lambda - \lambda_n) \qquad\qquad\]

\(N_{\mathrm{R}}\rm:\)
Number of "true" resonances
\((-1 < \lambda_n < 1)\)

Obtain ensemble average (arithmetic mean) of the DOR over many disorder realizations (ADOR).


R. Huang, G. Korniss, S. Nayak, 2009

In [14]:
w, h = 8, 5
In [15]:
# load data
# data is saved in the form data[N][p]
data = {}
filenames = glob.glob('/home/djung/work/condynet/dorswe/n*/p*h5')
for filename in filenames:
    f = h5obj.File(filename)
    try:
        reson = f['reson']
        dens = f['dens']
        stderr = f['stderr']
        param = f['param']
    except:
        continue
    size = param.size
    p = param.p
    if size not in data:
        data[size] = {}
    if p in data[size]:
        raise ValueError('data for N=%i and p=%02f already exists' % (size, p))
    data[size][p] = bundle.Bundle(reson=reson, dens=dens, stderr=stderr, param=param)
In [16]:
figsize(w, h)
start, stop, step = 4, None, 10
for (N, color, marker) in zip([100, 200, 400, 1000], ['red', 'green', 'blue', 'magenta'],
                              ['o', 's', 'D', '^']):
    for (p, facecolor) in zip([N/2, N-3], ['white', color]):
        d = data[N][p]
        errorbar(d.reson[start:stop:step], d.dens[start:stop:step], yerr=d.stderr[start:stop:step], label='N=%i, p=%g' % (N, p),
                 marker=marker, markeredgecolor=color, color=facecolor, markersize=12, linestyle='--')
axis([-1, 1, 0, None])
legend(loc=1)
title('density of resonances of the small world model with binary link disorder, \n' +
      'with system size $N$, shortcut density $p$ and composite ratio $q=1/2$', fontsize=12)
xlabel(r'$\lambda$', fontsize=24)
ylabel(r'$\rho(\lambda)$', fontsize=24)
savefig('pic/fig1c.svg', bbox_inches='tight')
close(gcf())
In [17]:
figsize(w, h)
start, stop, step = 4, None, 10
N = 1000
for (p, color, marker) in zip([.01, .1, 1, 2, 4, 10, 100], ['purple', 'red', 'blue', 'cyan', 'magenta', 'brown', 'lightgreen'], 
                              ['o', 's', 'D', '^', '<', 'v', '>']):
    d = data[N][p]
    plot(d.reson[start:stop:step], d.dens[start:stop:step], label='p=%g' % p,
         marker=marker, markeredgecolor=color, color=color, markersize=12, linestyle='--')
axis([-1, 1, 0, None])
legend(loc=1)
title('density of resonances of the small world model with binary link disorder, \n' +
      'with system size $N=1000$, shortcut density $p$ and composite ratio $q=1/2$', fontsize=12)
xlabel(r'$\lambda$', fontsize=24)
ylabel(r'$\rho(\lambda)$', fontsize=24)
savefig('pic/fig1b.svg', bbox_inches='tight')
close(gcf())

Resonance spectra II: Small World Model

Large-\(p\) limit

\(\Rightarrow\) Confirming results by Huang et al.


R. Huang, G. Korniss, S. Nayak, 2009

In [44]:
figsize(w, h)
start, stop, step = 4, None, 10
p = .01
for (N, color, marker) in zip([100, 200, 400, 1000], ['red', 'lightgreen', 'blue', 'magenta'], 
                              ['o', 's', 'D', '^']):
    d = data[N][p]
    plot(d.reson[start:stop:step], d.dens[start:stop:step], label='N=%i' % N,
         marker=marker, markeredgecolor=color, color=color, markersize=12, linestyle='--')
axis([-1, 1, 0, None])
legend(loc=2)
title('density of resonances of the small world model with binary link disorder, \n' +
      'with system size $N$, shortcut density $p=0.01$ and composite ratio $q=1/2$', fontsize=12)
xlabel(r'$\lambda$', fontsize=24)
ylabel(r'$\rho(\lambda)$', fontsize=24)
savefig('pic/fig2a.svg', bbox_inches='tight')
close(gcf())
In [45]:
figsize(w, h)
start, stop, step = 4, None, 10
p = .1
for (N, color, marker) in zip([100, 200, 400, 1000], ['red', 'lightgreen', 'blue', 'magenta'], 
                              ['o', 's', 'D', '^']):
    d = data[N][p]
    plot(d.reson[start:stop:step], d.dens[start:stop:step], label='N=%i' % N,
         marker=marker, markeredgecolor=color, color=color, markersize=12, linestyle='--')
axis([-1, 1, 0, None])
legend(loc=2)
title('density of resonances of the small world model with binary link disorder, \n' +
      'with system size $N$, shortcut density $p=0.1$ and composite ratio $q=1/2$', fontsize=12)
xlabel(r'$\lambda$', fontsize=24)
ylabel(r'$\rho(\lambda)$', fontsize=24)
savefig('pic/fig2b.svg', bbox_inches='tight')
close(gcf())

Resonance spectra III: Small World Model

Small-\(p\) limit

\(\Rightarrow\) Largely confirming results by Huang et al.

\(\Rightarrow\) Difference: Peak at \(\lambda = 0\).


R. Huang, G. Korniss, S. Nayak, 2009

In [20]:
figsize(w, h)
start, stop, step = 4, None, 10
for (N, color, marker) in zip([100, 200, 400, 1000], ['red', 'green', 'blue', 'magenta'],
                              ['o', 's', 'D', '^']):
    for (p, facecolor) in zip([N/2, N-3], ['white', color]):
        d = data[N][p]
        errorbar(d.reson[start:stop:step]*sqrt(p), d.dens[start:stop:step]/sqrt(p), yerr=d.stderr[start:stop:step],
                 label='N=%i, p=%g' % (N, p), marker=marker, markeredgecolor=color, color=facecolor,
                 markersize=12, linestyle='--')
axis([-6, 6, 0, None])
legend(loc=1)
title('density of resonances of the small world model with binary link disorder, \n' +
      'with system size $N$, shortcut density $p$ and composite ratio $q=1/2$', fontsize=12)
xlabel(r'$\lambda\sqrt{p}$', fontsize=24)
ylabel(r'$\rho(\lambda)/\sqrt{p}$', fontsize=24)
savefig('pic/fig1ci.svg', bbox_inches='tight')
close(gcf())

Resonance spectra IV: Scaling law

Test scaling law, valid in the large-\(p\) limit:

\[\rho(\lambda) = p^{\frac{1}{2}} \phi(p^{\frac{1}{2}} \lambda)\]


R. Huang, G. Korniss, S. Nayak, 2009

Summary & Outlook

Summary

  • Description of LC networks
    • simple graphs
    • binary distribution of edge impedances
  • Calculation of resonance frequencies
    and the density of resonances

Outlook

  • Different topologies (triangular grid, honeycomb grid, and realistic network topologies).
  • Other impedance distributions (also continuous distributions), also including ohmic resistances.
  • Beyond the resonance case (current and power flow calculations).
Thank you for your attention!

References

Jacobs University Bremen
School of Engineering and Science
Solid State Physics Group

Collective Nonlinear Dynamics of
Electricity Networks (CoNDyNet)