(文章复现)考虑分布式电源不确定性的配电网鲁棒动态重构

参考文献：

[1]徐俊俊,吴在军,周力,等.考虑分布式电源不确定性的配电网鲁棒动态重构[J].中国电机工程学报,2018,38(16):4715-4725+4976.

1.摘要

间歇性分布式电源并网使得配电网网络重构过程需要考虑更多的不确定因素。在利用仿射数对分布式电源出力的不确定性进行合理分析与建模基础上，建立以重构周期内开关动作耗费与网络有功损耗等综合成本最低为目标函数，以网络安全运行为约束条件的配电网鲁棒动态重构模型。为精确求解该数学模型，引入基于最佳等距思想的分段线性逼近方法将原目标函数松弛为线性可解形式，并根据对偶定理将模型进一步等效转化为双层混合整数线性规划问题；最后采用列约束生成算法对模型进行高效求解。修改的 PG&E 69节点系统测试分析结果表明，与现有的配电网确定性动态重构方法比较，所提鲁棒动态重构方法在抗系统不确定性扰动方面具有明显的优势。

2.原理介绍

2.1分布式电源出力区间预测

这部分是采用粒子群算法和神经网络对风电和光伏的输出功率进行预测，该部分内容和实现原理比较简单，而且不是文献的重点内容，这里不再过多介绍。这部分需要得到的结果就是下面两个公式和两个图。

2.2 配电网鲁棒动态重构模型

配电网鲁棒重构模型中所有节点注入功率不再用某一确定的预测值模糊表示，而是均以仿射数分别予以刻画，在给定 DG 和负荷不确定范围内搜索到最恶劣波动场景下的最优网损(第一阶段)以及制定出 DG 和负荷处于最恶劣波动场景下，满足网络经济运行的重构方案(第二阶段)。为此，建立如下式所示的配电网鲁棒动态重构数学模型：

1）目标函数。

由目标函数可知，第一阶段是以负荷需求和分布式电源出力的不确定扰动为决策变量，也即基于当前网络拓扑结构计算出不确定扰动最恶劣情形下的最低网损成本；第二阶段则以支路开关状态为决策变量，也即在所有网络可能存在的拓扑结构中，寻求出能够确保重构周期内开关动作耗费与网络有功损耗等综合成本最低的唯一网络拓扑结构。显然，第二阶段目标函数也即鲁棒重构总的目标函数。

2）约束条件。

综上所述，所建立的考虑节点注入功率不确定性的配电网鲁棒动态重构数学模型以式(8)为目标函数，以式(9)—(12)为约束条件。

2.3模型求解方法

3.文献中的问题分析

3.1 分段线性逼近

如文献中所述，对于有功功率和无功功率的平方项，可采用分段线性逼近的方法进行处理。由该方法的原理可知，如果分段数太少，线性化后的结果误差会很大，如果分段数太多，又会增加很多额外的变量，加大求解难度。原文3.3.1小节中设分段数为265段，由此增加了265×2×75×96≈384万个变量，模型求解将非常困难，如果设备不是非常好的话建议还是不要参考该方法(我用自己的垃圾电脑大概尝试过，Yalmip+gurobi跑一整天也才收敛到90%，要求收敛精度1%以内的话估计得跑好几个月，文中算例分析上写的只需要400多秒，不知道是怎么得到的)。这个文献对目标函数进行线性化其实就是为了可以将两阶段鲁棒优化子问题转为对偶问题。但实际上，没有线性化之前的子问题是一个二次规划问题，可以采用KKT条件进行转换，得到的结果更精确，也不需要线性化。

3.2 数学模型的细节问题

1）对于功率平衡方程（12），没有解释变量 $P_{s,it}$ 的的含义，根据分析可知该变量应为主电源节点的输出功率。

2）该文章标题是动态重构，但实际上设置的开关状态变量并未考虑时序性，也就是在长时间内都会保持一种运行方案，其实也就是静态重构，不知道他这个动态重构体现在哪里。

3）目标函数中，支路有功功率和无功功率使用的都是分段后的变量，但是其余约束中支路功率又都用了分段前的变量，上下文没有统一，且没有给出两者之间的关系。

4）不确定变量没有给定波动的范围(也就是96个时段中最多有多少个数据点可以取得波动上限)，那么最恶劣场景一定是所有DG出力取最小值，所有负荷需求取最大值，鲁棒优化结果过于保守。所以我在代码中加入了不确定预算，避免两阶段鲁棒优化的结果过于保守。新增约束公式如下：

3.3 部分参数没有提供

1）文中并没有提供一天之内的负荷预测值，因此代码是找了一个典型日负荷曲线带入。

2）文中没有提供支路的最大有功和无功功率，代码里是参考其他文献设置的功率上限。

3）文中没有提供新增联络线的电阻和电抗，代码中是参考其他文献进行设置。

4）文中并未提供动态重构时最大的动作次数，代码将其设置为8次。

3.4 算例分析结果的问题

原文献中表1的结果显示，确定性重构时一天内总有功损耗为481.844kWh，总成本为302.623元，但是式(7)后的解释将C1设置为0.2，C2设置为0.8，0.8×481.844，得到的结果和302.623完全不同，不知道这个结果是怎么得到的。

4.编程思路

4.1参数和变量定义

表1 相关参数

表2 决策变量

4.2编程思路

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

步骤1：输入所需数据

这一步比较简单。PG&E 69节点系统参数来源于matpower工具箱，部分未提供的参数需要自己假设，然后将所有需要的数据，按照表1的定义格式输入即可。需要注意的是，matpower中69节点系统编号和原文中不完全一致，为了编程更方便，代码中以matpower工具箱所提供的编号方式为准，新的编号见下图（其中红色虚线为文中新增的联络线）：

步骤2：建立确定性优化模型并求解

文中将两阶段鲁棒动态重构模型和确定性动态重构的结果进行对比，因此复现时还需要先求出确定性动态重构的结果，具体结果如下：

步骤3：建立两阶段鲁棒优化模型并求解

可以参考我之前写的博客对该问题的两阶段鲁棒优化形式进行分析(鲁棒优化入门(6)-CSDN博客和鲁棒优化入门(7)-CSDN博客)。标准的两阶段鲁棒优化问题的形式为：

可以采用Yalmip工具箱中的函数depends、getbase、getbasematrix、see写出约束矩阵取值，具体如何操作可以参考我之前的博客（Yalmip使用教程(6)-将约束条件写成矩阵形式-CSDN博客）。

4.部分Matlab代码

程序共有4个m文件和一个mat文件，其中case69.m是69节点系统的数据文件，main_do.m是确定性优化的主程序，运行这个代码即可得到确定性优化结果；main_ro.m是两阶段鲁棒优化的主程序，运行即可得到两阶段鲁棒动态重构的结果；Matrix.m是求系数矩阵的程序，运行即可得到系数矩阵的求解结果，并将结果存在Matrix.mat文件中方便读取，其中main_do.m的部分代码如下所示：

%% 1.确定性动态重构

%% 清除内存空间
clc
clear
close all
warning off

%% 系统参数
mpc = case69;
nb = length(mpc.bus(:,1));                          % 节点数
ns = 1;                                             % 主电源节点
nl = length(mpc.branch(:,1));                       % 支路数目
nT = 96;                                            % 调度时段数
Y_pv = [6,21,46];                                   % 光伏接入节点
Y_wt = [52 64];                                     % 风电接入节点
Data = xlsread('风光负荷数据.xlsx');                % 读取风光负荷数据
P_PV_max = Data(:,2)'/1000/mpc.baseMVA;             % 光伏出力上限
P_PV_min = Data(:,3)'/1000/mpc.baseMVA;             % 光伏出力下限
P_WT_max = Data(:,4)'/1000/mpc.baseMVA;             % 风电出力上限
P_WT_min = Data(:,5)'/1000/mpc.baseMVA;             % 风电出力下限
PL_curve = Data(:,6);                               % 负荷日变化曲线
P_PV0 = (P_PV_max + P_PV_min)/2;                    % 光伏出力均值
dP_PV0 = P_PV_max - P_PV0;                          % 光伏出力最大波动
P_WT0 = (P_WT_max + P_WT_min)/2;                    % 风电出力均值
dP_WT0 = P_WT_max - P_WT0;                          % 风电出力最大波动
phi = 0.85;                                         % DG的功率因数
P_L0 = mpc.bus(:,3)/mpc.baseMVA*PL_curve';          % 有功负荷
Q_L0 = mpc.bus(:,4)/mpc.baseMVA*PL_curve';          % 无功负荷
P_L0(P_L0 == 0) = 1e-6;                             % 加上一个很小的数防止0注入节点出现
Q_L0(Q_L0 == 0) = 1e-6;                             % 加上一个很小的数防止0注入节点出现
R_ik = mpc.branch(:,3);                             % 线路电阻
L_ik0 = mpc.branch(:,11);                           % 初始线路开断状态
C1 = 0.2;                                           % 支路开关动作一次所需要的成本系数
C2 = 0.8;                                           % 网络重构期间有功损耗所对应的成本系数
P_ik_max = 6;                                       % t时段支路 ik 上允许流过的最大有功功率
Q_ik_max = 5;                                       % t时段支路 ik 上允许流过的最大无功功率
Vi = 1;                                             % 根据文献式(6)后的解释将节点电压设为常数1
N = 8;                                              % 最大重构次数
Ps_max = 10;                                        % 上级电源输出有功功率最大值
Ps_min = 0;                                         % 上级电源输出有功功率最小值
Qs_max = 10;                                        % 上级电源输出无功功率最大值
Qs_min = -10;                                       % 上级电源输出无功功率最小值
branch_to_node = zeros(nb,nl);                      % 流入节点的支路
branch_from_node = zeros(nb,nl);                    % 流出节点的支路
for k = 1:nl
    branch_to_node(mpc.branch(k,2),k) = 1;
    branch_from_node(mpc.branch(k,1),k) = 1;
end

%% 决策变量
L_ik = binvar(nl,1);                                % 支路 ik 上开关的状态信息
u_ik = binvar(nb,nb,'full');                        % 表示节点关系

Pikt = sdpvar(nl,nT);                               % t时刻支路ik在l断面的有功功率
Qikt = sdpvar(nl,nT);                               % t时刻支路ik在j断面的无功功率
Ps_it = sdpvar(ns,nT);                              % t时刻主电源节点i的有功出力
Qs_it = sdpvar(ns,nT);                              % t时刻主电源节点i的无功出力

e_Git = zeros(5,nT);                                % t时段导致 DG 节点 i 注入功率不确定的扰动因子(包含光伏和风机)，确定性优化时取值为0
e_Lit = zeros(nb,nT);                               % t时段导致负荷节点i注入功率不确定的扰动因子，确定性优化时取值为0

%% 约束条件
Constraints = [];

%% 约束(10)
此处省略。。。。

%% 约束(11)
此处省略。。。。

%% 约束(12)
此处省略。。。。


%% 目标函数
此处省略。。。。

%% 设求解器
% gurobi求解器
ops = sdpsettings('verbose', 3, 'solver', 'gurobi','showprogress',1,'debug',1);
ops.gurobi.TimeLimit = 7200;                         % 运行时间限制
ops.gurobi.MIPGap = 0.01;                            % 收敛精度限制为0.01

% cplex求解器
% ops = sdpsettings('verbose', 3, 'solver', 'cplex','showprogress',1,'debug',1);
% ops.cplex.timelimit = 7200;                        % 运行时间限制
% 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 = 7200;         % 运行时间限制
% ops.mosek.MSK_DPAR_MIO_TOL_REL_GAP = 0.01;           % 收敛精度限制为0.01

sol = optimize(Constraints,objective,ops);


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

%% 结果
L_ik = value(L_ik);
u_ik = value(u_ik);
% disp('******************重构前******************')
% disp('开断支路为：')
% disp([num2str(mpc.branch(70,1)),'-',num2str(mpc.branch(70,2)),',',...
%       num2str(mpc.branch(71,1)),'-',num2str(mpc.branch(71,2)),',',...
%       num2str(mpc.branch(72,1)),'-',num2str(mpc.branch(72,2)),',',...
%       num2str(mpc.branch(73,1)),'-',num2str(mpc.branch(73,2)),',',...
%       num2str(mpc.branch(74,1)),'-',num2str(mpc.branch(74,2))])
% disp(['系统网损为：','36579.1335kW'])
disp('******************确定性优化重构结果******************')
open_branch = find(L_ik~=1)';
disp('开断支路为：')
disp([num2str(mpc.branch(open_branch(1),1)),'-',num2str(mpc.branch(open_branch(1),2)),',',...
      num2str(mpc.branch(open_branch(2),1)),'-',num2str(mpc.branch(open_branch(2),2)),',',...
      num2str(mpc.branch(open_branch(3),1)),'-',num2str(mpc.branch(open_branch(3),2)),',',...
      num2str(mpc.branch(open_branch(4),1)),'-',num2str(mpc.branch(open_branch(4),2)),',',...
      num2str(mpc.branch(open_branch(5),1)),'-',num2str(mpc.branch(open_branch(5),2))])
disp(['系统网损为：',num2str(value(objective2)*1000*mpc.baseMVA),'kW'])
disp(['开关动作次数为：',num2str(value(objective1)),'次'])
disp(['总运行成本为：',num2str(value(objective)),'元'])

figure
plot(P_PV_max*1000*mpc.baseMVA,'k-','linewidth',1)
hold on
plot(P_PV_min*1000*mpc.baseMVA,'k-','linewidth',1)
hold on
plot(P_PV0*1000*mpc.baseMVA,'r:','linewidth',1)
legend('光伏区间出力上限','光伏区间出力下限','光伏实际出力');
xlabel('时间')
ylabel('功率/kw')

figure
plot(P_WT_max*1000*mpc.baseMVA,'k-','linewidth',1)
hold on
plot(P_WT_min*1000*mpc.baseMVA,'k-','linewidth',1)
hold on
plot(P_WT0*1000*mpc.baseMVA,'r:','linewidth',1)
legend('风电区间出力上限','风电区间出力下限','风电实际出力');
xlabel('时间')
ylabel('功率/kw')

figure
plot(1.1*sum(P_L0)*1000*mpc.baseMVA,'k-','linewidth',1)
hold on
plot(0.9*sum(P_L0)*1000*mpc.baseMVA,'k-','linewidth',1)
hold on
plot(sum(P_L0)*1000*mpc.baseMVA,'r:','linewidth',1)
legend('负荷需求上限','负荷需求下限','负荷实际需求');
xlabel('时间')
ylabel('功率/kw')

经过测试，如果在main_ro.m中将代码81行的收敛精度ops.gurobi.MIPGap设置为0.05时，两阶段鲁棒优化大约需要10min即可收敛；如果将其设置为0.1时，5min左右即可收敛。大家可以根据自身需求对计算精度和运行时间的要求选择合适的收敛精度。