发布日期：2024-08-08 06:24 点击次数：156
该算法浮浅易懂草榴社区邀请码，很适合刚初学多意见算法的。且思改造其他单意见优化算法为多意见的，悉数可以在此算法框架上告成修改！今天从头推出一下。
多意见粒子群算法是期骗最世俗，亦然最经典的多意见寻优算法。各式硕士博士著作，王人将其期骗在各式万般的畛域。今天就为公共带来一期多意见粒子群算法。
与网上大多半多意见粒子群代码不同，本期给出的多意见粒子群优化算法，唯唯一个剧本和一个函数，很适合生手学习，何况出图精湛！
在经典的多意见测试函数“ZDT1”，“ZDT2”，“ZDT3”，“ZDT6”，“Kursawe”，“Schaffer”，“Poloni”，“Viennet2”，“Viennet3”中对多意见粒子群进行测试，限度如下：
其中绿色的线代表真确的Pareto前沿面，玄色的圆圈默示寻优获取的Pareto值，红色的圈默示其他粒子。
ZDT1
猎U者图片
ZDT2
图片
ZDT3
图片
ZDT6
图片
Kursawe
图片
Schaffer
图片
Viennet2
图片
Viennet3
图片
可以看到，在这几个经典函数中的测试，多意见粒子群的限度还诅咒常可以的，但也有可改造的空间。
接下来告成上代码！
当先是主函数：
clear all; clc;% Multi-objective function% MultiObjFnc = 'Schaffer';% MultiObjFnc = 'Kursawe';% MultiObjFnc = 'Poloni';% MultiObjFnc = 'Viennet2';% MultiObjFnc = 'Viennet3';% MultiObjFnc = 'ZDT1';% MultiObjFnc = 'ZDT2';% MultiObjFnc = 'ZDT3';% MultiObjFnc = 'ZDT6';switch MultiObjFnc case 'Schaffer' % Schaffer MultiObj.fun = @(x) [x(:).^2， (x(:)-2).^2]; MultiObj.nVar = 1; MultiObj.var_min = -5; MultiObj.var_max = 5; load('Schaffer.mat'); MultiObj.truePF = PF; case 'Kursawe' % Kursawe MultiObj.fun = @(x) [-10.*(exp(-0.2.*sqrt(x(:，1).^2+x(:，2).^2)) + exp(-0.2.*sqrt(x(:，2).^2+x(:，3).^2)))， ... sum(abs(x).^0.8 + 5.*sin(x.^3)，2)]; MultiObj.nVar = 3; MultiObj.var_min = -5.*ones(1，MultiObj.nVar); MultiObj.var_max = 5.*ones(1，MultiObj.nVar); load('Kursawe.mat'); MultiObj.truePF = PF; case 'Poloni' % Poloni's two-objective A1 = 0.5*sin(1)-2*cos(1)+sin(2)-1.5*cos(2); A2 = 1.5*sin(1)-cos(1)+2*sin(2)-0.5*cos(2); B1 = @(x，y) 0.5.*sin(x)-2.*cos(x)+sin(y)-1.5.*cos(y); B2 = @(x，y) 1.5.*sin(x)-cos(x)+2.*sin(y)-0.5.*cos(y); f1 = @(x，y) 1+(A1-B1(x，y)).^2+(A2-B2(x，y)).^2; f2 = @(x，y) (x+3).^2+(y+1).^2; MultiObj.fun = @(x) [f1(x(:，1)，x(:，2))， f2(x(:，1)，x(:，2))]; MultiObj.nVar = 2; MultiObj.var_min = -pi.*ones(1，MultiObj.nVar); MultiObj.var_max = pi.*ones(1，MultiObj.nVar); case 'Viennet2' % Viennet2 f1 = @(x，y) 0.5.*(x-2).^2+(1/13).*(y+1).^2+3; f2 = @(x，y) (1/36).*(x+y-3).^2+(1/8).*(-x+y+2).^2-17; f3 = @(x，y) (1/175).*(x+2.*y-1).^2+(1/17).*(2.*y-x).^2-13; MultiObj.fun = @(x) [f1(x(:，1)，x(:，2))， f2(x(:，1)，x(:，2))， f3(x(:，1)，x(:，2))]; MultiObj.nVar = 2; MultiObj.var_min = [-4， -4]; MultiObj.var_max = [4， 4]; load('Viennet2.mat'); MultiObj.truePF = PF; case 'Viennet3' % Viennet3 f1 = @(x，y) 0.5.*(x.^2+y.^2)+sin(x.^2+y.^2); f2 = @(x，y) (1/8).*(3.*x-2.*y+4).^2 + (1/27).*(x-y+1).^2 +15; f3 = @(x，y) (1./(x.^2+y.^2+1))-1.1.*exp(-(x.^2+y.^2)); MultiObj.fun = @(x) [f1(x(:，1)，x(:，2))， f2(x(:，1)，x(:，2))， f3(x(:，1)，x(:，2))]; MultiObj.nVar = 2; MultiObj.var_min = [-3， -10]; MultiObj.var_max = [10， 3]; load('Viennet3.mat'); MultiObj.truePF = PF; case 'ZDT1' % ZDT1 (convex) g = @(x) 1+9.*sum(x(:，2:end)，2)./(size(x，2)-1); MultiObj.fun = @(x) [x(:，1)， g(x).*(1-sqrt(x(:，1)./g(x)))]; MultiObj.nVar = 30; MultiObj.var_min = zeros(1，MultiObj.nVar); MultiObj.var_max = ones(1，MultiObj.nVar); load('ZDT1.mat'); MultiObj.truePF = PF; case 'ZDT2' % ZDT2 (non-convex) f = @(x) x(:，1); g = @(x) 1+9.*sum(x(:，2:end)，2)./(size(x，2)-1); h = @(x) 1-(f(x)./g(x)).^2; MultiObj.fun = @(x) [f(x)， g(x).*h(x)]; MultiObj.nVar = 30; MultiObj.var_min = zeros(1，MultiObj.nVar); MultiObj.var_max = ones(1，MultiObj.nVar); load('ZDT2.mat'); MultiObj.truePF = PF; case 'ZDT3' % ZDT3 (discrete) f = @(x) x(:，1); g = @(x) 1+(9/size(x，2)-1).*sum(x(:，2:end)，2); h = @(x) 1 - sqrt(f(x)./g(x)) - (f(x)./g(x)).*sin(10.*pi.*f(x)); MultiObj.fun = @(x) [f(x)， g(x).*h(x)]; MultiObj.nVar = 5; MultiObj.var_min = 0.*ones(1，MultiObj.nVar); MultiObj.var_max = 1.*ones(1，MultiObj.nVar); load('ZDT3.mat'); MultiObj.truePF = PF; case 'ZDT6' % ZDT6 (non-uniform) f = @(x) 1 - exp(-4.*x(:，1)).*sin(6.*pi.*x(:，1)); g = @(x) 1 + 9.*(sum(x(:，2:end)，2)./(size(x，2)-1)).^0.25; h = @(x) 1 - (f(x)./g(x)).^2; MultiObj.fun = @(x) [f(x)， g(x).*h(x)]; MultiObj.nVar = 10; MultiObj.var_min = 0.*ones(1，MultiObj.nVar); MultiObj.var_max = 1.*ones(1，MultiObj.nVar); load('ZDT6.mat'); MultiObj.truePF = PF;end% Parametersparams.Np = 200; % Population sizeparams.Nr = 200; % Repository sizeparams.maxgen = 100; % Maximum number of generationsparams.W = 0.4; % Inertia weightparams.C1 = 2; % Individual confidence factorparams.C2 = 2; % Swarm confidence factorparams.ngrid = 20; % Number of grids in each dimensionparams.maxvel = 5; % Maxmium vel in percentageparams.u_mut = 0.5; % Uniform mutation percentage% MOPSOREP = MOPSO(params，MultiObj);% Display infodisplay('Repository fitness values are stored in REP.pos_fit');display('Repository particles positions are store in REP.pos');然后是多意见粒子群函数代码：
function REP = MOPSO(params，MultiObj) % Parameters Np = params.Np; Nr = params.Nr; maxgen = params.maxgen; W = params.W; C1 = params.C1; C2 = params.C2; ngrid = params.ngrid; maxvel = params.maxvel; u_mut = params.u_mut; fun = MultiObj.fun; nVar = MultiObj.nVar; var_min = MultiObj.var_min(:); var_max = MultiObj.var_max(:); % Initialization POS = repmat((var_max-var_min)'，Np，1).*rand(Np，nVar) + repmat(var_min'，Np，1); VEL = zeros(Np，nVar); POS_fit = fun(POS); if size(POS，1) ~= size(POS_fit，1) warning(['The objective function is badly programmed. It is not returning' ... 'a value for each particle， please check it.']); end PBEST = POS; PBEST_fit= POS_fit; DOMINATED= checkDomination(POS_fit); REP.pos = POS(~DOMINATED，:); REP.pos_fit = POS_fit(~DOMINATED，:); REP = updateGrid(REP，ngrid); maxvel = (var_max-var_min).*maxvel./100; gen = 1; % Plotting and verbose if(size(POS_fit，2)==2) h_fig = figure(1); h_par = plot(POS_fit(:，1)，POS_fit(:，2)，'or'); hold on; h_rep = plot(REP.pos_fit(:，1)，REP.pos_fit(:，2)，'ok'); hold on; try set(gca，'xtick'，REP.hypercube_limits(:，1)'，'ytick'，REP.hypercube_limits(:，2)'); axis([min(REP.hypercube_limits(:，1)) max(REP.hypercube_limits(:，1)) ... min(REP.hypercube_limits(:，2)) max(REP.hypercube_limits(:，2))]); grid on; xlabel('f1'); ylabel('f2'); end drawnow; end if(size(POS_fit，2)==3) h_fig = figure(1); h_par = plot3(POS_fit(:，1)，POS_fit(:，2)，POS_fit(:，3)，'or'); hold on; h_rep = plot3(REP.pos_fit(:，1)，REP.pos_fit(:，2)，REP.pos_fit(:，3)，'ok'); hold on; try set(gca，'xtick'，REP.hypercube_limits(:，1)'，'ytick'，REP.hypercube_limits(:，2)'，'ztick'，REP.hypercube_limits(:，3)'); axis([min(REP.hypercube_limits(:，1)) max(REP.hypercube_limits(:，1)) ... min(REP.hypercube_limits(:，2)) max(REP.hypercube_limits(:，2))]); end grid on; xlabel('f1'); ylabel('f2'); zlabel('f3'); drawnow; axis square; end display(['Generation #0 - Repository size: ' num2str(size(REP.pos，1))]); % Main MPSO loop stopCondition = false; while ~stopCondition % Select leader h = selectLeader(REP); % Update speeds and positions VEL = W.*VEL + C1*rand(Np，nVar).*(PBEST-POS) ... + C2*rand(Np，nVar).*(repmat(REP.pos(h，:)，Np，1)-POS); POS = POS + VEL; % Perform mutation POS = mutation(POS，gen，maxgen，Np，var_max，var_min，nVar，u_mut); % Check boundaries [POS，VEL] = checkBoundaries(POS，VEL，maxvel，var_max，var_min); % Evaluate the population POS_fit = fun(POS); % Update the repository REP = updateRepository(REP，POS，POS_fit，ngrid); if(size(REP.pos，1)>Nr) REP = deleteFromRepository(REP，size(REP.pos，1)-Nr，ngrid); end % Update the best positions found so far for each particle pos_best = dominates(POS_fit， PBEST_fit); best_pos = ~dominates(PBEST_fit， POS_fit); best_pos(rand(Np，1)>=0.5) = 0; if(sum(pos_best)>1) PBEST_fit(pos_best，:) = POS_fit(pos_best，:); PBEST(pos_best，:) = POS(pos_best，:); end if(sum(best_pos)>1) PBEST_fit(best_pos，:) = POS_fit(best_pos，:); PBEST(best_pos，:) = POS(best_pos，:); end % Plotting and verbose if(size(POS_fit，2)==2) figure(h_fig); delete(h_par); delete(h_rep); h_par = plot(POS_fit(:，1)，POS_fit(:，2)，'or'); hold on; h_rep = plot(REP.pos_fit(:，1)，REP.pos_fit(:，2)，'ok'); hold on; try set(gca，'xtick'，REP.hypercube_limits(:，1)'，'ytick'，REP.hypercube_limits(:，2)'); axis([min(REP.hypercube_limits(:，1)) max(REP.hypercube_limits(:，1)) ... min(REP.hypercube_limits(:，2)) max(REP.hypercube_limits(:，2))]); end if(isfield(MultiObj，'truePF')) try delete(h_pf); end h_pf = plot(MultiObj.truePF(:，1)，MultiObj.truePF(:，2)，'.'，'color'，'g'); hold on; end grid on; xlabel('f1'); ylabel('f2'); drawnow; axis square; end if(size(POS_fit，2)==3) figure(h_fig); delete(h_par); delete(h_rep); h_par = plot3(POS_fit(:，1)，POS_fit(:，2)，POS_fit(:，3)，'or'); hold on; h_rep = plot3(REP.pos_fit(:，1)，REP.pos_fit(:，2)，REP.pos_fit(:，3)，'ok'); hold on; try set(gca，'xtick'，REP.hypercube_limits(:，1)'，'ytick'，REP.hypercube_limits(:，2)'，'ztick'，REP.hypercube_limits(:，3)'); axis([min(REP.hypercube_limits(:，1)) max(REP.hypercube_limits(:，1)) ... min(REP.hypercube_limits(:，2)) max(REP.hypercube_limits(:，2)) ... min(REP.hypercube_limits(:，3)) max(REP.hypercube_limits(:，3))]); end if(isfield(MultiObj，'truePF')) try delete(h_pf); end h_pf = plot3(MultiObj.truePF(:，1)，MultiObj.truePF(:，2)，MultiObj.truePF(:，3)，'.'，'color'，'g'); hold on; end grid on; xlabel('f1'); ylabel('f2'); zlabel('f3'); drawnow; axis square; end display(['Generation #' num2str(gen) ' - Repository size: ' num2str(size(REP.pos，1))]); % Update generation and check for termination gen = gen + 1; if(gen>maxgen)， stopCondition = true; end end hold off;end% Function that updates the repository given a new population and its% fitnessfunction REP = updateRepository(REP，POS，POS_fit，ngrid) % Domination between particles DOMINATED = checkDomination(POS_fit); REP.pos = [REP.pos; POS(~DOMINATED，:)]; REP.pos_fit= [REP.pos_fit; POS_fit(~DOMINATED，:)]; % Domination between nondominated particles and the last repository DOMINATED = checkDomination(REP.pos_fit); REP.pos_fit= REP.pos_fit(~DOMINATED，:); REP.pos = REP.pos(~DOMINATED，:); % Updating the grid REP = updateGrid(REP，ngrid);end% Function that corrects the positions and velocities of the particles that% exceed the boundariesfunction [POS，VEL] = checkBoundaries(POS，VEL，maxvel，var_max，var_min) % Useful matrices Np = size(POS，1); MAXLIM = repmat(var_max(:)'，Np，1); MINLIM = repmat(var_min(:)'，Np，1); MAXVEL = repmat(maxvel(:)'，Np，1); MINVEL = repmat(-maxvel(:)'，Np，1); % Correct positions and velocities VEL(VEL>MAXVEL) = MAXVEL(VEL>MAXVEL); VEL(VEL<MINVEL) = MINVEL(VEL<MINVEL); VEL(POS>MAXLIM) = (-1).*VEL(POS>MAXLIM); POS(POS>MAXLIM) = MAXLIM(POS>MAXLIM); VEL(POS<MINLIM) = (-1).*VEL(POS<MINLIM); POS(POS<MINLIM) = MINLIM(POS<MINLIM);end% Function for checking the domination between the population. It% returns a vector that indicates if each particle is dominated (1) or notfunction dom_vector = checkDomination(fitness) Np = size(fitness，1); dom_vector = zeros(Np，1); all_perm = nchoosek(1:Np，2); % Possible permutations all_perm = [all_perm; [all_perm(:，2) all_perm(:，1)]]; d = dominates(fitness(all_perm(:，1)，:)，fitness(all_perm(:，2)，:)); dominated_particles = unique(all_perm(d==1，2)); dom_vector(dominated_particles) = 1;end% Function that returns 1 if x dominates y and 0 otherwisefunction d = dominates(x，y) d = all(x<=y，2) & any(x<y，2);end% Function that updates the hypercube grid， the hypercube where belongs% each particle and its quality based on the number of particles inside itfunction REP = updateGrid(REP，ngrid) % Computing the limits of each hypercube ndim = size(REP.pos_fit，2); REP.hypercube_limits = zeros(ngrid+1，ndim); for dim = 1:1:ndim REP.hypercube_limits(:，dim) = linspace(min(REP.pos_fit(:，dim))，max(REP.pos_fit(:，dim))，ngrid+1)'; end % Computing where belongs each particle npar = size(REP.pos_fit，1); REP.grid_idx = zeros(npar，1); REP.grid_subidx = zeros(npar，ndim); for n = 1:1:npar idnames = []; for d = 1:1:ndim REP.grid_subidx(n，d) = find(REP.pos_fit(n，d)<=REP.hypercube_limits(:，d)'，1，'first')-1; if(REP.grid_subidx(n，d)==0)， REP.grid_subidx(n，d) = 1; end idnames = [idnames '，' num2str(REP.grid_subidx(n，d))]; end REP.grid_idx(n) = eval(['sub2ind(ngrid.*ones(1，ndim)' idnames ');']); end % Quality based on the number of particles in each hypercube REP.quality = zeros(ngrid，2); ids = unique(REP.grid_idx); for i = 1:length(ids) REP.quality(i，1) = ids(i); % First， the hypercube's identifier REP.quality(i，2) = 10/sum(REP.grid_idx==ids(i)); % Next， its quality endend% Function that selects the leader performing a roulette wheel selection% based on the quality of each hypercubefunction selected = selectLeader(REP) % Roulette wheel prob = cumsum(REP.quality(:，2)); % Cumulated probs sel_hyp = REP.quality(find(rand(1，1)*max(prob)<=prob，1，'first')，1); % Selected hypercube % Select the index leader as a random selection inside that hypercube idx = 1:1:length(REP.grid_idx); selected = idx(REP.grid_idx==sel_hyp); selected = selected(randi(length(selected)));end% Function that deletes an excess of particles inside the repository using% crowding distancesfunction REP = deleteFromRepository(REP，n_extra，ngrid) % Compute the crowding distances crowding = zeros(size(REP.pos，1)，1); for m = 1:1:size(REP.pos_fit，2) [m_fit，idx] = sort(REP.pos_fit(:，m)，'ascend'); m_up = [m_fit(2:end); Inf]; m_down = [Inf; m_fit(1:end-1)]; distance = (m_up-m_down)./(max(m_fit)-min(m_fit)); [~，idx] = sort(idx，'ascend'); crowding = crowding + distance(idx); end crowding(isnan(crowding)) = Inf; % Delete the extra particles with the smallest crowding distances [~，del_idx] = sort(crowding，'ascend'); del_idx = del_idx(1:n_extra); REP.pos(del_idx，:) = []; REP.pos_fit(del_idx，:) = []; REP = updateGrid(REP，ngrid); end% Function that performs the mutation of the particles depending on the% current generationfunction POS = mutation(POS，gen，maxgen，Np，var_max，var_min，nVar，u_mut) % Sub-divide the swarm in three parts [2] fract = Np/3 - floor(Np/3); if(fract<0.5)， sub_sizes =[ceil(Np/3) round(Np/3) round(Np/3)]; else sub_sizes =[round(Np/3) round(Np/3) floor(Np/3)]; end cum_sizes = cumsum(sub_sizes); % First part: no mutation % Second part: uniform mutation nmut = round(u_mut*sub_sizes(2)); if(nmut>0) idx = cum_sizes(1) + randperm(sub_sizes(2)，nmut); POS(idx，:) = repmat((var_max-var_min)'，nmut，1).*rand(nmut，nVar) + repmat(var_min'，nmut，1); end % Third part: non-uniform mutation per_mut = (1-gen/maxgen)^(5*nVar); % Percentage of mutation nmut = round(per_mut*sub_sizes(3)); if(nmut>0) idx = cum_sizes(2) + randperm(sub_sizes(3)，nmut); POS(idx，:) = repmat((var_max-var_min)'，nmut，1).*rand(nmut，nVar) + repmat(var_min'，nmut，1); endend代码中零落各个函数真确的Pareto前沿数据草榴社区邀请码，数据没法告成复制在文中。因此作家打包了一下。
本站仅提供存储处事，悉数本体均由用户发布，如发现存害或侵权本体，请点击举报。