一、引言
1.问题背景
配电网重构是指在满足配电网运行基本约束的前提下，通过改变配电网中一个或多个开关的状态对配电网中一个或多个指标进行优化。这里以一种改进二进制粒子群算法为例，进行配电网的重构研究。
2.二进制粒子群算法
2.1简介
粒子群的基本思想来自于对一群鸟类在觅食过程中行为的模拟，常用于连续空间的优化问题，而配电网中的开关只有断开或者闭合两种状态，属于离散空间的优化问题。而二进制粒子群算法可用于配网重构的优化问题。
2.2Sigmoid函数
Sigmoid函数是一个在生物学中常见的S型函数，也称为S型生长曲线。在其他领域，由于其单增以及反函数单增等性质，$Sigmoid$函数常被用渝将变量映射到$[0,1]$之间。其数学表达式如下：$$S(x)=\frac{1}{1+e^{-x}}\; (1)$$
2.3二进制粒子群算法
在二进制粒子群算法中，每个粒子拥有速度和位置两个变量，第$i$个粒子的位置向量和速度分别为$X_i$和$V_i$，粒子自身的最优位置为$P_i$，群体的最优位置为$P_g$，其中粒子的位置被限定在0或1两个值中，粒子的速度通过$Sigmoid$函数变换到[0,1]中，用于表示粒子的位置在更新时选择0或者1的概率，粒子的迭代规则如式(2)~(4)所示：$$S(V)=\frac{1}{1+e^{-V}}\; (2)$$ $$V_i^{k+1}=\omega V_i^k+C_1R_1(P_i^k-X_i^k)+C_2R_2(P_g^k-X_i^k)\;(3)$$
$$X_i^{k+1}=\left\{ \begin{array}{c}0,R_3式中，$\omega$为惯性因子，$C_1$和$C_2$分别为群体学习因子和个体学习因子，$R_1$~$R_3$为区间$[0,1]$上均匀分布的随机数。
2.4算法的改进
传统二进制粒子群算法容易陷入局部最优。针对这一问题，从两个方面进行改进：1.初始化和更新粒子时考虑配电网拓扑约束以缩小粒子搜索范围，增强算法收敛能力；2.加入变异的机制。
下文将结合配电网重构模型说明具体的改进措施。
二、配电网重构模型
1.目标函数
配电网重构的目标一般有降低系统网损、均衡负荷、提高电压稳定性等，这里选择目标函数为最小化系统网损：$$minP_{loss}=r_i\frac{P_i^2+Q_i^2}{V_i^2}\;(5)$$
2.约束条件
2.1拓扑约束
配电网重构之后，不能出现环路或者孤岛。如果在初始化和更新粒子时，采用完全随机的方式，会出现大量无效的配电网拓扑结构。如果可以将拓扑约束考虑到粒子群初始化和更新的过程中，就可以保证所有的粒子都处于有效的状态，同时也缩小了粒子群的搜索范围，减少局部收敛的可能性。假设先将配网中所有开关闭合，会形成一些回路，在每个回路中断开一个开关，就能避免出现环网结构；同时，在不同环路的公共支路中最多只有一个开关断开，在有n条回路汇集的节点周围最多只有n-1个开关断开，就能避免出现孤岛。
以下图为例，说明考虑拓扑约束是如何对算法进行改进的。
这是一个简单配电网，包含14个节点，13个常闭分段开关和2个常开联络开关。
在初始化粒子群时，如果不考虑拓扑约束，会出现大量的无效状态。在此例的配电网中共有15个可操作开关，则系统的状态一共有$2^{15}=32768$种；如果考虑拓扑约束，先确定系统中存在的两个回路，分别包含支路[3,4,5,7,8,9,14]和[2,3,4,5,6,10,11,12,13,15]，公共支路为[3,4,5]，说明有效状态为$2^{7+4}+2^3×(2^4+2^7)=3200$。由此可见，考虑拓扑约束，可以减少将近90%的无效状态，大大提升算法效率。
2.2其他约束
其他约束条件包括潮流约束，节点电压约束，支路电流约束，支路功率约束等等，网络上和各类文献都有非常详细的说明，这里不再赘述。
3.对算法的改进
为了提高算法效率，一共采取了两点改进措施，其中，对拓扑约束的考虑已有详细说明，下面重点介绍另外变异机制：
粒子群算法比较严重的问题就是容易陷入局部最优，在配电网重构问题中体现的尤其明显。当某一粒子找到一个相对较优的位置时，其他粒子会迅速向其靠拢，如果这个位置不是全局最优点，那么算法就停留在局部最优中了。
对于这一问题，如果在每次迭代的过程中，引入粒子变异的机制，可以得到一定的缓解：每次粒子更新之后，如果存在粒子处于上次迭代的最优位置，那么就令其有一定的概率发生变异，这样就可以在一定程度上避免所有粒子都向局部最优聚集。但这个概率的设置需要比较合理，太大了会降低收敛速度，太小了效果也不明显。
针对传统的二进制粒子群算法，考虑拓扑约束，并引入变异机制，就能缓解算法早熟收敛的问题。
三、算例分析与matlab代码
1.算例分析
算例采用IEEE33节点系统：
该系统包含了33个节点，32个分段开关和5个联络开关，不考虑拓扑约束时状态数量一共有$2^{37}=1.37×10^{11}$，考虑拓扑约束时，系统状态一共包括$2^{4+4+3+3+9}+2^3×(2^4+2^4)+2^1×(2^4+2^3)+2^2×(2^4+2^9)+2^4×(2^4+2^9)+2^3×(2^3+2^3)+2^1×(2^3+2^9)=8.4×10^6$，与原来状态相差5个数量级，大大减少了算法的搜索范围。
2.仿真结果
算法运行时重构过程通过动画进行动态显示：
重构前后节点电压得到了明显的改善：
命令行的输出结果：
3.matlab代码
3.1输入系统参数
系统参数存储在IEEE33.m文件中：
3.2列写存在的回路
首先需要将IEEE33节点系统中所有开关闭合后形成的回路与回路两两之间的公共支路列写出来：
% 将系统中所有开关闭合形成的环路
huan1=[2,3,4,5,6,7,18,19,20,33];
huan2=[3,4,5,22,23,24,25,26,27,28,37];
huan3=[8,9,10,11,21,33,35];
huan4=[9,10,11,12,13,14,34];
huan5=[6,7,8,15,16,17,25,26,27,28,29,30,31,32,34,36];
huan12=[3,4,5];
huan13=[33];
huan15=[6,7];
huan25=[25,26,27,28];
huan34=[9,10,11];
huan45=[34];
3.3初始化粒子群
粒子群算法的参数设置：
function [xlimit_max,xlimit_min,vlimit_max,vlimit_min] = set_up_2()%二进制粒子群的基本设置
xlimit_max=ones(1,37);
xlimit_min=zeros(1,37);
vlimit_max=4*ones(1,37);
vlimit_min=-1*vlimit_max;
end
粒子群初始化时，需考虑拓扑约束：
function [pop1,pop2,pop3,pop_v] = CreatPSO_2(n,vlimit_max,vlimit_min)%二进制粒子群的初始化函数
IEEE33;
kaiguan=ones(n,branch_num+5);
huan;
pop3=zeros(n,5);
i=1;
while i<=n
pop3(i,1)=randi([1,length(huan1)]);
pop3(i,2)=randi([1,length(huan2)]);
pop3(i,3)=randi([1,length(huan3)]);
pop3(i,4)=randi([1,length(huan4)]);
pop3(i,5)=randi([1,length(huan5)]);
duan(1)=huan1(pop3(i,1));
duan(2)=huan2(pop3(i,2));
duan(3)=huan3(pop3(i,3));
duan(4)=huan4(pop3(i,4));
duan(5)=huan5(pop3(i,5));
kaiguan(i,38:end)=duan;
kaiguan(i,duan)=0;
if(huanlu(duan))
i=i+1;
else
kaiguan(i,:)=ones(1,branch_num+5);
continue;
end
end
pop1=kaiguan(:,1:37);
pop2=kaiguan(:,38:end);
for i=1:37
for j=1:n
pop_v(j,i) = vlimit_min(i)+(vlimit_max(i) - vlimit_min(i))*rand; % 初始种群的速度
end
end
end
3.4 $Sigmoid$函数
function S=sigmoid(V)
S=1./(1+exp(-V));
end
3.5粒子群的更新
function new_pop=update_pop(pop_v)%粒子的更新
S=sigmoid(pop_v);
huan;
index=zeros(1,5);
[~,index(1)]=min(S(huan1));
[~,index(2)]=min(S(huan2));
[~,index(3)]=min(S(huan3));
[~,index(4)]=min(S(huan4));
[~,index(5)]=min(S(huan5));
duan(1)=huan1(index(1));
duan(2)=huan2(index(2));
duan(3)=huan3(index(3));
duan(4)=huan4(index(4));
duan(5)=huan5(index(5));
new_pop=ones(1,37);
new_pop(duan)=0;
end
3.6主函数
%% 清空环境
clc;clear;
%% 读取数据
IEEE33;
%% 设置种群参数
sizepop = 50; % 初始种群个数
dim = 37; % 空间维数
ger = 22; % 最大迭代次数
[xlimit_max,xlimit_min,vlimit_max,vlimit_min] = set_up_2();
c_1 = 0.8; % 惯性权重
c_2 = 1.5; % 自我学习因子
c_3 = 1.5; % 群体学习因子
%% 生成初始种群
% 首先随机生成初始种群位置
% 然后随机生成初始种群速度
% 然后初始化个体历史最佳位置，以及个体历史最佳适应度
% 然后初始化群体历史最佳位置，以及群体历史最佳适应度
% 设置禁忌对象和其对应的持续迭代次数
[pop1,pop2,pop3,pop_v] = CreatPSO_2(sizepop,vlimit_max,vlimit_min);
pop=[pop1,pop2,pop3];
gbest = pop1; % 每个个体的历史最佳位置
success = huanlu(pop2);
for j=1:sizepop
if success(j)
% [SAIFI,SAIDI,CAIDI,ASAI,Lambda_load,Gamma_load,U_load] = reliability(pop1(j,:));
% fitness_gbest(j) = SAIDI; % 每个个体的历史最佳适应度
[V,u,Loss] = loss(pop2(j,:));
fitness_gbest(j) = Loss;
else
fitness_gbest(j) = 10^10;
end
end
zbest = pop1(1,:); % 种群的历史最佳位置
forbidden=zbest; % 种群的禁忌对象
fitness_zbest = fitness_gbest(1); % 种群的历史最佳适应度
for j=1:sizepop
if fitness_gbest(j) < fitness_zbest % 如果求最小值，则为<; 如果求最大值，则为>;
zbest = pop1(j,:);
fitness_zbest=fitness_gbest(j);
end
end
%% 粒子群迭代
% 更新速度并对速度进行边界处理
% 更新位置并对位置进行边界处理
% 进行约束条件判断并计算新种群各个个体位置的适应度
% 新适应度与个体历史最佳适应度做比较
% 个体历史最佳适应度与种群历史最佳适应度做比较
% 再次循环或结束
iter = 1; %迭代次数
record = zeros(ger, 1); % 记录器
record_average=zeros(ger, 1);
%画图的一些准备
node_location=[(1:14)',2*ones(14,1);(15:18)',2.5*ones(4,1);(2:5)',ones(4,1);(3:5)',4*ones(3,1);(6:13)',3*ones(8,1)];%节点的位置
line_plot=branch(:,2:3);
figure(1)
while iter <= ger
......省略
end
%% 输出结果
disp(['最小功率损耗：',num2str(fitness_zbest),'kW'])
[MinV1,u1,Loss1] = loss(pop22);
[MinV2,u2,Loss2] = loss(33:37);
disp(['最低电压:',num2str(MinV1*12.66),'kV'])
disp(['断开的支路:',num2str(pop22)])
figure(1)
plot(1:ger,record)
title('算法收敛情况');
xlabel('迭代次数')
ylabel('种群最优适应度')
figure(2)
plot(1:33,u1,'-')
title('各节点电压');
xlabel('节点')
ylabel('电压幅值/kV')
hold on
plot(1:33,u2,'--')
legend('重构后的系统','初始系统')
