%% 该代码为基于PSO和BP网络的预测
%% 清空环境
clc
clear
%======原始数据输入========
p=[2.13 93.33 223.37 292.69 391.28 434.78 384.44 306.82 73.45 6.36 0 11 3 9 60;
0 3.57 16.05 29.35 32.79 24.41 31.7 48.76 46.28 7.74 0 7 -3 4 140;
0.29 62.58 226.56 415.97 521.5 526.53 478.68 369.35 212.08 41.56 0.7 7 -2 4.5 70;
5.32 105.55 250.57 356.98 270.46 253.85 125.34 201.06 67.63 16.2 0.07 7 -3 4 150;
1.7 68.61 242.08 423.18 510.63 555.61 525.2 389.15 214.82 39.83 0.81 8 2 7 80;
0.89 74.81 170.84 210.63 259.03 247 400.11 341.39 128.1 25.1 0.12 9 -4 4.5 130;
1.92 59.3 233.79 381.45 470.54 502.08 470.34 372.55 204.01 33.48 0.16 2 -7 -0.5 160;
1.38 62.47 284.79 427.32 550.2 596.79 511.49 400.28 256.57 37.52 0.33 3 -6 0.5 170;
2.49 74.25 340.92 497.27 553.45 617.26 588.96 436.19 258.96 39.5 0.68 5 -4 2.5 180;
1.57 63.88 266.29 404.16 503.21 489.28 388.83 271.36 139.09 17.1 0 6 1 5.5 90;
1.18 40.47 168.13 289.84 381.33 397.83 340.78 246.78 120.32 21.32 0.14 8 4 8 110;
1.23 50.97 221.57 350.71 401.97 447.42 384.44 367.57 226.1 49.36 0.1 6 3 6.5 50;
0.25 3.85 27.56 21.37 54.21 64.84 23.73 19.55 3.22 0 0 8 4 8 20;
0 0.83 2.76 16.63 13.01 24.87 28.71 39.54 15.55 1.1 0 9 4 8.5 40;
0 0.83 7.68 33.41 52.62 34.62 24.26 15.74 4.24 0.23 0 8 2 7 30;
0 0.58 7.75 14.36 21.16 32.29 29.41 29.97 16.23 2.22 0.7 7 1 6 10;
0 0 1.35 6.84 16.99 7.39 5.73 6.68 6.87 1.69 0 7 -3 4 180]';
%===========期望输出=======
t=[0 3.57 16.05 29.35 32.79 24.41 31.7 48.76 46.28 7.74 0;
0.29 62.58 226.56 415.97 521.5 526.53 478.68 369.35 212.08 41.56 0.7;
5.32 105.55 250.57 356.98 270.46 253.85 125.34 201.06 67.63 16.2 0.07;
1.7 68.61 242.08 423.18 510.63 555.61 525.2 389.15 214.82 39.83 0.81;
0.89 74.81 170.84 210.63 259.03 247 400.11 341.39 128.1 25.1 0.12;
1.92 59.3 233.79 381.45 470.54 502.08 470.34 372.55 204.01 33.48 0.16;
1.38 62.47 284.79 427.32 550.2 596.79 511.49 400.28 256.57 37.52 0.33;
2.49 74.25 340.92 497.27 553.45 617.26 588.96 436.19 258.96 39.5 0.68;
1.57 63.88 266.29 404.16 503.21 489.28 388.83 271.36 139.09 17.1 0;
1.18 40.47 168.13 289.84 381.33 397.83 340.78 246.78 120.32 21.32 0.14;
1.23 50.97 221.57 350.71 401.97 447.42 384.44 367.57 226.1 49.36 0.1;
0.25 3.85 27.56 21.37 54.21 64.84 23.73 19.55 3.22 0 0;
0 0.83 2.76 16.63 13.01 24.87 28.71 39.54 15.55 1.1 0;
0 0.83 7.68 33.41 52.62 34.62 24.26 15.74 4.24 0.23 0;
0 0.58 7.75 14.36 21.16 32.29 29.41 29.97 16.23 2.22 0.7;
0 0 1.35 6.84 16.99 7.39 5.73 6.68 6.87 1.69 0;
0.34 58.48 260.31 633.53 640.07 598.47 594.65 480.97 408.77 79.82 1.83]';
p_test=[0.34 58.48 260.31 633.53 640.07 598.47 594.65 480.97 408.77 79.82 1.83 3 -1 3 110]';
t_test=[2.62 67.85 230.24 406.66 568.79 545.62 250.34 155.36 56.44 0 0]';
%节点个数
inputnum=17;
hiddennum=15;
outputnum=17;
%训练数据和预测数据
input_train=p';
input_test=p_test';
output_train=t';
output_test=t_test';
%选连样本输入输出数据归一化
[inputn,inputps]=mapminmax(input_train);
[outputn,outputps]=mapminmax(output_train);
%构建网络
net=newff(inputn,outputn,hiddennum);
% 参数初始化
%粒子群算法中的两个参数
c1 = 1.49445;
c2 = 1.49445;
maxgen=2; % 进化次数
sizepop=20; %种群规模
Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
for i=1:sizepop
pop(i,:)=5*rands(1,21);
V(i,:)=rands(1,21);
fitness(i)=fun(pop(i,:),inputnum,hiddennum,outputnum,net,inputn,outputn);
end
% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:); %全局最佳
gbest=pop; %个体最佳
fitnessgbest=fitness; %个体最佳适应度值
fitnesszbest=bestfitness; %全局最佳适应度值
%% 迭代寻优
for i=1:maxgen
i
for j=1:sizepop
%速度更新
V(j,:) = V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));
V(j,find(V(j,:)>Vmax))=Vmax;
V(j,find(V(j,:)<Vmin))=Vmin;
%种群更新
pop(j,:)=pop(j,:)+0.2*V(j,:);
pop(j,find(pop(j,:)>popmax))=popmax;
pop(j,find(pop(j,:)<popmin))=popmin;
%自适应变异
pos=unidrnd(21);
if rand>0.95
pop(j,pos)=5*rands(1,1);
end
%适应度值
fitness(j)=fun(pop(j,:),inputnum,hiddennum,outputnum,net,inputn,outputn);
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(i)=fitnesszbest;
end
%% 结果分析
plot(yy)
title(['适应度曲线 ' '终止代数＝' num2str(maxgen)]);
xlabel('进化代数');ylabel('适应度');
x=zbest;
%% 把最优初始阀值权值赋予网络预测
% %用遗传算法优化的BP网络进行值预测
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=B2;
%% BP网络训练
%网络进化参数
net.trainParam.epochs=1000;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
%网络训练
[net,per2]=train(net,inputn,outputn);
%% BP网络预测
%数据归一化
inputn_test=mapminmax('apply',input_test,inputps);
an=sim(net,inputn_test);
test_simu=mapminmax('reverse',an,outputps);
error=test_simu-output_test;
figure(2)
plot(error)
结果一直显示
??? Index exceeds matrix dimensions.
Error in ==> fun at 13
w1=x(1:inputnum*hiddennum);
Error in ==> PSO at 78
fitness(i)=fun(pop(i,:),inputnum,hiddennum,outputnum,net,inputn,outputn);
矩阵哪里出现错误了
