Computing Wardrop Equilibria via Convex Programming

Input

Enter graph and cost function data in this section or use one of the existing examples

Braess Graph

ClearAll["Global`*"]
(* Graph *)
graph = Graph[{s->v,s->w,v->w,v->t,w->t}];
(* Slopes of affine linear cost functions *)
a={1,0,0,0,1};
(* Offsets of affine linear cost functions *)
b={0,1,0,1,0};
(* Specify source and sink vertices *)
source=s;
sink=t;

Program

(* Prepare Lists/Values associated with the graph *)
vertices=VertexList[graph];
edges=EdgeList[graph];
sourceIndex=Position[vertices,source][[1,1]];
sinkIndex=Position[vertices,sink][[1,1]];
n=VertexCount[graph];
m=EdgeCount[graph];

Γ=IncidenceMatrix[graph];

(* Create cost functions *)
c[e_,x_]:=a[[e]]*x+b[[e]];
(* Create itegrated cost functios *)
cInt[e_,x_]:=Integrate[c[e,ξ],ξ]/.ξ->x;

(* Print Graph with cost functions *)
Print["Input Graph with cost functions:"];
edgelabels=Table[ edges[[e]]->c[e,x], {e,1,m} ];
Graph[graph,VertexLabels->Automatic,EdgeLabels->edgelabels]

(* Create list of flow-variables for the convex program *)
flowvars = Table[f[e],{e,1,m}];

(* Create the objective function == Sum over all iterated cost functios *)
objective=Total[Table[ cInt[e,flowvars[[e]]], {e,1,m} ]];
(* Create flow conditions, i.e, excess = -r / +r / 0 for v = s / t / otherwise *)
(* -> Compute the excess vector containing the excess of all vertices depending on the flow definded by flowvars *)
excess=Γ.flowvars;
(* -> Create a list of flow conservation rules *)
flowconditions=Table[ excess[[v]]==If[v==sourceIndex,-r,If[v==sinkIndex,r,0]] ,{v,1,n}];

(* Create non-negativity constraints f[e]≥0 *)
nonnegconditions=Table[ flowvars[[e]]≥0, {e,1,m} ];

(* Print the conxex program *)
Print["Convex program:"];
Print[" min ",objective];
For[i=1,i≤n,i++,
praefix=If[i==1," s.t. ", "      "];
Print[praefix, flowconditions[[i]]];
];
Print["      ", flowvars, " ≥ 0"];

(* Compute and print the solution *)
Print["Solution:"];
Manipulate[
(* Create new flowconditions vector that depends on the "rate" variable of the Manipulate function *)
fc=flowconditions/.r->rate;
(* Execute the actual convex optimization *)
solution=Minimize[{objective}~Join~fc~Join~nonnegconditions,flowvars];
(* Extract optimal value and the optimal solution from the solution variable *)
optvalue=solution[[1]];
optflow=solution[[2]];
(* Create Arrays that describe the labels and the style of the edges *)
edgeweights=(Table[ edges[[e]]->f[e],{e,1,m}])/.optflow;
edgestyle=(Table[ edges[[e]]-> Directive[Thickness[(1/2+f[e])/50],If[f[e]==0,Gray,Blue]], {e,1,m} ])/.optflow;
(* Create the actual graph, by copying the input graph and equipping it with EdgeWeights and Labels *)
Graph[graph,EdgeWeight->edgeweights,EdgeLabels->"EdgeWeight", VertexLabels->Automatic, EdgeStyle->edgestyle]
,{rate,0,4,1/1000}]