📄 fuzzy_control_system_for_tanker_ship.m
字号:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -