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

简介: 配电网重构是指在满足配电网运行基本约束的前提下,通过改变配电网中一个或多个开关的状态对配电网中一个或多个指标进行优化。通过配电网重构,可以在不增加设备投资的情况下,充分发挥配电系统的潜力,提高系统的性能指标,具有较好的经济效益。配电网重构的算法有许多,包括数学规划方法,如分支定界法、0-1 整数规划、单纯形法等;启发式算法,如最优流算法、开关交换法等;智能优化算法,如模拟退火算法、遗传算法、蚁群算法、粒子群算法、禁忌搜索算法等。这篇博客将采用混合整数二阶锥(Mixed-integer second-order cone programming, MISOCP)实现配电网重构问题的求解。

 参考资料:主动配电网网络分析与运行调控 (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;

image.gif

二、基于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:nl
    branch_to_node(mpc0.branch(k,2),k)=1;
    branch_from_node(mpc0.branch(k,1),k)=1;
end

image.gif

       随后需要定义决策变量:

%% 优化变量
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);% 电源无功功率

image.gif

       对于优化问题的约束条件,可以直接使用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];

image.gif

这样,我们就将这两个约束条件表达出来,并存放在Constraints变量中。在之后的优化求解中,yalmip会自动地考虑这些限制条件,并找到满足这些条件的最优解。

三、配电网重构基本原理

       本博客将用数学规划的形式对网络重构问题进行描述,但由于严格的数学模型属于非线性非凸规划,理论上是一个NP难题而缺乏有效的求解方法。对此,可以根据需要引人二阶锥凸松弛技术和线性化方法,将原优化模型转化为具有凸可行域的形式以便求解。下面将针对网络重构问题,建立其对应的混合整数二阶锥规划模型。

1.1 目标函数

       主动配电网进行网络重构的目标有多种,主要包括:降低网络有功损耗,降低系统一定时间段内能量损耗,使线路负载均衡,提高系统供电可靠性,以及提高电压稳定性等。假设目标函数选择为降低系统网络损耗:

image.gif

       使用matlab代码可表示为:

objective = sum(I_ij.^2*r_ij);

image.gif

       或者:

objective =sum((P_ij.^2+Q_ij.^2)/V_i.^2.*r_ij);

image.gif

1.2拓扑约束

       为了保护整定和减小短路电流,一般要求配电网呈辐射状运行,即网络中不存在环状拓扑结构。可以表示为:

image.gif

       Matlab代码表示为:

Constraints = [Constraints , sum(z_ij) == nb-ns];

image.gif

1.3电压安全约束

image.gif

       Matlab代码表示为:

Constraints = [Constraints, Vmin <= V_i,V_i <= Vmax];

image.gif

1.4 支路容量约束

       配电网的每条线路都具有一定的传输容量,若实际传输电量过大,将使得电线发热量剧增,导致线路损耗增加且损伤输电线。线路容量可用视在功率的平方限值进行描述:

image.gif

       Matlab代码表示为:

Constraints = [Constraints, Sij_max'.^2.*z_ij >= P_ij.^2+Q_ij.^2];

image.gif

1.5 DistFlow潮流约束

image.gif

       对于网络重构而言,基本的 DistFlow潮流约束并不完善,仍存在以下三个问题需要解决。

       (1)防止零注入孤立节点的存在。

       零注入节点是指本身没有净流出或净注入功率的节点,即P_L和Q_L均为0。

image.gif

       图中浅灰色节点表示零注入节点。重构方案A和B均满足辐射状约束式(5-38)和潮流约束,但方案B中却出现了孤立节点和环状结构并存的情况,显然不符合要求。为了防止这样的零注入孤立节点的存在,将式(5-42)和式(5-43)分别改进为

image.gif

       倘若支路ij处于断开状态,约束条件式(5-41)会将Pij和Qij限制为0而此时式(5-44)将强制不相连支路两端的电压幅值相等这是不合理的。对此需要引入 big-M方法,将式(5-44)转变为如下形式

image.gif

       使用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];

image.gif

1.6 非凸约束的二阶锥松弛

       由于式(5-45)、式(5-46)和式(5-47)的强非凸形式,该网络重构模型是一个非凸规划,属于NP难题而不易求解。为此,需要引入二阶锥松弛技术,将原模型转化为一个MISOCP问题。

image.gif

       在模型中添加式(5-48)和式(5-49)这两个新约束,用新变量替代原目标函数和约束条件中的对应项,并将支路容量约束式(5-41)用电流幅值约束式(5-50)替代:

image.gif

       Matlab代码:

U_i=sdpvar(nb,1);% 电压的平方
L_ij=sdpvar(nl,1);% 电流的平方

image.gif

       对于式(5-49)的非凸形式,在满足目标函数是L_ij的严格增函数及节点负荷无上界等条件下,可作如下变形:

image.gif

在yalmip中,标准二阶锥可以使用cone函数进行表示:

for k=1:nl
    Constraints = [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

image.gif

       最终,网络重构的MISOCP模型如下所示。

image.gif编辑image.gif

四、运行结果

image.gif

image.gif

image.gif



相关文章
|
1月前
|
算法
基于粒子群算法的分布式电源配电网重构优化matlab仿真
本研究利用粒子群算法(PSO)优化分布式电源配电网重构,通过Matlab仿真验证优化效果,对比重构前后的节点电压、网损、负荷均衡度、电压偏离及线路传输功率,并记录开关状态变化。PSO算法通过迭代更新粒子位置寻找最优解,旨在最小化网络损耗并提升供电可靠性。仿真结果显示优化后各项指标均有显著改善。
|
3月前
|
安全
【2023高教社杯】D题 圈养湖羊的空间利用率 问题分析、数学模型及MATLAB代码
本文介绍了2023年高教社杯数学建模竞赛D题的圈养湖羊空间利用率问题,包括问题分析、数学模型建立和MATLAB代码实现,旨在优化养殖场的生产计划和空间利用效率。
199 6
【2023高教社杯】D题 圈养湖羊的空间利用率 问题分析、数学模型及MATLAB代码
|
3月前
|
存储 算法 搜索推荐
【2022年华为杯数学建模】B题 方形件组批优化问题 方案及MATLAB代码实现
本文提供了2022年华为杯数学建模竞赛B题的详细方案和MATLAB代码实现,包括方形件组批优化问题和排样优化问题,以及相关数学模型的建立和求解方法。
128 3
【2022年华为杯数学建模】B题 方形件组批优化问题 方案及MATLAB代码实现
|
3月前
|
数据采集 存储 移动开发
【2023五一杯数学建模】 B题 快递需求分析问题 建模方案及MATLAB实现代码
本文介绍了2023年五一杯数学建模竞赛B题的解题方法,详细阐述了如何通过数学建模和MATLAB编程来分析快递需求、预测运输数量、优化运输成本,并估计固定和非固定需求,提供了完整的建模方案和代码实现。
90 0
【2023五一杯数学建模】 B题 快递需求分析问题 建模方案及MATLAB实现代码
|
5月前
|
算法 安全 数据库
基于结点电压法的配电网状态估计算法matlab仿真
**摘要** 该程序实现了基于结点电压法的配电网状态估计算法,旨在提升数据的准确性和可靠性。在MATLAB2022a中运行,显示了状态估计过程中的电压和相位估计值,以及误差随迭代变化的图表。算法通过迭代计算雅可比矩阵,结合基尔霍夫定律解决线性方程组,估算网络节点电压。状态估计过程中应用了高斯-牛顿或莱文贝格-马夸尔特法,处理量测数据并考虑约束条件,以提高估计精度。程序结果以图形形式展示电压幅值和角度估计的比较,以及估计误差的演变,体现了算法在处理配电网状态估计问题的有效性。
|
6月前
|
数据安全/隐私保护
耐震时程曲线,matlab代码,自定义反应谱与地震波,优化源代码,地震波耐震时程曲线
地震波格式转换、时程转换、峰值调整、规范反应谱、计算反应谱、计算持时、生成人工波、时频域转换、数据滤波、基线校正、Arias截波、傅里叶变换、耐震时程曲线、脉冲波合成与提取、三联反应谱、地震动参数、延性反应谱、地震波缩尺、功率谱密度
基于混合整数规划的微网储能电池容量规划(matlab代码)
基于混合整数规划的微网储能电池容量规划(matlab代码)
|
6月前
|
算法 调度
含多微网租赁共享储能的配电网博弈优化调度(含matlab代码)
含多微网租赁共享储能的配电网博弈优化调度(含matlab代码)
|
6月前
|
Serverless
基于Logistic函数的负荷需求响应(matlab代码)
基于Logistic函数的负荷需求响应(matlab代码)
|
6月前
|
供应链 算法
基于分布式优化的多产消者非合作博弈能量共享(Matlab代码)
基于分布式优化的多产消者非合作博弈能量共享(Matlab代码)