(文章复现)考虑网络动态重构的分布式电源选址定容优化方法

参考文献：

[1]朱俊澎,顾伟,张韩旦,等.考虑网络动态重构的分布式电源选址定容优化方法[J].电力系统自动化,2018,42(05):111-119.

1.摘要

        以投资周期经济收益最高为目标，基于二阶锥规划提出了一种考虑网络动态重构的分布式电源选址定容优化方法。首先，针对闭环设计的配电网结构，提出计及联络线和分段开关状态的拓展DistFlow潮流模型，并基于此建立了考虑网络动态重构的分布式电源配置优化模型。引入“虚拟支路电压”的概念对模型进行线性化处理，提出了基于“有功流”的辐射形拓扑线性约束方法，同时配合电流、电压变量替换和二阶锥松弛，建立了统一优化模型的二阶锥形式。采用添加辅助电压约束的方法，解决了对含电压上限约束时二阶锥模型松弛不紧致的问题。在IEEE标准算例中测试了算法的有效性，结果表明，考虑网络动态重构可以提高电网对分布式电源的消纳能力，同时提升分布式电源投资周期内的总体经济效益。

2.原理介绍

2.1目标函数

        本节基于拓展 DistFlow潮流模型提出考虑网络动态重构的RDG选址定容优化模型。目标函数及控制变量表达式如下:

        各项成本的具体表达式见式(13)和式(14)。规划阶段的固定投资成本由RDG的容量决定,运行阶段的成本包含RDG运维成本Com 、配电网购电成本Cut和网络重构成本Crec 。

2.2 约束条件

2.3模型的简化

        在本文中,分布式电源投资周期设为15年﹐典型场景考虑了春、夏、秋、冬4个典型日内的负荷和RDG出力波动,时刻为一天内24个时刻,则总断面数为1 440个。考虑计算速度的要求和计算机内存容量的限制﹐做出如下简化。

        1)考虑负荷年增长率﹐将投资周期内的负荷简化为低、中、高3个负荷水平﹐分别对应1～5年,6～10年和 11～15年的平均负荷。

        2)每天24个时刻划分为12个时段﹐各时段负荷的取值为该时段内2个时刻的平均值。重构频率最高为每3个时段一次,即每天最多重构4次。

        3)网络重构起到平衡负荷、改善全网电压的作用,在负荷每年等比例增长时,各条线路负载的相对轻重状况和电压在全网的分布走势不发生变化,因此本文假定重构策略不随年份变化。

        经过简化,总断面数缩小为144个,二进制变量规模显著减小,使得算法满足计算效率要求。

3.编程思路

3.1参数和变量定义

表1 相关参数

表2 决策变量

3.2编程思路

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

步骤1：输入所需数据

        这一步比较简单，大部分所需数据文中都已给出，部分没有提供的数据可以自己假设一下。其中IEEE33节点配电网直接使用matpower中的数据。

步骤2：定义决策变量

        这一步比较简单，按照表2，初始化决策变量即可，同时每个决策变量的维度以及类型(sdpvar还是binvar)不要出错。需要注意的是，文中每个优化变量都有12个时刻，四种典型场景和三种负荷大小，共144个断面，定义时比较复杂。代码中为简化表示，将变量按照四种典型场景分别进行定义，每个典型场景对应的变量和约束均一致，命名时有所区分即可。

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

        这一步需要按照给定的数据和定义的变量，分别写出优化问题的目标函数和约束条件。由于优化问题比较复杂，写约束条件的时候很容易出错，需要非常仔细。

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

        文中构建的优化问题都是混合整数二阶锥规划形式，代码中是采用yalmip+gurobi进行求解。

3.3文献中存在的问题分析

1)RDG的接入位置前后不统一

        在文中3.1参数设置的第一段话中，RDG接入位置前后不统一，可能是存在笔误，代码中以17-18与31-33为准。

2)模型中没有考虑可控开关

        在文中3.1参数设置的第一段话中，提到了安装可控开关的支路，表示除了这些支路外其余支路都是不能动作的。但在模型中并未体现这一点，需要增加相关约束，公式如下：

3)模型中没有约束配电网从上级电网的购电量

        文中没有给变量Pi,uit设置上下限，在优化过程中可能会出现购电量为负数的情况，出现功率倒送问题，影响电网运行安全。因此需要给配电网购电量添加上下限约束：

4)光伏出力模型的细节问题

        原文献在附录表B2中给出了光伏四个典型场景的出力曲线，但DG运行模型中，仅用公式(23)约束了DG每个时段的出力不能超过其出力上限，并未考虑光伏出力曲线。此外，式23引入非常多的二次约束，会增加模型求解难度，可以通过给定的光伏出力曲线和光伏电站的容量确定每个时段的光伏出力。

4.部分Matlab代码

%% Case 3:考虑RDG的接入和重构,采用本文提出的方法优化求解。

%% 清除内存空间
clc
clear
close all
warning off
yalmip('clear')

%% 系统参数
mpc = IEEE33;
SB = mpc.baseMVA;                           % 基准功率,MVA
VB = mpc.bus(1,10);                         % 基准电压,kV
NL = length(mpc.branch(:,1));               % 线路总数
Nbus = length(mpc.bus(:,1));                % 节点总数
NT = 12;                                    % 时段数
r = 0.023;                                  % 通货膨胀率
Y = 15;                                     % 投资周期年数
Ocand = [17,18,31,32,33];                   % RDG接入候选节点
N_DG = length(Ocand);                       % RDG最多的连接节点数
O_autoL = [3,7,8,9,13,18,23,27,31,33:37];   % 安装可控开关的支路
N_autoL = length(O_autoL);                  % 安装可控开关的支路数
cinv_fix = 10000;                           % RDG固定投资成本
cinv_mar = 11;                              % RDG边际投资成本系数
com = 150;                                  % RDG运维成本系数
cuti = 0.4*ones(1,NT);                      % 谷时购电电价
cuti(5:10) = 0.7;                           % 峰时购电电价
crec = 0.5;                                 % 单次开关操作成本
S = 4;                                      % 场景总数
load('Load_PV_data.mat')                    % 读取典型负荷曲线和光伏出力曲线
Pi_L_sp = zeros(Nbus,NT);                   % 春季典型场景有功负荷
Qi_L_sp = zeros(Nbus,NT);                   % 春季典型场景无功负荷
Pi_L_su = zeros(Nbus,NT);                   % 夏季典型场景有功负荷
Qi_L_su = zeros(Nbus,NT);                   % 夏季典型场景无功负荷
Pi_L_au = zeros(Nbus,NT);                   % 秋季典型场景有功负荷
Qi_L_au = zeros(Nbus,NT);                   % 秋季典型场景无功负荷
Pi_L_wi = zeros(Nbus,NT);                   % 冬季典型场景有功负荷
Qi_L_wi = zeros(Nbus,NT);                   % 冬季典型场景无功负荷
PL_T = [1.0305, 1.1101, 1.1959];            % 低、中、高不同的负荷水平
Re_node = [6, 14];                          % 商业负荷接入节点
Co_node = setdiff(1:33,[6, 14]);            % 居民负荷接入节点
PL0 = PL0';
for t = 1:NT
    % 6和14节点接入商业负荷
    Pi_L_sp(Re_node,t) = mean(PL0(1, 2*t-1:2*t))*mpc.bus(Re_node,3)/SB;
    Qi_L_sp(Re_node,t) = mean(PL0(1, 2*t-1:2*t))*mpc.bus(Re_node,4)/SB;
    Pi_L_su(Re_node,t) = mean(PL0(3, 2*t-1:2*t))*mpc.bus(Re_node,3)/SB;
    Qi_L_su(Re_node,t) = mean(PL0(3, 2*t-1:2*t))*mpc.bus(Re_node,4)/SB;
    Pi_L_au(Re_node,t) = mean(PL0(5, 2*t-1:2*t))*mpc.bus(Re_node,3)/SB;
    Qi_L_au(Re_node,t) = mean(PL0(5, 2*t-1:2*t))*mpc.bus(Re_node,4)/SB;
    Pi_L_wi(Re_node,t) = mean(PL0(7, 2*t-1:2*t))*mpc.bus(Re_node,3)/SB;
    Qi_L_wi(Re_node,t) = mean(PL0(7, 2*t-1:2*t))*mpc.bus(Re_node,4)/SB;
    
    % 其余节点接入居民负荷
    Pi_L_sp(Co_node,t) = mean(PL0(2, 2*t-1:2*t))*mpc.bus(Co_node,3)/SB;
    Qi_L_sp(Co_node,t) = mean(PL0(2, 2*t-1:2*t))*mpc.bus(Co_node,4)/SB;
    Pi_L_su(Co_node,t) = mean(PL0(4, 2*t-1:2*t))*mpc.bus(Co_node,3)/SB;
    Qi_L_su(Co_node,t) = mean(PL0(4, 2*t-1:2*t))*mpc.bus(Co_node,4)/SB;
    Pi_L_au(Co_node,t) = mean(PL0(6, 2*t-1:2*t))*mpc.bus(Co_node,3)/SB;
    Qi_L_au(Co_node,t) = mean(PL0(6, 2*t-1:2*t))*mpc.bus(Co_node,4)/SB;
    Pi_L_wi(Co_node,t) = mean(PL0(8, 2*t-1:2*t))*mpc.bus(Co_node,3)/SB;
    Qi_L_wi(Co_node,t) = mean(PL0(8, 2*t-1:2*t))*mpc.bus(Co_node,4)/SB;
end
NDG_max = 3;                                    % 分布式电源个数上限
Ncap_max = [20;25;30;35;40];                    % 单个电站的机组上限
Sunit = 50/1e3/SB;                              % 分布式电源单位机组功率    
M = 100;                                        % 足够大的正数
rij = mpc.branch(:,3);                          % 支路电阻
xij = mpc.branch(:,4);                          % 支路电抗
Umax = 1.05;                                    % 节点电压上限
Umin = 0.95;                                    % 节点电压下限
Imax = 10;                                      % 支路电流上限
aij0 = mpc.branch(:,11);                        % 初始开关状态
Ns_d = [91 91 91 92];                           % 四个季节的天数
branch_to_node = zeros(Nbus,NL);                % 节点的上游支路
branch_from_node = zeros(Nbus,NL);              % 节点的下游支路
for k = 1:NL
    branch_to_node(mpc.branch(k,2),k) = 1;
    branch_from_node(mpc.branch(k,1),k) = 1;
end

%% RDG接入变量
xi_loc = binvar(N_DG,1);                        % 0-1决策变量，表示节点i是否接入RDG
Si = sdpvar(N_DG,1);                            % 节点i接入RDG的容量
Ni = intvar(N_DG,1);                            % 节点i分布式电站中机组的个数

%% 春季典型场景的决策变量
aij_sp = binvar(NL,NT/3);                       % 支路ij的开关状态
% RDG有功出力
Pi_DG_sp = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.8;   
% RDG无功出力
Qi_DG_sp = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.6;               
Pi_uit_sp = sdpvar(NT,3);                       % 上级电源有功出力
Qi_uit_sp = sdpvar(NT,3);                       % 上级电源无功出力
Ui_sp = sdpvar(Nbus,NT,3);                      % 节点i电压的平方
Iij_sp = sdpvar(NL,NT,3);                       % 支路电流的平方
Pij_sp = sdpvar(NL,NT,3);                       % 线路ij的有功功率
Qij_sp = sdpvar(NL,NT,3);                       % 线路ij的无功功率

%% 夏季典型场景的决策变量
aij_su = binvar(NL,NT/3);                       % 支路ij的开关状态
% RDG有功出力
Pi_DG_su = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.8;   
% RDG无功出力
Qi_DG_su = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.6;   
Pi_uit_su = sdpvar(NT,3);                       % 上级电源有功出力
Qi_uit_su = sdpvar(NT,3);                       % 上级电源无功出力
Ui_su = sdpvar(Nbus,NT,3);                      % 节点i电压的平方
Iij_su = sdpvar(NL,NT,3);                       % 支路电流的平方
Pij_su = sdpvar(NL,NT,3);                       % 线路ij的有功功率
Qij_su = sdpvar(NL,NT,3);                       % 线路ij的无功功率

%% 秋季典型场景的决策变量
aij_au = binvar(NL,NT/3);                       % 支路ij的开关状态
% RDG有功出力
Pi_DG_au = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.8;   
% RDG无功出力
Qi_DG_au = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.6;   
Pi_uit_au = sdpvar(NT,3);                       % 上级电源有功出力
Qi_uit_au = sdpvar(NT,3);                       % 上级电源无功出力
Ui_au = sdpvar(Nbus,NT,3);                      % 节点i电压的平方
Iij_au = sdpvar(NL,NT,3);                       % 支路电流的平方
Pij_au = sdpvar(NL,NT,3);                       % 线路ij的有功功率
Qij_au = sdpvar(NL,NT,3);                       % 线路ij的无功功率

%% 冬季典型场景的决策变量
aij_wi = binvar(NL,NT/3);                       % 支路ij的开关状态
% RDG有功出力
Pi_DG_wi = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.8;   
% RDG无功出力
Qi_DG_wi = Si*(Ppv0(1:2:23,1) + Ppv0(2:2:24,1))'/2*0.6;   
Pi_uit_wi = sdpvar(NT,3);                       % 上级电源有功出力
Qi_uit_wi = sdpvar(NT,3);                       % 上级电源无功出力
Ui_wi = sdpvar(Nbus,NT,3);                      % 节点i电压的平方
Iij_wi = sdpvar(NL,NT,3);                       % 支路电流的平方
Pij_wi = sdpvar(NL,NT,3);                       % 线路ij的有功功率
Qij_wi = sdpvar(NL,NT,3);                       % 线路ij的无功功率

%% DG配置与运行约束
省略

%% 春季典型场景约束条件
Constraints_sp = [];
省略

%% 夏季典型场景约束条件
Constraints_su = [];
省略

%% 秋季典型场景约束条件
Constraints_au = [];
省略

%% 冬季典型场景约束条件
Constraints_wi = [];
省略


%% 目标函数
% DG投资运维成本
Cinv = sum(xi_loc.*(cinv_fix + cinv_mar*Si*SB*1e6));                    % DG投资成本
Com = sum(com*Si*SB*1e3);                                               % DG运维成本

% 春季场景
Cuti_sp = sum(sum(Pi_uit_sp*Ns_d(1).*(cuti'*ones(1,3))))*SB*1e3;        % 购电成本
Nswitch_sp = sum(abs(aij_sp(:,1) - aij0));                              % 开关动作次数
for k = 2:4
    Nswitch_sp = Nswitch_sp + sum(abs(aij_sp(:,k) - aij_sp(:,k-1)));
end
Crec_sp = Ns_d(1)*crec*Nswitch_sp;

% 夏季场景
Cuti_su = sum(sum(Pi_uit_su*Ns_d(2).*(cuti'*ones(1,3))))*SB*1e3;        % 购电成本
Nswitch_su = sum(abs(aij_su(:,1) - aij0));                              % 开关动作次数
for k = 2:4
    Nswitch_su = Nswitch_su + sum(abs(aij_su(:,k) - aij_su(:,k-1)));
end
Crec_su = Ns_d(2)*crec*Nswitch_su;

% 秋季场景
Cuti_au = sum(sum(Pi_uit_au*Ns_d(3).*(cuti'*ones(1,3))))*SB*1e3;        % 购电成本
Nswitch_au = sum(abs(aij_au(:,1) - aij0));                              % 开关动作次数
for k = 2:4
    Nswitch_au = Nswitch_au + sum(abs(aij_au(:,k) - aij_au(:,k-1)));
end
Crec_au= Ns_d(1)*crec*Nswitch_au;

% 冬季场景
Cuti_wi = sum(sum(Pi_uit_wi*Ns_d(4).*(cuti'*ones(1,3))))*SB*1e3;        % 购电成本
Nswitch_wi = sum(abs(aij_wi(:,1) - aij0));                              % 开关动作次数
for k = 2:4
    Nswitch_wi = Nswitch_wi + sum(abs(aij_wi(:,k) - aij_wi(:,k-1)));
end
Crec_wi = Ns_d(2)*crec*Nswitch_wi;

Ctotal = Cinv + Com + 5*(Cuti_sp + Crec_sp + Cuti_su + Crec_su + Cuti_au + Crec_au) + 15*(Cuti_wi + Crec_wi);

%% 设求解器
% gurobi求解器
ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress', 1 , 'usex0',1,'debug',1, 'gurobi.Heuristics' , 0.2);
ops.gurobi.TimeLimit = 1200;                % 运行时间限制为30min
ops.gurobi.MIPGap = 0.005;                  % 收敛精度限制为0.005

sol = optimize([Constraints_DG, Constraints_sp, Constraints_su, Constraints_au, Constraints_wi], Ctotal, ops);


%% 分析错误标志
if sol.problem == 0
    disp('求解成功');
else
    disp('运行出错');
    yalmiperror(sol.problem)
end

%% 运行结果
show_result;

        以上仅为计算原文献场景3的部分代码，完整matlab代码可从以下链接获取(PS，代码中使用了gurobi求解器，没有的话需要提前安装好；另外，由于文中所提优化问题比较复杂，求解时很容易出现内存不足的情况。如果运行时出现内存不足的情况，一般可以通过修改gurobi最大求解时间和收敛精度来解决)：

(文章复现)考虑网络动态重构的分布式电源选址定容优化方法matlab代码资源-CSDN文库

5.代码运行结果

        原文中数据提供不全，所以代码复现结果和原文献相比会有偏差，但原理完全一样。

5.1 RDG选址定容、重构策略和经济性指标

5.2 配电网电压分布

5.3配电网损耗分布