clear		% Clear all variables in memory
eold=0; 	% Intial condition used to calculate c
rold=0; 	% Intial condition used to calculate r
yeold=0; 	% Intial condition used to calculate yc
ymold=0; 	% Initial condition for the first order reference model

% Next, initialize parameters for the fuzzy controller

nume=11; 	% Number of input membership functions for the e
			% universe of discourse
numc=11; 	% Number of input membership functions for the c
			% universe of discourse

ge=1/2;,gc=1/2;,gu=5;
		% Scaling gains for tuning membership functions for
		% universes of discourse for e, c and u respectively
		% These are tuned to improve the performance of the FMRLC
we=0.2*(1/ge);
	% we is half the width of the triangular input membership
	% function bases (note that if you change ge, the base width
	% will correspondingly change so that we always end
	% up with uniformly distributed input membership functions)
	% Note that if you change nume you will need to adjust the
	% "0.2" factor if you want membership functions that
	% overlap in the same way.
wc=0.2*(1/gc);
	% Similar to we but for the c universe of discourse
base=0.4*gu;
	% Base width of output membership fuctions of the fuzzy
	% controller

% Place centers of membership functions of the fuzzy controller:

%  Centers of input membership functions for the e universe of
% discourse for  of fuzzy controller (a vector of centers)
ce=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/ge);

% Centers of input membership functions for the c universe of
% discourse for  of fuzzy controller (a vector of centers)
cc=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/gc);
gf=0;

fuzzyrules=[-1  -1   -1   -1  -1    -1  -0.8 -0.6 -0.4 -0.2  0;
	    -1  -1   -1   -1  -1   -0.8 -0.6 -0.4 -0.2   0  0.2;
	    -1  -1   -1   -1  -0.8 -0.6 -0.4 -0.2   0   0.2 0.4;
	    -1  -1   -1  -0.8 -0.6 -0.4 -0.2   0   0.2  0.4 0.6;
	    -1  -1  -0.8 -0.6 -0.4 -0.2   0   0.2  0.4  0.6 0.8;
	    -1 -0.8 -0.6 -0.4 -0.2   0   0.2  0.4  0.6  0.8  1;
	  -0.8 -0.6 -0.4 -0.2   0   0.2  0.4  0.6  0.8   1   1;
	  -0.6 -0.4 -0.2   0   0.2  0.4  0.6  0.8   1    1   1;
	  -0.4 -0.2   0   0.2  0.4  0.6  0.8   1    1    1   1;
	  -0.2  0    0.2  0.4  0.6  0.8   1    1    1    1   1;
	   0   0.2   0.4  0.6  0.8   1    1    1    1    1   1]*gu*gf;

% Next, we define some parameters for the fuzzy inverse model

gye=1/2;,gyc=1/2;
	% Scaling gains for the error and change in error for
	% the inverse model
	% These are tuned to improve the performance of the FMRLC
gp=0.2;

numye=11; 	% Number of input membership functions for the ye
			% universe of discourse
numyc=11; 	% Number of input membership functions for the yc
			% universe of discourse

wye=0.2*(1/gye);	% Sets the width of the membership functions for
					% ye from center to extremes
wyc=0.2*(1/gyc);	% Sets the width of the membership functions for
					% yc from center to extremes
invbase=0.4*gp; % Sets the base of the output membership functions
				% for the inverse model

% Place centers of inverse model membership functions
% For error input for learning Mechanism
cye=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/gye);

% For change in error input for learning mechanism
cyc=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1]*(1/gyc);

% The next matrix contains the rule-base matrix for the fuzzy
% inverse model.  Notice that for simplicity we choose it to have
% the same structure as the rule-base for the fuzzy controller.
% While this will work for the control of the simple first order
% linear system for many nonlinear systems a different structure
% will be needed for the rule-base.  Again, the entries are
% the centers of the output membership functions, but now for
% the fuzzy inverse model.

inverrules=[-1  -1   -1   -1  -1    -1  -0.8 -0.6 -0.4 -0.2  0;
	    -1  -1   -1   -1  -1   -0.8 -0.6 -0.4 -0.2   0  0.2;
	    -1  -1   -1   -1  -0.8 -0.6 -0.4 -0.2   0   0.2 0.4;
	    -1  -1   -1  -0.8 -0.6 -0.4 -0.2   0   0.2  0.4 0.6;
	    -1  -1  -0.8 -0.6 -0.4 -0.2   0   0.2  0.4  0.6 0.8;
	    -1 -0.8 -0.6 -0.4 -0.2   0   0.2  0.4  0.6  0.8  1;
	  -0.8 -0.6 -0.4 -0.2   0   0.2  0.4  0.6  0.8   1   1;
	  -0.6 -0.4 -0.2   0   0.2  0.4  0.6  0.8   1    1   1;
	  -0.4 -0.2   0   0.2  0.4  0.6  0.8   1    1    1   1;
	  -0.2  0    0.2  0.4  0.6  0.8   1    1    1    1   1;
	   0   0.2   0.4  0.6  0.8   1    1    1    1    1   1]*gp;

% Next, we set up some parameters/variables for the
% knowledge-base modifier

d=1;
% This sets the number of steps the knowledge-base modifier looks
% back in time. For this program it must be an integer
% less than or equal to 10 (but this is easy to make larger)

% The next four vectors are used to store the information about
% which rules were on 1 step in the past, 2 steps in the past, ....,
% 10 steps in the past (so that picking 0<= d <= 10 can be used).

meme_int=[0 0 0 0 0 0 0 0 0 0];
	% sets up the vector to store up to 10 values of e_int
meme_count=[0 0 0 0 0 0 0 0 0 0];
	% sets up the vector to store up to 10 values of e_count
memc_int=[0 0 0 0 0 0 0 0 0 0];
	% sets up the vector to store up to 10 values of c_int
memc_count=[0 0 0 0 0 0 0 0 0 0];
	% sets up the vector to store up to 10 values of c_count

%
%  Next, we intialize the simulation of the closed-loop system.
%

k_p=1;	  % The numerator of the plant.  Change this value to study
		  % the ability of the FMRLC to control other plants.  Also,
		  % you can make this a time-varying parameter.
zeta_p=.707;
       % Damping ratio for the second order plant (could change this
	   % to see how the system will adapt to it)
w_p=1; % Undamped natural frequency for the plant (could change this
	   % to see how the system will adapt to it)
k_r=1;
 % The numerator of the reference model.  Change this value to study
 % the ability of the FMRLC to meet other performance specifications.
a_r=1;
    % The value of -a_r is the pole position for the reference model.
	% Change this value to study the ability of the FMRLC to meet
	% other performance specifications (e.g., a faster response).

t=0; 		% Reset time to zero
index=1;	% This is time's index (not time, its index).
tstop=64;	% Stopping time for the simulation (in seconds)
step=0.01;  % Integration step size
x=[0;0];	% Intial condition on state	of the plant

% Need a state space representation for the plant.  Since our
% plant is linear we use the standard form of xdot=Ax+Bu, y=Cx+Du
% Matrix A of state space representation of plant

A=[0 1;
   -w_p^2 -2*zeta_p*w_p];
B=[0; 1];	    % Matrix B of state space representation of plant
C=[k_p 0];	    % Matrix C of state space representation of plant


%
% Next, we start the simulation of the system.  This is the main
% loop for the simulation of the FMRLC.
%
while t <= tstop
	y(index)=C*x;     % Output of the plant

% Next, we define the reference input r as a sine wave

r(index)=sin(.6*t);



ym(index)=(1/(2+a_r*step))*((2-a_r*step)*ymold+...
                                    k_r*step*(r(index)+rold));

ymold=ym(index);
rold=r(index);
	% This saves the past value of the ym (r) so that we can use it
	% the next time around the loop

% Now that we have simulated the next step for the plant and reference
% model we will focus on the two fuzzy components.

% First, for the given fuzzy controller inputs we determine
% the extent at which the error membership functions
% of the fuzzy controller are on (this is the fuzzification part).

c_count=0;,e_count=0;   % These are used to count the number of
						% non-zero mf certainitie of e and c
e=r(index)-y(index);
			% Calculates the error input for the fuzzy controller
c=(e-eold)/step;
	% Calculates the change in error input for the fuzzy controller
eold=e;
% Saves the past value of e for use in the next time through the
% loop

% The following if-then structure fills the vector mfe
% with the certainty of each membership fucntion of e for the
% current input e

	if e<=ce(1)		% Takes care of saturation of the left-most
					% membership function
         mfe=[1 0 0 0 0 0 0 0 0 0 0]; % i.e., the only one on is the
         							  %left-most one
	 e_count=e_count+1;,e_int=1; 	  %  One mf on, it is the
	 								  %left-most one.
	elseif e>=ce(nume)				  % Takes care ofsaturation
									  %of the right-most mf
	 mfe=[0 0 0 0 0 0 0 0 0 0 1];
	 e_count=e_count+1;,e_int=nume; % One mf on, it is the
	 								%right-most one
	else      % In this case the input is on the middle part of the
			  % universe of discourse for e
			  % Next, we are going to cycle through the mfs to
			  % find all that are on
	   for i=1:nume
		 if e<=ce(i)
		  mfe(i)=max([0 1+(e-ce(i))/we]);
		  				% In this case the input isto the
		  				% left of the center ce(i)and we compute
						% the value of the mfcentered at ce(i)
						% for this input e
			if mfe(i)~=0
				% If the certainty is not equal to zerothen say
				% that have one mf on by incrementing our count
			 e_count=e_count+1;
			 e_int=i;	% This term holds the index last entry
			 			% with a non-zero term
			end
		 else
		  mfe(i)=max([0,1+(ce(i)-e)/we]);
		  						% In thiscase the input is to the
		  						% right ofthe center ce(i)
			if mfe(i)~=0
			 e_count=e_count+1;
			 e_int=i;  % This term holds the index of the
			 		   % last entry with a non-zero term