RBF神经网络

Halcom

RBF网络视频：链接：http://pan.baidu.com/s/1bpzy7Zx 密码：iios
RBF视频资料：链接：http://pan.baidu.com/s/1bOSOu6 密码：d8lb

RBF网络
         径向基函数(RBF-Radial Basis Function)神经网络是由J. Moody和C. Darken在80年代末提出的一种神经网络，它是具有单隐层的三层前馈网络。由于它模拟了人脑中局部调整、相互覆盖接收域（或称感受野-Receptive Field）的神经网络结构，因此，RBF网络是一种局部逼近网络，已证明它能任意精度逼近任意连续函数。
         RBF网络的学习过程与BP网络的学习过程类似，两者的主要区别在于各使用不同的作用函数。BP网络中隐层使用的是Sigmoid函数，其值在输入空间中无限大的范围内为非零值，因而是一种全局逼近的神经网络；而RBF网络中的作用函数是高斯基函数，其值在输入空间中有限范围内为非零值，因为RBF网络是局部逼近的神经网络。
        理论上，3层以上的BP网络能够逼近任何一个非线性函数，但由于BP网络是全局逼近网络，每一次样本学习都要重新调整网络的所有犬只，收敛速度慢，易于陷入局部最小，很难满足控制系统的高度实时性要求。RBF网络是一种3层前向网络，由输入到输出的映射是非线性的，而隐层空间到输出空间的映射是线性的，而且RBF网络局部逼近的神经网络，因而采用RBF网络大大加快学习速度并避免局部极小问题，适合于实时控制的要求。采用RBF网络构成神经网络控制方案，可有效提高系统的精度、鲁棒性和自适应性。
RBF网络特点
     (1) RBF网络的作用函数为高斯函数，是局部的，BP网络的作用函数为S函数，是全局的；
     (2) 如何确定RBF网络隐层节点的中心及基宽度参数是一个困难的问题；
     (3) 已证明RBF网络具有唯一最佳逼近的特性，且无局部极小
RBF网络结构
      RBF网络是一种三层前向网络，由于输入到输出的映射是非线性的，而隐含层空间到输出空间的映射是线性的，从而可以大大加快学习速度并避免局部极小问题。多输入单输出的RBF网络结构如图1所示。

具体原理参考

[1]     刘金琨.智能控制[M],北京：电子工业出版社,2008.

[2]     刘金琨.先进PID控制MATLAB仿真(第3版)[M],北京：电子工业出版社,2011.

[3]    余胜威. MATLAB优化算法案例分析与应用[M].清华大学出版社,2014.

[4]  rbf底层代码：http://cn.mathworks.com/matlabcentral/fileexchange/22173-radial-basis-function-network

Halcom

1. RBF网络
   
2. clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
   
3. warning off         % 消除警告
   
4. feature jit off     % 加速代码运行
   
5. format long         % 数据类型
   
6. tic                 % 运算计时
   
7. global wd sd fs md sj nh inputnum hiddennum outputnum net
   
8. load('librarydata.mat');         % 加载参数
   
9. wd = data(:,1);    % 温度
   
10. sd = data(:,2);    % 湿度
    
11. fs = data(:,3);    % 辐射
    
12. md = data(:,4);    % 密度
    
13. sj = data(:,5);    % 时间
    
14. nh = data(:,6);    % 能耗
    
15. % 训练数据和预测数据
    
16. n=size(data);         % 矩阵维数
    
17. input=data(:,1:5);    % 作为输入
    
18. output=[data(:,6)];   % 作为目标向量输出
    
19. input_train=input(1:n(1,1)-10,:)';      % 训练数据
    
20. input_test=input(n(1,1)-10+1:end,:)';   % 测试数据，6组测试数据
    
21. output_train=output(1:n(1,1)-10,:)';    % 训练数据输出量
    
22. output_test=output(n(1,1)-10+1:end,:)'; % 测试数据输出量
    
23. % 选连样本输入输出数据归一化
    
24. [inputn,inputps]=mapminmax(input_train);
    
25. [outputn,outputps]=mapminmax(output_train);
    
26. p = inputn;     % 网络输入
    
27. t = outputn;    % 网络输出
    
28. % 节点个数
    
29. inputnum = 5;             % 输入节点
    
30. hiddennum = size(p,2);    % 隐藏节点
    
31. outputnum = 1;            % 输出节点
    
32. %% RBF
    
33. % Dimensions
    
34. [R,Q] = size(p);
    
35. [S,Q] = size(t);
    
36. % Architecture
    
37. net = network(1,2,[1;1],[1;0],[0 0;1 0],[0 1]);
    
38. % Simulation
    
39. net.inputs{1}.size = R;
    
40. net.layers{1}.size = Q;
    
41. net.inputWeights{1,1}.weightFcn = 'dist';
    
42. net.layers{1}.netInputFcn = 'netprod';
    
43. net.layers{1}.transferFcn = 'radbas';
    
44. net.layers{2}.size = S;
    
45. net.outputs{2}.exampleOutput = t;
    
46. 
    
47. [r,q] = size(p);
    
48. [s2,q] = size(t);
    
49. 
    
50. % 权值和阈值
    
51. w1 = p';       % 权值
    
52. b1 = ones(q,1)*sqrt(-log(.5));     % 阈值
    
53. a1 = radbas(dist(w1,p).*(b1*ones(1,q)));
    
54. x = t/[a1; ones(1,q)];
    
55. w2 = x(:,1:q);
    
56. b2 = x(:,q+1);
    
57. 
    
58. net.b{1} = b1;
    
59. net.iw{1,1} = w1;
    
60. net.b{2} = b2;
    
61. net.lw{2,1} = w2;
    
62. 
    
63. %% 结果预测
    
64. 
    
65. % 优化的网络进行预测 -- 训练数据
    
66. ty=sim(net,p);            % 预测值
    
67. err = t-ty;               % 误差
    
68. figure(2)   % 误差曲线
    
69. plot(err,'bo-','linewidth',2);title('训练数据误差曲线');
    
70. figure(3)   % 训练数据 实际值与预测值
    
71. plot(output_train,'bo-','linewidth',2);
    
72. train_simu=mapminmax('reverse',ty,outputps);  % 训练数据预测值
    
73. hold on;plot(train_simu,'ro-','linewidth',2);
    
74. legend('训练数据实际值','训练数据预测值');
    
75. train_RMSE = mse(abs(err))    % RMSE 均方根误差
    
76. train_R2 = corr(t',ty').^2    % R2   样本决定系数

复制代码

Halcom

1. % 改进的粒子群算法
   
2. clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
   
3. warning off         % 消除警告
   
4. feature jit off     % 加速代码运行
   
5. format long         % 数据类型
   
6. tic                 % 运算计时
   
7. global ani bni hni dni eni fni gnj si Dij Wnmax Mnmax Fnmax e b r alpha
   
8. parameters;         % 加载参数
   
9. %% PSO算法参数初始化
   
10. % PSO粒子群算法中的两个参数
    
11. c1 = 2;                 % 学习因子1
    
12. c2 = 2;                 % 学习因子2
    
13. w0 = 0.5;               % 惯性权重
    
14. iter = 30;              % 进化次数  
    
15. sizepop = 1000;          % 精英种群
    
16. Vmax = 20;              % 粒子速度上限
    
17. Vmin = -20;             % 粒子速度下限
    
18. %              P                  O              H               v             Y
    
19. minbound = [zeros(1,4),       zeros(1,4),     zeros(1,4),    20*ones(1,6),   zeros(1,4)]; % 粒子取值范围下限
    
20. maxbound = [5000*ones(1,4), 2000*ones(1,4), 3000*ones(1,4), 120*ones(1,6),  ones(1,4)];  % 粒子取值范围上限
    
21. %% 产生初始粒子和速度
    
22. for i=1:sizepop
    
23.     % 随机产生一个种群
    
24.     for j=1:length(minbound)
    
25.         pop(i,j) = round( minbound(j) + (maxbound(j)-minbound(j))*rand );
    
26.     end
    
27.     [fitness(i),pop(i,:)]= PSO_fitness(pop(i,:));     % 适应度值
    
28. end
    
29. V = Vmax*rands(size(pop,1),size(pop,2));     % 初始化速度
    
30. %% 找最好的粒子种群
    
31. [bestfitness,best_index] = min(fitness);  % 最小适应度值
    
32. zbest = pop(best_index,:);   % 全局最佳
    
33. fitnesszbest = bestfitness;  % 全局最佳适应度值
    
34. gbest = pop;                 % 个体最佳
    
35. fitnessgbest = fitness;      % 个体最佳适应度值
    
36. %% 迭代寻优
    
37. wmax = w0;     % 最大权值
    
38. wmin = 0.05;   % 最小权值
    
39. w = wmin:(wmax-wmin)/(sizepop-1):wmax;  % 权重线性递增
    
40. for i=1:iter
    
41.     disp(['当前迭代次数：',num2str(i)]) % 迭代次数
    
42.     for j=1:sizepop
    
43.         % 速度更新
    
44.         V(j,:) = w(j)*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
    
45.         V(j,find(V(j,:)>Vmax))=Vmax;  % 上限
    
46.         V(j,find(V(j,:)<Vmin))=Vmin;  % 下限
    
47.         
    
48.         % 种群更新
    
49.         pop(j,:) = round( (pop(j,:)+0.5*V(j,:)) );    % 取整
    
50.         % pop个体取值范围约束  
    
51.         for np = 1:size(pop,2)
    
52.             if pop(j,np)>maxbound(1,np)     % 上限
    
53.                 pop(j,np)=maxbound(1,np);
    
54.             elseif pop(j,np)<minbound(1,np) % 下限
    
55.                 pop(j,np)=minbound(1,np);
    
56.             end
    
57.         end
    
58.         
    
59.         % 适应度值更新
    
60.         [fitness(j),pop(j,:)] = PSO_fitness(pop(j,:));      
    
61.         % 个体最优更新
    
62.         if fitness(j) < fitnessgbest(j)
    
63.             gbest(j,:) = pop(j,:);
    
64.             fitnessgbest(j) = fitness(j);
    
65.         end
    
66.         % 群体最优更新
    
67.         if fitness(j) < fitnesszbest
    
68.             zbest = pop(j,:);
    
69.             fitnesszbest = fitness(j);
    
70.         end
    
71.     end     
    
72.     M(i)=fitnesszbest;     % 最优适应度值
    
73. end
    
74. % 结果分析
    
75. [Pni,Oni,Hni,vnj,Yni]=funy(zbest)
    
76. 
    
77. figure('color',[1,1,1])
    
78. plot(M,'Linewidth',2)
    
79. title(['适应度曲线  ' '终止代数＝' num2str(iter)]);
    
80. grid on
    
81. xlabel('进化代数');ylabel('适应度值Fit');
    
82. % 结果
    
83. 
    
84. 
    
85. % RBF网络逼近函数程序
    
86. function [fitness]=PSO_RBF(p)
    
87. global cy mc mj ql jb
    
88. % 网络学习参数
    
89. alfa=0.08;  % 动量因子
    
90. xite=0.7;   % 学习因子   
    
91. x=[0,0,0,0]';
    
92. 
    
93. % 节点个数
    
94. inputnum = 4;    % 输入节点
    
95. hiddennum = 3;   % 隐藏节点
    
96. outputnum = 1;   % 输出节点
    
97. 
    
98. % b=1.5*ones(hiddennum,1);   
    
99. % c=0.5*ones(inputnum,hiddennum);   
    
100. % w=rands(hiddennum,outputnum);   
     
101. b=[p(1);p(2);p(3)];   
     
102. c=[p(4),p(5),p(6);
     
103.     p(7),p(8),p(9);
     
104.     p(10),p(11),p(12);
     
105.     p(13),p(14),p(15)];
     
106. w=[p(16);p(17);p(18)];  
     
107. 
     
108. w_1=w;  w_2=w_1;
     
109. c_1=c;  c_2=c_1;
     
110. b_1=b;  b_2=b_1;
     
111. y_1=0;
     
112. 
     
113. for k=1:1:length(jb)
     
114.    
     
115.     x(1)=cy(k);  % 4个输入值
     
116.     x(2)=mc(k);
     
117.     x(3)=mj(k);
     
118.     x(4)=ql(k);
     
119.    
     
120.     y(k)=jb(k);  % 输出值
     
121. 
     
122.     for j=1:1:hiddennum
     
123.         h(j)=exp(-norm(x-c(:,j))^2/(2*b(j)*b(j)));% 高斯基函数
     
124.     end
     
125.     ym(k)=w'*h';
     
126.     em(k)=y(k)-ym(k);
     
127.    
     
128.     dd_w=0*w; d_b=0*b; d_c=0*c;
     
129.     for j=1:1:hiddennum
     
130.         d_w(j,1)=xite*em(k)*h(j);
     
131.         d_b(j,1)=xite*em(k)*w(j)*h(j)*(b(j)^-3)*norm(x-c(:,j))^2;
     
132.         for i=1:1:inputnum
     
133.             d_c(i,j)=xite*em(k)*w(j)*h(j)*(x(i)-c(i,j))*(b(j)^-2);
     
134.         end
     
135.     end
     
136.     w=w_1+ d_w+alfa*(w_1-w_2);
     
137.     b=b_1+d_b+alfa*(b_1-b_2);
     
138.     c=c_1+d_c+alfa*(c_1-c_2);
     
139.    
     
140.     y_1=y(k);
     
141.     w_2=w_1;    w_1=w;
     
142.     c_2=c_1;    c_1=c;
     
143.     b_2=b_1;    b_1=b;
     
144. end
     
145. fitness = sum(sum( abs(em) ));
     

复制代码

Halcom

RBF也可以简单的采用MATLAB工具箱函数：

1. clc;
   
2. clear;
   
3. close all;
   
4. 
   
5. %???????, generate the learing data
   
6. ld=400; %???????
   
7. x=rand(2,ld); %0-1
   
8. x=(x-0.5)*1.5*2; %-1.5, 1.5
   
9. x1=x(1,:);
   
10. x2=x(2,:);
    
11. F=20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2);
    
12. 
    
13. 
    
14. %???????
    
15. net=newrb(x,F);
    
16. 
    
17. %generate the testing data
    
18. interval=0.1;
    
19. [i, j]=meshgrid(-1.5:interval:1.5);
    
20. row=size(i);
    
21. tx1=i(:);
    
22. tx1=tx1';
    
23. tx2=j(:);
    
24. tx2=tx2';
    
25. tx=[tx1;tx2];
    
26. 
    
27. %testing
    
28. ty=sim(net,tx);
    
29. 
    
30. v=reshape(ty,row);
    
31. figure
    
32. subplot(1,3,2)
    
33. mesh(i,j,v);
    
34. zlim([0,60])
    
35. 
    
36. %plot the original function
    
37. interval=0.1;
    
38. [x1, x2]=meshgrid(-1.5:interval:1.5);
    
39. F = 20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2);
    
40. subplot(1,3,1)
    
41. mesh(x1,x2,F);
    
42. zlim([0,60])
    
43. 
    
44. %plot the error
    
45. subplot(1,3,3)
    
46. mesh(x1,x2,F-v);
    
47. zlim([0,60])

复制代码

%Generate some training data
clc;
clear;
interval=0.01;
x1=-1.5:interval:1.5;
x2=-1.5:interval:1.5;
F = 20+x1.^2-10*cos(2*pi*x1)+x2.^2-10*cos(2*pi*x2);
net=newrbe([x1;x2],F)

ty=sim(net,[x1;x2]);
figure
plot3(x1,x2,F,'g');
figure
plot3(x1,x2,ty,'b');

Halcom

Neural Network RBF Regression：

1. clc,clear,close all;
   
2. warning off;
   
3. 
   
4. n1 = 1:2:200;
   
5. x1 = sin(n1*0.1);
   
6. 
   
7. n2 = 2:2:200;
   
8. x2 = sin(n2*0.1);
   
9. 
   
10. xn_train = n1;         
    
11. dn_train = x1;         
    
12. xn_test = n2;         
    
13. dn_test = x2;         
    
14. 
    
15. %---------------------------------------------------
    
16. 
    
17. switch 1
    
18. case 1
    
19.         
    
20. 
    
21. P = xn_train;
    
22. T = dn_train;
    
23. spread = 40;               
    
24. net = newrbe(P,T,spread);
    
25. 
    
26. case 2
    
27.    
    
28. P = xn_train;
    
29. T = dn_train;
    
30. goal = 1e-12;               
    
31. spread = 40;               
    
32. MN = size(xn_train,2);      
    
33. DF = 1;                     
    
34. net = newrb(P,T,goal,spread,MN,DF);
    
35. 
    
36. case 3
    
37.    
    
38. P = xn_train;
    
39. T = dn_train;
    
40. spread = 0.5;               
    
41. net = newgrnn(P,T,spread);
    
42.    
    
43. end
    
44. 
    
45. err1 = sum((dn_train-sim(net,xn_train)).^2)     
    
46. 
    
47. X = sim(net,xn_test);                          
    
48. err2 = sum((dn_test-X).^2)                     
    
49. 
    
50. %---------------------------------------------------
    
51. 
    
52. plot(1:length(n2),x2,'r+:',1:length(n2),X,'bo:')

复制代码