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. PCYiIiIiIiMiIiQiIiU= np:=nops(N); #number of players IiIl with(combinat):#We need this to find all the possible coalitions. L:=powerset(N):M:=convert(L,list);# All possible coalitions NzI8IjwmIiIiIiIjIiIkIiIlPCVGJkYnRig8JEYnRig8JUYlRidGKDwjRig8JEYlRig8JEYmRig8JUYlRiZGKDwjRiU8I0YmPCRGJUYmPCNGJzwkRiVGJzwkRiZGJzwlRiVGJkYn K:=nops(L);# Number of coalitions IiM7 M:=sort(M,length); NzI8IjwjIiIlPCMiIiI8IyIiIzwjIiIkPCRGK0YlPCRGJ0YlPCRGKUYlPCRGJ0YpPCRGJ0YrPCRGKUYrPCVGKUYrRiU8JUYnRitGJTwlRidGKUYlPCVGJ0YpRis8JkYnRilGK0Yl 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 Shapley shapleyval:=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); shapley = 10/3 shapley = 23/6 shapley = 23/6 shapley = 2 Part 2: Calculate the Nucleolus Step 1: For each coalition S=M[k] set up sum_{i in S}y_i for 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)); 0, y4, y1, y2, y3, y3+y4, y1+y4, y2+y4, y1+y2, y1+y3, y2+y3, y2+y3+y4, y1+y3+y4, y1+y2+y4, y1+y2+y3, 0 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)); 0, -y4, -y1, -y2, -y3, 2-y3-y4, 3-y1-y4, 2-y2-y4, 4-y1-y2, 4-y1-y3, 6-y2-y3, 8-y2-y3-y4, 7-y1-y3-y4, 7-y1-y2-y4, 10-y1-y2-y3, 13 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}; end if; C:=C union {add(y||i,i=1..nops(N))=v(N)}; end do: Find the smallest a so that C(a) is not empty: with(simplex): minimize(a,C); PCcvSSN5NEc2IiMiIiQiIiMvSSN5MUdGJSMiIihGKC9JI3kyR0YlIyIiKkYoL0kjeTNHRiVGKy8mSSJhR0YlNiMiIiIjISIkRig= We may replace a with -3/2 as the smallest epsilon, but not the y's: a:=-3/2; IyEiJCIiIw== y||np:=v(N)-add(y||j,j=1..nops(N)-1); LCoiIzgiIiJJI3kxRzYiISIiSSN5MkdGJkYnSSN5M0dGJkYn 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}; end if; end do;C; PDAxLChJI3kxRzYiIiIiSSN5MkdGJkYnSSN5M0dGJkYnIyIjQiIiIzEsJEYlISIiIyEiJEYsMSwkRihGL0YwMSwkRilGL0YwMSwmRiVGL0YoRi8jISM2RiwxLCZGJUYvRilGL0Y4MSwmRihGL0YpRi8jISM6RiwxRiUjIiIoRiwxRigjIiIqRiwxRilGRDEsKEYlRi9GKEYvRilGLyMhI0JGLDEsJkYlRidGKEYnIyIjPkYsMSwmRihGJ0YpRicjIiM8RiwxLCZGJUYnRilGJ0ZN Step 4: To find out if C has more than one allocation, Solve the inequalities : with(SolveTools:-Inequality): glc:=LinearMultivariateSystem(C,[y1,y2,y3]); PCc3JTwjL0kjeTFHNiIiIiQ8Iy9JI3kyR0YnIiIlPCMvSSN5M0dGJyMiIioiIiM3JTwkMUYmIyIiKEYyMkYoRiY8JDFGK0YwMiwmRjciIiJGJiEiIkYrPCMvRi8sKEYmRj5GK0Y+IyIjQkYyRj03JTwkRjgyRiZGNjwjL0YrRjxGLTclPCMvRiZGNjwjL0YrRjZGLTclRiQ8JEY6MkYsRis8Iy9GLywmIyIjPEYyRj1GK0Y+ The only set of inequalities that are all consistent with the others is the one requiring y3=-y1-y2+23/2, so we set y3:=23/2-y1-y2; LChJI3kxRzYiISIiSSN5MkdGJEYlIyIjQiIiIyIiIg== 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}; end if; end do: 3, "excess(", 3, ")=", -y1 4, "excess(", 4, ")=", -y2 5, "excess(", 5, ")=", y1+y2-23/2 6, "excess(", 6, ")=", -11+y1+y2 7, "excess(", 7, ")=", 3/2-y1 8, "excess(", 8, ")=", 1/2-y2 9, "excess(", 9, ")=", -y1-y2+4 10, "excess(", 10, ")=", -15/2+y2 11, "excess(", 11, ")=", -11/2+y1 12, "excess(", 12, ")=", -5+y1 13, "excess(", 13, ")=", -6+y2 14, "excess(", 14, ")=", 11/2-y1-y2 The next set for the next least core is: C; PC4xLCYjISM6IiIjIiIiSSN5Mkc2IkYoJkkiYUdGKjYjRicxLCYjISM2RidGKEkjeTFHRipGKEYrMSwkRjIhIiJGKzEsJEYpRjVGKzEsKEYyRihGKUYoIyEjQkYnRihGKzEsKEYxRihGMkYoRilGKEYrMSwmISImRihGMkYoRisxLCYhIidGKEYpRihGKzEsKCMiIzZGJ0YoRjJGNUYpRjVGKzEsJiMiIiRGJ0YoRjJGNUYrMSwmI0YoRidGKEYpRjVGKzEsKEYyRjVGKUY1IiIlRihGKw== minimize(a,C); PCUvSSN5MUc2IiMiIzgiIiUvJkkiYUdGJTYjIiIjIyEiKEYoL0kjeTJHRiVGKA== Now we know the second least core with a: a:=-7/4; C; IyEiKCIiJQ== PC4xSSN5MUc2IiMiIzoiIiUxRiQjIiM4RigxSSN5MkdGJSMiIzxGKDEsJkYkISIiRi1GMiMhI0hGKDEsJEYtRjIjISIoRigxLCZGJCIiIkYtRjsjIiNSRigxRjojIiNQRigxRi0jIiNCRigxLCRGJEYyRjcxRkUjISM4RigxRjYjISIqRigxRjEjISNCRig= Solve the ineqs to see that y1=13/4: glc:=LinearMultivariateSystem(C,[y1,y2]); PCM3JDwjL0kjeTFHNiIjIiM4IiIlPCQxIiIlSSN5Mkc2IjFJI3kyRzYiIyIjPCIiJQ== y1:=13/4; IyIjOCIiJQ== 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}; end if; end do: -y2 -33/4+y2 -31/4+y2 1/2-y2 3/4-y2 -15/2+y2 -6+y2 9/4-y2 Here is the 3rd least core: C; PCoxLCYjIiIiIiIjRiZJI3kyRzYiISIiJkkiYUdGKTYjIiIkMSwmI0YuIiIlRiZGKEYqRisxLCYjISM6RidGJkYoRiZGKzEsJiMiIipGMkYmRihGKkYrMSwmISInRiZGKEYmRisxLCRGKEYqRisxLCYjISNMRjJGJkYoRiZGKzEsJiMhI0pGMkYmRihGJkYr We look for the smallest a so this set is nonempty: minimize(a,C); PCQvJkkiYUc2IjYjIiIkIyEjOiIiKS9JI3kyR0YmIyIjTEYr a:=-15/8; IyEjOiIiKQ== glc:=LinearMultivariateSystem(C,[y2]); PCM3IzwjL0kjeTJHNiIjIiNMIiIp Finally we are left with one point: y2:=33/8; IyIjTCIiKQ== The Nucleolus!!!!!!!!!!!!!!! print(seq(y||i,i=1..4)); NiYjIiM4IiIlIyIjTCIiKUYmIyIiJCIiIw== The Shapley value!!!!!!!!!!!! print(seq(x[i],i=1..4)); NiYjIiM1IiIkIyIjQiIiJ0YmIiIj To set up the schedules for each player calculate the minutes each player saves under each allocation: First: Hours worked: h:=5;h:=7;h:=6;h:=3; IiIm IiIo IiIn IiIk 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; LyZJJm5taW5zRzYiNiMiIiIkIiQwIiIiIQ== JCIkKyIiIiE= LyZJJm5taW5zRzYiNiMiIiMkIisrKytEPCEiKA== JCIkIT4iIiE= LyZJJm5taW5zRzYiNiMiIiQkIisrKytENiEiKA== JCIkSSIiIiE= LyZJJm5taW5zRzYiNiMiIiUkIiMhKiIiIQ== JCIjZyIiIQ==