MATLAB实现NARX-ANFIS时间序列预测（完整程序和数据）

%% % By Piero Alessandro Riva Riquelme. % Undergraduate student. % Universidad de Concepci�n, Chile. % pieroriva@udec.cl. % November of 2020 %% Step 0: Cleaning clear all close all clc fprintf('ANFIS simulation started...'); tic %% Step 1: Process definition fprintf('\nStep 1: Process definition...'); % a) Time definition fprintf('\n\ta) Time definition'); ti = 0.01; step = 0.01; tf = 12; tt = ti:step:tf; nData = length(tt); % b) Process definition (DC motor) fprintf('\n\tb) Process definition...'); a1 = -0.680758164075105; a2 = -0.010439248103841; b0 = 9.217786558421309; theta = [a1 a2 b0]; % c) Generating input fprintf('\n\tc) Generating input...'); ut = zeros(nData,1); h = @(t) heaviside(t); r = @(t) heaviside(t).*(t); uSelect = 'NoisyRamp'; switch uSelect case 'NoisyRamp' u = @(t) -r(t)*4+r(t-3)*4+r(t-3)*4-r(t-9)*4-r(t-9)*4+r(t-12)*4; ut = u(tt); for i=1:nData ut(i) = ut(i)+(rand-0.5)*1.2; end end % d) Generating output fprintf('\n\td) Generating output...'); yt = zeros(nData,1); for i=3:nData yt(i) = -theta(1).*yt(i-1)-theta(2).*yt(i-2)+theta(3).*ut(i-1); end for i=2:nData yt(i) = yt(i)+0.2*(randn-0.5)*0; end % Transform umin = -12; % V umax = 12; % V ymin = -360; % rad/s ymax = 360; % rad/s MIN = 10; % % ? 4 bytes MAX = 90; % % ? 4 bytes for i=1:nData ut(i) = (ut(i)-umin)/(2*umax); yt(i) = (yt(i)-ymin)/(2*ymax); ut(i) = ut(i)*(MAX-MIN)+MIN; yt(i) = yt(i)*(MAX-MIN)+MIN; end % Deltas dut = zeros(nData,1); dyt = zeros(nData,1); for i=2:nData dut(i) = (ut(i)-ut(i-1))+50.0; dyt(i) = (yt(i)-yt(i-1))+50.0; end % f) Plotting fprintf('\n\tf) Plotting...'); figure(1); subplot(1,1,1); cla hold on; plot(tt, ut, '-', 'LineWidth', 1); plot(tt, yt, '-', 'LineWidth', 1); title('Proccess output given input'); xlabel('Time (s)'); ylabel('(%)') grid minor; legend('ut','yt'); figure(2); subplot(1,1,1); cla hold on; plot(tt, dut, '-', 'LineWidth', 1); plot(tt, dyt, '-', 'LineWidth', 1); title('Proccess output given input'); xlabel('Time (s)'); ylabel('(%)') grid minor; legend('dut','dyt'); %% Step 2: Training configuration and space generation fprintf('\nStep 2: Training configuration and space generation...'); % a) Training configuration fprintf('\n\ta) Training configuration...'); xx = 0:0.01:100; % For plotting MFs nTrainingData = nData; maxEpochs = 1000; % The more the merrier nInputs = 2; % Number of inputs, needs to be configured an ANFIS nFuzzy = 5; % Number of MFs per input nOutputs = 1; % Not changeable nRules = nFuzzy^nInputs; % Number of rules K = 0.02; % Initial K maxK = 0.25; % Maximum K, doesn't pass this number growthRate = 0.1; % Growth of K B = 15; % Backward values aIno = 20; % Initial a parameter of premise MFs aOuto = 0.01; % Initial a parameter of consequent MFs tol = 1e-6; % Tolerance for divisions by 0 % b) I/O OUTPUT = yt; % Target function INPUT1 = yt; % First input INPUT2 = ut; % Second input INPUT3 = yt; % Third input INPUT4 = yt; % Fourth input INPUT5 = ut; % Fifth input INPUT6 = ut; % Sixth input % c) Tools fprintf('\n\tb) Generating tools...'); gauss = @(x,a,c) exp(-(x-c).^2./(a.^2)); % Gauss MF for premise part sigmoid = @(x,a,c) 1./(1+exp(-(x-c)./a)); % Sigmoid MF for consequent par invSigmoid = @(x,a,c) c-a.*log(1./x-1); % Inverse of sigmoid function dGauss_da = @(x,a,c) (2.*exp(-(-c+x).*(-c+x)./(a.*a)).*(-c+x).*(-c+x))./(a.*a.*a); % Partial w/r to a dGauss_dc = @(x,a,c) (2.*exp(-(-c+x).*(-c+x)./(a.*a)).*(-c+x))./(a.*a); % Partial w/r to c dinvSigmoid_da = @(x,a,c) -log(1./x-1); % Partial of invSigmoid w/r to a dinvSigmoid_dc = @(x,a,c) 1; % Partial of invSigmoid w/r to c % d) Generating workspace fprintf('\n\tc) Generating workspace...'); % d.1) Initial fuzzy parameters aIne = zeros(nFuzzy,nInputs); cIne = zeros(nFuzzy,nInputs); aOute = zeros(nRules,nOutputs); cOute = zeros(nRules,nOutputs); aInfinal = aIne; cInfinal = cIne; aOutfinal = aOute; cOutfinal = cOute; for i=1:nInputs aIne(:,i) = aIno.*ones(nFuzzy,1); % Initial a for premise MFs cIne(:,i) = (0:(100/(nFuzzy-1)):100); % Initial c for premise MFs end for i=1:nOutputs aOute(:,i) = aOuto*ones(nRules,1); % Initial a for consequent MFs cOute(:,i) = (0:(100/(nRules-1)):100); % Initial c for consequent MFs end % d.2) Training workspace APE = zeros(maxEpochs,1); APEmin = 100000; epochFlag = 1; etaIna = zeros(nFuzzy,nInputs); etaInc = zeros(nFuzzy,nInputs); etaOuta = zeros(nRules,1); etaOutc = zeros(nRules,1); X = zeros(nInputs,1); O5 = zeros(nTrainingData,1); En = zeros(nTrainingData,1); muIn = zeros(nFuzzy,nInputs); w = zeros(nRules,1); wn = zeros(nRules,1); fi = zeros(nRules,1); fi_wn = zeros(nRules,1); dJn_dO5 = zeros(nTrainingData,1); dJn_dO2 = zeros(nRules,1); dO5_dfi = zeros(nRules,1); dO5_dO2 = zeros(1,nRules); dO2_dO1 = zeros(nRules,nFuzzy*nInputs); dfi_da = zeros(nRules,1); dfi_dc = zeros(nRules,1); dmu_daIn = zeros(nFuzzy,nInputs); dmu_dcIn = zeros(nFuzzy,nInputs); dJn_daOut = zeros(nRules,1); dJn_dcOut = zeros(nRules,1); dJp_daOut = zeros(nRules,1); dJp_dcOut = zeros(nRules,1); dJn_dmu = zeros(nFuzzy,nInputs); dJn_daIn = zeros(nFuzzy,nInputs); dJn_dcIn = zeros(nFuzzy,nInputs); dJp_daIn = zeros(nFuzzy,nInputs); dJp_dcIn = zeros(nFuzzy,nInputs); % d.3) Index table indexTable = zeros(nRules,nInputs); for k = 1:nInputs l=1; for j=1:nFuzzy^(k-1):nRules for i =1:nFuzzy^(k-1) indexTable(j+i-1,nInputs-k+1)=l; end l=l+1; if l>nFuzzy l=1; end end end % e) Initial fuzzy sets plotting figure(3) for i=1:nInputs subplot(nOutputs+1,nInputs,i); cla hold on for j=1:nFuzzy plot(xx,gauss(xx,aIne(j,i),cIne(j,i)),'DisplayName',[num2str(char(64+i)) num2str(j)],'LineWidth',0.5); end %legend; title(['Initial premise fuzzy set ' num2str(char(64+i))]); ylabel(['\mu(x_',num2str(i),')']); xlabel(['x_',num2str(i)]); grid minor; end for i=1:nOutputs subplot(nOutputs+1,nInputs,[(i*nInputs+1) (i*nInputs+nInputs)]); cla hold on; for j=1:nRules plot(xx,sigmoid(xx,aOute(j,i),cOute(j,i)),'LineWidth',0.1); end %legend; title(['Initial consequent fuzzy set Z_' num2str(i)]); ylabel(['\mu(z_',num2str(i),')']); xlabel(['z_',num2str(i)]); grid minor end %% Step 3: ANFIS training fprintf('\nStep 3: ANFIS training begin...'); for g=1:maxEpochs dJp_daOut = zeros(nRules,1); dJp_dcOut = zeros(nRules,1); dJp_daIn = zeros(nFuzzy,nInputs); dJp_dcIn = zeros(nFuzzy,nInputs); for i=1+B:nTrainingData %% a) ANFIS; % Prelayer: Selecting input variables X(1) = INPUT1(i-B); X(2) = INPUT2(i); % Layer 1: Input fuzzyfication for k=1:nInputs for j=1:nFuzzy