基于时间反向传播 (BPTT)递归神经网络实现非线性系统识别附matlab代码.zip

%% Non-Linear System Identification using Backpropagation Throught Time (BPTT) % Author: SHUJAAT KHAN (shujaat123@gmail.com) clc; clear all; close all; %% Training method BPTT_type='Epochwise'; % BPTT_type='Truncated'; %% system configuration ts=0.001; if (strcmp(BPTT_type,'Epochwise')) a=20; eta = 0.001; %%%%%% Learning Rate 0-1 eg .01, .05, .005 t = 0:ts:a; l=length(t); u = square(pi*t); else a=20; eta = 0.05; %%%%%% Learning Rate 0-1 eg .01, .05, .005 t = 0:ts:a; l=length(t); u = square(pi*t); end % u = awgn(u,30); plant=[1 -0.25 0.1045 0.0902]; % Minimum Phase Plant y(1:2)=0; for i=3:length(t) ph=[y(i-1) y(i-2) cos(2*u(i))+exp(-10*abs(u(i))) u(i-1)]; y(i)=ph*plant'+.001*randn(1,1); end %% RNN architecture inp=6; % Input Layer Neurons (2 nodes for u(i) and u(i-1), 2 for y(i-1) and y(i-2) for state 1 feedback + 1 bias N1=3; % Middle Layer Neurons N2=1; % Output Layer Neurons window_length=3; epoch=(length(t)-2)/window_length; w1=ones(N1,inp); w2=ones(N2,N1); y_NN = 0.*y; j=2; a1=zeros(N1,length(t)); e = ones(window_length,1); %% RNN configuration for k=1:epoch for i=1:window_length j=j+1; Input_Vec(:,i)=[-y_NN(j-1) -y_NN(j-2) u(j) u(j-1) a1(1,j-1) 1]'; n1 = w1*Input_Vec(:,i); a1(:,j)=tansig(n1); %%%% Hidden Layer Activation Function % tansig for [-1,+1] % logsig for [0,1] n2 = w2*a1(:,j); a2=tansig(n2); %%%% Output Layer Activation Function y_NN(j) = a2; e(i) = y(j) - y_NN(j); if(strcmp(BPTT_type,'Epochwise')) %% Epochwise BBTT Y2(:,i) = 2*dtansig(n2,a2)*e(i); %%%% For Output Layer Equation (15.40) for n=n1 Y1(:,i) = diag(dtansig(n1,a1(:,j)),0)*[e + w2'*Y2(:,i)]; %%%% For Hidden Layer elseif(strcmp(BPTT_type,'Truncated')) %% Truncated BPTT Y2(:,i) = 2*dtansig(n2,a2)*e(i); %%%% For Output Layer Y1(:,i) = diag(dtansig(n1,a1(:,j)),0)*w2'*Y2(:,i); %%%% For Hidden Layer end end MSE(k) = mse(e); w1 = w1 + eta*Y1*Input_Vec'; % Equation (15.41) and (15.44) w2 = w2 + eta*Y2(:,i)*a1(:,j)'; % w2 = w2 + eta*Y2*a1(:,j-i+1:j)'; end for n=1:length(MSE) MSE_trend(n)=mean(MSE(1:n)); end 10*log10(MSE_trend(n)) figure semilogy(MSE); hold on semilogy(MSE_trend,'r') title('Objective Function'); figure plot(t,y,'r') hold on plot(t,y_NN,'b') %% Testing a=10; ts=0.001; t = 0:ts:a; l=length(t); u = square(pi*t); u = awgn(u,30); plant=[1 -0.25 0.1045 0.0902]; % Minimum Phase Plant epoch=(length(t)-2)/window_length; y=0; y(1:2)=0; for i=3:length(t) ph=[y(i-1) y(i-2) cos(2*u(i))+exp(-10*abs(u(i))) u(i-1)]; y(i)=ph*plant'; end y_NN = 0.*y; j=2; a1=zeros(N1,length(t)); for k=1:epoch for i=1:window_length j=j+1; Input_Vec(:,i)=[-y_NN(j-1) -y_NN(j-2) u(j) u(j-1) a1(1,j-1) 1]'; n1 = w1*Input_Vec(:,i); a1(:,j)=tansig(n1); %%%% Hidden Layer Activation Function % tansig for [-1,+1] % logsig for [0,1] n2 = w2*a1(:,j); a2=tansig(n2); %%%% Output Layer Activation Function y_NN(j) = a2; e(i) = y(j) - y_NN(j); end MSE(k) = mse(e); end for n=1:length(MSE) MSE_trend(n)=mean(MSE(1:n)); end 10*log10(MSE_trend(n)) figure semilogy(MSE); hold on semilogy(MSE_trend,'r') title('Objective Function'); figure plot(t,y,'m') hold on plot(t,y_NN,'b') legend('Actual Output','RNN');