Solving for voltage potentials in an electrical network

Define a graph and gather some data

Define a graph, choose directions

For simplicity, we assume every edge has unit resistance.  Let’s start with a random graph.

n=10;
m=25;
g=RandomGraph[{n,m}]

We’ll want to assign directions to the edges in order to define the boundary matrix and solve for the currents.

directedg=Graph[EdgeList[g]/.UndirectedEdge[a_,b_]->DirectedEdge[a,b]]

It will be convenient to have lists of vertices and edges.

vertices=VertexList[directedg]
edges=EdgeList[directedg]

Notice that Mathematica does not store the vertex list in numerical order.  So use caution when interpreting the adjacency and incidence matrices!

The boundary matrix

boundarymatrix=Normal[IncidenceMatrix[directedg]];
MatrixForm[boundarymatrix]

The cycle matrix

Here, I illustrate how to form the cycle matrix.  However, in the next secion, we’ll use n absolute voltage potentials at each vertex, which will make our m edge potentials automatically satisfy Kirchoff’s potential law.  So, this section on the cycle matrix isn’t necessary for solving for the potentials.

We want a basis for the cycle space of the undirected graph, but we want the cycles to be oriented according to the orientations we have already specified.  Now, we could find a spanning tree and write a short program to find a basis for the cycle space but Mathematica has a built in function we can use.

cyclelist=FindFundamentalCycles[g]/.UndirectedEdge[a_,b_]->DirectedEdge[a,b]

Look, you can get nice pictures of the cycles:

cyclepictures=TableForm[Table[HighlightGraph[g,cycle],{cycle,FindFundamentalCycles[g]}]];

k=Length[cyclelist]

m-n+Length[ConnectedComponents[g]]

We want to express the cycles as elements of the edge space (Subscript[C, 1])(g), then collect them as rows of a matrix C.

cyclematrix=ConstantArray[0,{k,m}];

Do[Do[If[MemberQ[edges,edge],
cyclematrix[[j,Position[edges,edge][[1,1]]]]=1,
cyclematrix[[j,Position[edges,edge[[2]]\[DirectedEdge]edge[[1]]][[1,1]]]]=-1],{edge,cyclelist[[j]]}],{j,1,k}]

MatrixForm[cyclematrix]

Solve for voltage potentials using Kirchoff’s current law

Let’s collect some absolute voltage variables in a list.

vvect=Table[Subscript[v,i],{i,vertices}]

Then the voltage potentials of each edge are given by.

wvect=Table[Subscript[v,edge[[1]]]-Subscript[v,edge[[2]]],{edge,edges}]

edges

These voltage potentials also signify the current since I’m assuming the resistance of each edge is 1.  We should divide by the resistance of each edge if resistance differs from 1.

Notice that Kirchoff’s potential law is automatically satisfied for the potential differences on the edges when they’ve been defined via absolute potentials on the vertices.

cyclematrix.wvect

The current law, however, states that the current in to each vertex equals current out is a nontrivial condition.  Look at

boundarymatrix.wvect

Now, in order to have current flowing in the network, we specify a source and sink (you can do it randomly, but be sure they’re different!)

source=RandomChoice[vertices]
sink=RandomChoice[vertices]

vertexcurrents=ConstantArray[0,{n}];
vertexcurrents[[VertexIndex[directedg,source]]]=-c;
vertexcurrents[[VertexIndex[directedg,sink]]]=c;

vertexcurrents

Now, we can solve the equations in Kirchoff’s current law.

sols=Solve[{boundarymatrix.wvect==vertexcurrents,(Subscript[v, source])==1,(Subscript[v, sink])==0},Join[vvect,{c}]]

Let’s name the solution

voltages=(vvect/.sols)[[1]]

current=c/.sols [[1]]

Here are the currents along the edges

wvect/.sols[[1]]

Random Walks

The wander to the source game

For simplicity, we’ll weight each edge by 1.   Let’s program a little game on a graph.  Choose a source and sink.  Starting at an initial vertex, choose an adjacent vertex at random and continue on a random walk.  Reaching the source before landing on the sink is a win.  Reaching the sink before the source is a loss.  We’ll write the program to play the game a number of times and to output the percentage of wins.

randomwalkgame[graph_,source_,sink_,initial_,iterations_]:=

Module[
{currentvertex,newvertex,wins,losses}, (* here are the variables we will use in our program *)

wins=0;
losses=0;

If [initial==source, wins=iterations,If[initial==sink,losses=iterations]];
(* always win if initial is source, always lose if initial is sink *) 

currentvertex=initial;  (* set initial vertex *)

While[wins+losses<iterations,currentvertex=RandomChoice[AdjacencyList[graph,currentvertex]];

If [currentvertex==source,
wins++;
currentvertex=initial];
If[currentvertex==sink,
losses++;
currentvertex=initial];
];

wins/iterations(* output the proportion of visits to each webpage *)

]

expectedreturns=N[Table[randomwalkgame[g,source,sink,j,50000],{j,vertices}]]

As proved in Modern Graph Theory by Bela Bollobas (Theorem 8 on page 304) these expected returns are the voltage potentials.  Compare our numerical approximations of the expected returns to the actual voltage potentials computed earlier:

N[voltages]

The escape probability

In the previous game, choosing the source as the initial vertex always yields a win.  The “escape probability” is defined to the the probability that starting at the source, we reach the sink before returning to the source.

escapeprob[graph_,source_,sink_,iterations_]:=

Module[
{currentvertex,newvertex,wins,losses}, (* here are the variables we will use in our program *)

vertices=VertexList[graph]; (* get the vertices of the graph *) 
edges=EdgeList[graph];(* get the edges of the graph *) 

wins=0;
losses=0;

currentvertex=source;  (* set initial vertex to be the source *)

While[wins+losses<iterations,currentvertex=RandomChoice[AdjacencyList[graph,currentvertex]];

If [currentvertex==source,
losses++;
currentvertex=source];
If[currentvertex==sink,
wins++;
currentvertex=source];
];
wins/iterations(* output the proportion of visits to each webpage *)

]

ep=escapeprob[g,source,sink,50000]

According to the second part of Theorem 8, the escape probability times the degree of the source is the total current of the network

N[ep]*VertexDegree[g,source]

N[current]

Dirichlet’s Principle

In the context of electrical networks, Dirichlet’s principle is that absolute potentials that minimize energy produce electrical circuits satisfying Kirchoff’s current law.

Energy of a network

Again, for simplicity, we assume the resistance of every edge is one.  The energy of a point v in the vertex space (Subscript[C, 0])(g) is defined to be ,<v,Lv>.

adj=AdjacencyMatrix[directedg]+Transpose[AdjacencyMatrix[directedg]];
diag=DiagonalMatrix[Table[VertexDegree[g,v],{v,VertexList[directedg]}]];
lap=diag-adj;
MatrixForm[lap]

energy[v_]:=Dot[v,lap.v]

vvect

Notice that the energy is a quadratic function of the values {(Subscript[v, i])} assigned to the vertices.

energy[vvect]

For our purposes, we want the (Subscript[v, source])=1and (Subscript[v, sink])=0 and the remaining values will be variable

newvect=vvect/.{(Subscript[v, source])->1,(Subscript[v, sink])->0}

variables=Delete[newvect,{{VertexIndex[directedg,source]},{VertexIndex[directedg,sink]}}]

energy[newvect]

The theorem (Theorem 1 in Chapter IX of the book) is that energy is minimized when the values {(Subscript[v, i])} are the voltage potentials of an electrical network.  The energy is minimized where it’s derivative is zero:

Let’s name them dirichlet solutions

dirichletsols=Solve[Table[D[energy[newvect],x]==0,{x,variables}],variables][[1]];
dirichletvoltages=vvect/.Join[dirichletsols,{(Subscript[v, VertexIndex[directedg,source]])->1,(Subscript[v, VertexIndex[directedg,sink]])->0}]

Comparing them to the voltages obtained by solving the equations of Kirchoff’s current law, we see they are the same

voltages

You may be interested to notice that Mathematica has a built in Minimize function:

Minimize[energy[newvect],variables]

current

Exercises

Generalize the discussion above to include conductance.  Test your work on the network pictured below.  Use the vertex 6 (the red one) as the source and vertex 4 (the black one) as the sink.  Each each is labelled by its conductance (that will be more convenient than resistance).  So you don’t have to type all this data in, I’ve provided the adjacency matrix, weighted by the conductance of each edge.

(({{0, 5, 2, 0, 7, 0, 0, 0, 0, 0}, {5, 0, 9, 0, 2, 0, 1, 6, 0, 0}, {2, 9, 0, 2, 0, 0, 0, 4, 0, 0}, {0, 0, 2, 0, 0, 0, 0, 8, 5, 0}, {7, 2, 0, 0, 0, 2, 5, 0, 0, 0}, {0, 0, 0, 0, 2, 0, 1, 0, 0, 7}, {0, 1, 0, 0, 5, 1, 0, 0, 3, 9}, {0, 6, 4, 8, 0, 0, 0, 0, 8, 2}, {0, 0, 0, 5, 0, 0, 3, 8, 0, 8}, {0, 0, 0, 0, 0, 7, 9, 2, 8, 0}}))