Computing Equilibria, Optimal Flows, and the Price of Anarchy for Parallel Links

Input

Specify the cost functions of the two parallel links

Example 1

ClearAll["Global`*"]
(* Input of the cost functions of the two parallel links *)
c1[x_]:=2x;
c2[x_]:=x+2;

Example 2

ClearAll["Global`*"]
(* Input of the cost functions of the two parallel links *)
c1[x_]:=x+1;
c2[x_]:=2x

Example with one parameter

ClearAll["Global`*"]
(* Input of the cost functions of the two parallel links *)
c1[x_]:=a1*x;
c2[x_]:=1;

Fully Parametrized Example

ClearAll["Global`*"]
(* Input of the cost functions of the two parallel links *)
c1[x_]:=a1*x+b1;
c2[x_]:=a2*x+b2;

Pigou Example (with worst case PoA)

ClearAll["Global`*"]
(* Input of the cost functions of the two parallel links *)
c1[x_]:=x;
c2[x_]:=1;

Program

Compute Equilibrium and Optimal Flows

(* Define a functions "equilibrium" depending on cost functions cf1, cf2 that computes a Wardrop Equilibrium for two parallel links with these cost functions *)
equilibrium[cf1_,cf2_]:=
(* Define a Module with local variables c1, c2 that are initialized with the arguments of the "equilibrium"-function *)
Module[{c1=cf1,c2=cf2},
(* Solve the equilibrium conditions and take the first solution (assuming that there is a unique solution) *)
sol=Solve[ c1[f1] == c2[r-f1], f1][[1]];
(* Create a list of the two flows on both edges (f1, f2) where f2 = r - f1. Plug in the solution computed above *)
flows=FullSimplify[({f1,r-f1}/.sol)];

(* Cut off negative flows and flows > r *)
flows=(Max[0,Min[r,#]])&/@flows;
(* Return the resulting flow vector *)
Return[FullSimplify[PiecewiseExpand/@flows]];
];

(* Call the "equilibrium" function with normal cost functions *)
eqflows=equilibrium[c1,c2];
Print["Equilbrium Flows: x1=", eqflows[[1]], "   x2 = ",eqflows[[2]]];

(* Compute the total cost of the Equilibrium solution as (Subscript[x, 1])(Subscript[c, 1])((Subscript[x, 1])) + (Subscript[x, 2])(Subscript[c, 2])((Subscript[x, 2])) *)
eqcost = Simplify[flows[[1]]*c1[flows[[1]]]+flows[[2]]*c2[flows[[2]]]];
Print["Equilibrium Cost: ", eqcost]


(* Compute the marginal cost as c(x) + x c'[x] *)
margc1[x_]:=c1[x]+x*(D[c1[ξ],ξ]/.ξ->x);
margc2[x_]:=c2[x]+x*(D[c2[ξ],ξ]/.ξ->x);



(* Call the "equilibrium" function with normal cost functions *)
optflows=equilibrium[margc1,margc2];
Print["Optimal Flows: xOpt1=", optflows[[1]], "   x2 = ",optflows[[2]]];

(* Compute the total cost of the Equilibrium solution *)
optcost = Simplify[flows[[1]]*c1[optflows[[1]]]+flows[[2]]*c2[optflows[[2]]]];
Print["Optimal Cost: ", optcost]


(* Compute the price of Anarchy *)
(* Assumptions are only needed for parametrized example *)
poa=Simplify[PiecewiseExpand[Simplify[eqcost/optcost, Assumptions->{r>0,a1>0,a2>0,b2≥b1,b1>0}]], Assumptions->{r>0,a1>0,a2>0,b2≥b1,b1>0}];
Print["Price of Anarchy: ", poa];

(* Compute the maximal price of anarchy *)
poasol=Maximize[poa,r];
Print["Maximal PoA is ", poasol[[1]], " at r=", r/.poasol[[2]]];

Plot functions

(* Plot the Equilibrium flow functions *)
Plot[eqflows,{r,0,2},
 PlotLabel->"Wardrop Equilibrium Flow Functions",
PlotLegends->"Expressions",
PlotLabels->{"x1", "x2"}]

(* Plot the optimal flow functions *)
Plot[optflows,{r,0,2},
 PlotLabel->"Wardrop Equilibrium Flow Functions",
PlotLegends->"Expressions",
PlotLabels->{"xOpt1", "xOpt2"}]

Plot[poa,{r,0,2}, PlotLabel->"Price of Anarchy", PlotLegends->"Expressions"]

Compute Price of Anarchy

To use this computation, load the parametrized example first and compute the equilibrium and optimal functions

Maximize[{poa,a1≥0,a2>0,b2≥b1,b1≥0,r≥0},{a1,a2,b1,b2,r}]