基于粒子群算法的神经网络非线性函数拟合

基于粒子群算法的神经网络非线性函数拟合

文章初心

最近在学机器学习，自己的方向是智能算法，课程报告需要，于是试着把机器学习和粒子群算法相结合，写出来供大家参考，交流。

文末有这部分内容相关的代码，已开源。

摘要

BP神经网络的结构简单，可操作性强，但也存在收敛速度慢，易陷入局部最优值的缺点。随着科学技术的飞速发展，智能算法也逐步崭露头角，它是智能领域中的重要组成部分，现在已被广泛应用到了数据处理、模式识别、机器学习等多个领域。粒子群算法模拟鸟群觅食行为，以其简单、高效特点被广泛运用，但其也存在易陷入局部最优解、收敛精度低等缺点。本文基于基本粒子群算法的特点，提出改进的粒子群算法——模拟退火粒子群算法（Simulated Annealing Particle Swarm Optimization Algorithm, SAPSO）和混沌模拟退火粒子群算法（Chaotic Simulated Annealing Particle Swarm Optimization Algorithm, CSAPSO），并将其运用于BP神经网络非线性函数拟合中，同时分析了超参数对神经网络模型的影响。仿真实验验证了粒子群算法在BP神经网络运用的有效性。

关键词

粒子群算法；混沌；模拟退火；BP神经网络；函数拟合

BP神经网络

BP（Back-Propagation）神经网络是1986年由Rumelhart和McClelland为首的科学家提出的概念，是一种按照误差逆向传播算法训练的多层前馈神经网络，是应用最广泛的神经网络。BP神经网络是多层前馈型神经网络中的一种，属于人工神经网络的一类，可以对任何一种非线性输入输出关系进行模仿，因此被广泛应用在分类识别（classification）、回归（regression）、压缩（compression）、逼近（fitting）等领域。在工程应用中，大约百分之八十的神经网络模型都选择采用BP神经网络或者改进的BP神经网络。BP神经网络的结构简单，可操作性强，可以很好的处理在模拟线性与非线性的输入输出关系方面的问题。

1.1. BP神经网络结构

BP网络由输入层、隐层和输出层组成，隐层可以有一层或多层。BP网络的优势就是能学习和储存大量的输入输出的关系，而不用事先指出这种数学关系，同时它能以任意精度逼近任一非线性函数。

图中， $x$ 为输入数据， $y$ 为输出数据， $i$为输入层， $j$为隐藏层， $k$为输出层； $W_{ji}$为隐含层权值， $W_{kj}$ 为输出层权值， $f(x)$为神经元激励函数，BP算法常用 $sigmoid$ 函数作为其激励函数。

1.2. BP神经网络训练算法

（1）初始化：置所有权值为较小的随机数；
（2）提供训练集；
（3）计算实际输出：计算隐含层、输出层各神经元输出；
（4）计算目标值与实际输出的偏差MSE；
（5）计算 $\Delta W_{kj}$,$\Delta W_{ji}$；
（6）返回（2）重复计算，直到误差满足要求为止。

基本思想：利用输出后的误差来估计输出层前一层的误差，再用这层误差来估计更前一层误差，如此获取所有各层误差估计。这里的误差估计可以理解为某种偏导数，根据这种偏导数来调整各层的连接权值，再用调整后的连接权值重新计算输出误差。直到输出的误差达到符合的要求或者迭代次数。BP的传播对象就是“误差”，传播目的就是得到所有层的估计误差最小时系统的各参数值。

1.3. 本文的改进方案概述

针对BP神经网络特点，本文利用全局搜索能力较强的PSO算法优化该网络的网络参数。通过反向传播不断调整网络的权值和阈值，最后使全局误差系数最小。它的学习本质就是：对各连接权值的动态调整。同时针对PSO算法对离散的优化问题处理不佳并容易陷入局部最优的缺点，结合PSO算法的相应改进算法，对BP神经网络进行改进。对比分析不同算法对BP神经网络的改进效果，验证算法的有效性和性能。以下为BP的算法流程图：

2. 基本粒子群算法及其改进算法

2.1 基本粒子群算法（Particle Swarm Optimization, PSO）

2.1.1 PSO算法原理

粒子群优化算法是由Kennedy和Eberhart通过模拟鸟群的社会行为提出的一种进化计算方法[1]，其思想来源于人工生命和演化计算理论。粒子群算法的基本思想是：设定优化问题的每个潜在解为粒子，初始化阶段中随机选取一个粒子种群,每个粒子的性能由一个目标函数来衡量，通过此目标函数的适应值来衡量粒子所代表解的优劣，不断迭代寻优，从而得到更优的粒子种群，这个过程一直循环往复，直到达到满意的收敛条件[2]。

其中， $p_{pd}$,$p_{gd}$分别表示粒子个体最优位置和群体最优位置。总得来说，公式（1）和（2）决定了粒子速度、位置的更新方式，体现在迭代中，即表现为粒子不断逼近最优解的过程。由于粒子群算法操作简单，收敛速度快的优点，被广泛应用于复杂问题的优化。

2.1.2 PSO算法的缺陷

PSO收敛快,特别是在算法的早期,但也存在着精度较低,易发散等缺点。若加速系数、最大速度等参数太大,粒子群可能错过最优解,算法不收敛;而在收敛的情况下,由于所有的粒子都向最优解的方向飞去,所以粒子趋向同一化(失去了多样性),使得后期收敛速度明显变慢,同时算法收敛到一定精度时,无法继续优化,因此很多学者都致力于提高PSO算法的性能。

这里给出基本粒子群算法的MATLAB主程序（改进的PSO由于涉及部分实验室项目，因此不作公布，望理解），代码程序包见文末附件。

% 清空环境变量
clc
clear
% 生成训练数据，数量100
x1 = linspace(1,100,100);
x2 = linspace(1,100,100);
X = [x1;x2];
Y = zeros( 100,100);
for row = 1 : 1 : 100    
    for col = 1 : 1 : 100        
        Y( row,col) = sin(10*x1(row))-x2(col).^3+(x1(row).^2) .* x2(col);
    end
end

% 生成检验数据，数量100
xt1 = linspace(1,100,100);
xt2 = linspace(1,100,100);
XT = [xt1;xt2];
Y2 = zeros( 100,100);
for row = 1 : 1 : 100
    for col = 1 : 1 : 100
        Y2( row,col) = sin(10*xt1(row))-xt2(col).^3+(xt1(row).^2) .* xt2(col); 
    end
end
% 对样本输入X输出Y作归一化处理，数据范围限制在[-1,1]，归一化数据结构保存在ps
[Data_target,ps_output] = mapminmax(Y,-1,1);
[Data_input,ps_input] = mapminmax(X,-1,1);
% 对检验数据做归一化处理
Data_test = mapminmax('apply',XT,ps_input);

%节点个数
inputnum=size(Data_input,1);       % 输入层神经元个数 
outputnum=size(Data_target,1);     % 输出层神经元个数
hiddennum=10;
% 创建网络；
net1 = newff(Data_input,Data_target,hiddennum);
% net2 = newff(Data_input,Data_target,hiddennum);
% net3 = newff(Data_input,Data_target,hiddennum);
%节点总数 2*5 + 5 + 5 + 1 = 21 
numsum=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;

%% 粒子群算法求权值和阈值
%粒子群算法参数设置
N = 20;
c1 = 2;
c2 = 2;
w = 0.6;
M = 100;
D = numsum;
x = zeros(1,D);
% 调用粒子群算法函数
[xm1,fv1,Pbest1] = NNPSO(x,hiddennum,net1,Data_input,Data_target,N,w,c1,c2,M,D);
% [xm2,fv2,Pbest2] = NNSAPSO(x,hiddennum,net2,Data_input,Data_target,N,w,c1,c2,M,D);
% [xm3,fv3,Pbest3] = NNCSAPSO(x,hiddennum,net3,Data_input,Data_target,N,w,c1,c2,M,D);

%% 把最优初始阀值权值赋予网络预测
% 用粒子群算法优化的BP网络进行值预测
w1_1=xm1(1:inputnum*hiddennum);
B1_1=xm1(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2_1=xm1(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2_1=xm1(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);

net1.iw{1,1}=reshape(w1_1,hiddennum,inputnum);
net1.lw{2,1}=reshape(w2_1,outputnum,hiddennum);
net1.b{1}=reshape(B1_1,hiddennum,1);
net1.b{2}=reshape(B2_1,outputnum,1);

% % 用模拟退火粒子群算法优化的BP网络进行值预测
% w1_2=xm2(1:inputnum*hiddennum);
% B1_2=xm2(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
% w2_2=xm2(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
% B2_2=xm2(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);
% 
% net2.iw{1,1}=reshape(w1_2,hiddennum,inputnum);
% net2.lw{2,1}=reshape(w2_2,outputnum,hiddennum);
% net2.b{1}=reshape(B1_2,hiddennum,1);
% net2.b{2}=reshape(B2_2,outputnum,1);

% 用混沌模拟退火粒子群算法优化的BP网络进行值预测
% w1_3=xm3(1:inputnum*hiddennum);
% B1_3=xm3(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
% w2_3=xm3(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
% B2_3=xm3(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);
% 
% net3.iw{1,1}=reshape(w1_3,hiddennum,inputnum);
% net3.lw{2,1}=reshape(w2_3,outputnum,hiddennum);
% net3.b{1}=reshape(B1_3,hiddennum,1);
% net3.b{2}=reshape(B2_3,outputnum,1);

%% BP网络训练
%粒子群网络进化参数
net1.trainParam.epochs=100;
net1.trainParam.lr=0.1;
net1.trainParam.goal=1e-3;
% 
% %模拟退火粒子群网络进化参数
% net2.trainParam.epochs=100;
% net2.trainParam.lr=0.1;
% net2.trainParam.goal=1e-6;

%混沌模拟退火粒子群网络进化参数
% net3.trainParam.epochs=100;
% net3.trainParam.lr=0.1;
% net3.trainParam.goal=1e-3;

% 训练网络
net1 = train(net1,Data_input,Data_target); % 粒子群
% net2 = train(net2,Data_input,Data_target); % 模拟退火粒子群
% net3 = train(net3,Data_input,Data_target); % 混沌模拟退火粒子群

%% 仿真测试
test_sim1 = sim(net1,Data_test); % 粒子群
% test_sim2 = sim(net2,Data_test); % 模拟退火粒子群
% test_sim3 = sim(net3,Data_test); % 混沌模拟退火粒子群

% 输出数据反归一化，Test_sim为测试数据通过神经网络的预测输出值
Test_sim1 = mapminmax('reverse',test_sim1,ps_output); % 粒子群
% Test_sim2 = mapminmax('reverse',test_sim2,ps_output); % 模拟退火粒子群
% Test_sim3 = mapminmax('reverse',test_sim3,ps_output); % 混沌模拟退火粒子群

%% 算法结果分析 
% figure(1)
% t = 1:M;
% plot(t,Pbest1,'b',t,Pbest2,'g',t,Pbest3,'r');
% title('算法收敛过程');
% xlabel('进化代数');
% ylabel('最小均方误差值（MSE值）');
% legend('基本粒子群算法','模拟退火粒子群算法','混沌模拟退火粒子群算法');
% 
% figure(2)
% mesh(xt1,xt2,Y2);
% title('函数实际图形');

%% 拟合图形对比输出

%基础粒子群对比输出
figure(3)
subplot(1,2,1)
mesh(xt1,xt2,Y2);
title('函数实际图形');
xlabel('X1取值');ylabel('X2取值');zlabel('非线性函数输出值');

subplot(1,2,2)
mesh(xt1,xt2,Test_sim1);
title('基础粒子群算法拟合图形');
xlabel('X1取值');ylabel('X2取值');zlabel('非线性函数输出值');

% % 模拟退火粒子群对比输出
% figure(4)
% subplot(1,2,1)
% mesh(xt1,xt2,Y2);
% title('函数实际图形');
% xlabel('X1取值');ylabel('X2取值');zlabel('非线性函数输出值');
% 
% subplot(1,2,2)
% mesh(xt1,xt2,Test_sim2);
% title('模拟退火粒子群算法拟合图形');
% xlabel('X1取值');ylabel('X2取值');zlabel('非线性函数输出值');

%%混沌模拟退火粒子群拟合图形单独输出
% figure(5)
% subplot(1,2,1)
% mesh(xt1,xt2,Y2);
% title('函数实际图形');
% xlabel('X1取值');ylabel('X2取值');zlabel('非线性函数输出值');
% 
% subplot(1,2,2)
% mesh(xt1,xt2,Test_sim3);
% title('混沌模拟退火粒子群算法拟合图形');
% xlabel('X轴');ylabel('Y轴');zlabel('Z轴');

2.2 模拟退火PSO算法(Simulated Annealing PSO algorithm, SAPSO)原理

模拟退火算法是一种基于蒙特卡洛[3]思想设计的用于近似求解最优化问题的著名方法。该算法的基本思想主要是通过模拟物理退火[4]过程实现搜索最优解。

从物理过程来说，模拟退火过程一般指物理固体降温、内能逐步减小，最终趋于最低能量状态的结晶过程。而这一过程中，温度的变化应该是缓慢的。固体温度被升高后，其内部粒子的状态为无序且非稳定，随着退火过程的持续物体的粒子状态则逐步趋于有序且最终达到稳态的平衡，而这一过程，则类似算法中粒子不断寻求最优解的过程。

模拟退火算法的特点是在搜索过程中具有概率突跳的能力[5]，能够有效地避免搜索过程陷入局部最优解，退火过程中不但接受更优的解，而且还以一定概率 接受差的解，同时这种概率受到温度参数的控制，其大小随着温度的下降而减小。

在与粒子群算法结合的过程中，必须确定初始温度和退温方式，其公式如下：

2.3 混沌模拟退火PSO算法(Chaotic Simulated Annealing PSO algorithm, CSAPSO)原理

将混沌思想加入到PSO算法是利用混沌对与粒子速度更新相关的参数进行自适应调整[6],在产生混沌序列时,采用logistic模型：

其中， $x_i^t$为$x_i$在第$t$步混沌演变后的值。 $x_i\in [0,1],1\le u\le 4$ 当$u=4$且$x_i\notin 0.25,0.5,0.75$时，系统表现出完全混沌特性。

本文混沌模拟退火粒子群算法由三部分组成：

采用混沌思想对粒子群参数惯性权值$w$ 和随机数 $rand$进行调整；
运用基本粒子群算法公式(1)、(2)引领粒子群的全局搜索方向；
运用模拟退火算法公式(3)引导粒子群的局部搜索方向。
算法流程图如图所示：

3. 函数选择

拟合函数表达式：

训练集数据选择： ，等间隔取100个点； ，等间隔取100个点。
测试集数据选择： ，等间隔取100个点； ,等间隔取100个点。

4. 实验仿真

4.1 实验参数设置

4.2 智能算法函数拟合结果

如图可知，基础粒子群算法和改进的粒子群算法在神经网络的非线性函数拟合中效果都很不错。

4.3 算法稳定性分析

如图，混沌模拟退火粒子群算法的全局最优值为4.7266，基本粒子群算法和模拟退火粒子群算法则全局最优值分别为5.7571、5.8775。

4.4 神经网络超参数影响分析

仅调整BP网络训练目标误差，选定基础粒子群算法作为网络参数更新算法，观察网络拟合效果。

由图对比可知，当训练误差率设置值偏大时，拟合的精度会下降，反映出的拟合效果也会有所降低。
在以上实验基础上，仅调整BP网络隐层节点数，依然选定基础粒子群算法作为网络参数更新算法，观察网络拟合效果。

由图可知，增加隐层节点数目后，拟合效果反而下降了。因此，在神经网络模型优化中，并不是越多的层数越好，有可能层数越多反而拟合效果下降。要求在具体的实践中选择合适的神经元个数。

5. 总结与展望

本文对BP神经网络模型分析及其改进方法，基于基本粒子群算法，提出其改进算法，并将其应用于BP网络非线性函数拟合中，得到了效果不错的函数拟合结果。仿真结果验证了粒子群算法在神经网络改进中的有效性，同时验证了不同超参数设置对拟合效果的影响，为自己后期智能算法与深度学习算法的学习夯实了基础。

开源地址。

本人小白，如果代码中有不妥之处，还望各路前辈批评指正。

参考文献

[1] Kennedy J, Eberhart RC. Particle swarm optimization [C]//Proceedings of the 1995 IEEE International Conference on Neural Networks.1995:1942-1948．
[2] TAO Ai . Summary of Particle Swarm Optimization.[J] College of Science .2005.
[3] XU Tian-dong, SUN Li-jun, GENG Yuan-jing, etc. Appliance on the Estimation of LAN Transportation Information based on Montel Carlo Method.[J]. Computer Engineering & Appliance, 2008(15):206-208,234.
[4] WEI Ping. Analysis & Research of Simulated Annealing Algorithm[J]. Equipment Manufacturing Technology,2008(7):1-3.
[5] GONG Chun, WANG Zheng-lin. Master of Optimization Calculation on MATLAB[M]. Publishing House of Electronics Indus-try,2012:491-505.
[6] OTTE,GREBOGIC,YORKEJA. Controlling chaos[J]. Physical Review Letters,1990，64（11）：1196-1199.
[7] 尤晓东, 苏崇宇, 汪毓铎. BP 神经网络算法改进综述[J]. 民营科技, 2018, 4.