st1.m

Matlab程序用作混沌时间序列的预测。

% Numenclature
% ---------------
% st1 :singleton type-1 FLS
% mydata: MG time series sampled my data
% fc => ForeCast output

% TO RUN
% ---------
% mydata = mgts(3000);
% n = 0;
% RMSE = ttsf(mydata)


function fc = st1(mydata, np)
global TR;
% ------------------------------------------------
% Preparing 's' matrix of Training Data
% ------------------------------------------------
tgd = mydata(1:504);

k = 0;
for i = 1:500
    for j = 1:5
        s(i, j) = tgd(k+j);
    end
    k = k+1;
end

% ------------------------------------------------
% Preparing 'si' matrix of Test Data
% ------------------------------------------------
tsd = mydata(505:1000);

k = 0;
for i = 1:491
    for j = 1:5
        si(i, j) = tsd(k+j);
    end
    k = k+1;
end

% ------------------------------------------------
% Fuzzy Sets Locations
% ------------------------------------------------

x = 0:0.01:2;
fs = fsloc(x);

% -----------------------------------------------
% Rule Construction Segment
% -----------------------------------------------

N=500;           % No. of Rules
for l=1:N
   for j = 1:5
        for i = 1:201
            if x(i) >= (s(l, j) - 0.005) && x(i) < (s(l, j) + 0.005)
                rn(j) = i;
            end
        end
    end

    RDeg = 1;
    for j = 1:5                 
        for fn = 1:7
            if fs(rn(j), fn) == max(fs(rn(j), :))
                fsn(j) = fn;                                         % Fuzzy set no.
                mgr(j) = max(fs(rn(j), :));             % Membership Grade
            end
        end
        RDeg = RDeg*mgr(j);
    end
    Rule(l,:) = [fsn RDeg];
end

% ------------------------------------------------
% Rule Reduction Segment
% ------------------------------------------------

R=Rule;
for k = 1:N
    for j = 1:N
        if j~=k
            if Rule(k, 1:5)==Rule(j, 1:5) 
                if Rule(k, 6)>Rule(j, 6)
                    Rule(j, :) = zeros(1,6);                
                else
                    Rule(k, :) = zeros(1,6);
                end
            end
        end
    end
end

RRM = [];
for j=1:N
    if Rule(j, :)~=zeros(1, 6)
        RRM = [RRM; Rule(j, :)];
    end
end
[M, m]=size(RRM);

% ---------------------------------------------------------
% Passing Test Data by giving np
% ---------------------------------------------------------
% np=2;
ts = si(np, 1:4);

fsip = fsloc(ts);

fsc = [];
for j = 1:4
    k=0;
    for i = 1:7
        if fsip(j, i) ~= 0
            k=k+1;
            fsc(j, k) = i;
            mgc(j, k) = fsip(j, i);
        end
    end
end

FRM = [fsc(1,1) fsc(2,1) fsc(3,1) fsc(4,1);     % FRM: Fireable Rule Matrix
              fsc(1,1) fsc(2,1) fsc(3,1) fsc(4,2);
              fsc(1,1) fsc(2,1) fsc(3,2) fsc(4,1);
              fsc(1,1) fsc(2,1) fsc(3,2) fsc(4,2); 
              fsc(1,1) fsc(2,2) fsc(3,1) fsc(4,1);
              fsc(1,1) fsc(2,2) fsc(3,1) fsc(4,2);
              fsc(1,1) fsc(2,2) fsc(3,2) fsc(4,1);
              fsc(1,1) fsc(2,2) fsc(3,2) fsc(4,2); 
              fsc(1,2) fsc(2,1) fsc(3,1) fsc(4,1);
              fsc(1,2) fsc(2,1) fsc(3,1) fsc(4,2);
              fsc(1,2) fsc(2,1) fsc(3,2) fsc(4,1);
              fsc(1,2) fsc(2,1) fsc(3,2) fsc(4,2); 
              fsc(1,2) fsc(2,2) fsc(3,1) fsc(4,1);
              fsc(1,2) fsc(2,2) fsc(3,1) fsc(4,2);
              fsc(1,2) fsc(2,2) fsc(3,2) fsc(4,1);
              fsc(1,2) fsc(2,2) fsc(3,2) fsc(4,2)];

MGM = [mgc(1,1) mgc(2,1) mgc(3,1) mgc(4,1);     % MGM: Membership Grade Matrix
              mgc(1,1) mgc(2,1) mgc(3,1) mgc(4,2);
              mgc(1,1) mgc(2,1) mgc(3,2) mgc(4,1);
              mgc(1,1) mgc(2,1) mgc(3,2) mgc(4,2); 
              mgc(1,1) mgc(2,2) mgc(3,1) mgc(4,1);
              mgc(1,1) mgc(2,2) mgc(3,1) mgc(4,2);
              mgc(1,1) mgc(2,2) mgc(3,2) mgc(4,1);
              mgc(1,1) mgc(2,2) mgc(3,2) mgc(4,2); 
              mgc(1,2) mgc(2,1) mgc(3,1) mgc(4,1);
              mgc(1,2) mgc(2,1) mgc(3,1) mgc(4,2);
              mgc(1,2) mgc(2,1) mgc(3,2) mgc(4,1);
              mgc(1,2) mgc(2,1) mgc(3,2) mgc(4,2); 
              mgc(1,2) mgc(2,2) mgc(3,1) mgc(4,1);
              mgc(1,2) mgc(2,2) mgc(3,1) mgc(4,2);
              mgc(1,2) mgc(2,2) mgc(3,2) mgc(4,1);
              mgc(1,2) mgc(2,2) mgc(3,2) mgc(4,2)]; 
              
FR = [];                                % FR: Actually Fired Rule from Rulebase
for i=1:16
    for j=1:M
        if FRM(i, :) == RRM(j, 1:4)
            FR(i, :) = RRM(j, 1:5);
        end
    end
end

RFRM = [];                          % RFRM: Reduced Fire Rule Matrix after removing Zeros
RMGM = [];                          % RFRM: Reduced Mem. Grade Matrix after removing Zeros
[X, x]=size(FR);
for j=1:X
    if FR(j, :) ~= zeros(1, 5)
        RFRM = [RFRM; FR(j, :)];
        RMGM = [RMGM; MGM(j, :)];
    end
end
% RMGM = [RMGM prod(RMGM, 2)];

% -------------------------------------------------
% Height Defuzzification Method
% -------------------------------------------------

[MM mm] = size(RFRM);
for j=1:MM
    switch RFRM(j, 5)
        case 1
            y(j) = 0.2;
        case 2
            y(j) = 0.5;
        case 3
            y(j) = 0.8;
        case 4
            y(j) = 1.0;
        case 5
            y(j) = 1.2;
        case 6
            y(j) = 1.5;
        case 7
            y(j) = 1.8;        
        otherwise
            disp('No Rule');
    end
end

Mu = prod(RMGM, 2);

SNum = 0;
SDen = 0;
for j=1:MM
    SNum = SNum + y(j)*Mu(j);
    SDen = SDen + Mu(j);
end

HDfuzz = SNum/SDen;

fc = HDfuzz;

% -----------------------------------------
% For Better Display
% -----------------------------------------
TR = TR+MM;
AR = TR/np; 
dsp = [np   MM  AR];
disp('~~~~~~~~~~~~~~~~~~~~~~~~')
disp(' ')
disp('  FCast         FRules       ARules')
disp (dsp)

% -------------------------------------------------
% End of st1(mydata, np)
% -------------------------------------------------