%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都是怎么得到的
最新发布