(文章复现)基于断线解环思想的配电网辐射状拓扑约束建模方法-编程知识

参考文献：

[1]王颖,许寅,和敬涵,等.基于断线解环思想的配电网辐射状拓扑约束建模方法[J].中国电机工程学报, 2021, 41(7):9.

1.摘要

配电网采用闭环设计、开环运行的方式，因此在利用数学规划方法解决配电网中的某些优化问题时，需要建立可解析表达的辐射状拓扑约束。该文分析并证明目前使用较多的生成树约束是辐射状拓扑的必要不充分条件，即其无法确保生成的拓扑为辐射状，并给出使其成立的辅助条件。提出基于断线解环思想的辐射状拓扑约束，并证明该组约束是确保辐射状网络的充分必要条件。然后基于该文提出的约束，分别建立配电系统扩建规划和故障恢复问题的混合整数二阶锥优化模型。最后通过配电系统算例测试验证该文提出辐射状约束的有效性，并与常用的生成树约束和单商品流约束进行对比测试，验证该文提出的辐射状约束的优越性。

2.原理介绍

2.1辐射状拓扑约束建模原理

本文提出的辐射状拓扑约束建模基于断线解环思想，即若要确定辐射状拓扑，则原网络中所有环都不会存在，即所有环至少有 1 条线路需要断掉；此外最终确定的系统，所有根节点之间都不相连，因此所有根节点对之间的所有路径至少有 1 条线路需要断掉。因此辐射状拓扑约束可建模为：

式(4)表示对于网络中的所有环，最终的线路个数和小于环的线路个数，即所有环都不成立；式(5)表示对于所有根节点之间的路径，最终线路个数的和小于路径的线路个数，即所有根节点都不相连；式(6)表示最终线路个数和等于节点个数减去根节点个数，即节点–边数目关系。

2.2 配电系统扩建规划问题

1）目标函数

2）约束条件

约束条件(8)为决策变量相关约束，即 xij、yi 为0-1 变量；约束(9)为原有线路和新建线路的关系，若 i 和 j 之间原本有线路，则新建线路变量 xij 为 0，若 i 和 j 之间原本无线路，则 xij小于等于 1，约束(10) 同理；约束(11)和(12)为节点功率平衡约束；约束(13)和(14)为相邻两点电压关系，若节点 i 和 j 之间有线路，则电压关系满足vi - vj = 2(RijPij + XijQij) - (R2 ij + X2 ij)Lij，若无线路，则不作约束；式(15)为节点电压约束；式(16)为线路电流约束，若节点i 和 j 之间有线路，则线路电流约束起作用，若无线路，则两点间电流为 0；式(17)、(18)为线路功率约束，若节点 i 和 j 之间无线路，两点间功率为 0，若无线路，则不作约束；约束(19)为功率定义式的松弛形式；约束(20)为变电站流出功率约束，即流出功率小于变电站 i 相连的变压器的额定功率；式(21)为变电站连通约束，若变电站 n 和 m 同时存在，则约束两站间路径不连通，否则不作约束；式(22)、(23)为式(4)、(6)在本问题中的变形形式，与式(21) 共同组成辐射状拓扑约束。

综上，建立的配电系统扩建规划模型属于混合整数二阶锥规划，可利用成熟的商业软件进行求解。

2.3配电系统故障恢复问题

1）目标函数

式(24)为配电系统故障恢复问题的目标函数，前一项为归一化的最大化加权负荷恢复数目，后一项为最小化网损。

2）约束条件

约束(25)表示决策变量 xij和 zi为 0-1 变量；约束(26)、(27)表示节点功率平衡约束；约束(28)表示分布式电源有功无功功率容量约束；约束(29)、(30) 对应约束(13)、(14)，表示相邻两点电压关系；约束(31)—(33)对应约束(16)—(18)；约束(4)—(6)构成辐射状约束，若采用多源协同的恢复思路，则仅需要约束(4)、(6)，且(6)中 | R| = 1。

可以看出，与配电网扩建规划模型一样，配电网故障恢复问题亦为混合整数二阶锥规划模型，可以用成熟的商业求解器进行求解。

3.编程思路

3.1参数和变量定义

表1 配电网拓建规划的参数

表2 配电网拓建规划的决策变量

表3 配电网故障恢复的参数

表4 配电网故障恢复的决策变量

4.2编程思路

根据对文献内容的解读，可以设计下面的编程思路：

步骤1：输入所需数据

这一步比较简单，大部分所需数据文中都已给出，部分没有提供的数据可以自己假设一下。此外，文中用到的两个配电网数据分别来源于文献[9]和[20]，其中文献[20]是谷歌网盘，点开链接下载就可以得到对应数据；文献[9]是作者本文之前所发会议论文，所谓的改进32节点配电网其实就是我们常用的IEEE33节点配电网，为了方便，直接使用matpower工具箱中数据，并将编号方式和作者相对应，具体如下：

步骤2：定义决策变量

这一步比较简单，按照表2和表4，初始化决策变量即可，同时每个决策变量的维度以及类型(sdpvar还是binvar)不要出错。

步骤3：写目标函数和约束条件

这一步也比较简单，按照给定的数据和定义的变量，分别写出两个优化问题的目标函数和约束条件即可。其中需要重点分析的是约束条件(4)-(6)，也就是辐射状拓扑约束。需要采用图论的方法，遍历配电网中的所有环和根节点之间的所有路径。

步骤4：使用yalmip工具箱对优化问题进行求解

文中构建的两个优化问题都是混合整数二阶锥规划形式，采用mosek求解器的求解性能最佳，如果没有mosek求解器，使用cplex或者gurobi进行求解也可以。其中，main_plan.m文件为配电网拓展规划问题，main_restoration.m为配电网故障恢复问题。

如需切换求解器，在代码中将所需求解器下面三行代码取消注释，其余两个求解器下面的三行代码改为注释即可。

4.部分Matlab代码

程序共有6个m文件和2个mat文件，其中，main_plan.m文件为配电网拓展规划问题的主函数，运行即可得到24节点配电网拓展规划问题的结果，main_restoration.m为配电网故障恢复问题的主函数，运行即可得到改进32节点配电网故障恢复结果。case24.m和IEEE33节点分别为两个配电网的数据文件，main24_topology.m和main32_topology.m分别为采用图论方法遍历配电网中所有环路和根节点之间路径的方法，path_cycle_24.mat和path_cycle_32.mat分别为两个配电网存储的配电网环路和根节点之间的路径。

其中main_plan.m文件部分代码如下：

%% 清除内存空间
clc
clear
close all
warning off%% 系统参数
mpc = case24;
SB = mpc.baseMVA;                           % 基准功率,VA
VB = 20;                                    % 基准电压,kV
ZB = (VB*1e3)^2/(SB*1e6);                   % 基准电阻
NL = length(mpc.branch(:,1));               % 线路总数
NS = length(mpc.substation(:,1));           % 变电站总数
NB = length(mpc.bus(:,1));                  % 节点总数
kl = 0.071;                                 % 配电线路的投资回收率
ks = 0.071;                                 % 变电站的投资回收率
cij = 15020;                                % 配电线路单位长度的建设费用
cf_i = [0;0;750000;950000];                 % 变电站建设费用
l_ij = mpc.branch(:,4);                     % 线路ij所需建设长度
dl = 8760*0.03*0.08*0.8;                    % 配电线路有功损耗的单位电量年本息系数
ds = 8760*0.03*0.08;                        % 变电站的单位电量年本息系数
cv_i = 2800;                                % 变电站 i 单位能量平方的运行费用
Rij = l_ij*0.732/ZB;                        % 配电线路i-j 的电阻
Xij = l_ij*0.732/ZB;                        % 配电线路i-j 的电抗
xij0 = mpc.branch(:,5);                     % 表示为节点 i 和 j 之间原本有无线路，若有为 1，否则为 0
yi0 = [1;1;0;0];                            % 表示节点i处原本有无变电站，若有为1，否则为0
cos_phi = 0.9;                              % 功率因数
duration_time = [2000,5760,1000];           % 三种负荷水平各自的持续时间
wi = (ones(NB,1)*duration_time)/8760;       % 每个节点每种负荷水平的权重
pload_i = sum(mpc.bus(:,2:4).*wi,2)/SB;     % 节点 i 负荷有功功率需求
qload_i = pload_i*sqrt(1 - cos_phi^2);      % 节点 i 负荷无功功率需求
Vmax = 1.1^2;                               % 节点电压最大值的平方
Vmin = 0.9^2;                               % 节点电压最小值的平方
Lij_max = 10;                               % 线路 ij 的电流幅值上限
Srate_i = mpc.substation(:,3)/SB/cos_phi;   % 与变电站 i 相连的变压器的额定功率
load('path_cycle_24.mat')                   % 读取24节点配电网中所有环以及根节点之间的所有路径
branch_to_node = zeros(NB,NL);              % 节点的上游支路
branch_from_node = zeros(NB,NL);            % 节点的下游支路
for k = 1:NLbranch_to_node(mpc.branch(k,3),k) = 1;branch_from_node(mpc.branch(k,2),k) = 1;
end%% 优化变量
xij = binvar(NL,1);                         % 0-1 决策变量，表示线路ij是否建设
yi = binvar(NS,1);                          % 0-1 决策变量，表示变电站i是否建设
Lij = sdpvar(NL,1);                         % 线路ij电流的平方
ps_i = sdpvar(NB,1);                        % 变电站节点 i 的从上级电网得到有功功率
qs_i = sdpvar(NB,1);                        % 变电站节点 i 的从上级电网得到无功功率
vi = sdpvar(NB,1);                          % 节点i电压的平方
Pij = sdpvar(NL,1);                         % 线路ij的有功功率
Qij = sdpvar(NL,1);                         % 线路ij的无功功率%% 约束条件
Constraints = [];。。。省略%% 目标函数
省略%% 设求解器
% gurobi求解器
ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1);
ops.gurobi.TimeLimit = 600;                 % 运行时间限制为10min
ops.gurobi.MIPGap = 0.01;                   % 收敛精度限制为0.01% cplex求解器
% ops = sdpsettings('verbose', 3, 'solver', 'cplex','showprogress',1,'debug',1);
% ops.cplex.timelimit = 600;                  % 运行时间限制为10min
% ops.cplex.mip.tolerances.mipgap = 0.01;     % 收敛精度限制为0.01% mosek求解器
% ops=sdpsettings('verbose', 3, 'solver', 'MOSEK','cachesolvers',1);
% ops.mosek.MSK_DPAR_OPTIMIZER_MAX_TIME=600;% 运行时间限制为10min
% ops.mosek.MSK_DPAR_MIO_TOL_REL_GAP=0.01;  % 收敛精度限制为0.01
sol = optimize(Constraints, objective, ops);%% 分析错误标志
if sol.problem == 0disp('求解成功');
elsedisp('运行出错');yalmiperror(sol.problem)
end%% 运行结果
xij = value(xij);
yi = value(yi);
index_xij = find(xij ~= 0);
index_yi = find(yi ~= 0);
sub_location = 21:24;
objective = value(objective);
disp('************************************新建线路编号************************************')
for k = 1:length(index_xij) - 1fprintf('%d-%d , ', mpc.branch(index_xij(k),2) , mpc.branch(index_xij(k),3));
end
fprintf('%d-%d \n', mpc.branch(index_xij(end),2) , mpc.branch(index_xij(end),3));
disp('***********************************新建变电站位置***********************************')
for k = 1:length(index_yi) - 1fprintf('%d,', sub_location(index_yi(k)));
end
fprintf('%d \n', sub_location(index_yi(end)));
disp('***********************************拓展规划总成本***********************************')
fprintf('配电网拓展规划总成本为%f元 \n', objective);