声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3379|回复: 5

[控制理论] 动态模糊神经网络

[复制链接]
发表于 2012-6-6 20:59 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 牛小贱 于 2015-3-6 11:10 编辑
  1. function [CRBF width rule e RMSE ]=GDFNN(p,t,width0,parameter)
  2. %this is GD-FNN training program
  3. %input P inout data  r*q de matrix
  4. %      T outpue data s*q matrix
  5. %      q is the number of sample data
  6. %parameter(1)=kdmax  (max of accommodation criterion)
  7. %parameter(2)=kdmin  (min of accommodation criterion)
  8. %parameter(3)=gama  (decay constant)
  9. %parameter(4)=emax  (max of output error)
  10. %parameter(5)=emin  (min of output error)
  11. %parameter(6)=beta  (convergence constant)
  12. %parameter(7)=width0  (the width of the frist rule)
  13. %parameter(8)=k  (overlap factor of RBF units)
  14. %parameter(9)=ks  (width updating factor)
  15. %parameter(10)=kerr  (signnificance of a rule)
  16. %output:
  17. %CRBF is the centers of the RBF units which is u*r matrix
  18. %width is the widths of RBF units which is a r by u matrix
  19. %rule is the number of rules for each iteration
  20. %e is the output error for each iteration
  21. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化系统%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  22. if nargin<4 ;error('not enough input aruments');end
  23. if size(p,2)~=size(t,2);error('the input data are not correct');end
  24. [r,q]=size(p);[s2,q]=size(t); %三个输入3*199 个样本,1*199个输出
  25. pp=p.';p1=min(pp);p2=max(pp);rang=[p1' p2'];scope=(rang(:,2)-rang(:,1)); %系统的输入最大值与最小值。
  26. %width0=scope*ones(1,2)/2;
  27. kdmax=parameter(1);kdmin=parameter(2);gama=parameter(3);emax=parameter(4);%参考动态模糊神经网络这本书
  28. emin=parameter(5);beta=parameter(6);k=parameter(7);km=parameter(8);ks=parameter(9);kerr=parameter(10);
  29. ALLIN=[];ALLOUT=[];CRBF=[];width=[];
  30. %when frist sample data coming
  31. ALLIN=p(:,1);ALLOUT=t(:,1);

  32. %--------------seting up the intial FNN---------------------%
  33. diffc=abs(rang-p(:,1)*ones(1,2));
  34. for j=1:r
  35. cvar=abs(diffc(j,:));
  36. [cdmin,nmin]=min(cvar);
  37. if cdmin<=km
  38.     CRBF(j,1)=rang(j,nmin);
  39.     width(j,1)=width0(j,nmin);
  40. else
  41.    [cdb,nb]=max(cvar);
  42.    CRBF(j,1)=p(j,1);
  43.    width(j,1)=cdb/0.8;
  44. end
  45. end
  46. w1=CRBF.';rule(1)=1;
  47. %cacluating the first out error
  48. a0=exp(-mdist(ALLIN,CRBF,width));
  49. a01=[a0 p(:,1).'];
  50. w2=ALLOUT/a01.';
  51. a02=w2*a01.';
  52. sse(1)=sumsqr(ALLOUT-a02)/s2;rmse(1)=sqrt(sse(1));
  53. %%%%%%%%%%%%%%%%%%%%主循环%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  54. for i=2:q
  55. IN=p(:,i);OUT=t(:,i);ALLOUT=[ALLOUT OUT];ALLIN=[ALLIN IN];
  56. [r,N]=size(ALLIN);[r,s]=size(CRBF);
  57. dd=mdist(IN,CRBF,width);
  58. md=sqrt(dd);
  59. [d ind]=min(md);
  60. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  61. %--------------计算实际输出误差----------%
  62. ai=exp(-dd);
  63. ai1=transf(ai,IN);
  64. ai2=w2*ai1;
  65. errout=t(:,i)-ai2;
  66. errout1=errout.*errout;
  67. errout2=sum(errout1)/s2;
  68. e(i)=sqrt(errout2);
  69. %-------------参数选择-----------------%
  70. if i<=q/3
  71. ke=emax;
  72. kd=kdmax;
  73. ks=0.8;
  74. elseif i>q/3 && i<2*q/3
  75.     ke=max(emax*beta.^(i-q/3),emin);
  76.      kd=max(kdmax*gama.^(i-q/3),kdmin);
  77.      ks=0.9;
  78. else
  79.         ke=emin;kd=kdmin;ks=0.95;
  80. end
  81.   %%%%%%%%%%%%%%%%%%%%%%%%%%%马氏距离大于给定值的处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  82. if d>kd  %假如最小马氏距离的开根号,大于预先给定的Kd ,与模糊规则的完备性有关
  83.     if e(i)>ke
  84.         %计算宽度和中心值
  85.         extc=[CRBF rang];   %因半闭模糊集没有讨论在其论域的起点或终点是否存在两个隶属度函数
  86.         extw=[width width0];
  87.         diffc=extc-IN*ones(1,s+2); %计算输入数据域边界集的欧式距离
  88.         for J=1:r  %代表r
  89.             cvar=abs(diffc(J,:)); %计算欧距离
  90.             [cdmin nmin]=min(cvar);
  91.             if cdmin<=km %预定义的参数,该常数控制相邻隶属函数的相似性
  92.                 CRBF(J,s+1)=extc(J,nmin); %因而某个输入变量的隶属函数可能会重复
  93.                 width(J,s+1)=extw(J,nmin);
  94.             else
  95.                 CRBF(J,s+1)=IN(J);
  96.                 width(J,s+1)=k*cdmin; %给定新的宽度值
  97.             end
  98.         end
  99. %%
  100. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  101. %-----------------------完全可以删去----------------------------------%
  102. %但留下来,最初始阶段的近似有较好的效果
  103.     %比较compare the CRBF delect the repeating ones
  104.   oldCRBF=CRBF;oldwidth=width;
  105.   CCRBF=[];newwd=[];
  106.   [r,s1]=size(CRBF);
  107.   while s1>1
  108.       redc=CRBF;
  109.       redc(:,1)=[];
  110.       for j2=1:s1-1
  111.       dcc=CRBF(:,1)-redc(:,j2);
  112.       if all(dcc==0);break;end
  113.       end
  114.      if any(dcc~=0)
  115.      CCRBF=[CCRBF CRBF(:,1)];
  116.      newwd=[newwd width(:,1)];
  117.      end
  118.    CRBF(:,1)=[];width(:,1)=[];
  119.    s1=s1-1;
  120.   end
  121.   CCRBF=[CCRBF CRBF(:,1);];
  122.   newwd=[newwd width(:,1)];
  123.   CRBF=CCRBF;
  124.   width=newwd;
  125. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  126.   %%     
  127.          
  128.         
  129.    w1=CRBF.';[u r]=size(w1);
  130. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算所有规则的ERR%%%%%%%%%%%%%%%%%%%%%%%%%
  131.   A=exp(-mdist(ALLIN,CRBF,width));
  132.   A0=transf(A,ALLIN);
  133.   if u*(r+1)<=N
  134.   %-------------计算误差衰减率--%
  135.   err=cperr(A0,ALLOUT);
  136.   errT=err.';
  137.   err1=zeros(u,s2*(r+1));
  138.   err1(:)=errT;
  139.   err21=err1.'; err22=sum(err21.*err21)/(s2*(r+1));
  140.   err23=sqrt(err22); %err23 代表各条规则的最终误差的减少率
  141.   No=find(err23<kerr);
  142.   %%%%%%%%%%%%%%%%%%%%%%%%%规则是否需要剔除%%%%%%%%%%%%%%%%%%%%%%%%%
  143.   if ~isempty(No)
  144.    %   delc=CRBF(:,No);
  145.    %   delw=width(:,No);
  146.       CRBF(:,No)=[];
  147.       w1(No,:)=[];
  148.       width(:,No)=[];
  149.     %  rCRBF=CRBF;
  150.     %  rwidth=width;
  151.      % err21(:,No)=[];
  152.       [u1 r]=size(w1);
  153. % ------------结果参数调整与总体误差计算----------%
  154.        A=exp(-mdist(ALLIN,CRBF,width));
  155.      A0=transf(A,ALLIN);
  156.      w2=ALLOUT/A0;
  157.      A1=w2*A0;
  158.      derr=ALLOUT-A1;
  159.      sse=sumsqr(derr)/(i*s2);
  160.      rmse=sqrt(sse);RMSE(i)=rmse;
  161.      rule(i)=u1;
  162.   else
  163.   %-------在 最小马氏距离大于给定参数且误差大于也大于给定误差限,后经调整后。误差率值大于给定的值
  164.   %此时该规则不用提出。因而只需参数调整和计算误差。
  165.       w2=ALLOUT/A0; %结果参数调整
  166.      A1=w2*A0;
  167.      derr=ALLOUT-A1;
  168.      sse=sumsqr(derr)/(s2*i);
  169.      rmse=sqrt(sse);RMSE(i)=rmse;
  170.      rule(i)=u;
  171.   end
  172.   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  173.   else
  174.        w2=ALLOUT/A0;
  175.      A1=w2*A0;
  176.      derr=ALLOUT-A1;
  177.      sse=sumsqr(derr)/(s2*i);
  178.      rmse=sqrt(sse);RMSE(i)=rmse;
  179.      rule(i)=u;      
  180.   end  
  181.   %马氏距离大于给定值且误差小于给定值,直接调整结果参数
  182.     else
  183.         [w2 A1 A]=mweight(ALLIN,CRBF,width,ALLOUT);
  184.          derr=ALLOUT-A1;
  185.      sse=sumsqr(derr)/(i*s2);
  186.      rmse=sqrt(sse);RMSE(i)=rmse;
  187.      rule(i)=s;
  188.     end
  189.     %%
  190.     %%%%%%%%%%%%%%%%%%%%%%%%%%%小于最小马氏距离的处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  191.    %------------------或者小于最小马氏距离的开根号------------------------%
  192. else  
  193. %%
  194. %高斯宽度修正
  195.     if e(i)>ke   % %在最小马氏距离小于给定值情况下,如果系统误差大于系统给定的阀值
  196.         if s*(r+1)<=N %在规则数满足该条件下才运行本模块,N代表,此时具有的样本数,是一个不断递增的数
  197.          aa2=exp(-mdist(ALLIN,CRBF,width));
  198.      aa3=transf(aa2,ALLIN);
  199.      werr=cperr(aa3,ALLOUT) %计算误差下降率
  200.      werrT=werr.';
  201.      werr1=zeros(s,s2*(r+1));  %s stanf for the number of rule and se stand for output vatieble
  202.      werr1(:)=werrT;
  203.      werr21=werr1.';   %矩阵的全新排列有1* (r+1)*u==r+1 *u
  204.      err31=reshape(werr21(:,ind),r+1,s2);
  205.      err32=abs(err31.');
  206.      if s2>1
  207.      err33=sum(err32);
  208.      else
  209.          err33=err32;
  210.      end
  211.      err3=sum(err33);
  212.      err30=err33(1);
  213.      rate0=err30/err3;
  214.      err33(1)=[];  %在该条件下,虽然Xk 能被相邻的模糊规则包含,但该规则并没有重要到该RBF单元能
  215.      err3r=err33;  %很好的覆盖Xk 或者反映系统的输出在这个椭球区域有急剧变化。因此该椭球区域需要修正
  216.      rater=err3r./(err3-err30);
  217.      nm=find(rater<=1/r);
  218.      width(nm,ind)=ks.*width(nm,ind); %ind 代表,对一个新的样本,找到以马氏距离最靠近该样本的那条规则
  219.      [w2 A1 A]=mweight(ALLIN,CRBF,width,ALLOUT);
  220.      derr=ALLOUT-A1; %A1,代表神经网络的输出
  221.      sse=sumsqr(derr)/(i*s2); %输出项全部对应的误差的平方和,再除以总共的项
  222.      rmse=sqrt(sse);RMSE(i)=rmse;
  223.      rule(i)=s;
  224.         end
  225.     else %在最小马氏距离小于给定值情况下,如果系统误差小于系统给定的阀值
  226.          [w2 A1 A]=mweight(ALLIN,CRBF,width,ALLOUT);
  227.      derr=ALLOUT-A1;       %什么都不做,直接进去下一个循环。
  228.      sse=sumsqr(derr)/(i*s2); %输出项全部对应的误差的平方和,再除以总共的项
  229.      rmse=sqrt(sse);RMSE(i)=rmse;
  230.      rule(i)=s;
  231.     end
  232. %%
  233. end
  234. end
  235. end


  236. function y=mdist(p ,c ,b)
  237. %计算马氏距离,相当于第三层网络的输入。
  238. %p= r*q 矩阵 是采样输入矩阵
  239. %c=r*s 是中心矩阵
  240. %b =r*s 矩阵是基宽矩阵
  241. [r,q]=size(p); [r,s]=size(c);y=zeros(s,q);
  242. if r==1   
  243. for i=1:s %s 代表规则数,r 代表输入的量
  244. x=c(:,i)*ones(1,q);
  245. d=abs(p-x);
  246. xx=b(:,i)*ones(1,q);
  247. dd=d./xx;
  248. y(i,:)=dd.*dd; %同一条规则中计算每个样本的马氏距离
  249. end
  250. else
  251.     for i=1:s
  252.         x=c(:,i)*ones(1,q);  %将C的第i列复制q倍,即样本数
  253.     d=abs(p-x); %使每个输入都有相应的规则计算距离
  254. xx=b(:,i)*ones(1,q);
  255. dd=d./xx; %同一条规则中计算每个样本的马氏距离
  256. %假设  dd=[  0    0.0022    0.0044    0.0067
  257. %         0    0.0022    0.0044    0.0067]
  258. %代表有两个输入,四个样本
  259. y(i,:)=sum(dd.*dd);
  260.     end   
  261. end
  262. end

  263. function [w y a]=mweight(in,c,wb,out)
  264. %调整TSK FNNS 权值
  265. %[w2 A1 A]=mweight(p ,CRBF, width, t);
  266. %in 为采样数据样本,r个输入,n de sample
  267. %C the centered  with u rule
  268. %out s2 as output number
  269. [r,n]=size(in);[r u]=size(c);
  270. [s2 n]=size(out);
  271. a=exp(-mdist(in,c,wb));%输出结果为马氏距离,后成第三层的输出=规则数*样本数
  272. a1=transf(a,in); %输出结果为
  273. w=out/a1;
  274. y=w*a1;

  275. end

  276. function BB=transf(A ,P)
  277. %计算结果参数
  278. %A output of RBF
  279. %P sample input
  280. [u N]=size(A);[r N]=size(P);
  281. for j=1:N   %样本数
  282.   for i=1:r %输入的变量
  283.   PA((i-1)*u+1:i*u,j)=P(i,j)*A(:,j); %TSK模型,是输入变量的线性组合。各个输入变量需乘以各个规则中的响应的取值。
  284.                                       %每个输入样本对应的规则的数值应该是不一样的。
  285.   end
  286. end
  287. BB=[A;PA];
  288. end

  289. function y=cperr(a,b)
  290.    %  werr=cperr(aa3,ALLOUT);
  291.    %b 代表系统的期望输出误差
  292.    %a 代表TSK模型中,最后计算的动态变化部分,(包含了系统的输入量)
  293. %计算误差衰减率
  294. tT=b';PAT=a';
  295. [WW AW]=orthogonalize(PAT);
  296. WSSW=sum(WW.*WW).';
  297. WSStT=sum(tT.*tT).';
  298. y=((WW.'*tT).'.^2)./(WSStT*WSSW.'); %源于线性回归的计算。。需差线性回归的计算方法
  299. end

  300. function [w a]=orthogonalize(p)
  301. %该变换是每个W的基向量各自期望输出能量的贡献成为可能
  302. %p=W*a ,其中p and w have the same column,a is squared matrix
  303. [u v]=size(p);w(:,1)=p(:,1);a=eye(v);
  304. for k=2:v
  305.     b=zeros(u,1);
  306.     for i=1:k-1
  307.      a(i,k)=w(:,i).'*p(:,k)/(w(:,i).'*w(:,i));
  308.      b=b+a(i,k)*w(:,i);
  309.     end
  310.     w(:,k)=p(:,k)-b;
  311. end
  312. end
复制代码


点评

赞成: 5.0
以后不建议楼主这样向论坛里发帖大量的程序,望有自己的原创内容!  发表于 2015-3-6 11:11
赞成: 5
感谢分享  发表于 2015-2-2 11:31
回复
分享到:

使用道具 举报

发表于 2015-2-2 11:31 | 显示全部楼层
发表于 2016-9-18 11:24 | 显示全部楼层
发表于 2016-9-18 13:04 | 显示全部楼层
这个是干什么的?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 13:14 , Processed in 0.079773 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表