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

news/2024/2/21 4:20:50

1.引言

        配电网重构作为配电网优化运行的手段之一,通过改变配电网的拓扑结构,以达到降低网损、改善电压分布、提升系统的可靠性与经济性等目的。近年来,随着全球能源消耗快速增长以及环境的日趋恶化,清洁能源飞速发展,分布式电源(Distributed Generator, DG)大量接入配电网中。DG 因其随机性和波动性,大量接入给配电网带来巨大冲击,也给配电网重构带来严峻的挑战,因此有必要研究适用于高比例清洁能源接入下的配电网重构方法。本文在高比例清洁能源接入的背景下,提出计及需求响应的配电网重构模型,有效利用需求响应进一步降低配电网重构费用并减少弃风弃光率,提高配电网对清洁能源的消纳能力。在求解算法方面,本文基于混合整数二阶锥规划对配电网重构模型进行求解,针对配电网重构的非凸模型,通过引入中间变量并对配电网重构模型进行合理二阶锥松弛,获得混合整数凸规划模型并进行求解。

2.原理介绍

2.1目标函数

        以提高清洁能源的消纳率、减小配电网的运行成本为目标,本文在高比例清洁能源接入背景下提出了考虑需求响应的配电网重构模型。该模型的目标函数为配电网综合运行成本最小,其中考虑了网损成本、弃风弃光成本以及分段开关操作惩罚成本。目标函数的计算方法如式(1)所示。

2.2约束条件

①.潮流约束:

这篇文章采用了极坐标形式描述节点电压功率平衡约束,比较少见。而且功率平衡表达式存在问题,只有节点功率而没有将负荷与电源功率分开。因此我们采用通常的做法,是使用DistFlow模型约束系统潮流:

②.节点电压约束

在进行配电网重构时,配电网中节点电压需要满足上下限约束,即

③.支路电流约束

④.DG功率约束

⑤.网络结构约束

        原文中的公式(9)存在一定的问题,采用虚拟潮流约束还需加上虚拟电源,不然无法平衡。

⑥.储能约束

⑦.电容器约束

⑧.需求响应约束

3.文献中的问题分析

3.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;

3.2 潮流约束不正确

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

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

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

3.3 需求响应模型描述不清

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

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

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

3.4 算例信息描述不清

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

4.部分Matlab代码

        由于该问题的数学模型是一个大规模的混合整数二阶锥规划问题,我测试了cplex和gurobi,收敛精度为0.01的情况下,两个求解器分别需要6小时和3小时左右才能求出最优解。收敛精度同样是0.01,mosek求解器5s左右即可求出最优解。其中部分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:nlbranch_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 == 0disp('求解成功');
elsedisp('运行出错');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;

5.代码运行结果与原文献对比

        原文中数据提供不全,所以代码复现结果会有偏差,但原理完全一样,并修改了原文献中笔误的地方。以下所有结果图,均为复现结果在前,原文结果在后。

5.1 图2

5.2 图3

5.3 表1

5.4 图4

5.5 图5

5.6 图6

5.7 表2

5.8 图7

5.9 图8

6.完整代码获取

        想要完整代码可以戳这个链接(代码运行需要提前安装好matpower+yalmip+mosek工具箱):

https://download.csdn.net/download/weixin_44209907/87788838


http://www.ppmy.cn/news/71351.html

相关文章

Ubuntu安装时没注册root用户密码,怎么登录root

一、Ubuntu设置root密码: 1.桌面空白处右键菜单 点击“打开终端“。 2.输入命令 "sudo passwd root

内存一致性(Memory Consistency)模型简介

这里写自定义目录标题 1. 前言2 为什么需要内存一致性(Memory Consistency)模型3. 什么是内存一致性(Memory Consistency)模型4. 各种内存一致性(Memory Consistency)模型4.1 顺序一致性(SC: Sequential Consistency)模型4.2 完全存储定序(TSO: Total Store Order)模型4.3 部分…

【SQLServer】sqlserver数据库导入oracle

将sqlserver数据库导入到oracle 实用工具: SQL Server Management Studio 15.0.18424.0 SQL Server 管理对象 (SMO) 16.100.47021.07eef34a564af48c5b0cf0d617a65fd77f06c3eb1 Microsoft Analysis Services 客户端工具 15.0.19750.0 Microsoft 数据访问组件 (MDAC) …

SQL常用的内置函数

replace() replace(str,‘a’,‘b’) 将字符串str 中的a字符串替换为b。 regexp_replace() regexp_replace(str,‘正则表达式’,b) 将字符串str中正则匹配的地方替换为b。 translate() translate(str,‘str1’,‘str2’) 对str字符串和str1字符串进行比对,字符相同的…

从眼中窥视:Google AI 模型如何通过眼睛预测你的年龄

新的模型可以通过分析眼部照片揭示衰老的秘密 近年来,谷歌一直在研究各种人工智能模型,可以分析眼睛(内部和外部)的图像并监测某些参数。正如之前提到的,开发能够从眼睛中提取信息的 AI 模型意味着能够以经济高效和无创…

大疆无人机 MobileSDK(遥控器/手机端)开发 v5版<1>

文章目录 概要整体架构流程技术细节SDK 架构体系概述层级架构智能任务空白项目集成 MSDK新建空白项目新建 MyApplication.kt 文件修改 build.gradle(Module) 文件修改 AndroidManifest.xml 文件修改 MainActivity.kt 文件导入 UXSDK 开源框架4.X 和 5.X 版本差异说明DJIKey差异…

Cloud Studio 内核升级之触手可及

前言 Cloud Studio是基于浏览器的集成式开发环境(IDE),为开发者提供了一个永不间断的云端工作站。用户在使用 Cloud Studio 时无需安装,随时随地打开浏览器就能使用。云端开发体验与本地几乎一样,上手门槛更低&#x…

[C++]AVL树、红黑树以及map、set封装

目录 前言: 1 AVL树 1.1 AVL树的概念 1.2 AVL树结点的定义 1.3 AVL树插入 1.4 插入结点的调整 1.5 AVL树的旋转调整 1.5.1 右单旋 1.5.2 左单旋 1.5.3 左右双旋 1.5.4 右左双旋 1.5.4种旋转的判断方式 2 红黑树 2.1 红黑树概念 2.2 红黑树与AVL树的比…

轻松实现一个Python+Selenium的自动化测试框架

首先你得知道什么是Selenium? Selenium是一个基于浏览器的自动化测试工具,它提供了一种跨平台、跨浏览器的端到端的web自动化解决方案。Selenium主要包括三部分:Selenium IDE、Selenium WebDriver 和Selenium Grid。 Selenium IDE&#xff1…

Unity VR开发教程 OpenXR+XR Interaction Toolkit 番外(一)用 Grip 键, Trigger 键和摇杆控制手部动画

文章目录 📕制作手部动画📕设置 Animation Controller📕添加触摸摇杆的 Input Action📕代码部分 在大部分 VR 游戏中,手部的动画通常是由手柄的三个按键来控制的。比如 Grip 键控制中指、无名指、小拇指的弯曲&#xf…

【JS】1686- 重学 JavaScript API - Clipboard API

🏝 1. 什么是 Clipboard API 1.1 概念介绍 Clipboard API[1] 是一组 JavaScript API,用于在浏览器中操作剪贴板。通过 Clipboard API,开发者可以将文本、图片和其他数据复制到剪贴板,也可以从剪贴板中读取数据,实现复制…

Python语言基本控制结构

Python语言基本控制结构包括:条件语句:if、elif、else 循环语句:for、while 跳转语句:break、continue、return 下面是它们的基本用法: 条件语句 if condition1: statement1 elif condition2: statement2 else: stat…

【2023秋招】每日一题:P1087-美团3-18真题 + 题目思路 + 所有语言带注释

2023大厂笔试模拟练习网站(含题解) www.codefun2000.com 最近我们一直在将收集到的各种大厂笔试的解题思路还原成题目并制作数据,挂载到我们的OJ上,供大家学习交流,体会笔试难度。现已录入200道互联网大厂模拟练习题&…

Elasticsearch:了解和使用 match 查询

Match query 是针对多个用例的最常见和最强大的查询。 它是一个全文搜索查询,返回符合指定条件的文档。 match query 可以即兴使用来查询多个选项。在我之前的文章 “开始使用 Elasticsearch (2)” 对它有很多的描述。 Match 查询的格式 让我…

C++基础篇:01 简介

1 历史由来 本贾尼.斯特劳斯特卢普,于1979年在贝尔实验室分析UNIX系统分布式内核的流量时,特别希望有一种更加模块化的工具,于是在1979年10月时开始着手开发一块新的编程语言,在c语言的基础上增加了面向对象的机制,这就…

韩国访问学者签证D-2-5材料准备及签证流程

韩国的签证种类很多,对于申请访问学者签证来说,较常见的签证种类是D-2-5签证和E-3签证,本篇知识人网小编先介绍D-2-5签证。 签证的材料准备 根据韩国大使馆2023年4月12日最新发布的“签证申请与准备材料指导”内容, D-2-5签证的签发对象及准…

Mongodb Shell 常用操作命令

目录 一、启动与关闭mongodb服务 二、进入shell操作 三、常用shell命令 一、启动与关闭mongodb服务 启动:命令: ./mongod -config ../data/mongodb.conf 关闭命令: ./mongod -config ../data/mongodb.conf -shutdown 二、进入shell操作 命令:./mongo 三、常用shell命令 sh…

ES6中Iterator和for...of

1.Iterator 的遍历过程 说明: 遍历对象本质上就是一个指针对象。第一次调用指针对象的next方法,可以将指针指向数据结构的第一个成员。第二次调用指针对象的next方法,可以将指针指向第二个成员。依次进行。 2.默认Iterator接口 说明&…

【JavaEE】wait/notify方法 和 单例模型

目录 前言 1、 wait和notify 1.1、wait()方法 1.2、notify()方法 1.3、wait和sleep 的对比 2、单例模式 2.1、饿汉模式 2.2、懒汉模式 2.3、上述懒汉模式和饿汉模式在多线程情况下是否安全 2.3.1、解决懒汉模式多线程不安去问题 前言 这里补充…

AutoSizer.exe:自动调整窗口大小的便捷工具

AutoSizer.exe是一款实用的桌面应用程序,它旨在帮助用户自动调整窗口大小,提供更好的用户体验。无论您是在使用Windows操作系统进行日常工作还是进行多任务处理,AutoSizer.exe可以简化您的工作流程,提高效率。本文将介绍AutoSizer.exe的下载地址、功能介绍、使用方法以及其…
最新文章