This code calculates the Shapley value for an arbitrary n-person cooperative game in characteristic function form. The second part of the code calculates the Nucleolus for the game. The example we use is from problem 5.19 regarding the fair scheduling problem for doctors Moe-Larry-Curly and Shemp.restart:N:={1,2,3,4};#Set up the Players.np:=nops(N); #number of playerswith(combinat):#We need this to find all the possible coalitions.L:=powerset(N):M:=convert(L,list);# All possible coalitionsK:=nops(L);# Number of coalitionsM:=sort(M,length);The Characteristic function is defined here, coalition by coalition.for k from 1 to K do if nops(M[k])<=1 then v(M[k]):=0; end if;end do;v({1,2}):=4:v({1,3}):=4:v({2,3}):=6:v({1,4}):=3:v({2,4}):=2:v({3,4}):=2:
v({1,2,3}):=10:v({1,3,4}):=7:v({2,3,4}):=8:v({1,2,4}):=7:v({1,2,3,4}):=13:Procedure to calculate Shapleyshapleyval:=proc(v,x,N)
local i,k,shapley;
for i from 1 to nops(N)
do
x[i]:=0:
for k from 1 to K do
if member(i,M[k]) and nops(M[k])>=1 then
x[i]:=x[i]+(v(M[k])-v(M[k] minus {i}))*
((nops(M[k])-1)!*(nops(N)-nops(M[k]))!)/nops(N)!
end if;
end do;
lprint(shapley[i]=x[i]);
end do:
end proc:shapleyval(v,x,N);Part 2: Calculate the NucleolusStep 1: For each coalition S=M[k] set up sum_{i in S}y_ifor k from 1 to nops(M)
do
z[k]:=0:
for i from 1 to np do
if member(i,M[k]) and M[k]<>N and M[k]<>{} then
z[k]:=z[k]+y||i
end if;
end do;
end do:lprint(seq(z[i],i=1..16));Step 2: Now set up the excess functions for each coalition:excess:=k->v(M[k])-z[k]:lprint(seq(excess(k),k=1..K));Step 3: setup for the first least core:C:={}:for k from 1 to K do
if M[k]<>{} and M[k]<> N then C:=C union {excess(k)<=a[1]};
end if;
C:=C union {add(y||i,i=1..nops(N))=v(N)};
end do:Find the smallest a[1] so that C(a[1]) is not empty:with(simplex):minimize(a[1],C);We may replace a[1] with -3/2 as the smallest epsilon, but not the y's:a[1]:=-3/2;y||np:=v(N)-add(y||j,j=1..nops(N)-1);Here is the least core, but not in reduced form:C:={}:for k from 1 to K do
if M[k]<>{} and M[k]<> N then C:=C union {excess(k)<=a[1]};
end if;
end do;C;Step 4:To find out if C has more than one allocation, Solve the inequalities :with(SolveTools:-Inequality):glc:=LinearMultivariateSystem(C,[y1,y2,y3]);The only set of inequalities that are all consistent with the others is the one requiring y3=-y1-y2+23/2, so we sety3:=23/2-y1-y2;Step 5: Recalculate C and the excess functions, deleting the excess functions which are now constants.
C:={}:y||nops(N):=v(N)-add(y||j,j=1..nops(N)-1):for k from 1 to K do
if M[k]<>{} and M[k]<> N and
(y1 in indets(excess(k)) or y2 in indets(excess(k)))
then lprint(k,"excess(",k,")=",excess(k));
C:=C union {excess(k)<=a[2]};
end if;
end do:The next set for the next least core is:C;minimize(a[2],C);Now we know the second least core with a[2]:a[2]:=-7/4; C;Solve the ineqs to see that y1=13/4:glc:=LinearMultivariateSystem(C,[y1,y2]);y1:=13/4;Step 6: Since this set still has more than one point we set up the third least core:C:={}:for k from 1 to K do
if M[k]<>{} and M[k]<> N and (y2 in indets(excess(k)))
then lprint(excess(k));
C:=C union {excess(k)<=a[3]};
end if;
end do:Here is the 3rd least core:C;We look for the smallest a[3] so this set is nonempty:minimize(a[3],C);a[3]:=-15/8;glc:=LinearMultivariateSystem(C,[y2]);Finally we are left with one point:y2:=33/8;The Nucleolus!!!!!!!!!!!!!!!print(seq(y||i,i=1..4));The Shapley value!!!!!!!!!!!!print(seq(x[i],i=1..4));To set up the schedules for each player calculate the minutes each player saves under each allocation:First: Hours worked:h[1]:=5;h[2]:=7;h[3]:=6;h[4]:=3;The minutes each player must work is then the original time to work minus the amount of time saved under each allocation:for i from 1 to 4 do nmins[i]=evalf((h[i]-y||i)*60); smins[i]:=evalf((h[i]-x[i])*60); end do;