zdt2_mymopso.asv

硕士论文三个点中的一个

function [stPop, actEvals] = ZDT2_mymopso;
% leadCh:
% 指导元
% 2种选择
%   1) 随机一元 ronelead
%   2) 多元     mlead

% 测试适应度函数的作用
% 参数设置：
%   适应度：            基于共享法
%   变异概率mut_now:    0,无变异
%   惯性权值w_now:      常量＝0.2
%   替换差解:           有
%   指导元:             见上

n=30;
fCount=2;
o_share1=0.01;
o_share2=1;

gBCount=0;
fitInd = n + fCount + 1;
pInd = fitInd + 1;

popsize=100;
evaluations=200;
% % % w_start = 0.95;   %Initial inertia weight's value
% % % w_end = 0.15;       %Final inertia weight
% % % w_varyfor = floor(0.7*evaluations); %Weight change step. Defines total number of evaluations for which weight is changed.
% % % w_now = w_start;
w_now = 0.2;
% % % inertdec = (w_start-w_end)/w_varyfor; %Inertia weight's change per iteration
c1=2;
c2=2;
vMax=1;
replaceRate = 1;

% % % mut_start = 0.9; %最初的变异概率
% % % mut_end = 0;
% % % mut_varyfor = floor(0.7 * evaluations);
% % % mut_now = mut_start;
% % % mutdec = (mut_start - mut_end) / mut_varyfor;
mut_now = 0;



%初始化种群，并计算f1.f2的值存放于矩阵str的最后两列，solNum_tp存放各种f1,f2值各种组合下解的数目
initialPop=InitializePop(n,popsize);
for i=1:popsize
    for j=1:n
        Vel(i,j)=0;
    end
end

% %计算机种群适应值，并按适应值从高到低排序，得到新的种群stPop，f输出相应的f1,f2,F的值i_pop,i_n, i_fitInd, o_share1,o_share2)
% % % if (fitCh == 0)
	[fPop fD fN]=FitnessEvaluate(initialPop,n,fitInd,o_share1,o_share2);
	F=fPop(:,n+3);
	f=sort(F);
	stPop=sortD(fPop,F,f); 
	D=sortD(fD,F,f);
	N=sortD(fN,F,f);
% % % else
% % %     [fPop]=FitnessEvaluate_RandSum(initialPop,n,fitInd);
% % % 	F=fPop(:,n+3);
% % % 	f=sort(F);
% % % 	stPop=sortD(fPop,F,f);     
% % % end

for (i = 1:popsize)
    stPop(i, pInd) =i; %记录粒子编号，用于后面pBest排序
end
t_pop=strcat('pop=',num2str(popsize));
t_evals=strcat('evaluations=0');
output2figure(1,fPop,n,t_pop,t_evals);

% 更新全局最优(多个)和局部最优
[gBest gBCount]=GetStGBest(stPop, fitInd);
pBest=stPop;
%pBest
%gBest

% 开始进化计算    
tempPop=[];
tempPop=stPop;
allPareto = 0;
i = 1;
while (i <= evaluations) && (allPareto == 0)
   % i
% % %     % Update the value of the inertia weight w
% % %     if (i<=w_varyfor) & (i > 1)
% % %         w_now = w_now - inertdec; %Change inertia weight ＃＃ inertdec = (w_start-w_end)/w_varyfor;
% % %     end % rank the w

    %计算新的速度和位置
    for particle=1:popsize
        if (gBCount==1)
            gBi=1;
        else  
            gBi=ceil(rand(1)*(gBCount-1));
        end  
            
        for dem=1:n
            gBAdd = c2 *rand(1)*( gBest(gBi,dem)-tempPop(particle,dem) );
            newVel=w_now*Vel(particle,dem) + c1 *rand(1)*( pBest(particle,dem)-tempPop(particle,dem) ) + gBAdd; 
            if (newVel > vMax)
                newVel=vMax;
            elseif (newVel < -vMax)
                newVel=-vMax;
            end
            %newVel(1,1)
            Vel(particle,dem)=newVel;
            %更新位置
            newPox=tempPop(particle,dem)+newVel;
            if (newPox > 1)
                newPox = 1;
            elseif (newPox < 0)
                newPox = 0;
            end                      
            tempPop(particle,dem)=newPox;
        end
        
% % %         if (dem==1) %选一个变异的自变量xi
% % %             mutInd=1;
% % %         else  
% % %             mutInd=ceil(rand(1)*(dem-1));
% % %         end
% % %         % Update the value of the mutation rate
% % %             if (i<=mut_varyfor) & (i > 1)
% % %                 switch mutCh
% % %                     case 1  % 线性变小
% % %                         mut_now =  mut_now - mutdec; %Change mutation rate ＃＃ mutdec = (mut_start-mut_end)/mut_varyfor;
% % %                     case 2  % 固定概率
% % %                         mut_now =  mutFixed;
% % %                     case 3  % 无变异
% % %                         mut_now = 0;
% % %                 end
% % %             end % rank the w
% % %             mut = rand(1);
% % %             if mut <= mut_now
% % %                tempPop(particle,mutInd)= rand(1); %newPox *(1 + rand(1));
% % %             end
        
        for fk=(n+1):(n+3)
            tempPop(particle,fk)=0;
        end
        %tempPop(particle, :)
    end
    %更新粒子的函数值
    for j=1:popsize
        sumx=0;
        for k=2:n
            sumx=sumx+tempPop(j,k);
        end
        tempPop(j,n+1)=tempPop(j,1);
        gx=1+9*(sumx/(n-1));
        tempPop(j,n+2)=gx*(1-(tempPop(j,n+1)/gx)^2);
    end
    
    %更新粒子的适应度
    allPop=[tempPop; gBest]; %pop为当前的新种群和前一代的全局最优个体的并集
    allPop(:, fitInd)=0;
    
    [pop,aD,aN]=FitnessEvaluate(allPop,n,fitInd,o_share1,o_share2);    
   
    ftempPop=pop;
    %ftempPop
    %更新全局最优位置和局部最优位置 
    F=ftempPop(:,n+3);
    f=sort(F);
    sortedPop=sortD(ftempPop,F,f);
    % 更新全局最优(多个)
    gBCount_before= gBCount;
    [gBest gBCount]=GetStGBest(sortedPop, fitInd);
    
    % 将种群中适应度最小的gBCount_before个粒子替换成新的全局最优
    for (repPop = popsize+1 : popsize + gBCount_before)
        repRan = 1;%rand(1);
        if (repRan >= replaceRate) % 只有大于替换概率replaceRate时，才替换（为了最后的多样性，可进一步考虑变replaceRate）
            badInd = sortedPop(repPop, pInd);
            if (badInd >= 1) && (badInd <= popsize) % 在正常范围，参见GetStGBest中对全局最优解的粒子号定义（为-1）
                if (gBCount==1)
                    selGBInd=1;
                else  
                    selGBInd=ceil(rand(1)*(gBCount-1));
                end
                %strcat('Replace: ', num2str(badInd), ' to ', num2str(selGBInd))
                tempPop(badInd,1:n + 2) = gBest(selGBInd,1:n+2);  
            end
        end % end of if (repRan >= replaceRate)
    end
    
    %更新局部最优
    for j=1:popsize
         if (tempPop(j,n+1)<=stPop(j,n+1)) && (tempPop(j,n+2)<=stPop(j,n+2)) && ...
           ((tempPop(j,n+1)+tempPop(j,n+2))<(stPop(j,n+1)+stPop(j,n+2)))
            pBest(j,:)=tempPop(j,:);
         end
    end
    %重新计算粒子在当前种群中的适应度
    pop1=[];
    pop1=tempPop;
    
    [pop,D,N]=FitnessEvaluate(pop1,n,fitInd,o_share1,o_share2);    

    %更新局部最优
    %for j=1:popsize
    %     if (pop(j,n+1)<=stPop(j,n+1)) && (pop(j,n+2)<=stPop(j,n+2)) && ...
    %       ((pop(j,n+1)+pop(j,n+2))<(stPop(j,n+1)+stPop(j,n+2)))
    %        pBest(j,:)=pop(j,:);
    %    end
    %end
    stPop=pop;      
%         
    if (mod(i,20)==0)
        t_pop=strcat('pop=',num2str(popsize));
        t_evals=strcat('evaluations=',num2str(i));
        output2figure(i,stPop,n,t_pop,t_evals);
    end

    averD = mean(D);
    if (averD == 0) %&& (i > 200)
        allPareto = 1;
    end
    i = i + 1;
end
%averN = mean (N);
% 输出最终解的情况
t_pop=strcat('pop=',num2str(popsize));
actEvals = i - 1;
t_evals=strcat('evaluations=',num2str(actEvals));%, 'averN',num2str(averN));
output2figure(actEvals,stPop,n,t_pop,t_evals);%
% % % stPop
% % % D
% % % N