(文章复现)《高比例清洁能源接入下计及需求响应的配电网重构》（含matlab代码）_mosek求解器

作者：我家小花儿 | 2024-05-19 13:42:08

踩

mosek求解器

参考文献：[1]黄鸣宇,张庆平,张沈习等.高比例清洁能源接入下计及需求响应的配电网重构[J].电力系统保护与控制,2022,50(01):116-123.DOI:10.19783/j.cnki.pspc.210284.

1.原理介绍

该文献考虑到清洁能源消纳率的提升，将目标函数选取为配电网综合运行成本，包括了网损成本、弃风弃光成本和分段开关操作惩罚成本，考虑了潮流约束、节点电压约束、支路电流约束、DG功率约束、网络结构约束、储能约束、电容器约束、需求响应约束，具体内容这里不再赘述，可以自行下载文献查看。

2.文献中的问题分析

2.1 参数取值不全

文献中部分参数的取值没有提供，代码里参考其他文献给这些参数设定了值：

网络损耗成本系数c1为1元/kWh、弃风弃光惩罚系数1.2元/kWh，c2为以及分段开关操作惩罚成本系数c3为15元/次；

支路最大电流设定为2pu，节点电压上下限分别设定为1.06和0.94；

储能充放电功率上限值设定为0.2MW，容量上下限分别设定为0.15/0.8MWh，充放电效率分别为0.9/0.85，初始荷电状态设定为0.3MWh；

2.2 潮流约束不正确

文献中没有给出潮流平衡约束，式(5)中只有节点功率表达式，而且没有和节点负荷、支路功率以及电源输出功率相关联，无法实现功率平衡。为了模型求解，我在代码里使用的是最常用的二阶锥松弛后的DistFlow潮流约束。详细介绍可以参考我的之前的博客：

基于混合整数二阶锥(MISOCP)的配电网重构(附matlab代码)

另外，虚拟潮流约束中，漏了虚拟电源，同样无法实现功率平衡(只有负荷和支路功率，没有电源咋平衡)，在代码中同样需要加上。

2.3 需求响应模型描述不清

文献对需求响应约束进行描述时，没有解释变量ρi,t的含义，导致模型难以理解，对此我是这样理解的：

变量ρi,t代表没有实施需求响应之前的电价，是一个固定值。实施需求响应后，电价由固定值改为分时电价，负荷高峰期电价很高，平时段电价一般，谷时段电价较低，所以就可以引导一些负荷从高峰段转移到谷时段用电，达到降低负荷峰谷差的目的。初始固定电价没有提供，我在代码里假设为0.35元/kWh。

另外，式(15)也存在一定问题。假设调度时段为一天，T=24，这个公式的含义是一天内系统全部节点的在所有时间内的负荷总量保持不变。这么写似乎负荷不仅可以在时间上转移，还能在空间上转移，从一个节点转向另一个节点，显然是不太合理的。需求响应负荷应该只有在时间上转移的特性，所以公式(15)正确表达形式应该是：

2.4 算例信息描述不清

文献采用改进的IEEE33节点系统，但是光伏和风机的额定容量，电容器组的可接入位置都没说清楚。算例中提供了新能源的渗透率信息，我在代码中假设所有DG的额定容量都相同，那么就可以求出DG的额定容量=系统总负荷量×新能源渗透率/DG数量。算例图中好像在节点33处画了电容器，所以我假设电容器组接在节点33处，单个电容器无功补偿容量为0.15Mvar，节点最大可投切电容器数目为5。

3.部分Matlab代码

由于该问题的数学模型是一个大规模的混合整数二阶锥规划问题，我测试了cplex和gurobi，收敛精度为0.01的情况下，两个求解器分别需要6小时和3小时左右才能求出最优解。收敛精度同样是0.01，mosek求解器20分钟左右即可求出最优解。其中部分matlab代码如下：


%% 清除内存空间
clc
clear
close all
warning off
 
%% 系统参数
mpc = IEEE33;
% 风光负荷曲线
P_wind0=[0.21 0.07 0.11 0.21 0.38 0.42 0.12 0.19 0.22 0.47 0.55 0.71 0.80 0.99 0.89 0.99 0.99 0.98 0.99 0.99 0.98 0.77 0.61 0.19];
P_pv0=[0 0 0 0 0.17 0.24 0.40 0.54 0.60 0.51 0.35 0.29 0.27 0.25 0.18 0.10 0.06 0 0 0 0 0 0 0];
P_L0=[0.37 0.33 0.31 0.28 0.27 0.28 0.28 0.27 0.26 0.24 0.30 0.76 0.82 0.86 0.76 0.54 0.43 0.65 0.81 0.95 0.99 0.91 0.65 0.19];
nb=33;                                      % 节点数
ns=1;                                       % 电源节点数
nl=37;                                      % 支路数
n_pv=2;                                     % 光伏数
n_wind=3;                                   % 风机数
n_ess=2;                                    % 储能数
T=24;                                       % 调度时段总数
F=0.6;                                      % 渗透率
P_DG=sum(mpc.bus(:,3))*F/mpc.baseMVA/5;     % DG额定容量
P_wind_max=P_DG*P_wind0;                    % 风机最大有功
P_pv_max=P_DG*P_pv0;                        % 光伏最大有功
P_load=mpc.bus(:,3)/mpc.baseMVA*P_L0;       % 有功负荷
Q_load=mpc.bus(:,4)/mpc.baseMVA*P_L0;       % 无功负荷
r_ij=mpc.branch(:,3)*ones(1,T);             % 线路电阻
x_ij=mpc.branch(:,4)*ones(1,T);             % 线路电抗
wind=[9 25 32];                             % 风机接入位置
pv=[17 22];                                 % 光伏接入位置
ess=[7 25];                                 % 储能接入位置
Umax=[1;1.06*1.06*ones(32,1)];              % 电压上限的平方
Umin=[1;0.94*0.94*ones(32,1)];              % 电压下限的平方
I_max=2;                                    % 电流上限值
P_ch_max=0.2/mpc.baseMVA;                   % 充电功率上限0.2MW
P_dis_max=0.2/mpc.baseMVA;                  % 放电功率上限0.2MW
E_min=0.15/mpc.baseMVA;                     % 储能容量下限0.15MWh
E_max=0.8/mpc.baseMVA;                      % 储能容量上限0.8MWh
n_ch=0.9;                                   % 充电效率为0.9
n_dis=0.85;                                 % 放电效率为0.85
E0=0.3/mpc.baseMVA;                         % 初始荷电状态为0.3MWh
Q_CB_st=0.15/mpc.baseMVA;                   % 单个电容器无功补偿容量0.15Mvar
N_CB_max=5;                                 % 最大可投切电容器数目
ksai=0.5;                                   % 弹性系数
c1=1;                                       % 网络损耗成本系数1元/kWh
c2=1.2;                                     % 弃风弃光惩罚系数1.2元/kWh
c3=15;                                      % 分段开关操作惩罚成本系数15元/次
rho=zeros(1,24);                            % 分时电价
rho([12:15,19:23])=1.026;                   % 峰时电价
rho([7:11,16:18])=0.691;                    % 平时电价
rho([1:6,24])=0.2561;                       % 谷时电价
rho0=0.35;                                  % 初始节点电价为0.35元/kWh
M=1.1*1.1 - 0.9*0.9;                        % 中间变量                   
P_g_max=10/mpc.baseMVA;                     % 电源有功功率最大值
Q_g_max=10/mpc.baseMVA;                     % 电源无功功率最大值
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
 
%% 优化变量
alpha_ij=binvar(nl,1);                      % 支路开断情况
U_i=sdpvar(nb,T);                           % 电压的平方
I_ij=sdpvar(nl,T);                          % 电流的平方
P_ij=sdpvar(nl,T);                          % 线路有功功率
Q_ij=sdpvar(nl,T);                          % 线路无功功率
P_wind=sdpvar(n_wind,T);                    % 风机输出功率
P_pv=sdpvar(n_pv,T);                        % 光伏输出功率
Q_wind=sdpvar(n_wind,T);                    % 风机输出功率
Q_pv=sdpvar(n_pv,T);                        % 光伏输出功率
P_ch=sdpvar(n_ess,T);                       % 储能充电功率
P_dis=sdpvar(n_ess,T);                      % 储能充电功率
y_ch=binvar(n_ess,T);                       % 储能充电状态
y_dis=binvar(n_ess,T);                      % 储能放电状态
E_ESS=sdpvar(n_ess,T);                      % 储能荷电状态
N_CB=intvar(1);                             % 投切的电容器数量
P_cur=sdpvar(nb,T);                         % 需求响应后的负荷量
P_g=sdpvar(nb,T);                           % 节点注入有功
Q_g=sdpvar(nb,T);                           % 节点注入无功
P_g_dot=sdpvar(nb,1);                       % 虚拟电源
P_L_dot=ones(nb,1);                         % 虚拟负荷
P_ij_dot=sdpvar(nl,1);                      % 虚拟功率
 
%% 约束条件
Constraints = [];
。。。。。。省略
 
%% 目标函数
g1=c1*I_ij.*r_ij*10000;
g2=c2*(sum(ones(n_pv,1)*P_pv_max-P_pv)+sum(ones(n_wind,1)*P_wind_max-P_wind))*10000;
g3=c3*abs(alpha_ij-[ones(32,1);zeros(5,1)]);
objective=sum(g1(:))+sum(g2)+sum(g3);
 
%% 设求解器
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 == 0
    disp('求解成功');
else
    disp('运行出错');
    yalmiperror(sol.problem)
end
 
%% 各项费用计算
cost_all=value(objective); 
cost_pv=value(sum(c2*(sum(ones(n_pv,1)*P_pv_max-P_pv)))*10000);
cost_wind=value(sum(c2*(sum(ones(n_wind,1)*P_wind_max-P_wind)))*10000);
cost_loss=value(sum(g1(:)));
switch_num=value(sum(abs(alpha_ij-[ones(32,1);zeros(5,1)])));
U_i=sqrt(value(U_i));
 
%% 结果呈现
show_result;