47. Permutations II 【M】【61】

本文介绍了一种算法,用于找出包含重复数字的所有唯一排列。通过逐步添加数字并处理重复情况,确保结果中没有重复的排列。文章详细展示了算法的实现过程,并提供了Python代码示例。

Given a collection of numbers that might contain duplicates, return all possible unique permutations.

For example,
[1,1,2] have the following unique permutations:
[1,1,2][1,2,1], and [2,1,1].


Subscribe to see which companies asked this question



the basic idea is, to permute n numbers, we can add the nth number into the resulting List<List<Integer>> from the n-1 numbers, in every possible position. 

For example, if the input num[] is {1,2,3}: First, add 1 into the initial List<List<Integer>> (let's call it "answer"). 

Then, 2 can be added in front or after 1. So we have to copy the List in answer (it's just {1}), add 2 in position 0 of {1}, then copy the original {1} again, and add 2 in position 1. Now we have an answer of {{2,1},{1,2}}. There are 2 lists in the current answer.

Then we have to add 3. first copy {2,1} and {1,2}, add 3 in position 0; then copy {2,1} and {1,2}, and add 3 into position 1, then do the same thing for position 3. Finally we have 2*3=6 lists in answer, which is what we want.

就是说,用 [1,1,3] 举例子

每次把一个数字放进去,每个位置都放一遍。那处理重复怎么办呢?遇见重复就break就行了。。

[1]
[[1]]

[3, 1]
[1, 3]

[[3, 1], [1, 3]]
[3, 3, 1]
[[3, 1], [1, 3]]
[3, 1, 3]
[1, 3, 3]







class Solution(object):
    def permuteUnique(self, nums):
        ans = [[]]
        for n in nums:
            new_ans = []
            for l in ans:
                #print ans
                for i in xrange(len(l)+1):
                    #print l[:i]+[n]+l[i:]
                    new_ans.append(l[:i]+[n]+l[i:])
                    if i<len(l) and l[i]==n: #i < len(l) 保证数组不越界
                        break
                    #handles duplication
            ans = new_ans
        return ans
            
            
        '''
        p = itertools.permutations(nums)
        res = set()
        for i in p:
            res.add(i)
        ress = []
        for i in res:
            ress += [list(i)]
        return ress
        
        '''
        


%MOPSO:63、311;example:42 % ----------------------------------------------------------------------- % clear clc mpc=loadcase('case33bw'); n_Bus=mpc.bus(end,1); Load=[ 0.8, 0.805, 0.81, 0.818, 0.83, 0.91, 0.95, 1, 1.05, 1.08, 1.03, 1.01, 0.95, 0.9, 0.905, 0.91, 0.93, 1.02, 1.1, 1.18, 1.2, 1.08, 0.9, 0.84]/ 1.4; PV=[0,0,0,0,0,1,2.5,4,5,5.5,5.8,5.7,5.5,5.3,5.1,5,3.8,2.5,1.2,0,0,0,0,0]/5.8; Wind=[0.95,0.9,0.92,0.93,0.95,0.97,1,0.95,0.85,0.65,0.55,0.45,0.43,0.4,0.37,0.4,0.44,0.48,0.52,0.57,0.62,0.68,0.72,0.8]; % 多目标函数 MultiObjFnc = 'Kursawe'; switch MultiObjFnc case 'Kursawe' % Kursawe MultiObj.fun = @(x) funmm(x); MultiObj.nVar = 3+3+48; MultiObj.var_min(1:3) = 2; MultiObj.var_max(1:3) = 33; MultiObj.var_min(4:6) = 0; MultiObj.var_max(4:6) = 5000; MultiObj.var_min(7:54) = 0; MultiObj.var_max(7:54) = 300; case 'Viennet2' % Viennet2 MultiObj.fun = @(x) fun3(x); MultiObj.nVar = 31; MultiObj.var_min(1:4)=2;%位置 MultiObj.var_max(1:4)=n_Bus; MultiObj.var_min(5:7)=0; MultiObj.var_max(5:7)=30000; MultiObj.var_min(8:31)=0; MultiObj.var_max(8:31)=300; end %% MOPSO参数 params.Np = 200; % 种群规模 params.Nr = 200; % 存储库大小 params.maxgen = 2; % 最大迭代次数 params.W = 0.4; % 惯性权重 params.C1 = 1.4; % Individual confidence factor params.C2 = 1.4; % Swarm confidence factor params.ngrid = 10; % Number of grids in each dimension params.maxvel = 10; % Maxmium vel in percentage params.u_mut = 0.8; % Uniform mutation percentage % MOPSO REP = MOPSO(params,MultiObj); size(REP) REP.pos(:,1:6)=floor(REP.pos(:,1:6)); %结果 disp('目标函数值:REP.pos_fit'); disp('DGs接入位置、容量:REP.pos'); figure scatter(REP.pos_fit(:,1),REP.pos_fit(:,2),'r','g') xlabel('f1');ylabel('f2'); figure plot(Load) hold on plot(PV) hold on plot(Wind) % % for i=1:1:5 % figure % stairs(REP.pos(i,9:32),'ro-') % figure % stairs(REP.pos(i,33:56),'ro-') % end Figure(REP,12) 第二部分 function Figure=Figure(x,K) %x分别表示了DG的数量 位置 还有容i昂,为一个矩阵,k是x的总个数; define_constants; Load=[ 0.8, 0.805, 0.81, 0.818, 0.83, 0.91, 0.95, 1, 1.05, 1.08, 1.03, 1.01, 0.95, 0.9, 0.905, 0.91, 0.93, 1.02, 1.1, 1.18, 1.2, 1.08, 0.9, 0.84]/ 1.4; %这边查资料看负荷曲线,根据论文《考虑电动汽车时空特性分布》只得到了P的值,而且 Cap=[150 400]; %% 接入前后对比 for ii= 1: 24 %每一维进行潮流计算,计算网损等 mpc=loadcase('case33bw_60'); mpc.bus(:,PD)=mpc.bus(:,PD)*Load(ii);%负荷 results=runpf(mpc,mpoption('verbose',0,'pf.alg','NR','out.all',0)); Bus_Voltage=results.bus(:,8); withoutminV(ii)=min(Bus_Voltage); withoutV(:,ii)=Bus_Voltage; withoutPLoss(ii)=sum(real(get_losses(results))); withoutQLoss(ii)=sum(imag(get_losses(results))); withoutresults(:,ii)=real(get_losses(results)); %脆弱性 syms P Q; P=results.branch(:,1);Q=results.branch(:,1); PP=results.bus(P,3);QQ=results.bus(Q,4); V=results.bus(P,8);phase_angle=results.bus(P,9); A=1+(results.branch(:,3)+1i.*results.branch(:,4)).*1i.*results.branch(:,4)/2; AA=angle(A);A=abs(A); B=results.branch(:,3)+1i.*results.branch(:,4); BB=angle(B);B=abs(B); LCPI=4*(A.*cos(AA)).*(PP.*B.*cos(BB)+QQ.*B.*sin(BB))./V./cos(phase_angle)/1000; withoutLCPI(ii)=max(LCPI); %电压稳定性 L=0;withoutL=0; X=results.branch(:,3);R=results.branch(:,4); b=results.branch(:,2); Pb=results.bus(b,3);Qb=results.bus(b,4); L=4*(X.*Pb-R.*Qb).^2+X.*Qb+R.*Pb; withoutL=L+withoutL; end withoutL=withoutL/24; for ii= 1: 24 %每一维进行潮流计算,计算网损等 mpc=loadcase('case33bw'); mpc.bus(:,PD)=mpc.bus(:,PD)*Load(ii);%负荷 %分布式电源 mpc.bus(x.pos(K,1),PD)=mpc.bus(x.pos(K,1),PD)-x.pos(K,6)*Wind(ii)*dgPf/1000; %分布式电源 mpc.bus(x.pos(K,1),QD)=mpc.bus(x.pos(K,1),QD)-x.pos(K,6)*Wind(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x.pos(K,2),PD)=mpc.bus(x.pos(K,2),PD)-x.pos(K,7)*PV(ii)*dgPf/1000; mpc.bus(x.pos(K,2),QD)=mpc.bus(x.pos(K,2),QD)-x.pos(K,7)*PV(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x.pos(K,3),PD)=mpc.bus(x.pos(K,3),PD)-x.pos(K,8)*PV(ii)*dgPf/1000; mpc.bus(x.pos(K,3),QD)=mpc.bus(x.pos(K,3),QD)-x.pos(K,8)*PV(ii)*(sqrt(1-dgPf*dgPf))/1000; %无功补偿 mpc.bus(7,QD)=mpc.bus(7,QD)-x.pos(K,6+ii)/1000; mpc.bus(25,QD)=mpc.bus(25,QD)-x.pos(K,30+ii)/1000; results=runpf(mpc,mpoption('verbose',0,'pf.alg','NR','out.all',0)); Bus_Voltage=results.bus(:,8); withminV(ii)=min(Bus_Voltage); withV(:,ii)=Bus_Voltage; withPLoss(ii)=sum(real(get_losses(results))); withQLoss(ii)=sum(imag(get_losses(results))); withresults(:,ii)=real(get_losses(results)); %脆弱性 syms P Q; P=results.branch(:,1);Q=results.branch(:,1); PP=results.bus(P,3);QQ=results.bus(Q,4); V=results.bus(P,8);phase_angle=results.bus(P,9); A=1+(results.branch(:,3)+1i.*results.branch(:,4)).*1i.*results.branch(:,4)/2; AA=angle(A);A=abs(A); B=results.branch(:,3)+1i.*results.branch(:,4); BB=angle(B);B=abs(B); LCPI=4*(A.*cos(AA)).*(PP.*B.*cos(BB)+QQ.*B.*sin(BB))./V./cos(phase_angle)/10000; withLCPI(ii)=max(LCPI); %电压稳定性 L=0;withL=0; X=results.branch(:,3);R=results.branch(:,4); b=results.branch(:,2); Pb=results.bus(b,3);Qb=results.bus(b,4); L=4*(X.*Pb-R.*Qb).^2+X.*Qb+R.*Pb; withL=L+withL; end withL=withL/24/10; withoutV24=sum(withoutV,2)/24;%电压变化 withV24=sum(withV,2)/24; withoutPLoss24=sum(withoutPLoss);%总有功损耗 withPLoss24=sum(withPLoss); withoutQLoss24=sum(withoutQLoss);%总无功损耗 withQLoss24=sum(withQLoss); withoutLoss24=[withoutPLoss24;withPLoss24]; withLoss24=[withoutQLoss24;withQLoss24]; withoutresults=sum(withoutresults,2);%各支路损耗 withresults=sum(withresults,2); %% 作图 figure plot(Load,'LineWidth',1) hold on plot(PV,'LineWidth',1) hold on plot(Wind,'LineWidth',1) set(gca,'XTick',(1:2:24)) legend('负荷','光伏') figure [XX,YY] =meshgrid(1:24,1:33 ); mesh(XX,YY,withoutV); xlabel('时刻(h)'); ylabel('节点序号'); zlabel('电压幅值(pu)'); title('24小时节点电压图'); figure [XX,YY] =meshgrid(1:24,1:33 ); mesh(XX,YY,withV); xlabel('时刻(h)'); ylabel('节点序号'); zlabel('电压幅值(pu)'); title('24小时节点电压图'); figure; plot(withoutV24,'red'); hold on; plot(withV24,'blue'); hold off; title('系统电压配置文件'); xlabel('节点') ; ylabel('电压幅值[p.u]') ; legend('DGs接入前电压','DGs接入后电压'); figure; losses=[ withoutLoss24,withLoss24]; legends=categorical({'DGs接入前网损', 'DGs接入后网损'}); legends = reordercats(legends,{'DGs接入前网损', 'DGs接入后网损'}); bar(legends,losses); title('系统总有功、无功损耗'); legend('有功损耗','无功损耗'); ylabel('kW/kVar') ; figure losses=[withoutresults,withresults]; bar(losses); title('各线路损耗 (KW)'); set(gca,'XTick',(1:2:38)) legend('without DGs','with DGs'); xlabel('支路序号') ; ylabel('kW') ; figure plot(withPLoss,'LineWidth',1) hold on plot(withoutPLoss,'LineWidth',1) legend('接入补偿设备','未接入补偿设备'); title('网损对比') figure plot(withminV,'LineWidth',1) hold on plot(withoutminV,'LineWidth',1) legend('接入补偿设备','未接入补偿设备'); title('最小电压') figure plot(withoutLCPI,'LineWidth',1) hold on plot(withLCPI,'LineWidth',1) legend('未接入补偿设备','接入补偿设备'); xlabel('时刻'); ylabel('LCPI'); title('电压脆弱性') figure plot(withoutL,'LineWidth',1) hold on plot(withL,'LineWidth',1) legend('未接入补偿设备','接入补偿设备'); xlabel('时刻'); ylabel('L'); title('电压稳定性') figure stairs(x.pos(K,7:30),'ro-') title('无功补偿装置出力') xlabel('时刻'); ylabel('kW'); figure stairs(x.pos(K,31:54),'ro-') xlabel('时刻'); ylabel('kW'); title('无功补偿装置出力') end 第三 %% funm %function f=funmm(x) x(:,1:5)=floor(x(:,1:5)); n=size(x,1);%lizishu dgPf=0.9; define_constants; mpc=loadcase('case33bw'); nBus=mpc.bus(end,1); Load=[ 0.8, 0.805, 0.81, 0.818, 0.83, 0.91, 0.95, 1, 1.05, 1.08, 1.03, 1.01, 0.95, 0.9, 0.905, 0.91, 0.93, 1.02, 1.1, 1.18, 1.2, 1.08, 0.9, 0.84]/ 1.4; PV=[0,0,0,0,0,1,2.5,4,5,5.5,5.8,5.7,5.5,5.3,5.1,5,3.8,2.5,1.2,0,0,0,0,0]/5.8; Wind=[0.95,0.9,.92,0.93,0.95,0.97,1,0.95,0.85,0.65,0.55,0.45,0.43,0.4,0.37,0.4,0.44,0.48,0.52,0.57,0.62,0.68,0.72,0.8]; for i=1:n f11=0;f22=0; for ii= 1: 24 %每一维进行潮流计算,计算网损等 mpc=loadcase('case33bw'); mpc.bus(:,PD)=mpc.bus(:,PD)+mpc.bus(:,PD)*Load(ii);%负荷 %分布式电源 mpc.bus(x(i,1),PD)=mpc.bus(x(i,1),PD)-x(i,4)*Wind(ii)*dgPf/1000; %dgPf功率因数 mpc.bus(x(i,1),QD)=mpc.bus(x(i,1),QD)-x(i,4)*Wind(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x(i,2),PD)=mpc.bus(x(i,2),PD)-x(i,5)*PV(ii)*dgPf/1000; %dgPf功率因数 mpc.bus(x(i,2),QD)=mpc.bus(x(i,2),QD)-x(i,5)*PV(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x(i,3),PD)=mpc.bus(x(i,3),PD)-x(i,6)*PV(ii)*dgPf/1000; %dgPf功率因数 mpc.bus(x(i,3),QD)=mpc.bus(x(i,3),QD)-x(i,6)*PV(ii)*(sqrt(1-dgPf*dgPf))/1000; %无功补偿装置 mpc.bus(7,4)=mpc.bus(7,4)-x(i,6+ii)/1000; mpc.bus(25,4)=mpc.bus(25,4)-x(i,30+ii)/1000; results=runpf(mpc,mpoption('verbose',0,'pf.alg','NR','out.all',0)); Bus_Voltage=results.bus(:,8); f11=f11+sum(abs(Bus_Voltage-ones(33,1)*1)./(ones(33,1)*1)); f22=f22+sum(real(get_losses(results))); end f1(i,1)=f11; f2(i,1)=f22; end f=[f1,f2]; end 第四 %储能运行策略 function f=fun3(x) x(:,1:4)=floor(x(:,1:4)); n=size(x,1);%种群数 dgPf=0.9; define_constants; mpc=loadcase('case33bw'); nBus=mpc.bus(end,1); Load=[ 0.8, 0.805, 0.81, 0.818, 0.83, 0.91, 0.95, 1, 1.05, 1.08, 1.03, 1.01, 0.95, 0.9, 0.905, 0.91, 0.93, 1.02, 1.1, 1.18, 1.2, 1.08, 0.9, 0.84]/ 1.4; PV=[0,0,0,0,0,1,2.5,4,5,5.5,5.8,5.7,5.5,5.3,5.1,5,3.8,2.5,1.2,0,0,0,0,0]/5.8; Wind=[0.95,0.9,.92,0.93,0.95,0.97,1,0.95,0.85,0.65,0.55,0.45,0.43,0.4,0.37,0.4,0.44,0.48,0.52,0.57,0.62,0.68,0.72,0.8]; for i=1:n VDI=0;PLI=0; for ii= 1: 24 %每一维进行潮流计算,计算网损等 mpc=loadcase('case33bw'); mpc.bus(:,PD)=mpc.bus(:,PD)*Load(ii);%负荷 mpc.bus(x(i,1),PD)=mpc.bus(x(i,1),PD)-x(i,5)*PV(ii)*dgPf/1000; %dgPf功率因数 mpc.bus(x(i,1),QD)=mpc.bus(x(i,1),QD)-x(i,5)*PV(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x(i,2),PD)=mpc.bus(x(i,2),PD)-x(i,6)*PV(ii)*dgPf/1000; %dgPf功率因数 mpc.bus(x(i,2),QD)=mpc.bus(x(i,2),QD)-x(i,6)*PV(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x(i,3),PD)=mpc.bus(x(i,3),PD)-x(i,7)*Wind(ii)*dgPf/1000; %dgPf功率因数 mpc.bus(x(i,3),QD)=mpc.bus(x(i,3),QD)-x(i,7)*Wind(ii)*(sqrt(1-dgPf*dgPf))/1000; mpc.bus(x(i,4),QD)=mpc.bus(x(i,4),QD)-x(i,7+ii)/1000; results=runpf(mpc,mpoption('verbose',0,'pf.alg','NR','out.all',0)); Bus_Voltage=results.bus(:,8); VDI=VDI+sum(abs(Bus_Voltage-ones(33,1)*1)./(ones(33,1)*1)); PLI=PLI+sum(real(get_losses(results))); end f1(i,1)=VDI;%目标函数1 f1(i,2)=PLI; end %% 目标函数 f=f1; end 第五 %% ----------------------------------------------------------------------- % % Function MOPSO performs a Multi-Objective Particle Swarm Optimization % % over continous functions. % % % % Input parameters: % % - params: Struct that contains the customized parameters. % % * params.Np: Number of particles. % % * params.Nr: Repository size (in particles). % % * params.maxgen: Maximum number of generations. % % * params.W: Inertia coefficient. % % * params.C1: Personal confidence factor. % % * params.C2: Swarm confidence factor. % % * params.ngrid: Number of hypercubes in each dimension. % % * params.maxvel: Maximum velocity (search space percentage)% % * params.u_mut: Uniform mutation percentage. % % - MultiObj: Struct that contains the parameters relative to the % % optimization functions. % % * MultiObj.fun: Anonymous multi-obj function to minimize. % % * MultiObj.nVar: Number of variables. % % * MultiObj.var_min: Vector that indicates the minimum values % % of the search space in each dimension. % % * MultiObj.var_max: Same than 'var_min' with the maxima. % % ----------------------------------------------------------------------- % % For an example of use, run 'example.m'. % % ----------------------------------------------------------------------- % % Author: Victor Martinez Cagigal % % Date: 17/03/2017 % % E-mail: vicmarcag (at) gmail (dot) com % % Version: 1.1 % % Log: % % - 1.0: Initial version without mutation [1] (15/03/2017). % % - 1.1: Crowding and mutation are implemented [2]. % % ----------------------------------------------------------------------- % % References: % % [1]Coello, C. A. C., Pulido, G. T., & Lechuga, M. S. (2004). Handling% % multiple objectives with particle swarm optimization. IEEE Tran- % % sactions on evolutionary computation, 8(3), 256-279. % % % % [2]Sierra, M. R., & Coello, C. A. C. (2005, March). Improving PSO- % % based multi-objective optimization using crowding, mutation and ?-% % dominance. In International Conference on Evolutionary Multi-Crite% % rion Optimization (pp. 505-519). Springer Berlin Heidelberg. % %% ----------------------------------------------------------------------- % 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); POS(:,1:5)=floor(POS(:,1:5)); 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',0.8.*ones(1,3)); 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',0.8.*ones(1,3)); 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 % fitness function 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 boundaries function [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 not function 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 otherwise function 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 it function 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 end end % Function that selects the leader performing a roulette wheel selection % based on the quality of each hypercube function 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 distances function 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); REP.pos(:,1:5)=floor(REP.pos(:,1:5)); end % Function that performs the mutation of the particles depending on the % current generation function 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); end end 所有代码都在此,详细分析REP.pos与REP.pos_fit都是怎么得到的
最新发布
07-18
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值