%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % For training Takagi-Sugeno fuzzy systems using the % Levenberg-Marquardt method % % By: Kevin Passino % Version: 2/16/99 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear % Specify the maximum number of iterations: Nlm=15; % 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)+1*z(i); % Adds in the influence of the auxiliary variable (or can multiply by 0) end % Next, we provide another option for specifying the data set for M=13 %clear x G z Gz %x=-6: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)+1*z(i); % Adds in the influence of the auxiliary variable (or can multiply by 0) %end % Next, we need to initialize the appoximator c(1,:,1)=-5:1:5; % The first 1 is the index for the first input (n=1 case) % and the second one is for the time index R=length(squeeze(c(1,:,1))) sigma(1,:,1)=0.5*ones(1,R); %sigma(1,:,1)=0.2*ones(1,R); % Another way to initialize %sigma(1,:,1)=1*ones(1,R); % Another way to initialize % We initialize all the parameters that enter linearly to be zero thetalin(:,1)=0*ones(2*R,1); % With this the a_i,0 parameters are first in the thetalin vector, % followed by the a_i,1 parameters, from i=1,2,...,R (thetalin % is the vector of parameters that enter linearly) % Next, form the whole parameter vector is theta(:,1)=[c(1,:,1)'; sigma(1,:,1)'; thetalin(:,1)]; p=length(theta(:,1)); % Set the number of parameters for this case % Next, pick the parameter for the update method (remember you want lambda % bigger to make the updates *smaller*) lambda=0.5; Lambda=lambda*eye(p); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, we start the main loop: for j=1:Nlm j % Print to the screen the interation number % First, check if parameters are in range, and fix if not: for jj=1:R if sigma(1,jj,j)<0.1 sigma(1,jj,j)=0.1; end if sigma(1,jj,j)>1 sigma(1,jj,j)=1; end if c(1,jj,j)>6 c(1,jj,j)=6; end if c(1,jj,j)<-6 c(1,jj,j)=-6; end end % Next, we need to form the vector (epsilon) and matrix (del_epsilon) that are used in the % Levenberg-Marquardt update formula. for i=1:M % Define the phi vector, for each data pair in Gz for jj=1:R mu(jj,i,j)=exp(-0.5*((x(i)-c(1,jj,j))/sigma(1,jj,j))^2); end denominator(i)=sum(mu(:,i,j)); xi(:,i,j)=mu(:,i,j)*(1/denominator(i)); phi(:,i,j)=[xi(:,i,j)', x(i)*xi(:,i,j)']'; Fts(i)=thetalin(:,j)'*phi(:,i,j); epsilon(i,j)=Gz(i)-Fts(i); for jj=1:R g(jj,i,j)=thetalin(jj,j)+thetalin(jj+R,j)*x(i); % First, compute it for the centers, then the spreads, then the consequent parameters temp=(g(jj,i,j)-Fts(i))/denominator(i); del_epsilon(jj,i,j)=-temp*mu(jj,i,j)*(x(i)-c(1,jj,j))*(1/(sigma(1,jj,j)^2)); del_epsilon(jj+R,i,j)=-temp*mu(jj,i,j)*((x(i)-c(1,jj,j))^2)*(1/(sigma(1,jj,j)^3)); del_epsilon(jj+2*R,i,j)=-xi(jj,i,j); % For the a_i,0 del_epsilon(jj+3*R,i,j)=-x(i)*xi(jj,i,j); % For the a_i,1 end end % Next, we seek to generate the next estimate of the parameter vector theta for % the Levenberg-Marquardt update formula: theta(:,j+1)=theta(:,j)-... inv(del_epsilon(:,:,j)*del_epsilon(:,:,j)'+Lambda)*del_epsilon(:,:,j)*epsilon(:,j); % Load the components into the relevant vectors: c(1,:,j+1)=theta(1:R,j+1)'; sigma(1,:,j+1)=theta((R+1):(2*R),j+1)'; thetalin(:,j+1)=theta((2*R+1):(4*R),j+1); %end % End of main loop for Levenberg-Marquardt update formula %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, compute the approximator values (i.e., at a test set) % Define the points at which it will be tested xt=-6:0.01:6; Mt=length(xt); % First, define the phi vector, for each data pair in Gz % (note that j+1 is the last time at which we have generated a theta % and we use it here - first, we have to compute phi) for i=1:Mt for jj=1:R mut(jj,i)=exp(-0.5*((xt(i)-c(1,jj,j+1))/sigma(1,jj,j+1))^2); end denominatort(i)=sum(mut(:,i)); xit(:,i)=mut(:,i)*(1/denominatort(i)); phit(:,i)=[xit(:,i)', xt(i)*xit(:,i)']'; Fts(i)=thetalin(:,j+1)'*phit(:,i); % Just call the jth one the last one end % Next, plot the basis functions, data and the approximator to compare %figure(1) figure(j) clf plot(x,Gz,'ko',xt,xit,'k--',xt,Fts,'k-') xlabel('x') T=num2str(j); T=strcat('Basis functions, data, and approximator mapping, j=',T); ylabel(T) title('Takagi-Sugeno fuzzy system, 11 rules') grid axis([min(x) max(x) 0 max(G)]) % Could omit the "end" on the main loop and add the following to see the mapping % shape at each iteration (in this case change figure(1) above to figure(j) if you want) pause % Can add this if you want to view the mappings at each iteration while % the program is running end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%