%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 caseclear all% Define payoff (cost) functions for each playertheta1=-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 2for 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;;	endend% Next, compute the reaction curves for each player and will plot on top of J1 and J2for 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 theta2endfor 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 theta1end% Compute the family of Pareto-optimal strategies:p=0:0.01:1; % Used for defining the Pareto cost matricesfor 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));	endfigure(1)clfcolormap(jet)contour(theta2,theta1,J1,14)hold oncontour(theta2,theta1,J2,14)hold onplot(theta2,R1,'k-')hold onplot(R2,theta1,'k--')hold onplot(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 gameclear all% Set the number of different possible values of the decision variablesm=5; % If change will need to modify specific J1 value chosen belown=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 +5J2=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 matricesfor 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-pointendfigure(2)clfsubplot(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 costclear all% Define payoff (cost) functions for each playertheta1=-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 2for 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)));	endend% Next, compute the reaction curves for each player and will plot on top of J1 and J2for 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 theta2endfor 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 theta1end% Compute the family of Pareto-optimal strategies:p=0:0.01:1; % Used for defining the Pareto cost matricesfor 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));	endfigure(3)clfsurf(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)clfcontour(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)clfcolormap(jet)contour(theta2,theta1,J1,14)hold oncontour(theta2,theta1,J2,14)hold onplot(theta2,R1,'k-')hold onplot(R2,theta1,'k--')hold onplot(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%-------------------------------------