%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Example of finding Pareto-optimal solutions for a multiobjective % optimization problem. % % Example of computation of Pareto-optimal strategies for % bimatrix games. Ignores issues of multiple optima in payoff % matrices and as p varies. % % Example of problems found in defining the Pareto cost % % Author: K. Passino % Version: 4/5/02 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % First, example: multiobjective optimization - quadratic case clear all % Define payoff (cost) functions for each player theta1=-4:0.05:4; % Ok, while we think of it as an infinite game computationally % we of course only study a finite number of points. m=length(theta1); theta2=-5:0.05:5; n=length(theta2); % A set of cost functions, J1 for player 1 and J2 for player 2 for ii=1:length(theta1) for jj=1:length(theta2) J1(ii,jj)=(theta1(ii)-2)^2 + (theta2(jj)-3)^2; J2(ii,jj)=(theta1(ii)+2)^2 + (theta2(jj)+2)^2;; end end % Next, compute the reaction curves for each player and will plot on top of J1 and J2 for j=1:length(theta2) [temp,t1]=min(J1(:,j)); % Find the min point on J1 with theta2 fixed, for all theta2 R1(j)=theta1(t1); % Compute the theta1 value that is the best reaction to each theta2 end for i=1:length(theta1) [temp,t2]=min(J2(i,:)); % Find the min point on J2 with theta1 fixed R2(i)=theta2(t2); % Compute the theta2 value that is the best reaction to each theta1 end % Compute the family of Pareto-optimal strategies: p=0:0.01:1; % Used for defining the Pareto cost matrices for k=1:length(p) for i=1:m for j=1:n 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 % Next, find the theta1 and theta2 values for plotting below theta1k(k)=theta1(paretooptimalstrat1(k)); theta2k(k)=theta2(paretooptimalstrat2(k)); end figure(1) clf colormap(jet) contour(theta2,theta1,J1,14) hold on contour(theta2,theta1,J2,14) hold on plot(theta2,R1,'k-') hold on plot(R2,theta1,'k--') hold on plot(theta2k,theta1k,'x') xlabel('\theta^2') ylabel('\theta^1') title('J_1, J_2, R_1 (-), R_2 (--), "x" marks Pareto solution') hold off %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Do another example: a bimatrix game clear all % Set the number of different possible values of the decision variables m=5; % If change will need to modify specific J1 value chosen below n=3; % Set the payoff matrices J1(i,j) and J2(i,j): % For generating random bimatrix games: J1=round(10*rand(m,n)-5*ones(m,n)) % Make it random integers between -5 and +5 J2=round(10*rand(m,n)-5*ones(m,n)) % Make it random integers between -5 and +5 % These can be used to find other equilibria. For instance, if you use the following two % payoff matrices then you will get uniqueness in the response of P2. J1 =[4 1 -4; -2 5 3; -3 2 -1; 4 4 4; -3 -5 2] J2 =[2 -1 0; -2 4 0; -3 0 -1; -3 3 4; -3 0 -5] % One set of payoff matrices that for i=4 by P1 gives two possible responses from the follower P2 % that both achieve minimum loss for P2. J1 =[-1 5 -3; -2 5 1; 4 3 -2; -5 -1 5; 3 0 2] J2 =[-1 2 -3; 2 -3 1; -2 3 1; -1 1 -1; 4 -4 1] % Compute the family of Pareto-optimal strategies: p=0:0.01:1; % Used for defining the Pareto cost matrices for k=1:length(p) for i=1:m for j=1:n 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(2) 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 (*)') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Do another example: To show problems in defining the Pareto cost clear all % Define payoff (cost) functions for each player theta1=-4:0.05:4; % Ok, while we think of it as an infinite game computationally % we of course only study a finite number of points. m=length(theta1); theta2=-5:0.05:5; n=length(theta2); % A set of cost functions, J1 for player 1 and J2 for player 2 for ii=1:length(theta1) for jj=1:length(theta2) J1(ii,jj)=-2*(exp( (-(theta1(ii)-2)^2)/5 +(-4*(theta1(ii)*theta2(jj))/20) + ((-(theta2(jj)-3)^2)/2))); J2(ii,jj)=-1*(exp( (-(theta1(ii)-1)^2)/4 + (5*(theta1(ii)*theta2(jj))/10) + ((-(theta2(jj)+1)^2)/2))); end end % Next, compute the reaction curves for each player and will plot on top of J1 and J2 for j=1:length(theta2) [temp,t1]=min(J1(:,j)); % Find the min point on J1 with theta2 fixed, for all theta2 R1(j)=theta1(t1); % Compute the theta1 value that is the best reaction to each theta2 end for i=1:length(theta1) [temp,t2]=min(J2(i,:)); % Find the min point on J2 with theta1 fixed R2(i)=theta2(t2); % Compute the theta2 value that is the best reaction to each theta1 end % Compute the family of Pareto-optimal strategies: p=0:0.01:1; % Used for defining the Pareto cost matrices for k=1:length(p) for i=1:m for j=1:n 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 % Next, find the theta1 and theta2 values for plotting below theta1k(k)=theta1(paretooptimalstrat1(k)); theta2k(k)=theta2(paretooptimalstrat2(k)); end figure(3) clf surf(theta2,theta1,0.5*J1+0.5*J2) colormap(white) xlabel('\theta^2') ylabel('\theta^1') zlabel('J_p') title('J_p for p=0.5') figure(4) clf contour(theta2,theta1,0.5*J1+0.5*J2,10) colormap(jet) xlabel('\theta^2') ylabel('\theta^1') zlabel('J_p') title('J_p contour for p=0.5') figure(5) clf colormap(jet) contour(theta2,theta1,J1,14) hold on contour(theta2,theta1,J2,14) hold on plot(theta2,R1,'k-') hold on plot(R2,theta1,'k--') hold on plot(theta2k,theta1k,'x') xlabel('\theta^2') ylabel('\theta^1') title('J_1, J_2, R_1 (-), R_2 (--), "x" marks Pareto solution') hold off %------------------------------------- % End of program %-------------------------------------