径向配电网潮流的求解(Matlab代码实现)

目录

1 概述

2 Matlab代码实现

3 知识点

---

1 概述

潮流分析对于配电网的设计有着重要的作用。一个高效、良好的配电系统潮流汇聚不仅有利于获得网络的电压和功率损耗，而且对于分支导体的准确选择和规划的其他方面也是必要的。像 Newton-Raphson 这样的潮流方法，以及 快速解耦潮流方法可以有效地用于传输系统。但是对于具有高 R/X 比的系统，上述方法未能收敛。一些研究人员 试图修改传统的方法来解决配电网络。 使用了梯形网络理论。 Shirmohammadi 等人提出了一种基于补偿的方法来解决弱网格网络。分支编号方案的开发是为了提高数值性能。 Baran 和 Wu 提出了一种基于 Newton-Raphson 的方法，该方法需要计算 Jacobian 矩阵来求解径向分布网络。所以这种方法需要更多的计算工作。 Renato 使用双二次方程来关联发送端和接收端电压。但是他的方法没有考虑电压的角度。 Goswami 和 Basu  提出了一种解决径向网络的直接方法，但他们限制了没有节点可以是超过三个分支的连接点。 Jasmon 和 Lee 将网络简化为单线等效模型，并使用了涉及发送端电压的潮流方程。达斯等人。 [11] 为节点和分支提出了一种独特的编号方案。他们使用简单的代数方程来求解径向分布系统。 Haque [12] 提出了径向和网状配电网络的解决方法。后向和前向扫描技术用于最终解决方案。最初将变电站电压作为每个节点的电压，在每次连续迭代后更新。 Ghosh 和 Das 提出了一种使用正向反向技术解决潮流的方法，但只计算了每个节点的电压幅度。阿拉文达巴布等人。

一种新的有效的径向配电网潮流求解方法。简单的超越方程用于关联配电系统各支路的发送端电压、接收端电压和电压降。线路充电电容的影响已纳入潮流解决方案。一种计算机算法的开发方式使得不需要采用任何顺序节点编号方案来解决网络问题。接收端电压的角度也与电压的大小一起计算。它是一种迭代方法。考虑从变电站开始到每个终端节点的平坦电压（1pu）。在每次连续迭代后更新电压幅度和角度，然后使用新获得的电压幅度和角度值计算电压降。该方法与其他方法的速度和内存要求的比较已得到验证，以显示其效率。

2 Matlab代码实现

算例给出的是IEEE33节点和IEEE69节点。

clc;
clear all;
tic;
global linedata loaddata
%% IEEE33数据

%         No   Nd    Nd    电阻 电抗
linedata =[1    1      2     0.0922    0.0470   
          2    2      3     0.4930    0.2511    
          3    3      4     0.3660    0.1864  
          4    4      5     0.3811    0.1941    
          4    5      6     0.8190    0.7070    
          6    6      7     0.1872    0.6188   
          7    7      8     0.7114    0.2351 
          8    8      9     1.0300    0.7400   
          9    9      10    1.0440    0.7400   
         10    10     11    0.1966    0.0650   
         11    11     12    0.3744    0.1238   
         12    12     13    1.4680    1.1550   
         13    13     14    0.5416    0.7129   
         14    14     15    0.5910    0.5260    
         15    15     16    0.7463    0.5450   
         16    16     17    1.2890    1.7210   
         17    17     18    0.7320    0.5740  
         18     2     19    0.1640    0.1565    
         19    19     20    1.5042    1.3554  
         20    20     21    0.4095    0.4784   
         21    21     22    0.7089    0.9373    
         22     3     23    0.4512    0.3083    
         23    23     24    0.8980    0.7091  
         24    24     25    0.8960    0.7011   
         25     6     26    0.2030    0.1034   
         26    26     27    0.2842    0.1447    
         27    27     28    1.0590    0.9337    
         28    28     29    0.8042    0.7006   
         29    29     30    0.5075    0.2585  
         30    30     31    0.9744    0.9630   
         31    31     32    0.3105    0.3619   
         32    32     33    0.3410    0.5302];
     
%% IEEE69节点数据
%节点    -----负荷-----    
%No.     有功功率(kW)   无功功率(kVAr)  
loaddata=[  1	0         0  	  0
            2	100       60 	  0
            3	90        40      0  	
            4	120       80	  0  
            5	60        30	  0  
            6	60        20      0 
            7	200       100	  0
            8	200       100	  0 
            9	60        20      0	
            10	60        20 	  0
            11	45        30 	  0
            12	60        35      0 
            13	60        35      0 
            14	120       80      0 
            15	60        10      0 
            16	60        20      0 
            17	60        20      0 
            18	90        40      0 
            19	90        40      0 
            20	90        40      0 
            21	90        40      0
            22	90        40      0 
            23	90        50      0
            24	420       200	  0
            25	420       200	  0
            26	60        25      0  
            27	60        25      0 
            28	60        20      0
            29	120       70      0
            30	200       600	  0
            31	150       70      0  
            32	210       100	  0
            33	60        40      0];

nl=linedata(:,2);% 左节点
nr=linedata(:,3);% 右节点
br=length(linedata(:,1));% 支路数
no=length(loaddata(:,1));% 节点数
MVAb=100;%基准功率MVA
KVb=11;%基准电压 kV
Zb=(KVb^2)/MVAb;% 基准阻抗

%% 标幺值
for i=1:br
    R(i,1)=(linedata(i,4))/Zb;
    X(i,1)=(linedata(i,5))/Zb;
end
for i=1:no
    P(i,1)=((loaddata(i,2))/(1000*MVAb));
    Q(i,1)=((loaddata(i,3))/(1000*MVAb));
    Qsh(i,1)=((loaddata(i,4))/(1000*MVAb));   
end
    
R;% 线路电阻
X;%线路电抗
j=sqrt(-1);
Z=R+j*X;%线路阻抗
P;% 负荷有功功率
Q;% 负荷无功功率
Pt=sum(P);%总有功功率
Qt=sum(Q);% 总无功功率
%Pt
%Qt
C=zeros(br,no);%零矩阵的初始化
for i=1:br
    a=linedata(i,2);
    b=linedata(i,3);
    for j=1:no
        if a==j
            C(i,j)=-1;
        end
        if b==j
            C(i,j)=1;
        end
    end
end
C;
e=1;
for i=1:no
    d=0;
    for j=1:br
        if C(j,i)==-1
            d=1;
        end
    end
    if d==0
        endnode(e,1)=i;
        e=e+1;
    end
end
endnode;
h=length(endnode);
for j=1:h
    e=2;
    
    f=endnode(j,1);
   % while (f~=1)
   for s=1:no
     if (f~=1)
       k=1;  
       for i=1:br
           if ((C(i,f)==1)&&(k==1))
                f=i;
                k=2;
           end
       end
       k=1;
       for i=1:no
           if ((C(f,i)==-1)&&(k==1));
                f=i;
                g(j,e)=i;
                e=e+1;
                k=3;
           end            
       end
     end
   end
end
for i=1:h
    g(i,1)=endnode(i,1);
end
g;
w=length(g(1,:));
for i=1:h
    j=1;
    for k=1:no 
        for t=1:w
            if g(i,t)==k
                g(i,t)=g(i,j);
                g(i,j)=k;
                j=j+1;
             end
         end
    end
end
g;
for k=1:br
    e=1;
    for i=1:h
        for j=1:w-1
            if (g(i,j)==k) 
                if g(i,j+1)~=0
                    adjb(k,e)=g(i,j+1);            
                    e=e+1;
                else
                    adjb(k,1)=0;
                end
             end
        end
    end
end
adjb;
for i=1:br-1
    for j=h:-1:1
        for k=j:-1:2
            if adjb(i,j)==adjb(i,k-1)
                adjb(i,j)=0;
            end
        end
    end
end
adjb;
x=length(adjb(:,1));
ab=length(adjb(1,:));
for i=1:x
    for j=1:ab
        if adjb(i,j)==0 && j~=ab
            if adjb(i,j+1)~=0
                adjb(i,j)=adjb(i,j+1);
                adjb(i,j+1)=0;
            end
        end
        if adjb(i,j)~=0
            adjb(i,j)=adjb(i,j)-1;
        end
    end
end
adjb;
for i=1:x-1
    for j=1:ab
        adjcb(i,j)=adjb(i+1,j);
    end
end
b=length(adjcb);

%% 电压电流程序

for i=1:no
    vb(i,1)=1;
end
for s=1:10
for i=1:no
    nlc(i,1)=conj(complex(P(i,1),Q(i,1)))/(vb(i,1));
end
nlc;
for i=1:br
    Ibr(i,1)=nlc(i+1,1);
end
Ibr;
xy=length(adjcb(1,:));
for i=br-1:-1:1
    for k=1:xy
        if adjcb(i,k)~=0
            u=adjcb(i,k);
            %Ibr(i,1)=nlc(i+1,1)+Ibr(k,1);
            Ibr(i,1)=Ibr(i,1)+Ibr(u,1);
        end
    end      
end
Ibr;
for i=2:no
      g=0;
      for a=1:b 
          if xy>1
            if adjcb(a,2)==i-1 
                u=adjcb(a,1);
                vb(i,1)=((vb(u,1))-((Ibr(i-1,1))*(complex((R(i-1,1)),X(i-1,1)))));
                g=1;
            end
            if adjcb(a,3)==i-1 
                u=adjcb(a,1);
                vb(i,1)=((vb(u,1))-((Ibr(i-1,1))*(complex((R(i-1,1)),X(i-1,1)))));
                g=1;
            end
          end
        end
        if g==0
            vb(i,1)=((vb(i-1,1))-((Ibr(i-1,1))*(complex((R(i-1,1)),X(i-1,1)))));
        end
end
s=s+1;
end
nlc;
Ibr;
vb;
vbp=[abs(vb) angle(vb)*180/pi];

toc;
for i=1:no
    va(i,2:3)=vbp(i,1:2);
end
for i=1:no
    va(i,1)=i;
end
va;


Ibrp=[abs(Ibr) angle(Ibr)*180/pi];
PL(1,1)=0;
QL(1,1)=0;

%% 功率损耗
for f=1:br
    Pl(f,1)=(Ibrp(f,1)^2)*R(f,1);
    Ql(f,1)=X(f,1)*(Ibrp(f,1)^2);
    PL(1,1)=PL(1,1)+Pl(f,1);
    QL(1,1)=QL(1,1)+Ql(f,1);
end

Plosskw=(Pl)*100000;
Qlosskw=(Ql)*100000;
PL=(PL)*100000;
QL=(QL)*100000;

voltage = vbp(:,1);
angle = vbp(:,2)*(pi/180);

for i=1:no
    if i==1
        Pg(i,1)=(Pt*100000)+PL;
        Qg(i,1)=(Qt*100000)+QL;
    else
        Pg(i,1)=0;
        Qg(i,1)=0;
    end
end

%% 显示潮流计算结果

Vm=voltage';deltad=angle';Pgen=Pg';Qgen=Qg';Pd=P*100000';Qd=Q*100000';Qinj=Qsh';

Pdt=sum(Pd);Qdt=sum(Qd);Pgent=sum(Pg);Qgent=sum(Qg);Qinjt=sum(Qsh);

n=1:no;
%% 标题为显示结果

tech=('                   径向分布潮流解决方案                       ');
disp(tech)
disp('=============================================================================')
head =['    Bus  Voltage  Angle    ------Load------    ---Substation---   Injected'
       '    No.  Mag.     Degree     kW       kVAr       kW       kVAr       kvar '
       '                                                                          '];
disp(head)
disp('=============================================================================')

% bus information including voltage profile, load profile and substation
for n=1:no
      fprintf(' %5g', n), fprintf(' %7.4f', Vm(n)),
     fprintf(' %8.4f', deltad(n)), fprintf(' %9.4f', Pd(n)),
     fprintf(' %9.4f', Qd(n)),  fprintf(' %10.4f', Pgen(n)),
     fprintf(' %10.4f ', Qgen(n)), fprintf(' %9.4f\n', Qinj(n))
end
disp('=============================================================================') 
    fprintf('      \n'), fprintf('    Total              ')
    fprintf(' %9.4f', Pdt), fprintf(' %9.4f', Qdt),
    fprintf(' %10.4f', Pgent), fprintf(' %10.4f', Qgent), fprintf(' %9.4f\n\n', Qinjt)
disp('=============================================================================') 

%% 各支路潮流及损耗

SLT = 0;
fprintf('\n')
fprintf('                     线路功率和损耗 \n\n')
disp('=====================================================================') 
fprintf('     --Line--  Power at bus & line flow    --Line loss--\n')
fprintf('     from  to       kW      kVAr            kW      kVAr\n')
disp('=====================================================================') 
for n = 1:no
busprt = 0;
   for L = 1:br;
       if busprt == 0
       fprintf('   \n'), fprintf('%6g', n), fprintf('      %9.4f', Pd(n))
       fprintf('%9.4f\n', Qd(n))
       

       busprt = 1;
       else
       end
       if nl(L)==n     k = nr(L);
       In = (Vm(n) - Vm(k))/Z(L);
       Ik = (Vm(k) - Vm(n))/Z(L);
       Snk = Vm(n)*conj(In)*100000;
       Skn = Vm(k)*conj(Ik)*100000;
       SL  = Snk + Skn;
       SLT = SLT + SL;
       elseif nr(L)==n  k = nl(L);
       In = (Vm(n) - Vm(k))/Z(L);
       Ik = (Vm(k) - Vm(n))/Z(L);
       Snk = Vm(n)*conj(In)*100000;
       Skn = Vm(k)*conj(Ik)*100000;
       SL  = Snk + Skn;
       SLT = SLT + SL;
       else
       end
         if nl(L)==n || nr(L)==n
         fprintf('%12g', k),
         fprintf('%10.4f', real(Snk)), fprintf('%10.4f', imag(Snk)),
         fprintf('%12.4f', real(SL)),
         fprintf('%12.4f\n', imag(SL))
         else
         end
  end
end
disp('=====================================================================') 
In;
Ik;
SLT = SLT/2;
A=SLT;
fprintf('   \n'), fprintf('    Total loss                         ')
fprintf('%12.4f', real(SLT)), fprintf('%12.4f\n', imag(SLT))
clear Ik In SL SLT Skn Snk

3 知识点

详细知识点可以下载这篇文档：