unknownfunc.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For plotting training data from the "unknown function" G(x,z)
%
% By: Kevin Passino
% Version: 1/13/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% For the M=7 case
x7=-6:2:6;


M=length(x7)

for i=1:M,
	G7(i)=exp(-50*(x7(i)-1)^2)-0.5*exp(-100*(x7(i)-1.2)^2)+atan(2*x7(i))+2.15+...
	0.2*exp(-10*(x7(i)+1)^2)-0.25*exp(-20*(x7(i)+1.5)^2)+0.1*exp(-10*(x7(i)+2)^2)-0.2*exp(-10*(x7(i)+3)^2);
	if x7(i) >= 0
		G7(i)=G7(i)+0.1*(x7(i)-2)^2-0.4;
	end
end

% For the M=121 case
x121=-6:0.1:6;

M=length(x121)

for i=1:M,
	z(i)=0.15*(rand-0.5)*2;  % Define the auxiliary variable
	G121(i)=exp(-50*(x121(i)-1)^2)-0.5*exp(-100*(x121(i)-1.2)^2)+atan(2*x121(i))+2.15+...
	0.2*exp(-10*(x121(i)+1)^2)-0.25*exp(-20*(x121(i)+1.5)^2)+0.1*exp(-10*(x121(i)+2)^2)-0.2*exp(-10*(x121(i)+3)^2);
	if x121(i) >= 0
		G121(i)=G121(i)+0.1*(x121(i)-2)^2-0.4;
	end
	G121n(i)=G121(i)+z(i); % Adds in the influence of the auxiliary variable
%	fpoly(i)=0.6+0.1*x(i);
end


% Next, plot the functions:

figure(1)
plot(x7,G7,'ko')
xlabel('x(i)')
ylabel('y(i)=G(x(i))')
title('M=7')
grid
axis([min(x7) max(x7) 0 max(G7)])

figure(2)
plot(x121,G121,'ko')
xlabel('x(i)')
ylabel('y(i)=G(x(i))')
title('M=121')
grid
axis([min(x121) max(x121) 0 max(G121)])

figure(3)
plot(x121,G121n,'ko')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i))')
title('M=121')
grid
axis([min(x121) max(x121) 0 max(G121)])

% Next, plot some approximator fits to the function
% First, a line

for i=1:M,
	Flinear(i)=((4.5-0.5)/(6-(-6)))*x121(i)+2.5;
end

figure(4)
plot(x121,G121n,'ko',x121,Flinear,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and line')
title('Linear approximation')
grid
axis([min(x121) max(x121) 0 max(G121)])

% Next, a polynomial (parabola)

for i=1:M,
	Fpoly(i)=0.6+0.035*(x121(i)+6)^2;
end

figure(5)
plot(x121,G121n,'ko',x121,Fpoly,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and parabola')
title('Polynomial approximation')
grid
axis([min(x121) max(x121) 0 max(G121)])

% Next, a single hidden layer perceptron, n1=1

w1=[1.5]';
b1=0;

w=[3]'; 
b=0.6;

for i=1:M,
	phi1=inv(1+exp(-b1-w1*x121(i)));
	phi=[phi1]';
	Fmlp(i)=b+w'*phi;
end

figure(6)
plot(x121,G121n,'ko',x121,Fmlp,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and perceptron output')
title('Neural network approximation, one neuron')
grid
axis([min(x121) max(x121) 0 max(G121)])

% Next, a single hidden layer perceptron, n1=2

w1=[1.5]';
b1=0;

w2=[1.25]'; 
b2=-6;

w=[3 1]'; 
b=0.6;

for i=1:M,
	phi1=inv(1+exp(-b1-w1*x121(i)));
	phi2=inv(1+exp(-b2-w2*x121(i)));
	phi=[phi1 phi2]';
	Fmlp(i)=b+w'*phi;
end

figure(7)
plot(x121,G121n,'ko',x121,Fmlp,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and perceptron output')
title('Neural network approximation, two neurons')
grid
axis([min(x121) max(x121) 0 max(G121)])

% Next, we plot the premise membership functions for
% a Takagi-Sugeno fuzzy system approximator

c11=-3.5;
sigma11=0.8;

c12=-0.25;
sigma12=0.6;

c13=2;
sigma13=0.4;

c14=4.5;
sigma14=0.8; 

for i=1:M,
	mu1(i)=exp(-0.5*((x121(i)-c11)/sigma11)^2);
	mu2(i)=exp(-0.5*((x121(i)-c12)/sigma12)^2);
	mu3(i)=exp(-0.5*((x121(i)-c13)/sigma13)^2);
	mu4(i)=exp(-0.5*((x121(i)-c14)/sigma14)^2);
	denominator(i)=mu1(i)+mu2(i)+mu3(i)+mu4(i); % Used below
end

figure(8)
plot(x121,G121n,'ko',x121,mu1,'k',x121,mu2,'k',x121,mu3,'k',x121,mu4,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)) and premise membership function values')
title('Training data and premise membership functions')
grid

% Next, we plot the basis functions for
% a Takagi-Sugeno fuzzy system approximator

for i=1:M,
	xi1(i)=mu1(i)/denominator(i);
	xi2(i)=mu2(i)/denominator(i);
	xi3(i)=mu3(i)/denominator(i);
	xi4(i)=mu4(i)/denominator(i);
end

figure(9)
plot(x121,G121n,'ko',x121,xi1,'k',x121,xi2,'k',x121,xi3,'k',x121,xi4,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)) and basis function values')
title('Training data and basis functions')
grid

% Next, we plot the Takagi-Sugeno fuzzy system approximator

a10=1; 
a11=0.5/6;

a20=2.25;
a21=4.4/4;

a30=2.9;
a31=1/12;

a40=1.3; 
a41=4.8/8;

for i=1:M,
	g1(i)=a10+a11*x121(i);
	g2(i)=a20+a21*x121(i);
	g3(i)=a30+a31*x121(i);
	g4(i)=a40+a41*x121(i);
	numerator=g1(i)*mu1(i)+g2(i)*mu2(i)+g3(i)*mu3(i)+g4(i)*mu4(i);
	Fts(i)=numerator/(mu1(i)+mu2(i)+mu3(i)+mu4(i));
end

figure(10)
plot(x121,G121n,'ko',x121,Fts,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and fuzzy system output')
title('Takagi-Sugeno approximation')
grid
axis([min(x121) max(x121) 0 max(G121)])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%