%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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 nonlinearity bbar=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 controller nume=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 functions we=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 discourse base=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 model gye=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 FMRLC gp=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 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 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 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 % 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 zero index=1; % This is time's index (not time, its index). tstop=20000; % Stopping time for the simulation (in seconds) step=1; % Integration step size T=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 trajectory yeold=0; % Intial condition used to calculate yc ymold=0; % Initial condition for the first order reference model x=[0;0;0]; % First, set the state to be a vector x(1)=0; % Set the initial heading to be zero x(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 deg if t>=100, psi_r(index)=45*(pi/180); end % Request heading of 45 deg if t>=1500, psi_r(index)=0; end % Request heading of 0 deg if t>=3000, psi_r(index)=45*(pi/180); end % Request heading of -45 deg if t>=4500, psi_r(index)=0; end % Request heading of 0 deg if t>=6000, psi_r(index)=45*(pi/180); end % Request heading of 45 deg if t>=7500, psi_r(index)=0; end % Request heading of 0 deg if t>=9000, psi_r(index)=45*(pi/180); end % Request heading of 45 deg if t>=10500, psi_r(index)=0; end % Request heading of 0 deg if t>=12000, psi_r(index)=45*(pi/180); end % Request heading of -45 deg if t>=13500, psi_r(index)=0; end % Request heading of 0 deg if t>=15000, psi_r(index)=45*(pi/180); end % Request heading of 45 deg if t>=16500, psi_r(index)=0; end % Request heading of 0 deg if t>=18000, psi_r(index)=45*(pi/180); end % Request heading of 45 deg if 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 % FMRLC counter=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 % 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 certainities of e and c e(index)=psi_r(index)-psi(index); % Calculates the error input for the fuzzy controller c(index)=(e(index)-eold)/T; % Sets the value of c eold=e(index); % Save the past value of e for use in the above % computation the next time around 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. We use triangular membership functions. if e(index)<=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(index)>=ce(nume) % Takes care of saturation % 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(index)<=ce(i) mfe(i)=max([0 1+(e(index)-ce(i))/we]); % In this case the input is to the % left of the center ce(i) and we compute % the value of the mf centered at ce(i) % for this input e if mfe(i)~=0 % If the certainty is not equal to zero then 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(index))/we]); % In this case the input is to the % right of the 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 end end end end % The following if-then structure fills the vector mfc with the % certainty of each membership fucntion of the c % for its current value (to understand this part of the code see the above % similar code for computing mfe). Clearly, it could be more efficient to % make a subroutine that performs these computations for each of % the fuzzy system inputs. if c(index)<=cc(1) % Takes care of saturation of left-most mf mfc=[1 0 0 0 0 0 0 0 0 0 0]; c_count=c_count+1; c_int=1; elseif c(index)>=cc(numc) % Takes care of saturation of the right-most mf mfc=[0 0 0 0 0 0 0 0 0 0 1]; c_count=c_count+1; c_int=numc; else for i=1:numc if c(index)<=cc(i) mfc(i)=max([0,1+(c(index)-cc(i))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end else mfc(i)=max([0,1+(cc(i)-c(index))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end end end end % The next two loops calculate the crisp output using only the non- % zero premise of error,e, and c. This cuts down computation time % since we will only compute the contribution from the rules that % are on (i.e., a maximum of four rules for our case). The minimum % is used for the premise (and implication for the center-of-gravity % defuzzification case). num=0; den=0; for k=(e_int-e_count+1):e_int % Scan over e indices of mfs that are on for l=(c_int-c_count+1):c_int % Scan over c indices of mfs that are on prem=min([mfe(k) mfc(l)]); % Value of premise membership function % This next calculation of num adds up the numerator for the % center of gravity defuzzification formula. rules(k,l) is the output center % for the rule. base*(prem-(prem)^2/2 is the area of a symmetric % triangle with base width "base" and that is chopped off at % a height of prem (since we use minimum to represent the % implication). Computation of den is similar but without % rules(k,l). num=num+rules(k,l)*base*(prem-(prem)^2/2); den=den+base*(prem-(prem)^2/2); % To do the same computations, but for center-average defuzzification, % use the following lines of code rather than the two above (notice % that in this case we did not use any information about the output % membership function shapes, just their centers; also, note that % the computations are slightly simpler for the center-average defuzzificaton): % num=num+rules(k,l)*prem; % den=den+prem; end end delta(index)=num/den; % Crisp output of fuzzy controller that is the input % to the plant. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, we perform computations for the fuzzy inverse model. ye(index)=ym(index)-psi(index); % Calculates ye yc(index)=(ye(index)-yeold)/T; % Calculates yc yeold=ye(index); % Saves the value of ye for use the % next time ye_count=0;,yc_count=0; % Counts the non-zero certainties % of ye and yc similar to how we did % for the fuzzy controller % The following if-then structure fills the vector mfye with the % certainty of each membership fucntion of ye (similar to the % fuzzy controller). Notice that we use the same number of % input membership functions as in the fuzzy controller - % just for convenience - it does not have to be this way if ye(index)<=cye(1) % Takes care of saturation mfye=[1 0 0 0 0 0 0 0 0 0 0]; ye_count=ye_count+1;,ye_int=1; elseif ye(index)>=cye(numye) % Takes care of saturation mfye=[0 0 0 0 0 0 0 0 0 0 1]; ye_count=ye_count+1;,ye_int=numye; else for i=1:numye if ye(index)<=cye(i) mfye(i)=max([0 1+(ye(index)-cye(i))/wye]); if mfye(i)~=0 ye_count=ye_count+1; ye_int=i; % This term holds last entry with % a non-zero term end else mfye(i)=max([0,1+(cye(i)-ye(index))/wye]); if mfye(i)~=0 ye_count=ye_count+1; ye_int=i; % This term holds last entry with % a non-zero term end end end end % The following if-then structure fills the vector mfyc with the % certainty of each membership fucntion of yc if yc(index)<=cyc(1) mfyc=[1 0 0 0 0 0 0 0 0 0 0]; yc_count=yc_count+1;,yc_int=1; elseif yc(index)>=cyc(numyc) mfyc=[0 0 0 0 0 0 0 0 0 0 1]; yc_count=yc_count+1;,yc_int=numyc; else for i=1:numyc if yc(index)<=cyc(i) mfyc(i)=max([0 1+(yc(index)-cyc(i))/wyc]); if mfyc(i)~=0 yc_count=yc_count+1; yc_int=i; end else mfyc(i)=max([0,1+(cyc(i)-yc(index))/wyc]); if mfyc(i)~=0 yc_count=yc_count+1; yc_int=i; end end end end % These for loops calculate the crisp output of the inverse model % using only the non-zero premise of error,ye, and change in % error,yc. This cuts down computation time (similar to the % fuzzy controller). Minimum is used for the premise and the % implication. To understand this part of the code see the % similar code for the fuzzy controller (of course you could use % center-average defuzzification). invnum=0; invden=0; for k=(ye_int-ye_count+1):ye_int for l=(yc_int-yc_count+1):yc_int prem=min([mfye(k) mfyc(l)]); invnum=invnum+inverrules(k,l)*invbase*(prem-(prem)^2/2); invden=invden+invbase*(prem-(prem)^2/2); end end % Next we compute the output of the fuzzy inverse model. if gp==0, p(index)=0; % Allow us to turn off the learning else p(index)=invnum/invden; % Crisp output of inverse model if abs(p(index))<0.01*gp, % Robustification term which is used p(index)=0; % to stop the adaptation when the fuzzy end % controller is close to where it should % be. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, we modify the centers of ouput membership functions % associated with the rules that were on d steps ago % The indexing sheme is similar to the one used in the % computation of the outputs of the fuzzy controller and % fuzzy inverse model. Notice that in the early steps these % two loops would not be executed since we initialized % meme_int, meme_count, memc_int, and memc_count with all zeros. % This implies that there are no updates until d executions of % the controller (e.g., for our case it executes the controller % at the first step and then on the next controller execution which % occurs at 10 sec. the rules will be modified). for k=(meme_int(d)-meme_count(d)+1):meme_int(d) for l=(memc_int(d)-memc_count(d)+1):memc_int(d) rules(k,l)=rules(k,l)+p(index); end end % Note: Here, do not save the entire rule base for the past d steps % so this implies that, e.g., for a large d value, between d steps ago % and the present time the rules could be modified and so the way this % is coded it would modify the one that was modified in this case, less % than d steps ago. To avoid this you would have to save all the rules % over the past d steps and when you modify you may be undoing what was % done on a previous step - and it is not clear that would be good either. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Store values of pointers/counters for rules that were on at this step % Next we will save the number of mfs that are on and the pointer % e_int as to which rules were on. This vector of length % 10 saves the last 10 values of e_count and e_int as time % progresses (hence, it is a moving window). % These will be used by the FMRLC knowledge-base modifier. meme_count=[e_count meme_count(1:9)]; meme_int=[e_int meme_int(1:9)]; % Next we will save the number of mfs that are on and the pointer % c_int as to which rules were on. This vector of length 10 % saves the last 10 values of c_count and e_int as time progresses % (hence, it is a moving window). These will be used by the FMRLC % knowledge-base modifier. memc_count=[c_count memc_count(1:9)]; memc_int=[c_int memc_int(1:9)]; else % This goes with the "if" statement to check if the counter=10 % so the next lines up to the next "end" statement are executed % whenever counter is not equal to 10 % Now, even though we do not compute the FMRLC at each % time instant, we do want to save the data at its inputs and output at % each time instant for the sake of plotting. Hence, we need to % compute these here (note that we simply hold the values): e(index)=e(index-1); c(index)=c(index-1); delta(index)=delta(index-1); ye(index)=ye(index-1); yc(index)=yc(index-1); p(index)=p(index-1); ym(index)=ym(index-1); end % This is the end statement for the "if counter=10" statement %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, the Runge-Kutta equations are used to find the next state. % Clearly, it would be better to use a Matlab "function" for % F (but here we do not, so we can have only one program). % (This does not implement the exact equations in the book but % if you think about the problem a bit you will see that these % really should be accurate enough. Here, we are not calculating % intermediate values of the controller output, only % one per integration step; we do this simply to save some % computations. Similarly, for the reference input.) time(index)=t; % First, we define a wind disturbance against the body of the ship % that has the effect of pressing water against the rudder %w(index)=0.5*(pi/180)*sin(2*pi*0.001*t); % This is an additive sine disturbance to % the rudder input. It is of amplitude of % 0.5 deg. and its period is 1000sec. %delta(index)=delta(index)+w(index); % Next, implement the nonlinearity where the rudder angle is saturated % at +-80 degrees if delta(index) >= 80*(pi/180), delta(index)=80*(pi/180); end if delta(index) <= -80*(pi/180), delta(index)=-80*(pi/180); end % The next line is used in place of the line following it to % change the speed of the ship if t>=1000000, %if t>=9000, % This switches the ship speed (unrealistically fast) u=3; % A lower speed else u=5; end % Next, we change the parameters of the ship to tanker to reflect % changing loading conditions (note that we simulate as if % the ship is loaded while moving, but we only change the parameters % while the heading is zero so that it is then similar to re-running % the simulation, i.e., starting the tanker operation at different % times after loading/unloading has occurred). % The next line is used in place of the line following it to keep % "ballast" conditions throughout the simulation if t>=1000000, %if t>=9000, % This switches the parameters in the middle of the simulation K_0=0.83; % These are the parameters under "full" conditions tau_10=-2.88; tau_20=0.38; tau_30=1.07; else K_0=5.88; % These are the parameters under "ballast" conditions tau_10=-16.91; tau_20=0.45; tau_30=1.43; end % The following parameters are used in the definition of the tanker model: K=K_0*(u/ell); tau_1=tau_10*(ell/u); tau_2=tau_20*(ell/u); tau_3=tau_30*(ell/u); % Now, for the first step, we set the initial condition for the % third state x(3). if t==0, x(3)=-(K*tau_3/(tau_1*tau_2))*delta(index); end % Next, we use the formulas to implement the Runge-Kutta method F=[ x(2) ; x(3)+ (K*tau_3/(tau_1*tau_2))*delta(index) ; -((1/tau_1)+(1/tau_2))*(x(3)+ (K*tau_3/(tau_1*tau_2))*delta(index))-... (1/(tau_1*tau_2))*(abar*x(2)^3 + bbar*x(2)) + (K/(tau_1*tau_2))*delta(index) ]; k1=step*F; xnew=x+k1/2; F=[ xnew(2) ; xnew(3)+ (K*tau_3/(tau_1*tau_2))*delta(index) ; -((1/tau_1)+(1/tau_2))*(xnew(3)+ (K*tau_3/(tau_1*tau_2))*delta(index))-... (1/(tau_1*tau_2))*(abar*xnew(2)^3 + bbar*xnew(2)) + (K/(tau_1*tau_2))*delta(index) ]; k2=step*F; xnew=x+k2/2; F=[ xnew(2) ; xnew(3)+ (K*tau_3/(tau_1*tau_2))*delta(index) ; -((1/tau_1)+(1/tau_2))*(xnew(3)+ (K*tau_3/(tau_1*tau_2))*delta(index))-... (1/(tau_1*tau_2))*(abar*xnew(2)^3 + bbar*xnew(2)) + (K/(tau_1*tau_2))*delta(index) ]; k3=step*F; xnew=x+k3; F=[ xnew(2) ; xnew(3)+ (K*tau_3/(tau_1*tau_2))*delta(index) ; -((1/tau_1)+(1/tau_2))*(xnew(3)+ (K*tau_3/(tau_1*tau_2))*delta(index))-... (1/(tau_1*tau_2))*(abar*xnew(2)^3 + bbar*xnew(2)) + (K/(tau_1*tau_2))*delta(index) ]; k4=step*F; x=x+(1/6)*(k1+2*k2+2*k3+k4); % Calculated next state t=t+step; % Increments time index=index+1; % Increments the indexing term so that % index=1 corresponds to time t=0. counter=counter+1; % Indicates that we computed one more integration step end % This end statement goes with the first "while" statement % in the program so when this is complete the simulation is done. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, we provide plots of data from the simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % First, we convert from rad. to degrees psi_r=psi_r*(180/pi); psi=psi*(180/pi); delta=delta*(180/pi); e=e*(180/pi); c=c*(180/pi); ym=ym*(180/pi); ye=ye*(180/pi); % Next, we provide plots figure(1) clf subplot(311) plot(time,psi,'k-',time,ym,'k--',time,psi_r,'k-.') zoom grid on title('Ship heading (solid) and desired ship heading (dashed), deg.') subplot(312) plot(time,delta,'k-') zoom grid on title('Rudder angle, output of fuzzy controller (input to the ship), deg.') subplot(313) plot(time,p,'k-') zoom grid on xlabel('Time (sec)') title('Fuzzy inverse model output (nonzero values indicate adaptation)') figure(2) clf subplot(211) plot(time,e,'k-') zoom grid on title('Ship heading error between ship heading and desired heading, deg.') subplot(212) plot(time,c,'k-') zoom grid on xlabel('Time (sec)') title('Change in ship heading error, deg./sec') figure(3) clf subplot(211) plot(time,ye,'k-') zoom grid on title('Ship heading error between ship heading and reference model heading, deg.') subplot(212) plot(time,yc,'k-') zoom grid on xlabel('Time (sec)') title('Change in heading error between output and reference model, deg./sec') end % This ends the if statement (on flag1) on whether you want to do a simulation % or just see the control surface %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Next, provide a plot of the fuzzy controller surface: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Request input from the user to see if they want to see the tuned % controller mapping: flag2=input('Do you want to see the nonlinear \n mapping implemented by the tuned fuzzy \n controller? \n (type 1 for yes and 0 for no) '); if flag2==1, % First, compute vectors with points over the whole range of % the fuzzy controller inputs plus 20% over the end of the range % and put 100 points in each vector e_input=(-(1/g1)-0.2*(1/g1)):(1/100)*(((1/g1)+0.2*(1/g1))-(-(1/g1)-... 0.2*(1/g1))):((1/g1)+0.2*(1/g1)); ce_input=(-(1/g2)-0.2*(1/g2)):(1/100)*(((1/g2)+0.2*(1/g2))-(-(1/g2)-... 0.2*(1/g2))):((1/g2)+0.2*(1/g2)); % Next, compute the fuzzy controller output for all these inputs for jj=1:length(e_input) for ii=1:length(ce_input) c_count=0;,e_count=0; % These are used to count the number of % non-zero mf certainities of e and c % The following if-then structure fills the vector mfe % with the certainty of each membership fucntion of e for the % current input e. We use triangular membership functions. if e_input(jj)<=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_input(jj)>=ce(nume) % Takes care of saturation % 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_input(jj)<=ce(i) mfe(i)=max([0 1+(e_input(jj)-ce(i))/we]); % In this case the input is to the % left of the center ce(i) and we compute % the value of the mf centered at ce(i) % for this input e if mfe(i)~=0 % If the certainty is not equal to zero then 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_input(jj))/we]); % In this case the input is to the % right of the 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 end end end end % The following if-then structure fills the vector mfc with the % certainty of each membership fucntion of the c % for its current value. if ce_input(ii)<=cc(1) % Takes care of saturation of left-most mf mfc=[1 0 0 0 0 0 0 0 0 0 0]; c_count=c_count+1; c_int=1; elseif ce_input(ii)>=cc(numc) % Takes care of saturation of the right-most mf mfc=[0 0 0 0 0 0 0 0 0 0 1]; c_count=c_count+1; c_int=numc; else for i=1:numc if ce_input(ii)<=cc(i) mfc(i)=max([0,1+(ce_input(ii)-cc(i))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end else mfc(i)=max([0,1+(cc(i)-ce_input(ii))/wc]); if mfc(i)~=0 c_count=c_count+1; c_int=i; % This term holds last entry % with a non-zero term end end end end % The next loops calculate the crisp output using only the non- % zero premise of error,e, and c. num=0; den=0; for k=(e_int-e_count+1):e_int % Scan over e indices of mfs that are on for l=(c_int-c_count+1):c_int % Scan over c indices of mfs that are on prem=min([mfe(k) mfc(l)]); % Value of premise membership function % This next calculation of num adds up the numerator for the % center of gravity defuzzification formula. num=num+rules(k,l)*base*(prem-(prem)^2/2); den=den+base*(prem-(prem)^2/2); % To do the same computations, but for center-average defuzzification, % use the following lines of code rather than the two above: % num=num+rules(k,l)*prem; % den=den+prem; end end delta_output(ii,jj)=num/den; % Crisp output of fuzzy controller that is the input % to the plant. end end % Convert from radians to degrees: delta_output=delta_output*(180/pi); e_input=e_input*(180/pi); ce_input=ce_input*(180/pi); % Plot the data figure(4) clf surf(e_input,ce_input,delta_output); view(145,30); colormap(white); xlabel('Heading error (e), deg.'); ylabel('Change in heading error (c), deg.'); zlabel('Fuzzy controller output (\delta), deg.'); title('FMRLC-tuned fuzzy controller mapping between inputs and output'); % Next, print out to the screen the rule base (output centers) rules end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % End of program % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%