partetooptimalsolution.m

一个用MATLAB编写的优化控制工具箱

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%-------------------------------------