clc
clear

ness=2;
T=24;
Max_Dt=300;
D=26ness;%搜索空间维数（未知数个数）
N=100;%粒子数
Pessmax=ones(N,Tness);%充放电
Pessmin=-Pessmax;
ESSmax=2.5ones(N,ness);%容量
ESSmin=0.5ones(N,ness);
ESS=zeros(N,ness);
xsitemax=33ones(N,ness);
xsitemin=2ones(N,ness);
X_max=[Pessmax ESSmax xsitemax];
X_min=[Pessmin ESSmin xsitemin];
v_1max=0.5ones(N,Tness);
v_1min=-v_1max;
v_2max=0.5ones(N,ness);
v_2min=-v_2max;
v_3max=31ones(N,ness);
v_3min=-v_3max;
V_max=[v_1max v_2max v_3max];
V_min=[v_1min v_2min v_3min];


w_max=0.9;
w_min=0.4;


%%%%%粒子初始化%%%%%

x=Pessmin+rand(N,Tness).(Pessmax-Pessmin);
ESS=ESSmin+rand(N,ness).(ESSmax-ESSmin);
xsite = round(xsitemin+rand(N,ness).(xsitemax-xsitemin)); %初始种群的位置（节点位置）四舍五入取整
v1=v_1min+rand(N,Tness).(v_1max-v_1min);
v2=v_2min+rand(N,ness).(v_2max-v_2min);
v3=v_3min+rand(N,ness).(v_3max-v_3min);
X=[x ESS xsite];
V=[v1 v2 v3];



% 储能约束
SOC1=zeros(N,T+1); %SOC
SOC1(:,1)=0.5ESS(:,1); %初始容量设为50%总容量
SOC2=zeros(N,T+1);
SOC2(:,1)=0.5ESS(:,2);
for i=1:N %储能容量约束设为20%~90%
for t=1:T
if t==24
X(i,24)=-sum(X(i,1:23),2);%保证储能充放电平衡
X(i,48)=-sum(X(i,25:47),2);%保证储能充放电平衡
end
SOC1(i,t+1)=SOC1(i,t)-X(i,t);
if SOC1(i,t+1)<0.2X(i,49)
X(i,t)=SOC1(i,t)-0.2X(i,49);
SOC1(i,t+1)=SOC1(i,t)-X(i,t);
elseif SOC1(i,t+1)>0.9X(i,49)
X(i,t)=SOC1(i,t)-0.9X(i,49);
SOC1(i,t+1)=SOC1(i,t)-X(i,t);
end
SOC2(i,t+1)=SOC2(i,t)-X(i,t+T);
if SOC2(i,t+1)<0.2X(i,50)
X(i,t+T)=SOC2(i,t)-0.2X(i,50);
SOC2(i,t+1)=SOC2(i,t)-X(i,t+T);
elseif SOC2(i,t+1)>0.9X(i,50)
X(i,t+T)=SOC2(i,t)-0.9X(i,50);
SOC2(i,t+1)=SOC2(i,t)-X(i,t+T);
end
end
end


% 初始化pbest和gbest

for i=1:N
fitness(i)=evaluationl2(X(i,:),T);
end
[fitnessgbest, bestindex]=min(fitness);
gbest=X(bestindex,:);
pbest=X;
fitnesspbest=fitness;
soc1=SOC1(bestindex,:);
soc2=SOC2(bestindex,:);

% 迭代
for k=1:Max_Dt
for i=1:N
for j=1:D
w=w_max-(w_max-w_min)k/Max_Dt;%惯性权重更新
c1=(0.5-2.5)k/Max_Dt+2.5; %认知
c2=(2.5-0.5)k/Max_Dt+0.5; %社会认识
V(i,j)=wV(i,j)+c1rand(pbest(i,j)-X(i,j))+c2rand(gbest(1,j)-X(i,j));
if V(i,j)>V_max(i,j)
V(i,j)=V_max(i,j);
elseif V(i,j)<V_min(i,j)
V(i,j)=V_min(i,j);
end
end

X(i,:)=X(i,:)+V(i,:);
X(i,51:52)=round(X(i,51:52));
%对粒子边界处理*****************************
for n=1:D
if X(i,n)>X_max(i,n)
X(i,n)=X_max(i,n);
V(i,n)=-V(i,n);
elseif X(i,n)<X_min(i,n)
X(i,n)=X_min(i,n);
V(i,n)=-V(i,n);
end
end



for t=1:T
if t==24
X(i,24)=-sum(X(i,1:23),2);%保证储能充放电平衡
X(i,48)=-sum(X(i,25:47),2);%保证储能充放电平衡
end
SOC1(i,1)=0.5X(i,49);
SOC1(i,t+1)=SOC1(i,t)-X(i,t);
if SOC1(i,t+1)<0.2X(i,49)
X(i,t)=SOC1(i,t)-0.2X(i,49);
SOC1(i,t+1)=SOC1(i,t)-X(i,t);
elseif SOC1(i,t+1)>0.9X(i,49)
X(i,t)=SOC1(i,t)-0.9X(i,49);
SOC1(i,t+1)=SOC1(i,t)-X(i,t);
end
SOC2(i,1)=0.5X(i,50);
SOC2(i,t+1)=SOC2(i,t)-X(i,t+T);
if SOC2(i,t+1)<0.2X(i,50)
X(i,t+T)=SOC2(i,t)-0.2X(i,50);
SOC2(i,t+1)=SOC2(i,t)-X(i,t+T);
elseif SOC2(i,t+1)>0.9X(i,50)
X(i,t+T)=SOC2(i,t)-0.9X(i,50);
SOC2(i,t+1)=SOC2(i,t)-X(i,t+T);
end
end
%对粒子进行评价，寻找最优值******************

fitness(i)=evaluationl2(X(i,:),T);
if fitness(i)<fitnesspbest(i)
pbest(i,:)=X(i,:);
fitnesspbest(i)=fitness(i);
end
if fitness(i)<fitnessgbest
gbest=X(i,:);
fitnessgbest=fitness(i);
soc1=SOC1(i,:);
soc2=SOC2(i,:);
end
record(1,i)=fitnessgbest;
Xrecord(i,:)=X(i,:);

plot(X())
end

fprintf(1,‘迭代数: %d, 最好: %.4f\n’, k,fitnessgbest);
best(k)=fitnessgbest;
end
以上为粒子群算法进行配电网储能选址定容的代码，运行后两个储能接入位置相同，具体如何编写代码可以使两个储能的接入位置不同