%% 初始化clearclose allclcwarning off
%% 数据读取load('data.mat')%输入输出数据input=data(:,1:end-1);    %data的第一列-倒数第二列为特征指标output=data(:,end);  %data的最后面一列为输出的指标值testNum = 21; %测试样本量%% 划分训练集、测试集input_train = input';output_train =output';input_test =input';output_test =output';
%% 数据归一化[inputn,inputps]=mapminmax(input_train,0,1);[outputn,outputps]=mapminmax(output_train);inputn_test=mapminmax('apply',input_test,inputps);
%% 获取输入层节点、输出层节点个数inputnum=size(input,2);outputnum=size(output,2);disp('/////////////////////////////////')disp('神经网络结构...')disp(['输入层的节点数为：',num2str(inputnum)])disp(['输出层的节点数为：',num2str(outputnum)])disp(' ')disp('隐含层节点的确定过程...')
%确定隐含层节点个数%采用经验公式hiddennum=sqrt(m+n)+a，m为输入层节点个数，n为输出层节点个数，a一般取为1-10之间的整数MSE=1e+5; %初始化最小误差for hiddennum=fix(sqrt(inputnum+outputnum))+1:fix(sqrt(inputnum+outputnum))+10        %构建网络    net=newff(inputn,outputn,hiddennum);    % 网络参数    net.trainParam.epochs=1000;         % 训练次数    net.trainParam.lr=0.01;                   % 学习速率    net.trainParam.goal=0.000001;        % 训练目标最小误差    % 网络训练    net=train(net,inputn,outputn);    an0=sim(net,inputn);  %仿真结果    mse0=mse(outputn,an0);  %仿真的均方误差    disp(['隐含层节点数为',num2str(hiddennum),'时，训练集的均方误差为：',num2str(mse0)])        %更新最佳的隐含层节点    if mse0<MSE        MSE=mse0;        hiddennum_best=hiddennum;    endenddisp(['最佳的隐含层节点数为：',num2str(hiddennum_best),'，相应的均方误差为：',num2str(MSE)])
%% 构建最佳隐含层节点的BP神经网络disp(' ')disp('标准的BP神经网络：')net0=newff(inputn,outputn,hiddennum_best,{'tansig','purelin'},'trainlm');% 建立模型
%网络参数配置net0.trainParam.epochs=1000;         % 训练次数，这里设置为1000次net0.trainParam.lr=0.01;                   % 学习速率，这里设置为0.01net0.trainParam.goal=0.00001;                    % 训练目标最小误差，这里设置为0.0001net0.trainParam.show=25;                % 显示频率，这里设置为每训练25次显示一次net0.trainParam.mc=0.01;                 % 动量因子net0.trainParam.min_grad=1e-6;       % 最小性能梯度net0.trainParam.max_fail=6;               % 最高失败次数
%开始训练net0=train(net0,inputn,outputn);
%预测an0=sim(net0,inputn_test); %用训练好的模型进行仿真
%预测结果反归一化与误差计算test_simu0=mapminmax('reverse',an0,outputps); %把仿真得到的数据还原为原始的数量级%误差指标[mae0,mse0,rmse0,mape0,error0,errorPercent0]=calc_error(output_test,test_simu0);
%% 鲸鱼优化算法寻最优权值阈值disp(' ')disp('WOA优化BP神经网络：')net=newff(inputn,outputn,hiddennum_best,{'tansig','purelin'},'trainlm');% 建立模型
%网络参数配置net.trainParam.epochs=1000;         % 训练次数，这里设置为1000次net.trainParam.lr=0.01;                   % 学习速率，这里设置为0.01net.trainParam.goal=0.00001;                    % 训练目标最小误差，这里设置为0.0001net.trainParam.show=25;                % 显示频率，这里设置为每训练25次显示一次net.trainParam.mc=0.01;                 % 动量因子net.trainParam.min_grad=1e-6;       % 最小性能梯度net.trainParam.max_fail=6;               % 最高失败次数
%% 初始化WOA参数popsize=30;   %初始种群规模maxgen=50;   %最大进化代数dim=inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+outputnum;    %自变量个数lb=repmat(-2,1,dim);    %自变量下限ub=repmat(2,1,dim);   %自变量上限%初始化位置向量和领导者得分Leader_pos=zeros(1,dim);Leader_score=10^20;   
%% 初始化种群for i=1:dim    ub_i=ub(i);    lb_i=lb(i);   Positions(:,i)=rand(popsize,1).*(ub_i-lb_i)+lb_i;endcurve=zeros(maxgen,1);%初始化收敛曲线
%% 循环开始h0=waitbar(0,'WOA-BP optimization...');for t=1:maxgen    for i=1:size(Positions,1)%对每个个体一个一个检查是否越界        %对每个个体一个一个检查是否越界        % 返回超出搜索空间边界的搜索代理        Flag4ub=Positions(i,:)>ub;        Flag4lb=Positions(i,:)<lb;        Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;%超过最大值的设置成最大值，超过最小值的设置成最小值        %目标函数值的计算        fit(i)=fitness(Positions(i,:),inputnum,hiddennum_best,outputnum,net,inputn,outputn,output_train,inputn_test,outputps,output_test);                % 更新领导者位置        if fit(i)<Leader_score            Leader_score=fit(i);            Leader_pos=Positions(i,:);        end    end        a=2-t*((2)/maxgen);    a2=-1+t*((-1)/maxgen);    %参数更新    for i=1:size(Positions,1)        r1=rand();r2=rand();        A=2*a*r1-a;        C=2*r2;               b=1;        l=(a2-1)*rand+1;                p = rand();                for j=1:size(Positions,2)%对每一个个体地多维度进行循环运算            %收缩包围机制            if p<0.5                if abs(A)>=1                    rand_leader_index = floor(popsize*rand()+1);%floor将 X 的每个元素四舍五入到小于或等于该元素的最接近整数                    X_rand = Positions(rand_leader_index, :);                    D_X_rand=abs(C*X_rand(j)-Positions(i,j));                    Positions(i,j)=X_rand(j)-A*D_X_rand;                elseif abs(A)<1                    D_Leader=abs(C*Leader_pos(j)-Positions(i,j));                    Positions(i,j)=Leader_pos(j)-A*D_Leader;                end                %螺旋更新位置            elseif p>=0.5                distance2Leader=abs(Leader_pos(j)-Positions(i,j));                Positions(i,j)=distance2Leader*exp(b.*l).*cos(l.*2*pi)+Leader_pos(j);            end        end    end    curve(t)=Leader_score;    waitbar(t/maxgen,h0)endclose(h0)setdemorandstream(pi);
%% 绘制进化曲线figureplot(curve,'r-','linewidth',2)xlabel('进化代数')ylabel('均方误差')legend('最佳适应度')title('WOA的进化曲线')w1=Leader_pos(1:inputnum*hiddennum_best);         %输入层到中间层的权值B1=Leader_pos(inputnum*hiddennum_best+1:inputnum*hiddennum_best+hiddennum_best);   %中间各层神经元阈值w2=Leader_pos(inputnum*hiddennum_best+hiddennum_best+1:inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum);   %中间层到输出层的权值B2=Leader_pos(inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+1:inputnum*hiddennum_best+hiddennum_best+hiddennum_best*outputnum+outputnum);   %输出层各神经元阈值%矩阵重构net.iw{1,1}=reshape(w1,hiddennum_best,inputnum);net.lw{2,1}=reshape(w2,outputnum,hiddennum_best);net.b{1}=reshape(B1,hiddennum_best,1);net.b{2}=reshape(B2,outputnum,1);%% 优化后的神经网络训练net=train(net,inputn,outputn);%开始训练，其中inputn,outputn分别为输入输出样本%% 优化后的神经网络测试an1=sim(net,inputn_test);test_simu1=mapminmax('reverse',an1,outputps); %把仿真得到的数据还原为原始的数量级%误差指标[mae1,mse1,rmse1,mape1,error1,errorPercent1]=calc_error(output_test,test_simu1);
%% 作图figure% plot(test_simu0,'r-v','linewidth',0.1,'markerfacecolor','r')% hold onplot(test_simu1,'b-o','linewidth',0.1,'markerfacecolor','w')hold onscatter(1:testNum, output_test,'r','*')legend('WOA-BP预测值','真实值')xlabel('测试样本编号')ylabel('指标值')title('WOA-BP神经网络预测值和真实值对比图')
figureplot(error0,'rv-','markerfacecolor','r')hold onplot(error1,'ko-','markerfacecolor','k')legend('BP预测误差','WOA-BP预测误差')xlabel('测试样本编号')ylabel('预测偏差')title('WOA优化前后的BP神经网络预测值和真实值误差对比图')
% disp(' ')% disp('/////////////////////////////////')% disp('打印结果表格')% disp('样本序号     实际值      BP预测值  WOA-BP值   BP误差   WOA-BP误差')% for i=1:testNum%     disp([i output_test(i),test_simu0(i),test_simu1(i),error0(i),error1(i)])% end