基于Nonconvex规划的配电网重构研究（Matlab代码实现）

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

💥1 概述

本文基于Nonconvex规划的配电网重构研究。并用Matlab代码实现之。

基于Nonconvex规划的配电网重构研究是针对配电网优化问题的一种方法。传统的配电网通常是基于线性或者凸规划进行设计和运行，但是实际配电网系统的复杂性往往导致非线性和非凸问题的出现。为解决这些问题，基于Nonconvex规划的方法被提出来更好地优化配电网系统。

配电网重构是指通过变换网络拓扑结构和配置设备参数，以改善配电网的性能和可靠性。基于Nonconvex规划的配电网重构研究通常包括以下几个方面：

1. 非线性建模：将配电网系统建模为非线性的数学模型。这包括考虑网络拓扑结构、设备参数、电流限制、电压限制等的非线性方程和约束条件。

2. 问题定义：明确定义配电网重构的目标，例如最小化损耗、提高电压稳定性、降低网络阻塞等。同时，考虑到配电网的约束条件，例如设备的额定容量、电压限制、工作模式等。

3. Nonconvex规划建模：将配电网重构问题转化为Nonconvex规划问题。这可能涉及到非线性约束和非凸目标函数，并且由于配电网的复杂性，问题可能具有多个局部最优解。

4. 优化算法：针对Nonconvex规划问题，需要选择适当的优化算法来求解最优解。常见的算法包括非线性规划算法、启发式算法（如遗传算法、粒子群算法等）以及近似方法等。这些算法可以搜索到全局或者局部最优解。

5. 结果分析和评估：根据求解得到的最优解，分析和评估配电网的性能指标。这可能包括网络损耗、电压稳定性、负荷均衡等方面的评估。

6. 重构方案实施：根据优化结果，制定并实施配电网重构方案。这可能涉及到改变配电网的拓扑结构、设备配置、控制策略等。

基于Nonconvex规划的配电网重构研究能够更好地应对实际配电网系统的复杂性，帮助提高能源利用效率，降低电力系统的运行成本，并提高系统的可靠性和稳定性。然而，由于Nonconvex规划问题的复杂性，求解过程可能较为困难，需要综合考虑求解效率和解的质量。

📚2 运行结果

14147.3秒得到了全局最优，网损为1.7430，AA=[1;1;1;1;1;1;0;1;0;1;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;0]

部分代码：

%SB=100MVA,UB=12.66kV,IEEE-33 Bus,Distflow SOC-ACOPF
%考虑重构，全天拓扑不变
clear
clc
tic
%网络数据，标幺值
Pload=。。。
Line(:,3)=Line(:,3)*100/(12.66^2);
r=real(Line(:,3));
x=imag(Line(:,3));
father=zeros(33,37);son=zeros(33,37);
for i=1:32
father(i,i)=1;
end
father(20,33)=1;father(14,34)=1;father(21,35)=1;father(32,36)=1;father(28,37)=1;
for i=[1:16,18:20,22:23,25:31]
son(i,i+1)=1;
end
son(1,18)=1;son(2,22)=1;son(5,25)=1;son(33,1)=1;son(7,33)=1;son(8,34)=1;son(11,35)=1;son(17,36)=1;son(24,37)=1;
Umax=[1.07*1.07*ones(32,24);1.05*1.05*ones(1,24)];
Umin=[0.93*0.93*ones(32,24);1.05*1.05*ones(1,24)];
Pgmax=[zeros(32,24);ones(1,24)];
Qgmax=[zeros(32,24);ones(1,24)];
%定义变量
P=sdpvar(37,24);%线路有功
Q=sdpvar(37,24);%线路无功
U=sdpvar(33,24);%电压的平方
I=sdpvar(37,24);%电流的平方
Pg=[zeros(32,24);sdpvar(1,24)];%发电机有功
Qg=[zeros(32,24);sdpvar(1,24)];%发电机无功
AA=binvar(37,1);%网架结构
A0=[ones(32,1);zeros(5,1)];%初始拓扑
assign(AA,A0);
Pin=-father*P+father*(I.*(r*ones(1,24)))+son*P;%节点注入有功
Qin=-father*Q+father*(I.*(x*ones(1,24)))+son*Q;%节点注入无功
P_tree=sdpvar(37,1);%虚拟有功
Pin_tree=-father*P_tree+son*P_tree;%虚拟节点注入有功
Ploss_total=sum(sum(I.*(r*ones(1,24))));%目标函数，网损最小
%约束条件
C1=[sum(AA)==32,U>=Umin,U<=Umax,Pg>=-Pgmax,Pg<=Pgmax,Qg>=-Qgmax,Qg<=Qgmax,I>=0,I<=0.11*AA*ones(1,24),-AA<=P_tree<=AA,-0.11*AA*ones(1,24)<=P<=0.11*AA*ones(1,24),-0.11*AA*ones(1,24)<=Q<=0.11*AA*ones(1,24)];%边界约束
C2=[Pin+Pload-Pg==0,Pin_tree(1:32)+0.01==0];%有功KCL约束
C3=[Qin+Qload-Qg==0];%无功KCL约束
C4=[-(1.07*1.07-0.93*0.93)*(1-AA)*ones(1,24)<=-U(Line(:,2),:)+U(Line(:,1),:)-2*(r*ones(1,24)).*P-2*(x*ones(1,24)).*Q+((r.^2+x.^2)*ones(1,24)).*I<=(1.07*1.07-0.93*0.93)*(1-AA)*ones(1,24)];%电压降落约束
C=[C1,C2,C3,C4];
toc%建模时间
ops=sdpsettings('solver','gurobi','usex0',1);
[model,recoverymodel,diagnostic,internalmodel] = export(C,Ploss_total,ops);%得到除去P^2+Q^2=UI的约束
params.Nonconvex=2;%启动gurobi非线性求解器
params.FeasibilityTol=1e-9;%由于gurobi采取的是双层模型，因此可行性步长应尽可能小
params.IntFeasTol=1e-9;%由于gurobi对非线性模型采用的是外嵌分支定界算法，相当于求解MIP问题，因此整数可行性要足够精确
params.Threads=8;%并行计算，8线程
L=length(model.obj);%决策变量数
%下面定义P^2+Q^2=UI的约束,模型为sum(Qval*x(Qrow)*x(Qcol))+q*x=rhs
for t=1:24
for j=1:37
model.quadcon(37*t-37+j).Qrow=[37*t-37+j,37*t-37+j+37*24,33*t-33+Line(j,1)+37*24*2];%P,Q,Uj
model.quadcon(37*t-37+j).Qcol=[37*t-37+j,37*t-37+j+37*24,37*t-37+j+37*24*2+33*24];%P,Q,I
model.quadcon(37*t-37+j).Qval=[1,1,-1];%P^2+Q^2-UI
model.quadcon(37*t-37+j).q=sparse(L,1);
model.quadcon(37*t-37+j).rhs=0;
model.quadcon(37*t-37+j).sense='=';%严格等号

%定义变量
P=sdpvar(37,1);%线路有功
Q=sdpvar(37,1);%线路无功
U=sdpvar(33,1);%电压的平方
I=sdpvar(37,1);%电流的平方
Pg=[zeros(32,1);sdpvar];%发电机有功
Qg=[zeros(32,1);sdpvar];%发电机无功
AA=binvar(37,1);%网架结构
A0=[ones(32,1);zeros(5,1)];%初始拓扑
assign(AA,A0);
Pin=-father*P+father*(I.*r)+son*P;%节点注入有功
Qin=-father*Q+father*(I.*x)+son*Q;%节点注入无功
P_tree=sdpvar(37,1);%虚拟有功
Pin_tree=-father*P_tree+son*P_tree;%虚拟节点注入有功
Ploss_total=sum(I.*r);%目标函数，网损最小
%约束条件
C1=[sum(AA)==32,U>=Umin,U<=Umax,Pg>=-Pgmax,Pg<=Pgmax,Qg>=-Qgmax,Qg<=Qgmax,I>=0,I<=0.11*AA,-AA<=P_tree<=AA,-0.11*AA<=P<=0.11*AA,-0.11*AA<=Q<=0.11*AA];%边界约束
C2=[Pin+Pload-Pg==0,Pin_tree(1:32)+0.01==0];%有功KCL约束
C3=[Qin+Qload-Qg==0];%无功KCL约束
C4=[-(1.07*1.07-0.93*0.93)*(1-AA)<=-U(Line(:,2),:)+U(Line(:,1),:)-2*r.*P-2*x.*Q+(r.^2+x.^2).*I<=(1.07*1.07-0.93*0.93)*(1-AA)];%电压降落约束
C=[C1,C2,C3,C4];
toc%建模时间
ops=sdpsettings('solver','gurobi','usex0',1);
[model,recoverymodel,diagnostic,internalmodel] = export(C,Ploss_total,ops);%得到除去P^2+Q^2=UI的约束
params.Nonconvex=2;%启动gurobi非线性求解器
params.FeasibilityTol=1e-9;%由于gurobi采取的是双层模型，因此可行性步长应尽可能小
params.IntFeasTol=1e-9;%由于gurobi对非线性模型采用的是外嵌分支定界算法，相当于求解MIP问题，因此整数可行性要足够精确
params.Threads=8;%并行计算，8线程
L=length(model.obj);%决策变量数
%下面定义P^2+Q^2=UI的约束,模型为sum(Qval*x(Qrow)*x(Qcol))+q*x=rhs