elman identification.m

自己编写的改进ELMAN网络辨识程序.非常实用

%基于改进Elman网络的辨识
clear all;
close all;

xite=0.25;%η
alfa=0.05;%α动量因子
alfa2=0.05;%自反馈因子
belte=0.01;%
x=[0,0,0,0,0,0]';   %  承接层向量
r=[0,0,0]';   %  输入向量
xc=[0,0,0,0,0,0,0,0,0]';%高斯函数的输入

ci=30*ones(9,6);
bi=40*ones(6,1);
w1=10*ones(6,1); % 隐层到输出层权值
w2=10*ones(6,1); % 承接层到隐层的权值



h=[0,0,0,0,0,0]';  %隐层输出
   
ci_1=ci;ci_2=ci_1;ci_3=ci_2;
bi_1=bi;bi_2=bi_1;bi_3=bi_2;
w1_1=w1;w1_2=w1_1;w1_3=w1_2;
w2_1=w2;w2_2=w2_1;w2_3=w2_2;
h_1=h;
x_1=x;
x=alfa2*x_1+h_1;%初始化

u_1=0;y_1=0;%函数中间变量
pc=[0,0,0]';%PID三量状态值
error_1=0;error_2=0;error=0;%误差初始化

%PID参数初始化
kp0=0.03;        
ki0=0.01;     
kd0=0.03;

kp_1=kp0;
kd_1=kd0;
ki_1=ki0;  

xitekp=0.20;
xitekd=0.20;
xiteki=0.20;

ts=0.001;
for k=1:1:2000
   time(k)=k*ts;
   rin(k)=1.0*sign(sin(2*pi*k*ts));
   yout(k)=(-0.1*y_1+u_1)/(1+y_1^2);  %非线性对象
   
   for j=1:1:6
      h(j)=exp(-norm(xc-ci(:,j))^2/(2*bi(j)*bi(j)));
   end
   ymout(k)=w2'*h;         

   d_w1=0*w1;
   for j=1:1:6
       d_w1(j)=xite*(yout(k)-ymout(k))*h(j)*(1-h(j))*h_1(j)*w2(j);
   end
   w1=w1_1+d_w1+alfa*(w1_1-w2_2)+belte*(w1_2-w1_3);
   
   d_w2=0*w2;
   for j=1:1:6
      d_w2(j)=xite*(yout(k)-ymout(k))*h(j);
   end
   w2=w2_1+d_w2+alfa*(w2_1-w2_2)+belte*(w2_2-w2_3);
   
   d_bi=0*bi;
   for j=1:1:6
      d_bi(j)=xite*(yout(k)-ymout(k))*w2(j)*h(j)*(bi(j)^-3)*norm(xc-ci(:,j))^2;
   end
   bi=bi_1+ d_bi+alfa*(bi_1-bi_2)+belte*(bi_2-bi_3);
   for j=1:1:6
     for i=1:1:9
      d_ci(i,j)=xite*(yout(k)-ymout(k))*w2(j)*h(j)*(xc(i)-ci(i,j))*(bi(j)^-2);
     end
   end
   ci=ci_1+d_ci+alfa*(ci_1-ci_2)+belte*(ci_2-ci_3);
  
%%%%%%%%%%%%%%%%%%%%%%Jacobian%%%%%%%%%%%%%%%%%%%%%%%
  yu=0;
  for  j=1:1:6
      yu=yu+w2(j)*h(j)*(-x(1)+ci(1,j))/bi(j)^2; 
  end
  dyout(k)=yu;
%%%%%%%%%%%%%%%%%%%%%%Start of Control system%%%%%%%%%%%%%%%%%%
   error(k)=rin(k)-yout(k);
   kp(k)=kp_1+xitekp*error(k)*dyout(k)*pc(1);
   kd(k)=kd_1+xitekd*error(k)*dyout(k)*pc(2);
   ki(k)=ki_1+xiteki*error(k)*dyout(k)*pc(3);  
   if kp(k)<0
      kp(k)=0;
   end
   if kd(k)<0
      kd(k)=0;
   end
   if ki(k)<0
      ki(k)=0;
   end
   
   M=0;
   switch M
   case 0
   case 1  %Only PID Control
		kp(k)=kp0;
		ki(k)=ki0;     
		kd(k)=kd0;
   end
	du(k)=kp(k)*pc(1)+kd(k)*pc(2)+ki(k)*pc(3); 
   u(k)=u_1+du(k);

%Return of parameters
   r(1)=du(k);
   r(2)=yout(k);
   r(3)=y_1;

   
	u_1=u(k);
   y_1=yout(k);
   
   ci_3=ci_2;
   ci_2=ci_1;
   ci_1=ci;
   
   bi_3=bi_2;
   bi_2=bi_1;
   bi_1=bi;
   
   %   权值更新
   w1_3=w1_2;
   w1_2=w1_1;
   w1_1=w1;
   w2_3=w2_2;
   w2_2=w2_1;
   w2_1=w2;
   
   h_1=h;
   x=alfa2*x_1+h_1;
   xc=[x;r];
   
   
   pc(1)=error(k)-error_1;             %Calculating P
   pc(2)=error(k)-2*error_1+error_2;   %Calculating D
   pc(3)=error(k);                     %Calculating I
   
   error_2=error_1;
   error_1=error(k);
   
   kp_1=kp(k);
   kd_1=kd(k);
   ki_1=ki(k);  
end
figure(1);
plot(time,rin,'b',time,yout,'r');
xlabel('time(s)');ylabel('rin,yout');
figure(2);
plot(time,yout,'r',time,ymout,'b');
xlabel('time(s)');ylabel('yout,ymout');
figure(3);
plot(time,dyout);
xlabel('time(s)');ylabel('Jacobian value');
figure(4);
subplot(311);
plot(time,kp,'r');
xlabel('time(s)');ylabel('kp');
subplot(312);
plot(time,ki,'r');
xlabel('time(s)');ylabel('ki');
subplot(313);
plot(time,kd,'r');
xlabel('time(s)');ylabel('kd');