Session 2 : Exercises

Exercise 1

Given three parallel edges with cost functions (Subscript[c, 1])(x)=(x^2)+x, (Subscript[c, 2])(x)=(x^4)+(x^2), (Subscript[c, 3])(x)=x+1, compute a Wardrop Equilibrium for demand rate r=5.

Solve[ (* ... *) ]

Solution

Solve[{x^2+x==y^4+y^2,y^2==z+1,x+y+z==5, x≥0,y≥0,z≥0},{x,y,z}]

Exercise 2

We want to write a function that computes an equilibrium for two parallel links with affine cost functions in dependence of the demand rate r.

a)

Compute the solution by solving (Subscript[c, 1])(x)=(Subscript[c, 2])(y), x+y=r for x and y. This yields solutions that may be negative. Write a function that cuts off negative flows and flows greater than the demand (i.e. something like max(0,min(r,x))) and apply it to the solution.
You may want to use the function PiecewiseExpand[ ] and apply it to the solution to get a nice, piecewise representation of the function.

(* Some example cost functions *)
c1[x_]:=x+1;
c2[x_]:=2x
(* Solve the system *)
sol=First[Solve[(*...*)]];
flows = (* ... *) /. sol

The solution for the given Example function should be
{({{{, ({{(1/3) (-1+2 r), r>(1/2)}, {0, True}})}}),({{{, ({{r, 0<r≤(1/2)}, {((1+r)/3), r>(1/2)}, {0, True}})}})}

Solution

(* Some example cost functions *)
c1[x_]:=x+1;
c2[x_]:=2x
(* Solve the system *)
sol=First[Solve[c1[f1]==c2[r-f1],f1]];
(* Create a solution vector containing the flows x1, x2 *)
flows=FullSimplify[({f1,r-f1}/.sol)];

(* Cut off negative flows and flows > r *)
(* (Max[0,Min[r,#]])&  --- this creates an anonymous function where # represents the variable 
/@ --- Applies the preceding function componentwise to the vector
*)
flows=(Max[0,Min[r,#]])&/@flows;
(* Creates a simplified, piecewise function *)
FullSimplify[PiecewiseExpand/@flows]

b)

In order to use the same code to compute the equilibrium flows and the optimal flows (which can be obtained by computing an equilibrium wrt. marginal costs), we want to write a function that carries out the code from part a).
To do this, we use the Module[ ] function. This functions creates a scope for local variables and used as follows.

globalvar=42;
f[x0_]:=
(* Define a module with local variables x and y. x is initialized with x0 *)
Module[ {x=x0,y},
(* scope for the local code begins *)
(* some random computations... *)
If[x>1,
y=x^2;
,
y=-x];
(* Global variables can be used *)
y=y+x+globalvar;
(* A line without ; determines the returned value of the module.
You can also use the function Return[ ... ] *)
y
(* End of the module *)
];

Use the code from a) to create a function that computes the equilibrium depending on the two cost functions (Subscript[c, 1]) and (Subscript[c, 2]). The test code should return the same result as in a).

equilibrium[cf1_,cf2_]:=Module[ { (*local vars *) }, (* code *) ];

(* Test *)
c1[x_]:=x+1;
c2[x_]:=2x
equilibrium[c1, c2]

Solution

equilibrium[cf1_,cf2_]:=Module[{c1=cf1,c2=cf2,flows,sol},
sol=First[Solve[c1[f1]==c2[r-f1],f1]];
(* Create a solution vector containing the flows x1, x2 *)
flows=FullSimplify[({f1,r-f1}/.sol)];

(* Cut off negative flows and flows > r *)
(* (Max[0,Min[r,#]])&  --- this creates an anonymous function where # represents the variable 
/@ --- Applies the preceding function componentwise to the vector
*)
flows=(Max[0,Min[r,#]])&/@flows;
(* Creates a simplified, piecewise function and returns it *)
Return[FullSimplify[PiecewiseExpand/@flows]];
];

(* Test *)
c1[x_]:=x+1;
c2[x_]:=2x
equilibrium[c1, c2]

Exercise 3

Now create a notebook that, after input of two (affine linear) cost functions (Subscript[c, 1]),(Subscript[c, 2]), computes the 1) equilibrium flow with its total cost 2) the optimal flow with its cost and 3) the price of anarchy function, which is cost(equilibrium)/cost(optimal flow).
Entering the cost functions

c1[x_]:=x+1;
c2[x_]:=2x

and then executing your code should return the following solution

(you might need to apply PiecewiseExpand and Simplify to your solutions)

Solution

Compute Equilibrium and Optimal Flows

(* Cost function input *)
c1[x_]:=x+1;
c2[x_]:=2x

(* 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]]];