MATLAB

%%%%%%%% K-fold ANN %%%%%%%
 
%数据归一化
%输出层有sigmoid归一化函数导致其输出值必在-1到1之间，意味着预测的粗糙度必在给的训练数据之间，无法超出训练数据的粗糙度范围，所以归一化输出数据的时候得多缩一些。
[pre,outputps]=mapminmax(pre_map)
[inputn,inputps]=mapminmax(input);
outputn=mapminmax('apply',output,outputps);
%%%[Y,PS] = mapminmax(X,YMIN,YMAX)
%%%X是预处理的数据，Ymin和Ymax是期望的每一行的最小值与最大值
 
cvFold = crossvalind('Kfold',length(inputn),3);       %第一个位置用于选用验证方法，这里选的是k折验证法；第二个位置填数据的组数；准备将该数据组分为几组，第三个位置就填几；
    H_net_tmp = newff(minmax(inputn),[9,1],{'tansig','tansig','purelin'},'traingdx');%PH为输入数据矩阵（必须是标准化后的），中括号里填输入输出的维数。
    H_net_tmp.trainParam.lr = 0.01; %学习率
    H_net_tmp.trainParam.mc = 0.9; %动量
 H_net_tmp.trainParam.epochs = 100; %最大迭代次数，我设置的是100.
 H_net_tmp.trainParam.goal = 1e-3;
 H_net_tmp.divideFcn = 'divideind';
 H_net_tmp.trainParam.showWindow = false;
 epoch_N = [];%用于存储没戏最后的训练次数
 for n = 1 : 3
        testIdx = (cvFold == n);
        trainIdx = ~testIdx;
        trInd = find(trainIdx);
        tstInd = find(testIdx);
                        
        H_net_tmp.divideParam.trainInd = trInd;
        H_net_tmp.divideParam.testInd = tstInd;
                        
        %training Network
        %reconfigured assignment of random initial weights
        H_net_tmp = init(H_net_tmp);%每次都对权重进行初始化，很重要，网上的很多代码都没有说明这点
        [H_net_tmp,tr] = train(H_net_tmp,inputn,outputn,'useParallel','yes','useGPU','yes');
        [~,I] = min(tr.tperf);
        I = I(end);
        epoch_N(n) = I;%I就是滴k次训练最优的那个次数
    end
 epoch_final = ceil(mean(epoch_N));
 
%测试样本归一化
%inputn_test=mapminmax('apply',input_test,inputps);
%%%?Y = mapminmax('apply',X,PS)在模式识别或者统计学里，PS是训练样本的数据的映射，即PS中包含了训练数据的最大值和最小值，这里的X是测试样本，对于测试样本来说，预处理应该和训练样本一致即最大值和最小值应该是训练集的最大值与最小值。假设y是测试样本，一共两个测试样本，则代码如下:x=[2,3,4,5,6;7,8,9,10,11]; y=[2,3;4,5];[xx,ps]=mapminmax(x,0,1);mapminmax('apply',y,ps)
    
%bp神经网络预测
an=sim(H_net_tmp,inputn);
 
%预测结果反归一化 
BPoutput=mapminmax('reverse',an,outputps)
 
%网络存储
%save data net inputps outputps

 
 
%网络预测图形
figure(1)
plot(BPoutput,':og')
hold on
plot(output,'-*');
%legend('预测输出','期望输出',
%%%%%%%% K-fold ANN %%%%%%%

%%%%%%%% PSO with mutation %%%%%%%
 %速度更新参数
c1 = 1.5;
c2 = 1.5;
 
maxgen = 50;  %迭代次数
sizepop = 20;  %种群规模，取大值收敛快
 
%个体和速度最大最小值
popmax1 = 5.6;popmin1 = 2.4;
popmax2 = 80;popmin2 = 40;
popmax3 = 80;popmin3 = 40;
 
Vmax1 = 0.3;Vmin1= -0.3;
Vmax2 = 4;Vmin2= -4;
Vmax3 = 4;Vmin3= -4;
 
for i=1:sizepop
    
    %随机产生一个种群
    pop(i,1)=1.6*rands(1,1)+4;             %初始化压力
    pop(i,2)=20*rands(1,1)+60              %初始化流速
    pop(i,3)=20*rands(1,1)+60              %初始化转速
  
    V(i,1)=0.3*rands(1,1);               %初始化粒子速度
    V(i,2)=4*rands(1,1);                 %初始化粒子速度
    V(i,3)=4*rands(1,1);                 %初始化粒子速度
    
    %计算粒子适应度值
    pop_test(:,i)=mapminmax('apply',pop(i,:)',inputps) %先对位置信息转置再归一化
    fitness(i)=sim(H_net_tmp,pop_test(:,i)); %送入神经网络计算适应度
end
[bestfitness bestindex]=min(fitness);   %fitness是个行向量，bestfitness是最小的那个值，bestindex是最小值所在位置
zbest=pop(bestindex,:);      %群体极值位置
gbest=pop;                   %用于存放每一个个体极值位置，现在这个gbest就是个初始值，后面适应度高的位置信息会挤掉当前位置信息。
fitnessgbest=fitness;        %个体极值适应度值
fitnesszbest=bestfitness;    %群体极值适应度值
 
for i=1:maxgen
    %粒子位置和速度更新
    for j=1:sizepop
        
        %粒子速度更新
        V(j,:)=V(j,:)+c1*rand*(gbest(j,:)-pop(j,:))+c2*rand*(zbest-pop(j,:));
        
        V(j,find(V(j,1)>Vmax1))=Vmax1;%找出所有大于Vmax1的数，用Vmax1代替，下类似
        V(j,find(V(j,2)>Vmax2))=Vmax2;%找出所有大于Vmax2的数，用Vmax2代替，下类似
        V(j,find(V(j,3)>Vmax3))=Vmax3;%找出所有大于Vmax3的数，用Vmax3代替，下类似
        V(j,find(V(j,1)<Vmin1))=Vmin1;%找出所有小于Vmin1的数，用Vmin1代替，下类似
        V(j,find(V(j,2)<Vmin2))=Vmin2;%找出所有小于Vmin2的数，用Vmin2代替，下类似
        V(j,find(V(j,3)<Vmin3))=Vmin3;%找出所有小于Vmin3的数，用Vmin3代替，下类似
        
        %粒子位置更新
        pop(j,:)=pop(j,:)+0.5*V(j,:);
        pop(j,find(pop(j,1)>popmax1))=popmax1;
        pop(j,find(pop(j,1)<popmin1))=popmin1;
        pop(j,find(pop(j,2)>popmax2))=popmax2;
        pop(j,find(pop(j,2)<popmin2))=popmin2;
        pop(j,find(pop(j,3)>popmax3))=popmax3;
        pop(j,find(pop(j,3)<popmin3))=popmin3;
        
        if rand>0.96     
            pop(j,:)=rands(1,3);
        end
        
        %新粒子适应度值
        pop_test(:,j)=mapminmax('apply',pop(j,:)',inputps) %先对位置信息转置再归一化
        fitness(j)=sim(H_net_tmp, pop_test(:,j));
    end
    
   %个体极值和群体极值更新
   for j=1:sizepop
       
       %个体极值更新
       if fitness(j)<fitnessgbest(j)
           gbest(j,:)=pop(j,:);
           fitnessgbest(j)=fitness(j);
       end
       
       %群体极值更新
       if fitness(j)<fitnesszbest
           zbest=pop(j,:);%找最优值的数据时就找这个，会有最优输入组合
           fitnesszbest=fitness(j);
       end
   end
   
   %每代最优值记录到yy数组中
   yy(i)=fitnesszbest;
   %反归一化
   yyy=mapminmax('reverse',yy,outputps)
end
 
%画出每代最优个体适应度值
plot(yyy)
title('最优个体适应度值','fontsize',12);
xlabel('进化次数','fontsize',12);ylabel('适应度','fontsize',12)

%%%%%%%% PSO with mutation %%%%%%%