基于混合整数二阶锥(MISOCP)的配电网重构(附matlab代码)
参考资料：主动配电网网络分析与运行调控 (sciencereading.cn)
配电网重构是指在满足配电网运行基本约束的前提下，通过改变配电网中一个或多个开关的状态对配电网中一个或多个指标进行优化。通过配电网重构，可以在不增加设备投资的情况下，充分发挥配电系统的潜力，提高系统的性能指标，具有较好的经济效益。配电网重构的算法有许多，包括数学规划方法，如分支定界法、0-1 整数规划、单纯形法等；启发式算法，如最优流算法、开关交换法等；智能优化算法，如模拟退火算法、遗传算法、蚁群算法、粒子群算法、禁忌搜索算法等。这篇博客将采用混合整数二阶锥(Mixed-integer second-order cone programming, MISOCP)实现配电网重构问题的求解。
一、算例信息
以IEEE33节点测试系统作为算例，数据来源于matpower工具箱，数据文件如下：
function mpc = IEEE33%% MATPOWER Case Format : Version 2
mpc.version = '2';%%----- Power Flow Data -----%%
%% system MVA base
mpc.baseMVA = 10;%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
mpc.bus = [ %% (Pd and Qd are specified in kW & kVAr here, converted to MW & MVAr below)1 3 0 0 0 0 1 1 0 12.66 1 1 1;2 1 100 60 0 0 1 1 0 12.66 1 1.1 0.9;3 1 90 40 0 0 1 1 0 12.66 1 1.1 0.9;4 1 120 80 0 0 1 1 0 12.66 1 1.1 0.9;5 1 60 30 0 0 1 1 0 12.66 1 1.1 0.9;6 1 60 20 0 0 1 1 0 12.66 1 1.1 0.9;7 1 200 100 0 0 1 1 0 12.66 1 1.1 0.9;8 1 200 100 0 0 1 1 0 12.66 1 1.1 0.9;9 1 60 20 0 0 1 1 0 12.66 1 1.1 0.9;10 1 60 20 0 0 1 1 0 12.66 1 1.1 0.9;11 1 45 30 0 0 1 1 0 12.66 1 1.1 0.9;12 1 60 35 0 0 1 1 0 12.66 1 1.1 0.9;13 1 60 35 0 0 1 1 0 12.66 1 1.1 0.9;14 1 120 80 0 0 1 1 0 12.66 1 1.1 0.9;15 1 60 10 0 0 1 1 0 12.66 1 1.1 0.9;16 1 60 20 0 0 1 1 0 12.66 1 1.1 0.9;17 1 60 20 0 0 1 1 0 12.66 1 1.1 0.9;18 1 90 40 0 0 1 1 0 12.66 1 1.1 0.9;19 1 90 40 0 0 1 1 0 12.66 1 1.1 0.9;20 1 90 40 0 0 1 1 0 12.66 1 1.1 0.9;21 1 90 40 0 0 1 1 0 12.66 1 1.1 0.9;22 1 90 40 0 0 1 1 0 12.66 1 1.1 0.9;23 1 90 50 0 0 1 1 0 12.66 1 1.1 0.9;24 1 420 200 0 0 1 1 0 12.66 1 1.1 0.9;25 1 420 200 0 0 1 1 0 12.66 1 1.1 0.9;26 1 60 25 0 0 1 1 0 12.66 1 1.1 0.9;27 1 60 25 0 0 1 1 0 12.66 1 1.1 0.9;28 1 60 20 0 0 1 1 0 12.66 1 1.1 0.9;29 1 120 70 0 0 1 1 0 12.66 1 1.1 0.9;30 1 200 600 0 0 1 1 0 12.66 1 1.1 0.9;31 1 150 70 0 0 1 1 0 12.66 1 1.1 0.9;32 1 210 100 0 0 1 1 0 12.66 1 1.1 0.9;33 1 60 40 0 0 1 1 0 12.66 1 1.1 0.9;
];%% generator data
% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf
mpc.gen = [1 0 0 10 -10 1 100 1 10 0 0 0 0 0 0 0 0 0 0 0 0;
];%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax
mpc.branch = [ %% (r and x specified in ohms here, converted to p.u. below)1 2 0.0922 0.0470 0 0 0 0 0 0 1 -360 360;2 3 0.4930 0.2511 0 0 0 0 0 0 1 -360 360;3 4 0.3660 0.1864 0 0 0 0 0 0 1 -360 360;4 5 0.3811 0.1941 0 0 0 0 0 0 1 -360 360;5 6 0.8190 0.7070 0 0 0 0 0 0 1 -360 360;6 7 0.1872 0.6188 0 0 0 0 0 0 1 -360 360;7 8 0.7114 0.2351 0 0 0 0 0 0 1 -360 360;8 9 1.0300 0.7400 0 0 0 0 0 0 1 -360 360;9 10 1.0440 0.7400 0 0 0 0 0 0 1 -360 360;10 11 0.1966 0.0650 0 0 0 0 0 0 1 -360 360;11 12 0.3744 0.1238 0 0 0 0 0 0 1 -360 360;12 13 1.4680 1.1550 0 0 0 0 0 0 1 -360 360;13 14 0.5416 0.7129 0 0 0 0 0 0 1 -360 360;14 15 0.5910 0.5260 0 0 0 0 0 0 1 -360 360;15 16 0.7463 0.5450 0 0 0 0 0 0 1 -360 360;16 17 1.2890 1.7210 0 0 0 0 0 0 1 -360 360;17 18 0.7320 0.5740 0 0 0 0 0 0 1 -360 360;2 19 0.1640 0.1565 0 0 0 0 0 0 1 -360 360;19 20 1.5042 1.3554 0 0 0 0 0 0 1 -360 360;20 21 0.4095 0.4784 0 0 0 0 0 0 1 -360 360;21 22 0.7089 0.9373 0 0 0 0 0 0 1 -360 360;3 23 0.4512 0.3083 0 0 0 0 0 0 1 -360 360;23 24 0.8980 0.7091 0 0 0 0 0 0 1 -360 360;24 25 0.8960 0.7011 0 0 0 0 0 0 1 -360 360;6 26 0.2030 0.1034 0 0 0 0 0 0 1 -360 360;26 27 0.2842 0.1447 0 0 0 0 0 0 1 -360 360;27 28 1.0590 0.9337 0 0 0 0 0 0 1 -360 360;28 29 0.8042 0.7006 0 0 0 0 0 0 1 -360 360;29 30 0.5075 0.2585 0 0 0 0 0 0 1 -360 360;30 31 0.9744 0.9630 0 0 0 0 0 0 1 -360 360;31 32 0.3105 0.3619 0 0 0 0 0 0 1 -360 360;32 33 0.3410 0.5302 0 0 0 0 0 0 1 -360 360;21 8 2.0000 2.0000 0 0 0 0 0 0 0 -360 360;9 15 2.0000 2.0000 0 0 0 0 0 0 0 -360 360;12 22 2.0000 2.0000 0 0 0 0 0 0 0 -360 360;18 33 0.5000 0.5000 0 0 0 0 0 0 0 -360 360;25 29 0.5000 0.5000 0 0 0 0 0 0 0 -360 360;
];%%----- OPF Data -----%%
%% generator cost data
% 1 startup shutdown n x1 y1 ... xn yn
% 2 startup shutdown n c(n-1) ... c0
mpc.gencost = [2 0 0 3 0 20 0;
];%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3; %% in Volts
Sbase = mpc.baseMVA * 1e6; %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);%% convert loads from kW to MW
mpc.bus(:, [PD, QD]) = mpc.bus(:, [PD, QD]) / 1e3;二、基于yalmip构建数学模型
YALMIP是一个开源的 MATLAB工具箱，用于建模与求解优化问题，包括线性规划、混合整数线性规划、二次规划、半正定规划、凸优化等。该工具箱提供了高层次的优化接口，支持使用 MATLAB 的熟悉语法来描述优化问题，使得用户能够快速地尝试不同的建模和算法选择，从而便于对优化问题进行建模、求解和验证。此外，YALMIP工具箱还支持多种开源的优化求解器，包括 MOSEK、CPLEX、GLPK、SDPT3等，使得用户在对优化问题求解时具有多种选择。
使用yalmip工具箱求解优化问题的第一步就是需要定义决策变量，配电网重构问题中，决策变量包括支路状态，节点电压、支路有功/无功功率、电源有功/无功出力等。
首先需要设置系统相关的参数：
%% 系统参数
mpc0 = IEEE33;
nb=33; % 节点数
ns=1; % 电源节点数
nl=37; % 支路数
P_load=mpc0.bus(:,3)/mpc0.baseMVA; % 有功负荷
Q_load=mpc0.bus(:,4)/mpc0.baseMVA; % 无功负荷
r_ij=mpc0.branch(:,3); % 线路电阻
x_ij=mpc0.branch(:,4); % 线路电抗
M=1.06*1.06 - 0.94*0.94;
% 电源功率最大值
P_g_max=zeros(nb,1);
P_g_max(1)=10/mpc0.baseMVA;
Q_g_max=zeros(nb,1);
Q_g_max(1)=10/mpc0.baseMVA;
% 电源功率最小值
P_g_min=zeros(nb,1);
Q_g_min=[-10/mpc0.baseMVA;zeros(nb-1,1)];
% 支路功率最大值
Sij_max=15/mpc0.baseMVA;
% 最大电压
Vmax=[1;1.1*ones(32,1)];
Vmin=[1;0.9*ones(32,1)];
% 流入节点的支路
branch_to_node=zeros(nb,nl);
% 流出节点的支路
branch_from_node=zeros(nb,nl);
for k=1:nlbranch_to_node(mpc0.branch(k,2),k)=1;branch_from_node(mpc0.branch(k,1),k)=1;
end随后需要定义决策变量：
%% 优化变量
z_ij=binvar(nl,1);% 支路开断情况
V_i=sdpvar(nb,1);% 节点电压
I_ij=sdpvar(nl,1);% 支路电流
P_ij=sdpvar(nl,1);% 线路有功功率
Q_ij=sdpvar(nl,1);% 线路无功功率
P_g=sdpvar(nb,1);% 电源有功出力
Q_g=sdpvar(nb,1);% 电源无功功率对于优化问题的约束条件，可以直接使用matlab语句添加到yalmip工具箱中，举个例子，如果我们要对一组变量x1, x2和x3进行优化，我们可能会需要满足以下两个约束条件：
1.x1 + x2 <= 10
2.2x2 + 3x3 >= 20
在yalmip中，我们可以使用 <= 或 >= 符号来表示这些限制条件，具体实现方法如下：
% 定义变量
x = sdpvar(3, 1);
% 定义约束条件
Constraints = [x(1) + x(2) <= 10, 2*x(2) + 3*x(3) >= 20];这样，我们就将这两个约束条件表达出来，并存放在Constraints变量中。在之后的优化求解中，yalmip会自动地考虑这些限制条件，并找到满足这些条件的最优解。
三、配电网重构基本原理
本博客将用数学规划的形式对网络重构问题进行描述,但由于严格的数学模型属于非线性非凸规划,理论上是一个NP难题而缺乏有效的求解方法。对此,可以根据需要引人二阶锥凸松弛技术和线性化方法,将原优化模型转化为具有凸可行域的形式以便求解。下面将针对网络重构问题,建立其对应的混合整数二阶锥规划模型。
1.1 目标函数
主动配电网进行网络重构的目标有多种,主要包括:降低网络有功损耗，降低系统一定时间段内能量损耗，使线路负载均衡，提高系统供电可靠性，以及提高电压稳定性等。假设目标函数选择为降低系统网络损耗：
使用matlab代码可表示为：
objective = sum(I_ij.^2*r_ij);或者：
objective =sum((P_ij.^2+Q_ij.^2)/V_i.^2.*r_ij);1.2拓扑约束
为了保护整定和减小短路电流,一般要求配电网呈辐射状运行,即网络中不存在环状拓扑结构。可以表示为：
Matlab代码表示为：
Constraints = [Constraints , sum(z_ij) == nb-ns];1.3电压安全约束
Matlab代码表示为：
Constraints = [Constraints, Vmin <= V_i,V_i <= Vmax];1.4 支路容量约束
配电网的每条线路都具有一定的传输容量,若实际传输电量过大,将使得电线发热量剧增,导致线路损耗增加且损伤输电线。线路容量可用视在功率的平方限值进行描述：
Matlab代码表示为：
Constraints = [Constraints, Sij_max'.^2.*z_ij >= P_ij.^2+Q_ij.^2];1.5 DistFlow潮流约束
对于网络重构而言,基本的 DistFlow潮流约束并不完善,仍存在以下三个问题需要解决。
(1)防止零注入孤立节点的存在。
零注入节点是指本身没有净流出或净注入功率的节点，即P_L和Q_L均为0。
图中浅灰色节点表示零注入节点。重构方案A和B均满足辐射状约束式(5-38)和潮流约束,但方案B中却出现了孤立节点和环状结构并存的情况,显然不符合要求。为了防止这样的零注入孤立节点的存在,将式(5-42)和式(5-43)分别改进为
倘若支路ij处于断开状态,约束条件式(5-41)会将Pij和Qij限制为0，而此时式(5-44)将强制不相连支路两端的电压幅值相等，这是不合理的。对此，需要引入 big-M方法,将式(5-44)转变为如下形式：
使用matlab代码表示为：
m_ij=(1-z_ij)*M;
Constraints = [Constraints, P_g-P_load+branch_to_node*P_ij-branch_to_node*(I_ij.*r_ij)-branch_from_node*P_ij == 0];
Constraints = [Constraints, Q_g-Q_load+branch_to_node*Q_ij-branch_to_node*(I_ij.*x_ij)-branch_from_node*Q_ij == 0];
Constraints = [Constraints,U_i(mpc0.branch(:,1))-U_i(mpc0.branch(:,2))<= m_ij + 2*r_ij.*P_ij + 2*x_ij.*Q_ij - ((r_ij.^2 + x_ij.^2)).*I_ij];
Constraints = [Constraints,U_i(mpc0.branch(:,1))-U_i(mpc0.branch(:,2))>= -m_ij + 2*r_ij.*P_ij + 2*x_ij.*Q_ij - ((r_ij.^2 + x_ij.^2)).*I_ij];
Constraints = [Constraints, U_i(mpc0.branch(:,1),:).*I_ij == P_ij.^2+Q_ij.^2];1.6 非凸约束的二阶锥松弛
由于式(5-45)、式(5-46)和式(5-47)的强非凸形式,该网络重构模型是一个非凸规划,属于NP难题而不易求解。为此,需要引入二阶锥松弛技术,将原模型转化为一个MISOCP问题。
在模型中添加式(5-48)和式(5-49)这两个新约束,用新变量替代原目标函数和约束条件中的对应项,并将支路容量约束式(5-41)用电流幅值约束式(5-50)替代:
Matlab代码：
U_i=sdpvar(nb,1);% 电压的平方
L_ij=sdpvar(nl,1);% 电流的平方对于式(5-49)的非凸形式,在满足目标函数是L_ij的严格增函数及节点负荷无上界等条件下,可作如下变形:
在yalmip中，标准二阶锥可以使用cone函数进行表示：
for k=1:nlConstraints = [Constraints, cone([2*P_ij(k) 2*Q_ij(k) L_ij(k)-U_i(mpc0.branch(k,1))],L_ij(k)+U_i(mpc0.branch(k,1)))];
end最终,网络重构的MISOCP模型如下所示。
四、运行结果
五、完整matlab代码
这个链接可以获取完整matlab代码：
基于混合整数二阶锥(MISOCP)的配电网重构的matlab代码
另外，简单提一下，求解MISOCP问题最好的求解器是mosek，如果需要做MISOCP问题，最好还是用mosek求解器，求解性能远远超过cplex和gurobi。
基于混合整数二阶锥(MISOCP)的配电网重构(附matlab代码)相关推荐
- 【配电网重构】基于yalmip求解含sop+二阶锥配电网重构附matlab代码
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
- 【配电网重构】基于粒子群算法实现最小化功率损耗的配电网重构附matlab代码
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
- 【配电网】基于遗传算法实现三相单目标配电网重构附matlab代码
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
- 【配电网重构】基于遗传算法实现配电网重构附matlab代码
✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信.
- 基于改进二进制粒子群算法的配电网重构（matlab实现）
目录 一.引言 1.问题背景 2.二进制粒子群算法 2.1简介 2.2 S i g m o i d Sigmoid Sigmoid函数 2.3二进制粒子群算法 2.4算法的改进 二.配电网重构模型 1 ...
- 混合整数二阶锥规划公式中具有不同动力学的定价惯性和频率响应
混合整数二阶锥规划公式中具有不同动力学的定价惯性和频率响应 测试环境:MATLAB yalmip 关键词:凸优化,频率响应,边际定价,可再生能源 能源和频率服务的联合市场清算,包括定价方法. 频率服务 ...
- 【图像分割】基于计算机视觉实现视网膜图像中的血管分割附matlab代码
1 简介 视网膜图像里的血管是可以被观察到的一类微血管,并且它是无创伤的,而其分布位置也属于深度部位[5].其分布.结构和形态特征的变化能在一定程度上反映病变的程度.而白血病.糖尿病以及高血压等疾病都 ...
- 【图像增强】基于DEHAZENET和HWD的水下去散射图像增强附matlab代码
1 内容介绍 去散射和边缘增强是解决水下图像的对比度严重衰减.颜色偏差和边缘模糊等问题的关键步骤.这篇论文提出了一种较好的水下图像增强的方法.首先使用经过端到端训练的卷积神经网络去测量输入图片,同时以 ...
- 【预测模型-ELM分类】基于鲸鱼算法优化核极限学习机实现数据分类附matlab代码
1 内容介绍 极限学习机(extreme learning machine,ELM)作为一种新兴的机器学习方法,已经成为了一个热门的研究方向. ELM 随机确定单隐含层网络的输入权值和隐含层节点偏置, ...
最新文章
- JSBridge的思考
- java面试基础_Java面试基础
- 【ARM】ARM其它指令
- 8 Regular Expressions You Should Know
- 前端电子表数字字体_爬虫：如何优雅应对字体反爬
- java part part,java泛型指导手册（part1），指导手册part1
- Oracle的重做日志
- WSL2扩展虚拟硬件磁盘( VHD大于256GB报错问题)
- TASKCTL-函数表达式分类
- cmos逻辑门传输延迟时间_组合逻辑电路详解、实现及其应用
- 2019数据安装勾选_勾选方式认证发票你不会？整理了最全的操作流程，快来看看吧...
- java 基础的数组 添加
- 网络中数据如何保证数据的安全性？
- Python: numpy tile()函数 可实现ndarray的横向纵向复制
- Swift语言入门学习
- 游戏BI数据统计分析相关
- android 自定义桌面,无需 root，3 个 App 帮你构建高效的 Android 桌面
- c语言双精度小数点后取几位_c语言float类型小数点后位数
- ios安装并信任证书
- RARP协议的使用场景