我根据论坛里找到的粒子群程序代码，试着修改了一下来求解一个六维函数最小值。在程序位置更新过程中，我也把各个粒子位置的约束条件写了上去，可是运行出来后发现粒子位置的取值一直会超出约束条件，导致结果一直错误，我是个菜鸟，写程序的基本功不行，琢磨了好久也没琢磨出程序的错误在什么地方，不知道有没有大神能帮我看一看？麻烦了!
适应度是：F=50*x(1)+14.5*x(2)+0.03*x(2)^2+4.2*x(3)+0.03*x(3)^2+0.031*x(2)*x(3)+36*x(4) +0.0435*x(4)^2+0.6*x(5)+0.02*x(5)^2+0.011*x(4)*x(5)+23.4*x(6)+4000，
约束条件是：0≤x(1)≤150，
81≤x(2)≤247，
0≤x(3)≤180，
40≤x(4)≤128.5，
0≤x(5)≤136.6，
0≤x(6)≤2695.2
x(2)+x(1)+x(4)=200
x(3)+x(5)+x(6)=115
然后根据条件写的粒子初始位置为，X(1)=150*rand，X(2)=166*rand+81，X(3)=180*rand，，X(4)=200-X(1)-X(2)，X(5)=136.6*rand，X(6)=115-X(3)-X(5)。
程序：
- %清空運行環境
- clc;
- clear;
- %參數初始化
- %速度更新參數
- c1=2;
- c2=2;
- maxgen=100; %迭代次數(進化次數)
- sizepop=100; %種群規模
- popmax1=150;
- popmin1=0;
- popmax2=247;
- popmin2=82;
- popmax3=180;
- popmin3=0;
- popmax4=119;
- popmin4=40;
- popmax5=136.6;
- popmin5=0;
- popmax6=115;
- popmin6=0;
- %產生初始粒子
- pop = zeros(sizepop,6);
- V = zeros(sizepop,6);
- fitness = zeros(sizepop,1);
- for i=1:sizepop
- %隨機產生一個種群
- pop(i,1)=150*rand; %初始化粒子
- pop(i,2)=166*rand+81;
- pop(i,3)=180*rand;
- pop(i,4)=200-pop(i,1)-pop(i,2);
- pop(i,5)=136.6*rand;
- pop(i,6)=115-pop(i,3)-pop(i,5);
- V(i,
 =1*rand(1,6); %初始化速度 - %計算粒子適應度
- fitness(i)=fun(pop(i,);%粒子的適應度
- end
- %找出最好的粒子
- [bestfitness, bestindex]=min(fitness);
- gbest = pop; %個體最佳解
- fitnessgbest=fitness; %個體最佳適應度值
- zbest = pop(bestindex,; %全局最佳解
- fitnesszbest=bestfitness;%全局最佳適應度
- %迭代尋優
- yy = zeros(maxgen,1);
- for i=1:maxgen
- for j=1:sizepop
- %速度更新
- V(j, = V(j,+c1*rand*(gbest(j,-pop(j,)+c2*rand*(zbest-pop(j,);
- %種群更新
- pop(j,=pop(j,+V(j,:);
- pop(j,pop(j,1)>popmax1)=popmax1;
- pop(j,pop(j,2)>popmax2)=popmax2;
- pop(j,pop(j,3)>popmax3)=popmax3;
- pop(j,pop(j,4)>popmax4)=popmax4;
- pop(j,pop(j,5)>popmax5)=popmax5;
- pop(j,pop(j,6)>popmax6)=popmax6;
- pop(j,pop(j,1)<popmin1)=popmin1;
- pop(j,pop(j,2)<popmin2)=popmin2;
- pop(j,pop(j,3)<popmin3)=popmin3;
- pop(j,pop(j,4)<popmin4)=popmin4;
- pop(j,pop(j,5)<popmin5)=popmin5;
- pop(j,pop(j,6)<popmin6)=popmin6;
- %新適應度值
- fitness(j)=fun(pop(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)%#ok
- end
- end
- %%結果分析
- yy(i)=fitnesszbest;
- end
- plot(yy)
- title(['適應度曲線' '終止代數=' num2str(maxgen)]);
- xlabel('進化次數');
- ylabel('適應度');
- zbest%#ok
复制代码
- function y=fun(x)
- y=50*x(1)+15.5*x(2)+0.03*x(2)^2+4.2*x(3)+0.03*x(3)^2+0.031*x(2)*x(3)+36*x(4)+0.0435*x(4)^2+0.6*x(5)+0.02*x(5)^2+0.011*x(4)*x(5)+23.4*x(6)+4000;
复制代码
|
|
发表于 2018-11-20 10:31:52
楼主