【配电网重构】基于粒子群算法实现实现直流潮流计算配电网网架结构重构附matlab代码.zip

clc; clear; close all %%%%%%%%%%%%%%%%%%% 主函数调用的主要分函数 %%%%%%%%%%%%%%%%%%%%%%%%%% % -f03 : 粒子解值集间的减法算子 % -f05 : 粒子解值集与速度值集间的加法算子 % -Tree : 构造初始辐射网络粒子 % -changes : 网络结构分层显示 % -cthd01 : 网络前推后代潮流计算 % -fuzhi ; 新粒子赋值判断 %%%%%%%%%%%%%%%%%%% 主函数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global popsize; %种群规模 global M; %节点原始参数，节点坐标 global N; %节点原始参数，节点功率 global inx ind; global gen; %迭代次数 global k Z; global c1 a; %个体最优导向系数 global c2 i; %全局最优导向系数 global best_fitness; %最优解 global best_in_history; %最优解变化轨迹 global gbest_x; %全局最优解矩阵 global Sb Ub Zb; %基准值 global exetime; %当前迭代次数 gen=200; Ci=0.155; % 线路年运行维护系数 Tm=4000;dj=0.05; % 负荷最大利用小时数和电价 popsize=10; %% 原始参数输入，程序初始化 %% Z=[];Zk=0.2*(0.85+j*0.385); M=[4.05 3.25 4.05 5.4 4.05 2 1.2;5.7 6.2 6.7 4.2 3.05 6.7 5.7]; N=[1 1.50+j*1.13 0.4+j*0.3 0.75+j*0.55 1.2+j*0.9 0.6+j*0.45 0.55+j*0.41] Sb=10; Ub=10.5; Zb=Ub^2/Sb; Te= Tm*dj; %% 负荷最大利用小时*电价 %% 种群初始化 %% [E,results]=Tree(10); %%生成初始辐射网络粒子群 for i=1:10 c(:,:,i,1) = E(7*(i-1)+1:7*i,:); end b(:,:,:,1)=randint(70,7,[-4,4]); %%生成粒子初始速度 for exetime=1:gen %%%%整个粒子群迭代次数 for i=1:10 f(:,:,i,exetime)=c(:,:,i,exetime); end for i=1:10 %%%--------导入网络拓扑矩阵，第一列存支路号，第二列存首节点号，第三列存尾节点号，第四列存支路自阻抗,第五列存尾节点给定功率，第六列存支路长度，第七列存支路投资费用----------%%% Z(:,1,i,exetime)=[1;2;3;4;5;6]; Z(:,2,i,exetime)=results(:,1,i,exetime); Z(:,3,i,exetime)=results(:,2,i,exetime); for j=1:6 Z(j,4,i,exetime)=(Zk*sqrt((M(1,results(j,1,i,exetime))-M(1,results(j,2,i,exetime)))^2+(M(2,results(j,1,i,exetime))-M(2,results(j,2,i,exetime)))^2))/Zb ; Z(j,5,i,exetime)=N(1,results(j,2,i,exetime))/Sb; Z(j,6,i,exetime)=sqrt((M(1,results(j,1,i,exetime))-M(1,results(j,2,i,exetime)))^2+(M(2,results(j,1,i,exetime))-M(2,results(j,2,i,exetime)))^2); Z(j,7,i,exetime)=18.5*sqrt((M(1,results(j,1,i,exetime))-M(1,results(j,2,i,exetime)))^2+(M(2,results(j,1,i,exetime))-M(2,results(j,2,i,exetime)))^2); end end A=[]; for i=1:10 %A1 至 A10 eval(['A',num2str(i),'=','c(:,:,i,exetime)']); %A1 至 A10 都是 7 行 7 列的矩阵 eval(['A=[A;A',num2str(i),'];']); %组合 A=[c1;c2;c3;...] end %% 求适应度函数 %% k=0; for i=1:10 [Im(:,:,i,exetime),Vm(:,:,i,exetime)]=qthd(Z(:,:,i,exetime)) for j=1:6 k=k+Te*Z(j,4,i,exetime)*Zb*(Im(j,:,i,exetime))^2+Ci*Z(j,7,i,exetime) end p(:,1,i,exetime)=k; end c1=2.0; %学习因子参数 c2=2.0; %学习因子参数 gbest_x(:,:,1,exetime)=c(:,:,1,exetime); %全局最优初始值为种群第一个粒子的位置 %% 计算适应值并赋最优值 %% for i=1:popsize q(:,1,i,exetime+1)=p(:,1,i,exetime); if exetime>1 if q(:,1,i,exetime)>p(:,1,i,exetime) %若当前适应值优于个体最优值，则进行个体最优信息的更新 q(:,1,i,exetime)=p(:,1,i,exetime); %适值更新 f(:,:,i,exetime)=c(:,:,i,exetime); %位置坐标更新 else f(:,:,:,exetime)=f(:,:,:,exetime-1); end else q(:,1,i,exetime)=p(:,1,i,exetime); end end C=[]; for i=1:10 %A1 至 A10 eval(['C',num2str(i),'=','f(:,:,i,exetime)']); eval(['C=[C;C',num2str(i),'];']); %组合 C=[C1;C2;C3;...] end %% 计算完适应值后寻找当前全局最优位置并记录其坐标 %% if nnz(p(:,1,:,exetime)==min(p(:,1,:,exetime)))>1 d=find(p(:,1,:,exetime)==min(p(:,1,:,exetime))); r(:,:,:,exetime)=f(:,:,d(1,1),exetime) else r(:,:,:,exetime)=f(:,:,find(p(:,1,:,exetime)==min(p(:,1,:,exetime))),exetime); %全局最优粒子的位置 end if exetime>1 if best_fitness>min(p(:,1,:,exetime)) best_fitness=min(p(:,1,:,exetime)); %全局最优值 else r(:,:,:,exetime)=r(:,:,:,exetime-1); end else best_fitness=min(p(1,1,1,exetime)); end for i=1:popsize gbest_x(:,:,i,exetime)=r(:,:,1,exetime); end best_in_history(exetime)=best_fitness; %记录当前全局最优 G=[]; for i=1:10 %A1 至 A10 eval(['G',num2str(i),'=','gbest_x(:,:,i,exetime)']); eval(['G=[G;G',num2str(i),'];']); %组合 G=[G1;G2;G3;...] end %% 粒子群速度与位置更新 %% c0=rand(1); y1=f03(A,C); %矩阵组合减法规则； y2=f03(G,A); w=0.9-(exetime/gen)*0.5; t1=w*b(:,:,:,exetime); t2=(c1*rand(1))*y1; t3=(c2*rand(1))*y2; b(:,:,:,exetime+1)=t1+t2+t3; A=f05(A,b(:,:,:,exetime+1)) %粒子更新位置 for i=1:10 c(:,:,i,exetime+1) = A(7*(i-1)+1:7*i,:); end for i=1:10 D(:,:,i,exetime+1)=fuzhi(c(:,:,i,exetime+1)) if norm(D(:,:,i,exetime+1))==0 c(:,:,i,exetime+1)=Tree(1); %% 辐射网判断 end end for i=1:10 c(:,:,i,exetime+1)=tril(c(:,:,i,exetime+1),-1); if length(find(c(:,:,i,exetime+1)>0))<6 c(:,:,i,exetime+1)=Tree(1); c(:,:,i,exetime+1)=tril(c(:,:,i,exetime+1),-1); end results(:,:,i,exetime+1)=changes(c(:,:,i,exetime+1),1); end %%% 实时输出结果 %%% if exetime-1>0 % exetime 当前迭代次数 line([exetime-1,exetime],[best_in_history(exetime-1),best_fitness]); title('粒子群最优解变化轨迹'); xlabel('迭代次数 x'); ylabel('最优解 y'); grid on; hold on; end end % 迭代次数结束；