基于改进粒子群算法的微电网多目标优化调度
参考文献《基于改进粒子群算法的微电网多目标优化调度》
摘要：微电网优化调度作为智能电网优化的重要组成部分，对降低能耗、环境污染具有重要意义。微电网的发展目标既要满足电力供应的基本需求，又要提高经济效益和环境保护。对此，提出了一种综合考虑微电网系统运行成本和环境保护成本的并网模式下微电网多目标优化调度模型。同时采用改进的粒子群算法对优化模型进行求解。仿真结果表明，该模型可以有效降低用户的用电成本和环境污染，促进微电网的优化运行，并验证了改进的粒子群算法的优越性能。
关键词：微电网；多目标；改进粒子群算法；优化调度
1 引言
多目标粒子群算法（MOPSO）是一种用于解决多目标优化问题的进化算法。它基于粒子群算法（PSO），通过模拟鸟群中鸟类的群体行为来寻找最优解。MOPSO算法的目标是找到一组解决方案，这些解决方案能够在多个目标函数下达到最优值，而不是单一目标函数的最优解。在实际应用中，许多问题都涉及到多个相互矛盾的目标，因此MOPSO算法具有广泛的应用前景。
然而，传统的MOPSO算法在处理多目标优化问题时存在一些问题，例如收敛速度慢、局部最优解问题以及对多样性的维护不足。因此，对MOPSO算法进行改进是非常必要的。首先，可以引入多目标优化算法的权衡因子，用于平衡不同目标之间的重要性。这样可以确保算法更好地探索和利用多目标空间，避免过度偏向某个目标而忽视其他目标。其次，设计更有效的粒子更新策略是改进MOPSO算法的关键。通过优化粒子的更新策略，可以提高算法的收敛速度和全局搜索能力，从而更快地找到最优解。另外，引入多样性维持机制也是一种重要的改进方法。多样性维持机制可以避免算法陷入局部最优解，同时保持种群的多样性，从而更好地探索整个解空间。
此外，结合进化算法或遗传算法也是一种改进MOPSO算法的有效途径。通过结合不同的优化算法，可以增加算法的搜索能力和收敛性，从而更好地解决多目标优化问题。另外，引入自适应参数调整机制也是一种重要的改进方法。通过根据问题的特性动态调整算法参数，可以提高算法的适应性和鲁棒性，从而更好地适应不同的多目标优化问题。最后，结合机器学习技术和多种启发式搜索策略，以增加算法的多样性和全局搜索能力，也是改进MOPSO算法的有效途径。通过引入更多的智能技术，可以提高算法的预测能力和决策能力，从而更好地探索和利用多目标空间。
总之，改进多目标粒子群算法是一个非常重要的研究课题，通过引入新的思想和技术，可以提高算法的性能和适用性，从而更好地解决实际的多目标优化问题。希望未来能够有更多的研究者投入到这个领域，共同推动多目标优化算法的发展。
2 微电网多目标优化模型
(1) 微电网的运行成本
在并网模式下的目标是最小化微电网的运行成本。
(2) 微电网的环境保护成本。
3 算例
本文的微电网系统包含各种分布式电源，包括 PV、WT、DE、MT 和储能。微网中各个 DG的运行参数与成本见表 1,各个 DG 污染物排放系数及成本见表 2，储能参数见表 3。
4 程序运行结果
1）帕累托前沿
2）总成本最低时光伏消纳
3）总成本最低时风机消纳
4）总成本最低时各机组出力
5 程序
1）主函数
%基于改进粒子群算法的微电网多目标优化调度
clear;
clc;
close all;
%定义全局变量
global P_load; %电负荷
global WT;%风电
global PV;%光伏
global buy_price;
global sell_price;
%数据
data=[52 0 2 0.38 0.36
50 0 2 0.38 0.36
50 0 2 0.38 0.36
51 0 2 0.38 0.36
55 0 2 0.38 0.36
63 0 1 0.38 0.36
70 0 2 0.82 0.36
74 0 1 0.82 0.36
75 4 2 0.82 0.36
78 6 3 1.35 0.36
76 10 8 1.35 0.36
62 12 10 1.35 0.36
60 23 4 1.35 0.36
60 20 3 1.35 0.36
68 6 2 0.82 0.36
78 4 1 0.82 0.36
80 1 2 0.82 0.36
85 0 2 1.35 0.36
90 0 2 1.35 0.36
88 0 3 1.35 0.36
70 0 2.5 1.35 0.36
65 0 2.5 1.35 0.36
63 0 2 0.38 0.36
55 0 1 0.38 0.36
];
P_load=data(:,1);%电负荷（MW）
PV=data(:,2);%光伏（MW）
WT=data(:,3);%风电(MW)
buy_price=data(:,4);%向电网买电价
sell_price=data(:,5);%向电网卖电价
%蓄电池最大放电功率（正表示为电负荷供电，即放电）
BESSMax_dischar=30;
%蓄电池最大充电功率
BESSMax_char=-30;
%柴油机最大发电功率
DEMax=30;
%柴油机最小发电功率
DEMin=6;
%燃气轮机最大发电功率
MTMax=30;
%燃气轮机最小发电功率
MTMin=3;
%主网交互最大功率(正表示为电负荷供电)
GridMax=30;
%主网交互最小功率
GridMin=-30;
%% 调用mopso函数
mm=mopso; %调用mopso函数
nn=length(mm.swarm); %非支配解数目
%% 比较不同目标函数寻优对调度结果的影响
%% 第1种.将两个目标函数值归一化相加，取相加后最小的目标值的粒子，即寻找折衷解并画图
%将非支配解中的运行成本和环境保护成本分别赋值给yyy,xxx
% for i=1:nn
% yyy(i)= mm.swarm(1,i).cost(1);
% xxx(i)= mm.swarm(1,i).cost(2);
% end
% m1=max(yyy);
% m2=max(xxx);
%
% for i=1:nn
% object(i)= mm.swarm(1,i).cost(1)./m1+ mm.swarm(1,i).cost(2)./m2;
% f1(i)=mm.swarm(1,i).cost(1)./m1;
% f2(i)=mm.swarm(1,i).cost(2)./m2;
% end
% [m,p]=min(object); %得到有着最小目标值的微粒所在的行数P
% pg=mm.swarm(1,p).x; %pg为折衷解
% Title = sprintf('折衷解情况下');
%% 第2种寻找总成本最低时的解并画图
for i=1:nn
object(i)= mm.swarm(1,i).cost(1)+ mm.swarm(1,i).cost(2);
end
[m,p]=min(object); %得到有着最小目标值的微粒所在的行数P
pg=mm.swarm(1,p).x; %pg为总成本最低时的解
Title = sprintf('总成本最低情况下');
%% 第3种寻找运行成本最低时的解并画图
% for i=1:nn
% object(i)= mm.swarm(1,i).cost(1);
% end
% [m,p]=min(object); %得到有着最小目标值的微粒所在的行数P
% pg=mm.swarm(1,p).x; %pg为运行成本最低时的解
% Title = sprintf('运行成本最低情况下');
%% 第4种寻找环境保护成本最低时的解并画图
% for i=1:nn
% object(i)= mm.swarm(1,i).cost(2);
% end
% [m,p]=min(object); %得到有着最小目标值的微粒所在的行数P
% pg=mm.swarm(1,p).x; %pg为环境保护成本最低时的解
% Title = sprintf('环境保护成本最低情况下');
%% 不同情况下的解赋值
for i=1:24
pg_PV(i)=pg(i);
end
for m=25:48
pg_WT(m-24)=pg(m);
end
for m=49:72
pg_BESS(m-48)=pg(m);
end
for m=73:96
pg_DE(m-72)=pg(m);
end
for m=97:120
pg_MT(m-96)=pg(m);
end
for m=121:144
pg_grid(m-120)=pg(m);
end
%% 画图
figure(2)
plot(pg_PV,'-*')
hold on;
plot(PV,'-r*')
xlim([1 24])
grid
legend('实际消纳光伏功率','预测光伏功率');
xlabel('时间');
ylabel('功率');
title('光伏发电机')
figure(3)
plot(pg_WT,'-*')
hold on;
plot(WT,'-r*')
xlim([1 24])
grid
legend('实际消纳风电功率','预测风电功率');
xlabel('时间');
ylabel('功率');
title('风力发电机')
figure(4)
plot(pg_DE,'-k*')
hold on;
xlim([1 24])
grid
xlabel('时间');
ylabel('功率');
title('柴油发电机')
figure(5)
plot(pg_BESS,'-*')
xlim([1 24])
grid
xlabel('时间');
ylabel('功率');
title('蓄电池')
figure(6)
plot(P_load,'-c*')
xlim([1 24])
grid
xlabel('时间');
ylabel('功率');
title('微电网负荷')
figure(7)
plot(pg_MT,'-m*')
grid
xlabel('时间');
ylabel('功率');
title('燃气轮机出力')
figure(8)
plot(pg_grid,'-g*')
xlim([1 24])
grid
xlabel('时间');
ylabel('功率');
title('主网交互')
figure(9)
plot(pg_BESS);
hold on
plot(pg_PV,'-d')
grid
plot(pg_WT,'-d');
plot( pg_DE,'-d')
plot(pg_grid,'-d')
plot(pg_MT,'-d')
legend('BESS','PV','WT','DE','grid','MT');
xlabel('时间/h')
ylabel('功率/kw')
title(Title)
2）子函数
%求解多目标粒子群的适应度
function [c,result]=fitness(x)
global P_load;
global buy_price;
global sell_price;
%% 约束条件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DE爬坡约束处理
DE_sum=0;
for i=73:95
DE_temp(i)=abs(x(i+1)-x(i));
DE_sum=DE_sum+DE_temp(i);
end
f_DE=0;
%阶梯式惩罚系数
if(DE_sum<=345)
f_DE=0;
elseif(DE_sum>345&&DE_sum<=400)
f_DE=2;%%%%%迭代次数
elseif(DE_sum>445&&DE_sum<=500)
f_DE=5;
else
f_DE=20;
end
%% MT爬坡约束处理
MT_sum=0;
for i=97:119
MT_temp(i)=abs(x(i+1)-x(i));
MT_sum=MT_sum+MT_temp(i);
end
f_MT=0;
%MT爬坡阶梯惩罚系数(未满足爬坡约束惩罚)
if(MT_sum<=345)
f_MT=0;
elseif(MT_sum>345&&MT_sum<=400)
f_MT=2;%%%%%迭代次数
elseif(MT_sum>445&&MT_sum<=500)
f_MT=5;
else
f_MT=20;
end
%% 储能SOC约束
%储能荷电状态
SOCMax=150;
SOCMin=5;
SOC_sum=50; %初始容量
SOC_sum_delt=0;
n=0.9;%充放电效率
for i=49:72
if x(i)>0 %充放电
n=-1/0.9;
else
n=0.9;
end
SOC_sum=SOC_sum+n*x(i);
if SOC_sum>SOCMax
SOC_sum_delt= SOC_sum_delt+max(0,SOC_sum-SOCMax);
end
if SOC_sum<SOCMin
SOC_sum_delt= SOC_sum_delt+abs(SOC_sum-SOCMin);
end
end
f_SOC=0.05;
%SOC容量阶梯惩罚系数(未满足SOC约束惩罚)
if(SOC_sum_delt<=0)
f_SOC=0;
elseif(SOC_sum_delt>0&&SOC_sum_delt<=10)
f_SOC=1;%%%%%迭代次数
elseif(SOC_sum_delt>10&&SOC_sum_delt<50)
f_SOC=2;
elseif(SOC_sum_delt>50&&SOC_sum_delt<=100)
f_SOC=5;
else
f_SOC=20;
end
%% 电功率平衡约束处理
ele_sum=0;
for i=1:24
ele_temp(i)=abs(x(i)+x(i+24)+x(i+48)+x(i+72)+x(i+96)+x(i+120)-P_load(i));
ele_sum=ele_sum+ele_temp(i);
end
f_ele=0;
%电平衡阶梯惩罚系数(未满足电平衡约束惩罚)
if(ele_sum==0)
f_ele=0.0;
elseif(ele_sum>0&&ele_sum<=100)
f_ele=1;
elseif(ele_sum>100&&ele_sum<=500)
f_ele=5;
elseif(ele_sum>500&&ele_sum<=800)
f_ele=10;
else
f_ele=50;
end
%% 判断是否为可行解
if ele_sum>4500
c=1;
else
c=0;
end
%% 目标函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%目标函数1：运行成本
%% BESS&&DE&&MT (运维成本 && 燃料成本)
C_DE=0;C_BESS=0;C_MT=0;
for i=1:144
if i>48&&i<73
C_BESS=C_BESS+(0.026)*abs(x(i));%运维成本
elseif i>72&&i<97
C_DE=C_DE+(0.128*x(i))+(0.00011*x(i)^2+0.180*x(i)+6); %运维成本 && 燃料成本
elseif i>96&&i<121
C_MT=C_MT+(0.0293*x(i))+2.55/9.8*x(i)/(0.0753*x(i)^3-0.3095*x(i)^2+0.4174*x(i)+0.1068); %运维成本 && 燃料成本
end
end
C_OM_F= C_DE+ C_MT+ C_BESS;
%% Grid成本
C_grid=0;
for i=121:144
if x(i)>0
C_grid=C_grid+buy_price(i-120)*x(i);
else
C_grid=C_grid-sell_price(i-120)*x(i);
end
end
result=C_grid+C_OM_F+f_DE*DE_sum+f_MT*MT_sum+f_SOC*SOC_sum_delt+f_ele*ele_sum;
function [REP]= mopso(c,iw,max_iter,lower_bound,upper_bound,swarm_size,rep_size,grid_size,alpha,beta,gamma,mu,problem)
%mopso 是多目标粒子群优化的实现
% 最小化问题的技术
%% 初始化参数
global PV;
global WT;
%蓄电池最大放电功率（正表示为电负荷供电，即放电）
BESSMax_dischar=30;
%蓄电池最大充电功率
BESSMax_char=-30;
%柴油机最大发电功率
DEMax=30;
%柴油机最小发电功率
DEMin=6;
%燃气轮机最大发电功率
MTMax=30;
%燃气轮机最小发电功率
MTMin=3;
%主网交互最大功率(正表示为电负荷供电)
GridMax=30;
%主网交互最小功率
GridMin=-30;
%% 种群初始化
if nargin==0 %nargin是判断输入变量个数的函数
c = [0.1,0.2]; % 加速因子
iw = [0.5 0.001]; % 惯性因子
max_iter =100; % 最大迭代次数
%各设备出力约束
for n=1:144 %粒子长度为144（光伏，风电，储能，柴油,燃气轮机，主网的6*24个小时出力）
if n<25
lower_bound(n)=0;
upper_bound(n) =PV(n);
end
if n>24&&n<49
lower_bound(n)=0;
upper_bound(n) =WT(n-24);
end
if n>48&&n<73
lower_bound(n)=BESSMax_char;
upper_bound(n) =BESSMax_dischar;
end
if n>72&&n<97
lower_bound(n)=DEMin;
upper_bound(n) =DEMax;
end
if n>96&&n<121
lower_bound(n)=MTMin;
upper_bound(n) =MTMax;
end
if n>120
lower_bound(n)=GridMin;
upper_bound(n) =GridMax;
end
end
swarm_size=100; % 种群个数
rep_size=100; % 存档库大小
grid_size=7; % 每个维度的网格数
alpha=0.1; % 通货膨胀率
beta=2; % 领导人选择压力
gamma=2; % 删除选择压力
mu=0.1; % 变异速率
problem=@prob; % 创建函数句柄为problem，函数为pro，可以简单理解为调用
end
%% 初始化粒子
fprintf('初始化种群中
')
w = @(it) ((max_iter - it) - (iw(1) - iw(2)))/max_iter + iw(2); %更新惯性因子--改进粒子群算法
pm = @(it) (1-(it-1)/(max_iter-1))^(1/mu); %类比遗传算法引入变异操作，更新变异速率，在particle函数的78-84行
swarm(1,swarm_size) = Particle(); %调用Particle函数，从obj中得到swarm_size
for i = 1:swarm_size
swarm(i)=Particle(lower_bound,upper_bound,problem);%调用Particle函数
retry = 0;
while swarm(i).infeasablity > 0 && retry < 100 %循环条件为：无不可行解且次数低于100
swarm(i)=Particle(lower_bound,upper_bound,problem);%调用Particle函数
retry = retry + 1;
end
end
REP = Repository(swarm,rep_size,grid_size,alpha,beta,gamma); %调用Repository函数
%% 算法循环
fprintf('优化算法开始循环中
')
for it=1:max_iter
leader = REP.SelectLeader(); %选择领导者
wc = w(it); %目前的惯性因子
pc=pm(it); %目前的变异因子
for i =1:swarm_size %更新种群
swarm(i)=swarm(i).update(wc,c,pc,leader,problem);
end
REP = REP.update(swarm);
Title = sprintf('迭代第 %d 次 , 存档库内非支配解个数 = %d',it,length(REP.swarm));
PlotCosts(swarm,REP.swarm,Title) %调用下面的PlotCosts函数
disp(Title);
end
function PlotCosts(swarm,rep,Title)
%画出粒子群的动态
figure(1)
feasable_swarm = swarm([swarm.infeasablity]==0); %可行解
infeasable_swarm = swarm([swarm.infeasablity]>0); %非可行解
LEG = {};
if ~isempty(feasable_swarm)
swarm_costs=vertcat(feasable_swarm.cost);
plot(swarm_costs(:,1),swarm_costs(:,2),'go')
hold on
LEG = {'mospo的可行解'};
Title = sprintf([Title '
可行解的个数=%d'],length(feasable_swarm));
end
if ~isempty(infeasable_swarm)
swarm_costs=vertcat(infeasable_swarm.cost);
plot(swarm_costs(:,1),swarm_costs(:,2),'ro')
hold on
LEG = [LEG, 'mopso的非可行解'];
if contains(Title,newline)
Title = sprintf([Title ', 非可行解的个数 =%d'],length(infeasable_swarm));
else
Title = sprintf([Title '
可行解的个数=%d'],length(infeasable_swarm));
end
end
rep_costs=vertcat(rep.cost);
plot(rep_costs(:,1),rep_costs(:,2),'b*')
xlabel('目标函数1：运行成本')
ylabel('目标函数2：环境保护成本')
grid on
hold off
title(Title)
legend([LEG ,'存档库内非占优解'],'location','best')
drawnow
figure(2)
plot(rep_costs(:,1),rep_costs(:,2),'m*')
xlabel('运行成本')
ylabel('环境保护成本')
grid on
hold off
title('pareto前沿解集')
drawnow
end
end
%多目标粒子这里的程序尽量不要动，直接套用即可
classdef Particle
properties %定义与obj相关的类变量
x
l
u
v
cost
infeasablity
pBest
pBestCost
pBestinfeasablity
GridIndex
isDominated
end
methods
function obj = Particle(lower,upper,problem) %得到粒子
if nargin > 0
obj.GridIndex = 0;
obj.isDominated = false;
obj.x = unifrnd(lower,upper);
obj.l = lower;
obj.u = upper;
obj.v = zeros(1,max(length(lower),length(upper)));
[obj.cost, obj.infeasablity] = problem(obj.x);
obj.pBest = obj.x;
obj.pBestCost = obj.cost;
obj.pBestinfeasablity = obj.infeasablity;
end
end
function obj = update(obj,w,c,pm,gBest,problem) %为粒子更新作准备
obj = obj.updateV(w,c,gBest);
obj = obj.updateX();
[obj.cost, obj.infeasablity] = problem(obj.x);
obj = obj.applyMutatation(pm,problem);
obj = obj.updatePbest();
end
function obj = updateV(obj,w,c,gBest) %更新粒子速度
obj.v = w.*obj.v + c(1).*rand.*(obj.pBest-obj.x) + c(2).*rand.*(gBest.x-obj.x);
end
function obj = updateX(obj) %更新粒子位置
i=find(or(obj.x+obj.v>obj.u,obj.x+obj.v<obj.l));
obj.v(i) = -obj.v(i);
obj.x = max(min(obj.x+obj.v,obj.u),obj.l);
end
function obj = updatePbest(obj) %更新个体最优
if obj.infeasablity == 0
if obj.pBestinfeasablity > 0
obj.pBest = obj.x;
obj.pBestCost = obj.cost;
obj.pBestinfeasablity = obj.infeasablity;
elseif all(obj.pBestCost >= obj.cost) && any(obj.pBestCost > obj.cost)
obj.pBest = obj.x;
obj.pBestCost = obj.cost;
obj.pBestinfeasablity = obj.infeasablity;
end
else
if obj.pBestinfeasablity >= obj.infeasablity
obj.pBest = obj.x;
obj.pBestCost = obj.cost;
obj.pBestinfeasablity = obj.infeasablity;
end
end
end
function obj = applyMutatation(obj,pm,problem) %粒子变异约束
if rand<pm
X=obj.Mutate(pm);
[X.cost,X.infeasablity]=problem(X.x);
if X.dominates(obj)
obj=X;
elseif ~obj.dominates(X)
if rand<0.5
obj=X;
end
end
end
end
function obj=Mutate(obj,pm) %粒子变异--改进粒子群算法
nVar=numel(obj.x);
j=randi([1 nVar]);
dx=pm*(obj.u(j)-obj.l(j));
lb=max(obj.x(j)-dx,obj.l(j));
ub=min(obj.x(j)+dx,obj.u(j));
obj.x(j)=unifrnd(lb,ub);
end
function d = dominates(obj,obj1) %得到支配个体数
if obj.infeasablity == 0
if obj1.infeasablity == 0
if all(obj.cost <= obj1.cost) && any(obj.cost < obj1.cost)
d = true;
else
d = false;
end
else
d = true;
end
elseif obj1.infeasablity == 0
d = false;
elseif obj.infeasablity < obj1.infeasablity
d = true;
else
d = false;
end
end
function obj=updateGridIndex(obj,Grid) %更新网格数
nObj=length(obj.cost);
nGrid=length(Grid(1).LB);
GridSubIndex=zeros(1,nObj);
for j=1:nObj
GridSubIndex(j)=find(obj.cost(j)<Grid(j).UB,1,'first');
end
obj.GridIndex=GridSubIndex(1);
for j=2:nObj
obj.GridIndex=obj.GridIndex-1;
obj.GridIndex=nGrid*obj.GridIndex;
obj.GridIndex=obj.GridIndex+GridSubIndex(j);
end
end
end
methods (Static)
function swarm = updateDomination(swarm) %更新支配个体
for index = 1:length(swarm)
swarm(index).isDominated = false;
for i = 1:length(swarm)
if i == index
continue
end
if swarm(i).dominates(swarm(index))
swarm(index).isDominated = true;
break
end
end
end
end
end
end
%得到多目标问题的解
function [y,c] = prob(x) %c=1则x为非可行解
[c,y(1)] = fitness(x); %得出粒子的适应度 %%运行成本赋值给y（1）
%% 目标函数2:环境保护成本
C_DE_en=0;C_grid_en=0;C_MT_en=0;
for i=1:144
if i>72&&i<97
C_DE_en=C_DE_en+(0.023*680+6*0.306+8*10.09)*x(i);
elseif i>96&&i<121
C_MT_en=C_MT_en+((0.023*889+6*1.8+8*1.6)*6*x(i));
end
end
for i=121:144
if x(i)>0
C_grid_en=C_grid_en+(0.023*724+6*0.0036+8*0.2)*x(i);
else
C_grid_en=C_grid_en+0;
end
end
C_en= C_DE_en+ C_MT_en+ C_grid_en;
y(2) =C_en; %环境保护成本赋值给y（2）
end
%%存档库（储藏室）作用：存档，用来储存过去产生的非支配解
%这里的程序尽量不要动，直接套用即可
classdef Repository %储藏室
properties
swarm
rep_size
Grid
grid_size
alpha
beta
gamma
end
methods
function obj = Repository(swarm,rep_size,grid_size,alpha,beta,gamma) %存档库
if nargin>0
obj.rep_size = rep_size;
swarm = Particle.updateDomination(swarm);
obj.swarm = swarm(~[swarm.isDominated]);
obj.grid_size=grid_size;
obj.alpha=alpha;
obj.beta = beta;
obj.gamma = gamma;
obj.Grid=obj.grid();
for i = 1:length(obj.swarm)
obj.swarm(i) = obj.swarm(i).updateGridIndex(obj.Grid);
end
end
end
function Grid = grid(obj) %更新网格
C = vertcat(obj.swarm.cost);
cmin = min(C,[],1);
cmax = max(C,[],1);
dc = cmax - cmin;
cmin = cmin - obj.alpha * dc;
cmax = cmax + obj.alpha * dc;
nObj = size(C,2);
empty_grid.LB = [];
empty_grid.UB = [];
Grid = repmat(empty_grid,nObj,1);
for j = 1:nObj
cj = linspace(cmin(j),cmax(j),obj.grid_size+1);
Grid(j).LB = [-inf, cj];
Grid(j).UB = [cj, +inf];
end
end
function leader = SelectLeader(obj) %选择领导者
GI = [obj.swarm.GridIndex];
OC = unique(GI);
N = zeros(size(OC));
for k = 1:length(OC)
N(k) = length(find(GI==OC(k)));
end
P = exp(-obj.beta*N);
P = P/sum(P);
sci = Repository.RouletteWheelSelection(P);
sc = OC(sci);
SCM = find(GI==sc);
smi = randi([1 length(SCM)]);
sm = SCM(smi);
leader = obj.swarm(sm);
end
function obj = DeleteOneRepMemebr(obj) %删除一个存档库的成员
GI=[obj.swarm.GridIndex];
OC=unique(GI);
N=zeros(size(OC));
for k=1:length(OC)
N(k)=length(find(GI==OC(k)));
end
P=exp(obj.gamma*N);
P=P/sum(P);
sci=Repository.RouletteWheelSelection(P);
sc=OC(sci);
SCM=find(GI==sc);
smi=randi([1 length(SCM)]);
sm=SCM(smi);
obj.swarm(sm)=[];
end
function obj = update(obj,swarm) %更新
swarm = Particle.updateDomination(swarm);
obj.swarm = [obj.swarm,swarm(~[swarm.isDominated])];
obj.swarm = Particle.updateDomination(obj.swarm);
obj.swarm = obj.swarm(~[obj.swarm.isDominated]);
obj.Grid=obj.grid();
for i = 1:length(obj.swarm)
obj.swarm(i) = obj.swarm(i).updateGridIndex(obj.Grid);
end
Extra=length(obj.swarm)-obj.rep_size;
for e=1:Extra
obj=obj.DeleteOneRepMemebr();
end
end
end
methods (Static)
function i = RouletteWheelSelection(P) %轮盘赌筛选进行删除
i = find(rand<=cumsum(P),1,'first');
end
end
end