fuzzcalc.m

用于汽车巡航控制系统的模糊控制算法

function [y, dxhatdc, dxhatdb, dxhatdgamma, dxhatdbeta] = ...
   FuzzCalc(nin, nclassin, nout, nclassout, b, c, bout, cout, rule, x, PartialFlag)

% function [y, dxhatdc, dxhatdb] = ...
%    FuzzCalc(nin, nclassin, nout, nclassout, b, c, bout, cout, rule, x, PartialFlag)
% Compute the defuzzified outputs corresponding to the fuzzy inputs.a
% Also compute the derivatives of the outputs with respect to the membership
% function parameters (if PartialFlag = 1).  This routine uses triangular
% membership functions.

% INPUTS
%   nin = number of fuzzy input variables.
%   nclassin = nin-element vector containing number of classes for each
%              fuzzy input variable.
%   nout = number of fuzzy output variables.
%   nclassout = nout-element vector containing number of classes for each
%               fuzzy output variable.
%   b = (nin x max(nclassin) x 2)-array containing the lower and upper
%       half widths of the input membership functions.
%       The first subscript specifies the input variable.
%       The second subscript specifies the fuzzy class number.
%       The third subscript specifies the lower half width or upper half width.
%   c = (nin x max(nclassin))-array containing the centroids of
%       the input membership functions.
%       The first subscript specifies the input variable.
%       The second subscript specifies the fuzzy class number.
%   bout = (nout x max(nclassout) x 2)-array containing the lower and upper
%          half widths of the output membership functions.
%   cout = (nout x max(nclassout))-array containing the centroids of
%          the output membership functions.
%   rule = fuzzy rules; an [nclassin(1) x nclassin(2) x nout] dimensional array
%          specififying which class the specified fuzzy output should be in when
%          input 1 and input 2 are in the specified classes.  For example, 
%          "rule(2,3,4)=5" means that if input 1 is in fuzzy class # 2 and input 2 is
%          in fuzzy class # 3, then output 4 will be in fuzzy class # 5.
%   x = inputs; an nin-element vector of crisp numbers.
%   PartialFlag = a flag (0 or 1) specifying whether this routine should
%                 calculate the derivatives of the fuzzy output with respect
%                 to the fuzzy membership function parameters.
% OUTPUTS
%   y = nout-element vector of outputs of the fuzzy system.
%   dxhatdc = derivaties of the fuzzy outputs with respect 
%             to the fuzzy input variable centroids.
%   dxhatdb = derivaties of the fuzzy outputs with respect 
%             to the fuzzy input variable half widths.
%   dxhatdgamma = derivaties of the fuzzy outputs with respect 
%                 to the fuzzy output variable centroids.
%   dxhatdbeta = derivaties of the fuzzy outputs with respect 
%                to the fuzzy output variable half widths.

NumRules = nclassin(1) * nclassin(2);

w = zeros(NumRules,1);
coutAll = zeros(NumRules,1);
boutAll = zeros(NumRules,2);
areaAll = zeros(NumRules,1);

p = 1;
y = zeros(nout);
m = zeros(nout, max(nclassout));
% The area array is an (nout x max(nclassout)) array that contains the 
% areas of each fuzzy output membership function.
area = (bout(:, :, 1) + bout(:, :, 2)) / 2;
for h = 1 : nout
   for g = 1 : nin
      % grade = fuzzy input membership grades; an [nin x max(nclassin)] dimensional array
      for i = 1 : nclassin(g)
         grade(g, i) = GradeCalc(x(g), b(g, i, 1), b(g, i, 2), c(g, i));
      end
      % Check for inputs outside the membership function ranges.
      if x(g) < c(g, 1) 
         grade(g, 1) = 1;
      elseif x(g) > c(g, nclassin(g)) 
         grade(g, nclassin(g)) = 1;
      end
   end
   % Compute the fuzzy output functions.
   % Note that this logic can handle only 1, 2, or 3 inputs right now.
   if nin == 1
      for i = 1 : nclassin(1)
         if rule(i, h) ~= 0
            m(h, rule(i, h)) = m(h, rule(i, h)) + grade(1, i);
         end
      end
   elseif nin == 2
      for k = 1 : nclassin(2)
         for i = 1 : nclassin(1)
            if rule(i, k, h) ~= 0
               % Compute w(p), the activation level of the p-th rule.
               w(p) = min(grade(1, i), grade(2, k));
               m(h, rule(i, k, h)) = m(h, rule(i, k, h)) + w(p);
               % Set the modal point and half-widths of the p-th rule.
               coutAll(p) = cout(h, rule(i, k, h));
               boutAll(p, 1) = bout(h, rule(i, k, h), 1);
               boutAll(p, 2) = bout(h, rule(i, k, h), 2);
               areaAll(p) = (boutAll(p, 1) + boutAll(p, 2)) / 2;
               p = p + 1;
            end
         end
      end
   elseif nin == 3
      for l = 1 : nclassin(3)
         for k = 1 : nclassin(2)
            for i = 1 : nclassin(1)
               if rule(i, k, l, h) ~= 0
                  m(h, rule(i, k, l, h)) = m(h, rule(i, k, l, h)) + ...
                     min(grade(1, i), grade(2, k), grade(3, l));
               end
            end
         end
      end
   end
   % Defuzzify using the center-of-area method.
   num = 0;
   den(h) = 0;
   for j = 1 : NumRules
      centroid(j) = boutAll(j, 2) * (3 * coutAll(j) + boutAll(j, 2));
      centroid(j) = centroid(j) + boutAll(j, 1) * (3 * coutAll(j) - boutAll(j, 1));
      centroid(j) = centroid(j) / 3 / (boutAll(j, 1) + boutAll(j, 2));
      num = num + w(j) * centroid(j) * areaAll(j);
      den = den + w(j) * areaAll(j);
   end
   if den(h) > 0
      y(h) = num / den(h);
   elseif num > 0
      y(h) = cout(h, nclassout(h));
   elseif num < 0
      y(h) = cout(h, 1);
   else
      y(h) = 0;
   end
end

dxhatdc = 0;
dxhatdb = 0;
dxhatdgamma = 0;
dxhatdbeta = 0;
if ~PartialFlag
   return;
end

% Compute the partial derivatives.
% This info will be used to update the partial derivatives of the RMS error 
% with respect to
% the centroids and widths of the input and output fuzzy sets.  This will
% then be used to optimize the centroids and widths.
% Restricted to a fuzzy system with 2 inputs and 1 output.

% Compute the partials of the fuzzy output w/r to the fuzzy output modal points.
dxhatdgamma = zeros(nout, max(nclassout));
for j = 1 : nclassout(1)
   for p = 1 : NumRules
      dxhatdgammaAll(p) = w(p) * areaAll(p) / den;
   end
end

% Compute the partials of the fuzzy output w/r to the fuzzy output half-widths.
dxhatdbeta = zeros(nout, max(nclassout), 2);
dxhatdbetaAll = zeros(NumRules, 2);
for j = 1 : NumRules
   temp1 = (centroid(j) - y) / 2;
   temp2 = areaAll(j) / 3 / (boutAll(j,1) + boutAll(j,2));
   dxhatdbetaAll(j, 1) = w(j) * (temp1 - temp2) / den;
   dxhatdbetaAll(j, 2) = w(j) * (temp1 + temp2) / den;
end

% Compute the partials of the fuzzy output w/r to the fuzzy input centroids.
dxhatdc = zeros(nin, max(nclassin));
dxhatdb = zeros(nin, max(nclassin), 2);

% Compute the r array. 
% r(i,1,k) = 1 if {input 1 is in fuzzy class i} is
% a premise of the k-th rule, and w(k) = grade(1, i).
% r(j,2,k) = 1 if (input 2 is in fuzzy class j) is
% a premise of the k-th rule, and w(k) = grade(2, j).
r = zeros(max(nclassin), 2, NumRules);
k = 1;
for j = 1 : nclassin(2)
   for i = 1 : nclassin(1)
      if abs(w(k) - grade(1,i)) <= eps
         r(i, 1, k) = 1;
      end
      if abs(w(k) - grade(2,j)) <= eps
         r(j, 2, k) = 1;
      end      
      k = k + 1;
   end
end

for j = 1 : nin
   for i = 1 : nclassin(j)
      for p = 1 : NumRules
         if r(i, j, p) == 1
            temp = areaAll(p) * (centroid(p) - y);
            if ((c(j, i) - b(j, i, 1)) <= x(j)) & (x(j) <= c(j, i))
               dfdc = -1 / b(j, i, 1);
            elseif (c(j, i) <= x(j)) & (x(j) <= (c(j, i) + b(j, i, 2)))
               dfdc = 1 / b(j, i, 2);
            else
               dfdc = 0;
            end
            dxhatdc(j, i) = dxhatdc(j, i) + temp * dfdc;
         end
      end
   end
end
dxhatdc = dxhatdc / den;

% Compute the partials of the fuzzy output w/r to the fuzzy input half-widths.
for j = 1 : nin
   for i = 1 : nclassin(j)
      for p = 1 : NumRules
         if r(i, j, p) == 1
            temp = areaAll(p) * (centroid(p) - y);
            if ((c(j, i) - b(j, i, 1)) <= x(j)) & (x(j) <= c(j, i))
               dfdb = (1 - grade(j, i)) / b(j, i, 1);
               dxhatdb(j, i, 1) = dxhatdb(j, i, 1) + temp * dfdb;
            elseif (c(j, i) <= x(j)) & (x(j) <= (c(j, i) + b(j, i, 2)))
               dfdb = (1 - grade(j, i)) / b(j, i, 2);
               dxhatdb(j, i, 2) = dxhatdb(j, i, 2) + temp * dfdb;
            end
         end
      end
   end
end
dxhatdb = dxhatdb / den;

% Now translate the derivatives from the number of rules to the number
% of output sets.
for i = 1 : nout
   p = 1;
   num = zeros(nclassin(1), nclassin(2));
   for k = 1 : nclassin(2)
      for l = 1 : nclassin(1)
         j = rule(l, k, 1); % j = output class of p-th rule
         dxhatdgamma(i, j) = dxhatdgamma(i, j) + dxhatdgammaAll(p);
         dxhatdbeta(i, j, 1) = dxhatdbeta(i, j, 1) + dxhatdbetaAll(p, 1);
         dxhatdbeta(i, j, 2) = dxhatdbeta(i, j, 2) + dxhatdbetaAll(p, 2);
         num(i,j) = num(i,j) + 1;
         p = p + 1;
      end
   end
end
dxhatdgamma(i, j) = dxhatdgamma(i, j) / num(i,j);
dxhatdbeta(i, j, 1) = dxhatdbeta(i, j, 1) / num(i,j);
dxhatdbeta(i, j, 2) = dxhatdbeta(i, j, 2) / num(i,j);