基于二阶锥规划（SOCP）松弛和线性离流的配电网规划（DNP）方法示例（Matlab代码实现）

    目录

💥1 概述

📚2 运行结果

🎉3 参考文献

👨‍💻4 Matlab代码

💥1 概述

配电网最优潮流 Optimal Power Flow, OPF) 问题是指在满足一定约束条件的情况下，通过控制配电网中的可控变量，使配电网达到优化运行的目的。由于OPF问题约束条件的特点，导致其为难以求解的非凸规划问题。目前OPF求解方法主要分为经典数学规划算法和智能优化算法两种。

近年来，很多学者不断探索高效求解OPF 问题的方法，随着研究的不断深入，二阶锥松弛（Second Order Cone Relaxation, SOCR）技术被逐步运用于求解OPF问题。有学者建立了以支路潮流计算为基础的OPF模型，针对OPF中的非凸性约束，采用SOCR技术将其松弛为二阶锥约束，整个 OPF模型则被转化为二阶锥规划（Second Order Cone Programming，SOCP）问题，对其求解可以得到全局最优解。还有学者在主动配电网最优潮流计算中采取SOCR技术处理非凸性约束，将优化模型转化为SOCP问题，得到了很好的求解效果，并对产生的松弛间误差进行分析，结果表明松弛误差满足计算准确度。

​本文包括两个基于二阶锥规划（SOCP）松弛和线性离流的配电网规划（DNP）方法示例。目标是最大限度地降低配电线路的投资成本和唯一的运营成本，即负载损失值（VOLL）。

建模由YALMIP完成，YALMIP是MATLAB中用于进行优化建模的工具箱，并由GUROBI解决。因此，您应该在PC中安装MATLAB，YALMIP和GUROBI，以使代码正常工作。

📚2 运行结果

主函数部分代码：

%% With/without s \in S_volt
WithSvolt=0; % 1 / 0 means with / without s\in S_volt
%% System Parameters
Ubase=12.35e3; % unit:V
Sbase=1e7;   % unit:VA
Ibase=Sbase/sqrt(3)/Ubase; % unit:A
Zbase=Ubase/Ibase/sqrt(3);  % unit: 
netpara=xlsread('SCE47','网络参数');
loadpoint=xlsread('SCE47','节点负荷');
L=size(netpara,1); % number of lines
r=netpara(:,4)/Zbase;% resistance, unit: p.u.
x=netpara(:,5)/Zbase;% reactance, unit: p.u.
So=loadpoint(:,3)*1e6/Sbase;% node power unit: MVA
judge=loadpoint(:,4); % node type
I_max=560.98/Ibase; % Ubase=12.35e3, assuming P_max is 12MW for each line, then I_max=12e6/Ubase/sqrt(3)= 560.98 A,
v_max=1.1^2;
v_min=0.9^2;
​
%% Plot the distribution network
I=netpara(:,2);
J=netpara(:,3);
G=graph(I,J);
h1=figure;
p=plot(G);
highlight(p,find(judge==1),'Marker','s','NodeColor','c','Markersize',10);  % Capacitator
highlight(p,find(judge==0),'Marker','v','NodeColor','y','Markersize',10);  % PV panels
In=myincidence(I,J); % node-branch incidence matrix
Inn=In;
Inn(Inn>0)=0;    % Inn is the negative part of I, showing the lines starting from nodes
%% Number of nodes
N=47;
%m=4;
%% Decision variables
P_ij=sdpvar(L,1); % Pij=sdpvar(numnodes,numnodes,'full');% Active power from node i to node j
Q_ij=sdpvar(L,1); % Qij=sdpvar(numnodes,numnodes,'full');% Reactive power from node i to node j
u=sdpvar(N,1); % Voltage
v=sdpvar(N,1); % u^2
l_ij=sdpvar(L,1) ;% Currents' magnitudes' square, I^2
Pi=sdpvar(N,1); % Active power injection of nodes
Qi=sdpvar(N,1); % Reactive power injection of nodes
​
%% Resistantce vector
w=sqrt(-1);
z=r+x*w;
z_c=conj(z); %conjugate of resistance
​
Cons=[];
​
%% Node types are different
Cons_Load=[];
Eta=2;  % Scale up the PV penetration rate, 5 means 5 times than Table I in [1]
Pi_max=zeros(N,1);
Qi_max=zeros(N,1);
for i=1:N
    if judge(i)==1     % Capacitator
        Cons_Load=[Cons_Load,Pi(i)==0];
        Qi_max(i)=So(i)*Eta;
        Cons_Load=[Cons_Load,0<=Qi(i)<=Qi_max(i)];   %st=[pl(i)==0,0<=ql(i)<=So(i)];
    elseif judge(i)==2  % Load node
        Pi_max(i)=-So(i)*0.9;
        Qi_max(i)=-So(i)*sqrt(1-0.9^2);
        Cons_Load=[Cons_Load,Pi(i)==Pi_max(i)];  % cos(0.9)???
        Cons_Load=[Cons_Load,Qi(i)==Qi_max(i)];  % st=[pl(i)==-So(i)*cos(0.9);ql(i)==-So(i)*sin(0.9)];
    elseif judge(i)==0  %PV panel
        Pi_max(i)=So(i)*Eta;
        Cons_Load=[Cons_Load,0<=Pi(i)<=Pi_max(i)];
        Cons_Load=[Cons_Load,0==Qi(i)];
    elseif judge(i)==3  %slack node (substation node)
        Cons_Load=[Cons_Load,-So(i)<=Pi(i)<=So(i)];
        Cons_Load=[Cons_Load,-So(i)*0.5<=Qi(i)<=So(i)*0.5];
    end
end
Cons=[Cons,Cons_Load];
display('Constraints on load type of nodes completed!')
% Cons_Load
​
%% Check whether C1 holds for the network
De = degree(G); % degree of each node
LeafNodes=find(De==1);
LeafNodes=LeafNodes(2:end); % delete node 1
Pij_hat=abs(inv(In(2:end,:))*Pi_max(2:end));
Qij_hat=abs(inv(In(2:end,:))*Qi_max(2:end));
IJ=[I J];
A_l=cell(length(LeafNodes),1);
Path_L=cell(length(LeafNodes),1);
Al_uij=zeros(length(LeafNodes),2);
for i=1:length(LeafNodes)
    lt=LeafNodes(i);
    Path = shortestpath(G,lt,1);
    A_l{i}=zeros(2,2,length(Path)-1);
    for j=1:length(Path)-1 % for all node in Path from lt1, calculate A_l
%         nl=[]; % the set of 1,2,...,nl in C1
        for k=1:L
            if (IJ(k,1)==Path(j) && IJ(k,2)==Path(j+1))||(IJ(k,2)==Path(j) && IJ(k,1)==Path(j+1))
                Path_L{i}=[Path_L{i},Path(j+1)];
                A_l{i}(:,:,j)=diag([1 1])-2/v_min*mtimes([r(k);x(k)],[Pij_hat(k),Qij_hat(k)]);
            end
        end
    end
​

🎉3 参考文献

[1]范金月, 白晓清. 基于扩展二次锥规划的最优潮流模型研究[J]. 电网与清洁能源, 2014(3):7.​

部分理论引用网络文献，若有侵权联系博主删除。