%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For training Takagi-Sugeno fuzzy systems using batch least squares
%
% By: Kevin Passino
% Version: 1/26/99
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% First, generate the training data, G

% For the M=121 case

x=-6:0.1:6;

M=length(x)

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Consider the R=20 case

% First, form the vector Y

Y=Gz';

% Next, find Phi, which involves processing x through phi

R=20;

c(1,:)=-5.4:0.6:6;
sigma(1,:)=0.1*ones(1,R);

% Initialize Phi

	for j=1:R
		mu(j,1)=exp(-0.5*((x(1)-c(1,j))/sigma(1,j))^2);
	end
	
	denominator(1)=sum(mu(:,1)); 
	
	for j=1:R
		xi(j,1)=mu(j,1)/denominator(1);
	end

Phi=[xi(:,1)', x(1)*xi(:,1)'];

% Form the rest of Phi

for i=2:M
	
	for j=1:R
		mu(j,i)=exp(-0.5*((x(i)-c(1,j))/sigma(1,j))^2);
	end
	
	denominator(i)=sum(mu(:,i)); 
	
	for j=1:R
		xi(j,i)=mu(j,i)/denominator(i);
	end

	Phi=[Phi; xi(:,i)', x(i)*xi(:,i)'];
	
end

% Find the least squares estimate

theta=Phi\Y
theta40=theta;  % This is for saving it for later use in other programs
save variables theta40

% Note that we tested the result by plotting the resulting approximator
% mapping and it produces a reasonable result.  It is for this reason
% that we trust the numerical computations, and do not seek to 
% use other methods for the computation of the estimate.

% Next, compute the approximator values

for i=1:M,
	phi=[xi(:,i)', x(i)*xi(:,i)']';
	Fts(i)=theta'*phi;
end

% Next, plot the basis functions, data and the approximator to compare

figure(1)
plot(x,xi,'k')
xlabel('x')
ylabel('Basis function values')
title('Takagi-Sugeno fuzzy system, 20 rules, basis functions')
grid

figure(2)
plot(x,Gz,'ko',x,Fts,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and fuzzy system output')
title('Takagi-Sugeno fuzzy system, 20 rules')
grid
axis([min(x) max(x) 0 max(G)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Next, consider a different value for sigma

Y=Gz';

% Next, find Phi, which involves processing x through phi

clear Phi

sigma(1,:)=1*ones(1,R);

% Initialize Phi

	for j=1:R
		mu(j,1)=exp(-0.5*((x(1)-c(1,j))/sigma(1,j))^2);
	end
	
	denominator(1)=sum(mu(:,1)); 
	
	for j=1:R
		xi(j,1)=mu(j,1)/denominator(1);
	end

Phi=[xi(:,1)', x(1)*xi(:,1)'];

% Form the rest of Phi

for i=2:M
	
	for j=1:R
		mu(j,i)=exp(-0.5*((x(i)-c(1,j))/sigma(1,j))^2);
	end
	
	denominator(i)=sum(mu(:,i)); 
	
	for j=1:R
		xi(j,i)=mu(j,i)/denominator(i);
	end

	Phi=[Phi; xi(:,i)', x(i)*xi(:,i)'];
	
end

% Find the least squares estimate

theta=Phi\Y

% Next, compute the approximator values

for i=1:M,
	phi=[xi(:,i)', x(i)*xi(:,i)']';
	Fts(i)=theta'*phi;
end

% Next, plot the basis functions, data and the approximator to compare

figure(3)
plot(x,xi,'k')
xlabel('x')
ylabel('Basis function values')
title('Takagi-Sugeno fuzzy system, 20 rules, basis functions')
grid

figure(4)
plot(x,Gz,'ko',x,Fts,'k')
xlabel('x(i)')
ylabel('y(i)=G(x(i),z(i)), and fuzzy system output')
title('Takagi-Sugeno fuzzy system, 20 rules')
grid
axis([min(x) max(x) 0 max(G)])

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