%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Competitive and Cooperative Foraging Games: Static One-Stage Case
% By: K. Passino
% Version: 6/17/02
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
N=2; % Number of players (can't change)
M=2; % The number of different types of resources (e.g., nutrients)
% Set the number of cells (locations) in the space (i.e., bins along a line)
Q=21;
% Set the number of different possible values of the decision variables for the players
D1=Q; % Sets that player 1 can make decisions i=1,2,...,D1. Decision corresponds to moving the
% player to position i.
D2=D1; % Sets that player 2 can make decisions j=1,2,...,D2, meaning placing player 2 at position j.
theta1=1:D1;
theta2=1:D2;
% Set consumption amount. Assume that this same amount is consumed of each type of
% resource in each cell.
z1=1; % Sets the amount of resources consumed on one visit to a cell by player 1
z2=1; % Same, but for player 2
% Define the initial resource map for each type of resource:
for m=1:M
% r(:,m)=rand(Q,1); % Just a way to initialize
r(:,m)=[1 2 3 4 5 6 7 8 9 10 11 10 9 8 7 6 5 4 3 2 1]'; % An example profile (for testing no energy case)
% r(:,m)=[1 2 3 4 5 4 3 2 1 2 3 4 5 6 6 6 5 4 3 2 1]'; % An example profile (for testing energy case, two humps, L=1)
if m==1
r(:,m)=[1 2 3 4 5 4 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0]'; % First resource is close to 0
end
if m==2
r(:,m)=[0 0 0 0 0 0 0 0 1 2 3 4 5 6 6 6 5 4 3 2 1]'; % Second resoure is further away
end
end
% Set depletion rate for resource m (independent of cell)
alpha=1*ones(M,1);
% Initialize consumption amounts
C=0*ones(N,M,D1,D2);
% Set resource priorities for each player
p1=ones(M,1);
p2=p1;
p2=[1 2]'; % To make it so that player 2 prefers resource 2
% Next set how much it costs in energy to go get a resource (line thinking of players starting at 0)
% (scales its importance relative to value of resources (so scaling must be correct or it won't move))
we1=0.1;
we2=0.1;
%we1=0; % To remove the effect of energy
%we2=0;
% Next, compute how much is consumed of each resource m,
% if players 1 and 2 make all possible moves (decisions, i and j)
for m=1:M
for i=1:D1
for j=1:D2
if i ~= j
C(1,m,i,j)=r(i,m)*(1-exp(-alpha(m)*z1)); % Set consumption for forager 1 when forager 2 not there
C(2,m,i,j)=r(j,m)*(1-exp(-alpha(m)*z2)); % Set consumption for forager 2 when forager 1 not there
end
if i==j
C(1,m,i,j)=0.5*r(i,m)*(1-exp(-alpha(m)*(z1+z2))); % Split return when both present
C(2,m,i,j)=0.5*r(j,m)*(1-exp(-alpha(m)*(z1+z2)));
end
end
end
end
% Assume a movement cost that corresponds to each player being at position 0 and having to move to i (j)
% Set weighting matrices on cost of movement (amount of energy to forage):
for i=1:D1
J1e(i)=we1*i; % So moving one unit costs we1
end
for j=1:D2
J2e(j)=we2*j;
end
% Set the cost matrices J1(i,j) and J2(i,j):
for i=1:D1
for j=1:D2
temp1=0;
temp2=0;
for m=1:M % Sum up resources consumed for play i,j at location q by both players
temp1=temp1+p1(m)*C(1,m,i,j);
temp2=temp2+p2(m)*C(2,m,i,j);
end
J1(i,j)=-sum(temp1)+J1e(i); % Sums over all cells the consumption for play i,j and the energy
% expended to get it. Note - sign.
J2(i,j)=-sum(temp2)+J2e(j);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
clf
bar(1:Q,r,'stacked')
colormap(cool)
xlabel('Cell number')
ylabel('Resource amount')
title('Initial resource distribution')
figure(2)
clf
surf(1:D2,1:D1,J1), colormap(white)
xlabel('Player 2 decision, j')
ylabel('Player 1 decision, i')
zlabel('J_1')
title('Cost function for player 1')
figure(3)
clf
surf(1:D2,1:D1,J2), colormap(white)
xlabel('Player 2 decision, j')
ylabel('Player 1 decision, i')
zlabel('J_2')
title('Cost function for player 2')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Nash equilibria:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
flag=0; % Flag for saying if there is no Nash equilibria
for i=1:D1
for j=1:D2
if J1(i,j)<=min(J1(:,j)) & J2(i,j)<=min(J2(i,:)), % Conduct two inequality tests
display('Nash equilibrium and outcome:') % If satisfied, then diplay solution
i
j
J1(i,j)
J2(i,j)
flag=1; % Indicates that there was one Nash equilibrium (or more)
end
end
end
if flag==0
display('There were no Nash equilibria')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Minimax strategy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the minimax strategy (the name "minimax" clearly arises from the computations involved in
% each player trying to minimize its maximum loss):
% First, compute the security strategy for player 1
[maxvals1,indexmax1]=max(J1'); % This is the max value for each row (note transpose)
[minimaxvalue1,secstratP1]=min(maxvals1) % Display the security value of P1 and its security strategy
% (P1 loses no more than minimaxvalue1, no matter what strategy
% P2 uses)
% Second, compute the security strategy for player 2 (note difference from computation of security
% strategies for matrix games since both payoff matrices are "loss" matrices, and note that computation of the
% minimax strategies is completely independent)
[maxvals2,indexmax2]=max(J2); % This is the min value for each column (note no transpose)
[minimaxvalue2,secstratP2]=min(maxvals2) % Display the security value of P2 and its security strategy
% The outcome of the game will be (with a minimax strategy)
outcome1=J1(secstratP1,secstratP2)
outcome2=J2(secstratP1,secstratP2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stackelberg strategy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the Stackelberg strategy:
% First, find follower reactions (note that this just finds the min loss by P2, so it can
% be that more than one P2 strategy results in this loss):
for i=1:D1
minJ2(i)=min(J2(i,:)); % Finds the minimum loss of P2 given the leader P1 chooses i
end
P1loss=-inf*ones(1,D1); % Initialization so that it will pick the first min value above it
jstar=0*ones(1,D1); % Initialization to nonvalid values (valid ones picked next)
for i=1:D1
for j=1:D2
if J2(i,j)==minJ2(i) % Test that j corresponds to a min point for P2, for a given i
if J1(i,j)>=P1loss(i) % Tests if it is above any previously stored value for the loss
% (this finds the security value against all possible P2 rational reactions)
P1loss(i)=J1(i,j); % Keeps the value to compare to any other possible P2 reactions later
jstar(i)=j; % Save the index of the reaction of P2 that results in worst loss for P1 if it uses i
end
end
end
end
% Specify the Stackelberg strategy:
[stackelbergcost,istar]=min(P1loss) % Display the P1 Stackelberg strategy and cost
jstar(istar) % Display the P2 follower strategy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pareto solutions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First, start with the Pareto cost defined via scalarization
p=0:0.01:1; % Used for defining the Pareto cost matrices
for k=1:length(p)
for i=1:D1
for j=1:D2
Jp(i,j)=p(k)*J1(i,j)+(1-p(k))*J2(i,j); % Form the Pareto cost matrix
end
end
[temp,indicesmin]=min(Jp); % Find the min elements of each column of Jp, and their indices
[paretooptimalval(k),paretooptimalstrat2(k)]=min(temp); % Gives min value and its column index
paretooptimalstrat1(k)=indicesmin(paretooptimalstrat2(k)); % Gives the row index
J1opt(k)=J1(paretooptimalstrat1(k),paretooptimalstrat2(k)); % Find the cost to J1 for the Pareto-point
J2opt(k)=J2(paretooptimalstrat1(k),paretooptimalstrat2(k)); % Find the cost to J2 for the Pareto-point
end
figure(4)
clf
subplot(311)
plot(p,paretooptimalval,'ko',p,J1opt,'k-',p,J2opt,'k--')
ylabel('Costs')
title('J_p (o), J_1 (-), and J_2 (--) for Pareto points')
subplot(312)
plot(p,paretooptimalstrat1,'kx')
ylabel('i^*')
title('Decision of player 1 (x)')
subplot(313)
plot(p,paretooptimalstrat2,'k*')
xlabel('Pareto parameter, p')
ylabel('j^*')
title('Decision of player 2 (*)')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute all Pareto solutions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PP=ones(size(J1)); % Initialize a matrix that will hold flags indicating if a point is a Pareto Point
% (initially it indicates that they are all Pareto points)
% Compute the family of Pareto-optimal strategies:
for ii=1:length(theta1) % These are the loops for the test points theta^*
for jj=1:length(theta2)
for iii=1:length(theta1) % These are the loops for the points theta
for jjj=1:length(theta2)
% Perform tests to determine if (ii,jj) is a Pareto point
if (iii ~=ii & jjj~=jj) &...
((J1(iii,jjj) <= J1(ii,jj)) & (J2(iii,jjj) <= J2(ii,jj))) &...
((J1(iii,jjj) < J1(ii,jj)) | (J2(iii,jjj) < J2(ii,jj)))
PP(ii,jj)=0; % If find one such time that the conditions hold then it is not a Pareto point
end
end
end
end
end
%%%%%%%%%%%%%%%%%%% Plot the Pareto points
figure(5)
clf
colormap(jet)
contour(theta2,theta1,J1,14)
hold on
contour(theta2,theta1,J2,14)
hold on
for ii=1:length(theta1) % These are the loops for the test points theta^*
for jj=1:length(theta2)
if PP(ii,jj)==1
plot(theta2(jj),theta1(ii),'kx')
hold on
end
end
end
xlabel('j')
ylabel('i')
title('J_1, J_2, "x" marks a Pareto solution')
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%