seqfis.asv

建立Fuzzy模型－序列预测matlab代码

% 序列预测建立Fuzzy模型
clc;
clear;
% 产生序列,序列为正值,否则ceil函数要变
k       = 0:.1:100;
y       = [0 0];
g       = 0.6.*sin(k.*pi)+0.3.*sin(3.*pi.*k)+0.1.*sin(5.*pi.*k);
for i = 3: 101
    y = [y 0.3*y(i-1)+0.6*y(i-2)+g(i-2)];
end
g       = 2+cos(k);

SLen    = [fix(length(k)/2)] ;



% 窗口长4,表示四个输入变量,输出变量1个
n       = 4;
% 序列的边界
gMin    = ceil(min(g));
gMax    = ceil(max(g));

% 输入的标准模糊集设15个,输出的标准模糊集设7个
InNum   = 15;
InStep  = (gMax-gMin)/(InNum-1);
OutNum  = 10;
OutStep = (gMax-gMin)/(OutNum-1);

% 三角隶属度函数参数表

InMFP   = [gMin,gMin,gMin+InStep];
for i = 1:InNum-2
    InMFP = [InMFP;gMin+(i-1)*InStep,gMin+i*InStep,gMin+(i+1)*InStep];
end
InMFP   = [InMFP;gMax-InStep,gMax,gMax];

OutMFP  = [gMin,gMin,gMin+OutStep];
for i = 1:OutNum-2
    OutMFP = [OutMFP;gMin+(i-1)*OutStep,gMin+i*OutStep,gMin+(i+1)*OutStep];
end
OutMFP  = [OutMFP;gMax-OutStep,gMax,gMax];

% 显示序列
figure(1);
subplot(3,1,1);
plot(k,g);
xlabel('产生的序列(x)');
ylabel('g(x)');
% 显示输入隶属度函数
xx      = linspace(gMin,gMax,1000);
subplot(3,1,2);
hold on;
for i =1:InNum
   yy    = trimf(xx,InMFP(i,:));
   plot(xx,yy);
end
xlabel('输入隶属度函数(xx)');
% 显示输出隶属度函数
subplot(3,1,3);
hold on;
for i =1:OutNum
   yy    = trimf(xx,OutMFP(i,:));
   plot(xx,yy);
end
xlabel('输出隶属度函数(xx)');

% 取50个点,建立规则库
% 规则库变量,最后列为强度
Rule    = [];
for i = 1:SLen-n+1
    % 对一对输入数据求最大隶属度的模糊集
    In0    = zeros(1,n);
    % 初始选择的输入数据对对应第一个模糊集
    InMF   = ones(1,n);

    Out0    = zeros(1);
    % 初始选择的输出数据对对应第一个模糊集
    OutMF   = ones(1);
    for j = 1:InNum
        % 得到隶属度
        In   = [evalmf(g(i),InMFP(j,:),'trimf'),evalmf(g(i+1),InMFP(j,:),'trimf'),...
            evalmf(g(i+2),InMFP(j,:),'trimf'),evalmf(g(i+3),InMFP(j,:),'trimf')];
        % 取最大隶属度的那个模糊集
        for ii = 1:n
            if In(ii)>In0(ii)
                InMF(ii) = j;
                In0(ii)  = In(ii);
            end
        end
    end
     % 对一对输出数据求最大隶属度的模糊集
    for j = 1:OutNum
        % 得到隶属度
        Out    = evalmf(g(i+4),OutMFP(j,:),'trimf');
        % 取最大隶属度的那个模糊集
        if Out>Out0
            OutMF = j;
            Out0  = Out;
        end
    end
    % 强度
    Dens = prod(In0)*Out0;
    % 消除冲突规则 
    if isempty(Rule)
        Rule = [Rule;InMF,OutMF,Dens];
    else
        % 冲突序号
        Coll = 0;
        for jj=1:size(Rule,1)
         % 冲突
          if (Rule(jj,1:n+1)==[InMF,OutMF])
              Coll = jj;
              break;
          end
        end
        if  (Coll>0) 
           if (Rule(Coll,end)<Dens)
               Rule(Coll,end) = Dens;
           end
        else 
           Rule = [Rule;InMF,OutMF,Dens];
        end 
    end 
 end
% 调用预测文件
 Predict;