%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fuzzy Model Reference Learning Control (FMRLC) System for a Tanker Ship%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% By: Kevin Passino % Version: 4/26/01%% Notes: This program has evolved over time and uses programming % ideas of Andrew Kwong, Jeff Layne, and Brian Klinehoffer.%% This program simulates an FMRLC for a tanker% ship.  It has a fuzzy controller with two inputs, the error% in the ship heading (e) and the change in that error (c).  The output% of the fuzzy controller is the rudder input (delta).  The FMRLC % adjusts the fuzzy controller to try to get tanker ship heading (psi) % to track the output of a "reference model" (psi_m) that has as an % input the reference input heading (psi_r). We simulate the tanker % as a continuous time system that is controlled by an FMRLC that % is implemented on a digital computer with a sampling interval of T=1 sec.  %% This program can be used to illustrate:%  - How to code an FMRLC (for two inputs and one output, %    for the controller and "fuzzy inverse model").%  - How to tune the input and output gains of an FMRLC.%  - How changes in plant conditions ("ballast" and "full" and speed) %    can affect performance and how the FMRLC can adapt to improve%    performance when there are such changes.%  - How the effects of sensor noise (heading sensor noise) and plant %    disturbances (wind hitting the side of the ship) can be reduced%    over the case of non-adaptive fuzzy control.%  - The shape of the nonlinearity synthesized by the FMRLC%    by plotting the input-output map of the fuzzy controller at the%    end of the simulation (and providing the output membership function%    centers).%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear		% Clear all variables in memory% Initialize ship parameters ell=350;			% Length of the ship (in meters)abar=1;             % Parameters for nonlinearitybbar=1;% Define the reference model (we use a first order transfer function % k_r/(s+a_r)):a_r=1/150;k_r=1/150;% Initialize parameters for the fuzzy controllernume=11; 	% Number of input membership functions for the e 			% universe of discourse (can change this but must also			% change some variables below if you make such a change)numc=11; 	% Number of input membership functions for the c 			% universe of discourse (can change this but must also			% change some variables below if you make such a change)% Next, we define the scaling gains for tuning membership functions for % universes of discourse for e, change in e (what we call c) and % delta.  These are g1, g2, and g0, respectively% These can be tuned to try to improve the performance. % First guess: g1=1/pi;,g2=100;,g0=8*pi/18; % Chosen since: 		% g1: The heading error is at most 180 deg (pi rad) 		% g2: Just a guess - that ship heading will change at most		%     by 0.01 rad/sec (0.57 deg/sec)		% g0: Since the rudder is constrained to move between +-80 deg% "Good" tuned values:g1=2/pi;,g2=250;,g0=8*pi/18; % Next, define some parameters for the membership functionswe=0.2*(1/g1);	  	% we is half the width of the triangular input membership 	% function bases (note that if you change g0, 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/g2);        	% Similar to we but for the c universe of discoursebase=0.4*g0;      	% 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 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/g1);% Centers of input membership functions for the c universe of% discourse 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/g2);% This next matrix specifies the rules of the fuzzy controller.  % The entries are the centers of the output membership functions.  % This choice represents just one guess on how to synthesize % the fuzzy controller.  Notice the regularity % of the pattern of rules (basiscally we are using a type of % saturated index adding).  Notice that it is scaled by g0, the % output scaling factor, since it is a normalized rule base.% The rule base can be tuned to try to improve performance.% The gain gf is a gain that can be set to one to initialize % the rule base with an initial guess at the fuzzy controller to be% synthesized or to zero to initialize the rule base to all zeros% (i.e., all centers at zero so the rules all say to put zero% into the plant).gf=1;rules=[1    1     1     1    1     1    0.8   0.6   0.3    0.1     0;  	   1    1     1     1    1   0.8    0.6   0.3   0.1    0     -0.1;	   1    1     1     1   0.8  0.6    0.3   0.1   0     -0.1   -0.3;	   1    1     1   0.8   0.6  0.3    0.1    0    -0.1   -0.3  -0.6;	   1    1   0.8   0.6   0.3  0.1     0   -0.1   -0.3   -0.6  -0.8;	   1  0.8   0.6   0.3   0.1   0    -0.1  -0.3   -0.6   -0.8   -1;	 0.8  0.6   0.3   0.1    0   -0.1  -0.3  -0.6   -0.8   -1     -1;	 0.6  0.3   0.1    0   -0.1  -0.3  -0.6  -0.8   -1     -1     -1;	 0.3  0.1    0   -0.1  -0.3  -0.6  -0.8   -1    -1     -1     -1;	 0.1   0   -0.1  -0.3  -0.6  -0.8   -1    -1    -1     -1     -1;	   0 -0.1  -0.3  -0.6  -0.8  -1     -1    -1    -1     -1     -1]*gf*g0;% Next, we define some parameters for the fuzzy inverse modelgye=2/pi;,gyc=10;    % Scaling gains for the error and change in error for 					% the inverse model					% These are tuned to improve the performance of the FMRLCgp=0.4;	% Scaling gain for the output of inverse model.  If you let gp=0 then% you turn off the learning mechanism and hence the FMRLC reduces to % a direct fuzzy controller that is not tuned.  Hence, from this % program you can also see how to code a non-adaptive fuzzy controller.% If gp is large then generally large updates will be made to the % fuzzy controller.  You should think of gp as an "adaptation gain."numye=11; 	% Number of input membership functions for the ye 			% universe of discoursenumyc=11; 	% Number of input membership functions for the yc 			% universe of discoursewye=0.2*(1/gye);	% Sets the width of the membership functions for 					% ye from center to extremeswyc=0.2*(1/gyc);	% Sets the width of the membership functions for 					% yc from center to extremesinvbase=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 mechanismcye=[-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 mechanismcyc=[-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 tanker % 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 modifierd=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_intmeme_count=[0 0 0 0 0 0 0 0 0 0];  	% sets up the vector to store up to 10 values of e_countmemc_int=[0 0 0 0 0 0 0 0 0 0];  	% sets up the vector to store up to 10 values of c_intmemc_count=[0 0 0 0 0 0 0 0 0 0];  	% sets up the vector to store up to 10 values of c_count% Now, you can proceed to do the simulation or simply view the nonlinear% surface generated by the fuzzy controller that is now fully defined.flag1=input('Do you want to simulate the \n FMRLC for the tanker?  \n (type 1 for yes and 0 for no) ');if flag1==1, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next, we initialize the simulation:t=0; 		% Reset time to zeroindex=1;	% This is time's index (not time, its index).  tstop=20000;	% Stopping time for the simulation (in seconds)step=1;     % Integration step sizeT=10;		% The controller is implemented in discrete time and			% this is the sampling time for the controller.			% Note that the integration step size and the sampling			% time are not the same.  In this way we seek to simulate			% the continuous time system via the Runge-Kutta method and			% the discrete time fuzzy controller as if it were			% implemented by a digital computer.  Hence, we sample			% the plant output every T seconds and at that time			% output a new value of the controller output.counter=10;	% This counter will be used to count the number of integration			% steps that have been taken in the current sampling interval.			% Set it to 10 to begin so that it will compute a fuzzy controller			% output at the first step.			% For our example, when 10 integration steps have been			% taken we will then we will sample the ship heading			% and the reference heading and compute a new output			% for the fuzzy controller.  eold=0;     % Initialize the past value of the error (for use            % in computing the change of the error, c).  Notice            % that this is somewhat of an arbitrary choice since             % there is no last time step.  The same problem is            % encountered in implementation.  			psi_r_old=0; % Initialize the reference trajectoryyeold=0; 	 % Intial condition used to calculate ycymold=0; 	 % Initial condition for the first order reference modelx=[0;0;0];	% First, set the state to be a vector            x(1)=0;		% Set the initial heading to be zerox(2)=0;		% Set the initial heading rate to be zero.  			% We would also like to set x(3) initially but this			% must be done after we have computed the output			% of the fuzzy controller.  In this case, by			% choosing the reference trajectory to be 			% zero at the beginning and the other initial conditions			% as they are, and the fuzzy controller as designed			% we will know that the output of the fuzzy controller			% will start out at zero so we could have set 			% x(3)=0 here.  To keep things more general, however, 			% we set the intial condition immediately after 			% we compute the first controller output in the 			% loop below.			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next, we start the simulation of the system.  This is the main % loop for the simulation of fuzzy control system.while t <= tstop% First, we define the reference input psi_r (desired heading).if t>=0, psi_r(index)=0; end			    % Request heading of 0 degif t>=100, psi_r(index)=45*(pi/180); end     % Request heading of 45 degif t>=1500, psi_r(index)=0; end    			% Request heading of 0 degif t>=3000, psi_r(index)=45*(pi/180); end    % Request heading of -45 degif t>=4500, psi_r(index)=0; end    			% Request heading of 0 degif t>=6000, psi_r(index)=45*(pi/180); end     % Request heading of 45 degif t>=7500, psi_r(index)=0; end    			% Request heading of 0 degif t>=9000, psi_r(index)=45*(pi/180); end     % Request heading of 45 degif t>=10500, psi_r(index)=0; end    			% Request heading of 0 degif t>=12000, psi_r(index)=45*(pi/180); end    % Request heading of -45 degif t>=13500, psi_r(index)=0; end    			% Request heading of 0 degif t>=15000, psi_r(index)=45*(pi/180); end     % Request heading of 45 degif t>=16500, psi_r(index)=0; end    			% Request heading of 0 degif t>=18000, psi_r(index)=45*(pi/180); end     % Request heading of 45 degif t>=19500, psi_r(index)=0; end    			% Request heading of 0 deg% Next, suppose that there is sensor noise for the heading sensor with that is% additive, with a uniform distribution on +- 0.01 deg.%s(index)=0.01*(pi/180)*(2*rand-1);s(index)=0;					  % This allows us to remove the noise.psi(index)=x(1)+s(index);     % Heading of the ship (possibly with sensor noise).if counter == 10,  % When the counter reaches 10 then execute the 				   % FMRLCcounter=0; 			% First, reset the counter% Reference model calculations:% The reference model is part of the controller and to simulate it% we take the discrete equivalent of the% reference model to compute psi_m from psi_r (if you use% a continuous-time reference model you will have to augment % the state of the closed-loop system with the state(s) of the % reference model and hence update the state in the Runge-Kutta % equations).%% For the reference model we use a first order transfer function % k_r/(s+a_r) but we use the bilinear transformation where we % replace s by (2/step)(z-1)/(z+1), then find the z-domain % representation of the reference model, then convert this % to a difference equation:ym(index)=(1/(2+a_r*T))*((2-a_r*T)*ymold+...                                    k_r*T*(psi_r(index)+psi_r_old));ymold=ym(index);  psi_r_old=psi_r(index);	% This saves the past value of the ym and psi_r so that we can use it	% the next time around the loop	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Fuzzy controller calculations:% First, for the given fuzzy controller inputs we determine% the extent at which the error membership functions