fc_tanker.m

一个用MATLAB编写的优化控制工具箱

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fuzzy Control System for a Tanker Ship
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% By: Kevin Passino 
% Version: 4/15/99
%
% Notes: This program has evolved over time and uses programming 
% ideas of Andrew Kwong, Scott Brown, and Brian Klinehoffer.
%
% This program simulates a fuzzy control system 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).  We want the 
% tanker ship heading (psi) to track the reference input heading
% (psi_r).  We simulate the tanker as a continuous time system 
% that is controlled by a fuzzy controller that is implemented on
% a digital computer with a sampling interval of T.  
%
% This program can be used to illustrate:
%  - How to code a fuzzy controller (for two inputs and one output, 
%    illustrating some approaches to simplify the computations, for
%    triangular membership functions, and either center-of-gravity or
%    center-average defuzzification).
%  - How to tune the input and output gains of a fuzzy controller.
%  - How changes in plant conditions ("ballast" and "full") 
%    can affect performance.
%  - How sensor noise (heading sensor noise), plant disturbances
%    (wind hitting the side of the ship), and plant operating
%    conditions (ship speed) can affect performance.
%  - How improper choice of the scaling gains can result in 
%    oscillations (limit cycles).
%  - How an improper choice of the scaling gains (or rule base) can
%    result in an unstable system.
%  - The shape of the nonlinearity implemented by the fuzzy controller
%    by plotting the input-output map of the fuzzy controller.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear		% Clear all variables in memory

% Initialize ship parameters 
% (can test two conditions, "ballast" or "full"):

ell=350;			% Length of the ship (in meters)
u=5;				% Nominal speed (in meters/sec)
%u=3;               % A lower speed where the ship is more difficult to control
abar=1;             % Parameters for nonlinearity
bbar=1;

% The parameters for the tanker under "ballast" conditions 
% (a heavy ship) are:

K_0=5.88;
tau_10=-16.91;
tau_20=0.45;
tau_30=1.43;

% The parameters for the tanker under "full" conditions (a ship
% that weighs less than one under "ballast" conditions) are:

%K_0=0.83;
%tau_10=-2.88;
%tau_20=0.38;
%tau_30=1.07;

% Some other parameters are:

K=K_0*(u/ell);
tau_1=tau_10*(ell/u);
tau_2=tau_20*(ell/u);
tau_3=tau_30*(ell/u);

% 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
% Tuning:
g1=1/pi;,g2=200;,g0=8*pi/18; % Try to reduce the overshoot
g1=2/pi;,g2=250;,g0=8*pi/18; % Try to speed up the response a bit but to do this 
							 % have to raise g2 a bit to avoid overshoot. Take these
							 % as "good" tuned values.  
%g1=2/pi;,g2=0.000001;,g0=2000*pi/18; % Values tuned to get an oscillation (limit
									  % cycle) for COG, ballast,
									  % and nominal speed with no sensor 
									  % noise or rudder disturbance):
									  % g1: Leave as before
								      % g2: Essentially turn off the derivative gain
									  %     since this help induce an oscillation
									  % g0: Make this big to force the limit cycle
									  % In this case simulate for 16,000 sec.
%g1=2/pi;,g2=250;,g0=-8*pi/18;       % Values tuned to get an instability
									% g0: Make this negative so that when there
									% is an error the rudder will drive the 
									% heading in the direction to increase the error

% 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.  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.

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]*g0;

% Now, you can proceed to do the simulation or simply view the nonlinear
% surface generated by the fuzzy controller. 

flag1=input('\n Do you want to simulate the \n fuzzy control system \n 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=4000;	% 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.  
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.
psi_r_old=0; % Initialize the reference trajectory


% 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<100, 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>2000, psi_r(index)=0; end      			% Then request heading of 0 deg
%if t>4000, psi_r(index)=45*(pi/180); end     % Then request heading of 45 deg
%if t>6000, psi_r(index)=0; end      			% Then request heading of 0 deg
%if t>8000, psi_r(index)=45*(pi/180); end     % Then request heading of 45 deg
%if t>10000, psi_r(index)=0; end      			% Then request heading of 0 deg
%if t>12000, psi_r(index)=45*(pi/180); end     % Then request heading of 45 deg

% Next, suppose that there is sensor noise for the heading sensor with that is
% additive, with a uniform distribution on [- 0.01,+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 
				   % fuzzy controller

counter=0; 			% First, reset the counter

% 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