NSPSOforZDT1.rar_Dominated_NSPSO_ZDT1

clear clc close all %% Set Parameter Npl=5; Ndc=5; Nsupplier=3; Ncustomer=4; No=2; %number of objective Nsol = 30; %%particle number Nitr = 300; %%number of iteration Nvar = Npl+Ndc; %% number of potential PL+DC Nrun =1; %% number of runs Cg = 0.5; Cp = 0.8; Cw = 0.95; X = zeros(Nsol,Nvar); Xt= zeros(Nsol,Nvar); F = zeros(Nsol,2);%% Fit_1,Fit_2 Ft= zeros(Nsol,2); w=[1/6 4/6 1/6]; %% load data load capacity_PL load capacity_DC load demand totaldemand=0; for i=1:Ncustomer for j=1:3 totaldemand=totaldemand+w(j)*demand(i,j); end end totaldemand=ceil(totaldemand); %% initializing swarm and velocities and position %for run = 1:Nrun %% 跑10次 start_time = cputime; Xij={}; Yjk={}; Zkl={}; X=round(rand(Nsol,Nvar)); %初始化 Npl+Ndc for i = 1:Nsol %% 跑30個particle %% check the number of weapons X(i,:) = CapOpened(X(i,:),Nvar,Npl,capacity_PL,capacity_DC,w,totaldemand); %Check Plant/DC capacity是否滿足demand %% distribution by minmum cost flow problem [Xij Yjk Zkl transportationcost] = Distribution(X(i,:),Npl,Nsupplier,Ncustomer,Ndc,w); %% Evaluate fitness function [F(i,:)] = MOCFLP_fit(X(i,:),Nvar,Npl,Ncustomer,Ndc,Yjk,Zkl,transportationcost); %% 挑出上個滿足限制式的解計算fitness function end %% 跑出30個Particle的初始解 X_p = X; %% 30*12 (初始 Pbest) F_p = F; scatter(F(:,1),F(:,2),'filled'); for iter = 1:Nitr %% 迭帶次數:300 (t=1) particle_init=[X F]; [particle front]=NDS_CD_cons(particle_init,Nvar,No); %population is nonDomSSOlist,Calculate Rank & crowded distance [nonDomSSOlist] = Resort(particle,Nvar,No,Nsol); % Resort according Rank & crowded distance nextPopList=[]; %% Start to update for i = 1:Nsol %% run 30 particles g=randperm(ceil(0.1*Nsol)); %random select a global best from specified top part. X_g= nonDomSSOlist(g(1),1:Nvar); f_g= nonDomSSOlist(g(1),Nvar+1:Nvar+No); for j = 1:Nvar r1 = rand(); if r1 < Cg Xt(i,j) = X_g(j); %% 新世代的解用Gbest取代 elseif r1 < Cp Xt(i,j) = X_p(i,j); %% 新世代的解用Pbest取代 elseif r1 < Cw Xt(i,j) = X(i,j); %% 沿用上一世代的值 else %%random movement to escape from local optimum Xt(i,j) = round(rand()); %% 偶數迴圈idx = 2;奇數迴圈idx = 1，且用upper/2為隨機值當作新世代的解 end end %% check the number of weapons Xt(i,:) = CapOpened(X(i,:),Nvar,Npl,capacity_PL,capacity_DC,w,totaldemand); %% Evaluate fitness function [Ft(i,:) Xij{i} Yjk{i} Zkl{i}] = MOCFLP_fit(X(i,:),Nvar,Npl,Ncustomer,Nsupplier,Ndc,w); nextPopList(i,:)=X(i,:); nextPopList(Nsol+i,:)=Xt(i,:); %Xt:該iteration中解 nextF(i,:)=F(i,:); nextF(Nsol+i,:)=Ft(i,:); end %% 更新前解與更新後形成2N 找pbest nextPopList=[nextPopList nextF]; [particle_inter_sorted front]=NDS_CD_cons(nextPopList,Nvar,No); %population is nonDomSSOlist %% Replacement - N nextSSOList=replacement(particle_inter_sorted,front,Nsol); X_p=nextSSOList(:,1:Nvar); %pbest F_p=nextSSOList(:,Nvar+1:Nvar+2); X=Xt; F=Ft; end hold on [nextSSOList] = Resort(nextSSOList,Nvar,No,Nsol); e= find(nextSSOList(:,Nvar+No+2)==1); scatter(nextSSOList(1:length(e),Nvar+1),nextSSOList(1:length(e),Nvar+2),'r','filled'); end_time=cputime; runtime = end_time - start_time; end %run end %---------------------------------------------- Q=[]; nextSSOList=sortrows(nextSSOList,Nvar+1); paretoset(run).trial=nextSSOList(:,1:Nvar+No+1); Q = [Q; paretoset(run).trial]; % Combining Pareto solutions obtained in each run end_time = cputime; runtime= end_time - start_time; hold on if run==1 plot(nextSSOList(:,Nvar+1),nextSSOList(:,Nvar+2),'*') else [pareto_filter front]=NDS_CD_cons(Q,Nvar,No); % Applying non domination sorting on the combined Pareto solution set rank1_index=find(pareto_filter(:,Nvar+No+2)==1); % Filtering the best solutions of rank 1 Pareto pareto_rank1=pareto_filter(rank1_index,1:Nvar+No); plot(pareto_rank1(:,Nvar+1),pareto_rank1(:,Nvar+2),'*') ; % Final Pareto plot end xlabel('objective function 1') ylabel('objective function 2')