基于粒子群PSO优化灰色神经网络的鞋销量预测

目标

背影
灰色模型原理
神经网络原理
灰色神经网络原理
粒子群算法的原理
粒子群改进灰色神经网络原理
基于粒子群PSO优化灰色神经网络的鞋销量预测
效果图
结果分析
展望
参考

背影
提前预测销量，然后再进货，有利于实时调整库存，提供资金利用率，本文基于粒子群PSO优化灰色神经网络的鞋销量预测

灰色模型

基本思想是用原始数据组成原始序列(0)，经累加生成法生成序列(1)，它可以弱化原始数据的随机性，使其呈现出较为明显的特征规律。对生成变换后的序列(1) 建立微分方程型的模型即GM模型。GM(1,1) 模型表示1阶的、1个变量的微分方程模型。GM(1,1) 模型群中，新陈代谢模型是最理想的模型。这是因为任何一个灰色系统在发展过程中，随着时间的推移，将会不断地有一些随即扰动和驱动因素进入系统，使系统的发展相继地受其影响。用GM(1,1) 模型进行预测，精度较高的仅仅是原点数据(0)(n) 以后的1到2个数据，即预测时刻越远预测的意义越弱[3]。而新陈代谢GM(1,1)模型的基本思想为越接近的数据，对未来的影响越大。也就是说，在不断补充新信息的同时，去掉意义不大的老信息，这样的建模序列更能动态地反映系统最新的特征，这实际上是一种动态预测模型。
有如下有点
1、不需要大量样本。
2、样本不需要有规律性分布。
3、计算工作量小。
4、定量分析结果与定性分析结果不会不一致。
5、可用于Recent、短期、中长期预测。
6、灰色预测准确度高。

神经网络的原理

BP神经网络的定义

人工神经网络无需事先确定输入输出之间映射关系的数学方程，仅通过自身的训练，学习某种规则，在给定输入值时得到最接近期望输出值的结果。作为一种智能信息处理系统，人工神经网络实现其功能的核心是算法。BP神经网络是一种按误差反向传播(简称误差反传)训练的多层前馈网络，其算法称为BP算法，它的基本思想是梯度下降法，利用梯度搜索技术，以期使网络的实际输出值和期望输出值的误差均方差为最小。

BP神经网络的基本结构

基本BP算法包括信号的前向传播和误差的反向传播两个过程。即计算误差输出时按从输入到输出的方向进行，而调整权值和阈值则从输出到输入的方向进行。正向传播时，输入信号通过隐含层作用于输出节点，经过非线性变换，产生输出信号，若实际输出与期望输出不相符，则转入误差的反向传播过程。误差反传是将输出误差通过隐含层向输入层逐层反传，并将误差分摊给各层所有单元，以从各层获得的误差信号作为调整各单元权值的依据。通过调整输入节点与隐层节点的联接强度和隐层节点与输出节点的联接强度以及阈值，使误差沿梯度方向下降，经过反复学习训练，确定与最小误差相对应的网络参数(权值和阈值)，训练即告停止。此时经过训练的神经网络即能对类似样本的输入信息，自行处理输出误差最小的经过非线形转换的信息。

bp神经网络的神经元

神经网络是一种模仿动物神经网络行为特征，进行分布式并行信息处理的算法数学模型。这种网络依靠系统的复杂程度，通过调整内部大量节点之间相互连接的关系，从而达到处理信息的目的。
神经网络由多个神经元构成，下图就是单个神经元的图1所示：

。。。。。。。。。。。。。。。。。。。。。。。。图1 ，神经元模型

bp神经网络激活函数及公式

![在这里插入图片描述](https://img-blog.csdnimg.cn/29edde342c3945939ad5945145ca8509.png

BP神经网络传递函数及公式

图2是Sigmoid函数和双极S函数的图像，其中Sigmoid函数的图像区域是0到1，双极S函数的区间是正负1，归一化的时候要和传递函数的区域相对应，不然，可能效果不好
神经网络就是将许多个单一的神经元联结在一起，这样，一个神经元的输出就可以是另一个神经元的输入。
例如，下图就是一个简单的神经网络：

粒子群算法的原理

粒子群优化算法（PSO）又翻译为粒子群算法、微粒群算法、或微粒群优化算法，PSO是由Kennedy和Eberhart共同提出，最初用于模拟社会行为，作为鸟群中有机体运动的形式化表示。自然界中各种生物体均具有一定的群体行为，Kennedy和Eberhart的主要研究方向之一是探索自然界生物的群体行为，从而在计算机上构建其群体模型。PSO是一种启发式算法，因为它很少或没有对被优化的问题作出假设，并且能够对非常大候选解决方案空间进行搜索。PSO算法初始化为一群随机的粒子，然后通过多次迭代找到最优解。每一次的迭代过程中，粒子通过本身所找到的最优解（被成为个体极值）和整个种群目前找到的最优解（被称为全局极值）来更新自己。也可以使用粒子本身的邻居（被称为局部极值）的极值来更新自己。

粒子群算法的主要参数

一、种群个数popsize，既算法中粒子的个数；
二、最大迭代次数gen，既算法迭代gen次后停止迭代；
三、种群维度dim，既需要优化的自变量个数；
四、种群位置pop，既每个粒子群的对应的自变量的值，一个粒子对应一组自变量，相当于一个解；
五、种群速度v，既粒子群每次迭代更新的飞行速度，粒子群位置更新的步长；
六、种群全局最优值gbest，既迭代过程中曾经出现的最优解，包括最优位置和对应的目标函数值；
七、个体最优，既每个粒子迭代过程中单个体曾经出现的个体最优解，，包括个体最优位置和对应的目标函数值；
八、个体学习因子c1，既个体最优解对粒子群飞行的影响能力；
九、全局学习因子c2，既全局最优值对粒子群飞行的影响能力；
十、惯性权重w，既个体位置所占的权重，权重越大，粒子群收敛越慢，全局搜索能力越强；

粒子群算法流程图

灰色神经网络模型

首先通过灰色预测建立预测模型，然后以灰色模型的输入为神经网络的输入，以灰色模型的误差为神经网络的输出，进行训练建模，然后通过训练好的BP神经网络，校正灰色模型的误差

基于粒子群改进的灰色神经网络的鞋销量预测

主要代码如下：
lc
clear
close all
num1 = xlsread(‘鞋.xlsx’);
num2 = num1(end,:);
mm = 1:2:19;
x = num2(mm)‘;
y = x(:,1);
num = [];
for ii = 1:8
z = y(ii:ii+2,1);
p(ii)=huise(z,1);
num =[ num z];
end
% y = [y;p(end)];
for ii = 9:10
z = y(end-2:end,1);
p(ii)=huise(z,1);
num =[ num z];
y = [y;p(end)];
end
p = p’;
p = [y(1:2);p];

ERR = p(3:end)-y(3:end);

pingjunjueduiwucha = sum(abs(ERR))/length§;
pingjunbaifenbiwucha = sum(abs(p-y))./length(y);
junfanggenwucha = sqrt(sum((p-y).^2)/length§);

[psobpoutput,net] = psobp(ERR,x(7:end,:));
n2 = [7 8 15 20 23];
zhenshi = y(n2);
huis = p(n2) ;
figure (2)
plot(zhenshi,‘k-*’)
hold on
plot(huis-psobpoutput’,‘r–o’)
hold on
plot(huis,‘b–^’)
hold off

legend(‘真实值’,‘灰色神经网络’,‘灰色预测’)

figure (3)
plot(zhenshi-huis,‘k-*’)
hold on
plot(zhenshi-(huis-psobpoutput’),‘r–o’)

hold off

legend(‘灰色预测误差’,‘灰色神经网络预测误差’)

%依次是优化后的隐含阈值；输入层到隐含层权值；输出层阈值；隐含层到输出层权值
aa1 = net.b{1}
aa2 = net.iw{1,1}
aa3 = net.b{2}
aa4 = net.lw{2,1}
function [EvaSamOut,net] = psobp(ERR,num)
n1 = randperm(19);
n2 = [7 8 15 20 23]-6;
input_train=num(n1,2:17)‘;
output_train=ERR(n1,1)’;
input_test=num(n2,2:17)‘;
output_test=ERR(n2,1)’;
[BPoutput1,error1,net] = bpp(ERR,num);
[AllSamInn,minAllSamIn,maxAllSamIn,AllSamOutn,minAllSamOut,maxAllSamOut]=premnmx(input_train,output_train);

EvaSamIn=input_test;
EvaSamInn=tramnmx(EvaSamIn,minAllSamIn,maxAllSamIn); % preprocessing
Ptrain = AllSamInn;
Ttrain = AllSamOutn;

indim=16;
hiddennum=20;
outdim=1;

% Initialize PSO
vmax=0.0151; % Maximum velocity
minerr=0.001; % Minimum error
wmax=0.90;
wmin=0.30;
% global itmax; %Maximum iteration number
itmax=100;
c1=2;
c2=2;
for iter=1:itmax
W(iter)=wmax-((wmax-wmin)/itmax)*iter; % weight declining linearly
end

%Between (m,n), (which can also be started from zero)
m=-1;
n=1;
% global N; % number of particles
N=20;
% global D; % length of particle
D=(indim+1)*hiddennum+(hiddennum+1)*outdim;
gbests = [reshape(net.IW{1,1},1,320) reshape(net.LW{2,1},1,20) reshape(net.b{1},1,20) reshape(net.b{2},1,1)] ;
% particles are initialized between (a,b) randomly
a=abs(gbests)*0.5+gbests;
b=-abs(gbests)*0.5+gbests;

% Initialize positions of particles
% rand(‘state’,sum(100*clock));
X = [];
for ii = 1:N
X =[X;a+(b-a).rand(1,D,1)]; %取值范围[-1,1] rand * 2 - 1 ,rand 产生[0,1]之间的随机数
end
%Initialize velocities of particles
V=0.2(m+(n-m)*rand(N,D,1));
%
% global fvrec;
MinFit=[];
BestFit=[];
net=newff(minmax(Ptrain),[hiddennum,outdim],{‘tansig’,‘tansig’},‘traingdx’);
fitness=fitcal(X,net,indim,hiddennum,outdim,D,Ptrain,Ttrain,minAllSamOut,maxAllSamOut);
fvrec(:,1,1)=fitness(:,1,1);
[C,I]=min(fitness(:,1,1));
MinFit=[MinFit C];
BestFit=[BestFit C];
L(:,1,1)=fitness(:,1,1); %record the fitness of particle of every iterations
B(1,1,1)=C; %record the minimum fitness of particle
gbest(1,:,1)=X(I,:,1); %the global best x in population

%Matrix composed of gbest vector
for p=1:N
G(p,:,1)=gbest(1,:);
end
for ii=1:N;
pbest(ii,:,1)=X(ii,:);
end
V(:,:,2)=W(1)V(:,:,1)+c1rand*(pbest(:,:,1)-X(:,:,1))+c2rand(G(:,:,1)-X(:,:,1));

for ni=1:N
for di=1:D
if V(ni,di,2)>vmax
V(ni,di,2)=vmax;
elseif V(ni,di,2)<-vmax
V(ni,di,2)=-vmax;
else
V(ni,di,2)=V(ni,di,2);
end
end
end
X(:,:,2)=X(:,:,1)+V(:,:,2);
for ni=1:N
for di=1:D
if X(ni,di,2)>1
X(ni,di,2)=1;
elseif X(ni,di,2)<-1
X(ni,di,2)=-1;
else
X(ni,di,2)=X(ni,di,2);
end
end
end

%******************************************************
for jj=2:itmax
disp(‘Iteration and Current Best Fitness’)
disp(jj-1)
disp(B(1,1,jj-1))
reset =1; % reset = 1时设置为粒子群过分收敛时将其打散，如果＝1则不打散
if reset1
bit = 1;
for k=1:N
bit = bit&(range(X(k,:))<0.02);
end
if bit1 % bit=1时对粒子位置及速度进行随机重置
for ik = 1:N
X(ik,:) = funx; % present 当前位置,随机初始化
X(ik,:) = [0.02rand()-0.01 0.02rand()-0.01]; % 速度初始化
end
end
end

% Calculation of new positions 
fitness=fitcal(X,net,indim,hiddennum,outdim,D,Ptrain,Ttrain,minAllSamOut,maxAllSamOut);

[C,I]=min(fitness(:,1,jj));
MinFit=[MinFit C];
BestFit=[BestFit min(MinFit)];
L(:,1,jj)=fitness(:,1,jj);
B(1,1,jj)=C;
gbest(1,:,jj)=X(I,:,jj);
[C,I]=min(B(1,1,:));
% keep gbest is the best particle of all have occured
if B(1,1,jj)<=C
gbest(1,:,jj)=gbest(1,:,jj);
else
gbest(1,:,jj)=gbest(1,:,I);
end
if C<=minerr
break
end
%Matrix composed of gbest vector
if jj>=itmax
break
end
for p=1:N
G(p,:,jj)=gbest(1,:,jj);
end
for ii=1:N;
[C,I]=min(L(ii,1,:));
if L(ii,1,jj)<=C
pbest(ii,:,jj)=X(ii,:,jj);
else
pbest(ii,:,jj)=X(ii,:,I);
end
end

V(:,:,jj+1)=W(jj)*V(:,:,jj)+c1*rand*(pbest(:,:,jj)-X(:,:,jj))+c2*rand*(G(:,:,jj)-X(:,:,jj));
for ni=1:N
    for di=1:D
        if  V(ni,di,jj+1)>vmax
            V(ni,di,jj+1)=vmax;
        elseif V(ni,di,jj+1)<-vmax
            V(ni,di,jj+1)=-vmax;
        else
            V(ni,di,jj+1)=V(ni,di,jj+1);
        end
    end
end 
X(:,:,jj+1)=X(:,:,jj)+V(:,:,jj+1);
for ni=1:N
    for di=1:D
        if  X(ni,di,jj+1)>1
            X(ni,di,jj+1)=1;
        elseif X(ni,di,jj+1)<-1
            X(ni,di,jj+1)=-1;
        else
            X(ni,di,jj+1)=X(ni,di,jj+1);
        end
    end
end

end
disp(‘Iteration and Current Best Fitness’)
disp(jj)
disp(B(1,1,jj))
disp(‘Global Best Fitness and Occurred Iteration’)
[C,I]=min(B(1,1,:));

% simulation network 网络拟合
for t=1:hiddennum
x2iw(t,:)=gbest(1,((t-1)indim+1):tindim,jj);
end
for r=1:outdim
x2lw(r,:)=gbest(1,(indimhiddennum+1):(indimhiddennum+hiddennum),jj);
end
x2b=gbest(1,((indim+1)*hiddennum+1)😄,jj);
x2b1=x2b(1:hiddennum).‘;
x2b2=x2b(hiddennum+1:hiddennum+outdim).’;
net=newff(minmax(AllSamInn),[20,1],{‘logsig’,‘tansig’},‘trainlm’);
net.IW{1,1}=x2iw;
net.LW{2,1}=x2lw;
net.b{1}=x2b1;
net.b{2}=x2b2;

%% BP网络训练
%网络进化参数
net.trainParam.epochs=5000;
net.trainParam.lr=0.1;
net.trainParam.goal=0.01;
% net.trainParam.show=100;
% net.trainParam.showWindow=1;
tic
%网络训练
net=train(net,AllSamInn,AllSamOutn);
toc
EvaSamOutn = sim(net,EvaSamInn);
EvaSamOut = postmnmx(EvaSamOutn,minAllSamOut,maxAllSamOut);%反归一化
error=EvaSamOut-output_test;
errormape=(EvaSamOut-output_test)./output_test;
p1 = sum(abs(error))/133;
MSE= mse(error)

figure(1)
grid
hold on
plot((BestFit),‘r’);
title(['粒子群算法优化bp ’ ‘最优代数＝’ I]);
xlabel(‘进化代数’);ylabel(‘误差’);
disp(‘适应度变量’);

% figure(2)
% grid
% plot(EvaSamOut,‘-^g’)
% hold on
% plot(BPoutput1,‘-ob’)
% hold on
% plot(output_test,‘-*r’);
%
% legend(‘粒子群优化BP预测输出’,‘BP预测输出’,‘期望输出’)
% title(‘粒子群优化BP网络预测输出’,‘fontsize’,12)
% ylabel(‘函数输出’,‘fontsize’,12)
% xlabel(‘样本’,‘fontsize’,12)

figure(4)
grid
plot(EvaSamOut,‘:og’)
hold on

hold on
plot(output_test,‘-*r’);

legend(‘粒子群优化BP预测输出’,‘真实误差’)
title(‘粒子群优化BP网络预测输出’,‘fontsize’,12)
ylabel(‘灰色预测的真实误差’,‘fontsize’,12)
xlabel(‘样本’,‘fontsize’,12)

%
% % figure(5)
% % plot(BPoutput1,‘:og’)
% % hold on
% % plot(output_test,‘-‘);
% % legend(‘预测输出’,‘期望输出’)
% % title(’·BP神经网络’,‘fontsize’,12)
% % ylabel(‘函数输出’,‘fontsize’,12)
% % xlabel(‘样本’,‘fontsize’,12)
% %预测误差
% error1=BPoutput1-output_test;
% junfanggen = mse(BPoutput1-output_test);
%
% figure(6)
% plot(error1,'-k’)
% title(‘BP网络预测误差’,‘fontsize’,12)
% ylabel(‘误差’,‘fontsize’,12)
% xlabel(‘样本’,‘fontsize’,12)
% %axis([1 2500 -0.5 0.5])
% figure(7)
% plot(error,‘ok-’);
% title(‘粒子群优化神经网络预测误差’)
% % axis([1 2500 -0.2 0.2])
%

% figure(3)
% plot(error1./output_test,‘*g-’);
% hold on
% plot((EvaSamOut-output_test)./output_test/2,‘*r-’)
% hold off
% title(‘粒子群优化BP神经网络预测误差百分比’)
% legend(‘BP预测输出’,‘粒子群优化BP预测输出’)
% figure(4)
% hist(errormape);
% title(‘PSO-BP神经网络预测误差频率分布直方图’);
% ylabel(‘频率（次）’,‘fontsize’,12)
% xlabel(‘相对误差’,‘fontsize’,12)

% MAE=(sum(abs(errormape)))/24 %绝对平均误差
% RMSE=sqrt((sum(errormape.^2))/24)%RMSE 均方根误差公式

end

结果图

结果分析

从结果图中可以看出，改进的灰色神经网络，比灰色模型预测准确率高很多，改进效果比较明显

展望

按照BP神经网络改进灰色模型的思路，可以用其他神经网络改进灰色模型，提高灰色模型的准确度，比如常见的SVM,CNN,RBF,ELMAN,RVM,DBN，LSTM等

参考

百科