show下我自己用matlab编写的潮流程序

小草1987 · 发表于 2009-11-18 22:00:45

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有账号？立即加入

×

本文的的案例是在福州大学编写的《电力系统计算程序及其实现》书中的第二页，所使用的方法也是本书中的第1，2，3章，考虑了节点优化和稀疏导纳矩阵储存问题，以下各子程序都是存在M文件中，只需把各个子程序拷贝在M文件中，分别调用即可。

LGP子程序（导入数据）：

SJ1=xlsread('branch.xls')
IZA=SJ1(:,1)    %存支路状态数%
IZ1=SJ1(:,2)    %存支路一端的节点号%
IZ2=SJ1(:,3)    %存支路另一端的节点号%
Z1=SJ1(:,4)       %存支路正序电阻%
Z2=SJ1(:,5)       %存支路正序电抗%
Z3=SJ1(:,6)       %存支路正序电纳或非标准变比%
N=SJ1(1,7)       %存网络节点数%
M=SJ1(1,8)       %存网络支路数%
SJ2=xlsread('generate.xls')
IWGA=SJ2(:,1)    %存发电机状态数%
IWG=SJ2(:,2)    %存发电机节点数%
W1=SJ2(:,3)       %存发电机有功功率%
W2=SJ2(:,4)       %存发电机无功功率%
IQ=SJ2(1,5)       %存发电机台数%
SJ3=xlsread('load.xls')
ILP=SJ3(:,1)    %存负荷状态数%
ILD=SJ3(:,2)    %存负荷节点号%
WL1=SJ3(:,3)    %存负荷有功功率%
WL2=SJ3(:,4)    %存负荷无功功率%
IP=SJ3(1,5)       %存负荷个数%
SJ4=xlsread('jd.xls')
N0=SJ4(1,1)       %存平衡节点的节点号%
U0=SJ4(1,2)       %存平衡节点的给定电压值%
IPV=SJ4(:,3)    %存PV节点的节点号%
PV=SJ4(:,4)       %存PV节点的给定电压值%
N1=SJ4(1,5)       %存PV节点数%
UP=SJ4(1,6)       %存PQ节点的电压初值%

LOP子程序（节点优化编号）：

for I=1:N             %寄存各节点连接支路数的ID数组要预先充零%
ID(I)=0;
end
for K=1:M
if IZA(K)==0|IZA(K)==4  %停运支路IZA(K)=0和对地支路IZA(K)=4不属于统计范围之内%
   continue
end
I=IZ1(K);          %先将当前支路两端的节点号从储存它们的IZ1和IZ2数组取出，存于I,J单元%
J=IZ2(K);
I1=ID(I);          %将I,J节点已连接的支路数从ID数组取出，存于I1,J1单元%
J1=ID(J);
X=0;
Y=0;
for K1=1:I1          %判别当前支路是否与前面支路并联%
   if IX(I,K1)==J
      X=1;
      break
   end
end
if X~= 0
   for K1=1:J1
      if IX(J,K1)==I
         Y=1;
      break
      end
   end
end
if X==0|Y==0          %如果是非并联支路，先将I,J的支路数各自增一，然后各自寄存对端的节点号与IX数组%
I1=I1+1;
J1=J1+1;
ID(I)=I1;
ID(J)=J1;
IX(I,I1)=J;
IX(J,J1)=I;
end
end
for I5=1:N                %在尚未编号的网络中查找连接支路最少的节点%
I1=1000;
for I=1:N
      I2=ID(I);
      if I2<=1
         I1=I2;
         J=I;
         break
      end
      if I2<I1
         I1=I2;
         J=I;
      end
end
INA(I5)=J;
INB(J)=I5;
ID(J)=1000;
for J1=1:I1                      %挨个取出一个与J节点相连接的节点号，与J相连接的支路数有I1条，即与J节点相连接的节点数有I1个%
      I3=IX(J,J1);                      %先将IX数组取出一个与J节点相连接的节点号，存于I3单元%
      J2=ID(I3);                         %再从ID数组取出I3节点连接的支路数存于J2中%
      for J3=1:J2                      %挨个取出一个与I3节点相连接的节点号%
         J4=IX(I3,J3);
         if J4==J                         %判断J4是否等于J%
            break
         end
      end
      IX(I3,J3)=IX(I3,J2);             %去掉与I3节点相连接的节点号J%
      ID(I3)=J2-1;
end
X=0;
for J1=1:I1-1                      %去掉J节点后，使原于J节点相连接的所有节点每两个之间增加一条新支路%
      K1=IX(J,J1);                      %挨个从IX数组取出一个原与J节点相连接的节点，存于K1单元%
      K3=ID(K1);                         %从ID数组取出K1节点所连接的支路数，，存于K3单元%
      for J2=J1+1:I1                      %挨个取出原与J节点相连接的K1节点之后的节点，存于K2%
         K2=IX(J,J2);
         for J3=1:K3
            if IX(K1,J3)==K2                   %从IX数组挨个取出与K1节点相连接的节点号IX(K1,J3)，与K2进行比较，如果均不等于K2，则K1与K2节点之间无支路关系%
               X=1;                               %如果IX(K1,J3)等于K2,则表明K1与K2节点之间已有支路连接，不必增加新支路%
               break
            end
         end
         if X==0                            %K1与K2节点之间增加一条新支路%
            K3=K3+1;                            %K1节点连接支路数增加一%
            IX(K1,K3)=K2;                      %寄存对端节点%
            ID(K1)=K3;
            K4=ID(K2)+1;                      %K2节点连接支路数增加一%
            IX(K2,K4)=K1;                      %寄存对端的节点%
            ID(K2)=K4;
         end
      end
end
end

LKP子程序（优化编号后的新节点储存）：

for K=1:M          %将支路原有旧节点号换成新的节点号%
I=IZ1(K);
J=IZ2(K);
IZ1(K)=INB(I);
if J==0
      continue
end
IZ2(K)=INB(J);
end
for K=1:IQ          %将发电机旧节点号换成新的节点号%
I=IWG(K);
IWG(K)=INB(I);
end
for K=1:IP          %将负荷旧的节点号换成新的节点号%
I=ILD(K);
ILD(K)=INB(I);
end
for K=1:N1          %将平衡节点的旧节点号换成新的节点号%
I=IPV(K);
IPV(K)=INB(I);
end
N0=INB(N0);          %将平衡节点旧的换成新的%

LDP子程序（形成导纳矩阵）：

for I=1:N                %储存自导纳的D11,D12数组要预先置零%
D11(I)=0;
D12(I)=0;
end
L=0;                   %非零互导纳元素的计数单元，开始置零%
for K=1:M
IG=IZA(K);          %IG是当前支路状态数的临时寄存单元%
if IG==0
      continue
end
I=IZ1(K);
J=IZ2(K);
R=Z1(K);
X=Z2(K);
B=Z3(K);
A=R*R+X*X;
if A==0
      continue
end
GIJ=R/A;             %计算支路导纳%
BIJ=-X/A;
Y=0;
if IG==1
   Y=1;
   GI=GIJ;          %计算输电线支路I节点自导纳，J节点自导纳和它们之间的互导纳%
   GJ=GIJ;
   BI=BIJ+B/2;
   BJ=BI;
end
if IG==2|IG==3
   Y=1;
   GJ=GIJ;          %计算变压器支路I节点自导纳，J节点自导纳和它们之间的互导纳%
   BJ=BIJ;
   GIJ=GIJ/B;
   BIJ=BIJ/B;
   GI=GIJ/B;
   BI=BIJ/B;
end
if Y==0
   D11(I)=D11(I)+GIJ; %对地支路时，只将当前支路的导纳累加到I节点的自导纳上%
   D12(I)=D12(I)+BIJ;
   continue
end
if Y==1
   D11(I)=D11(I)+GI;  %非对地支路时，将当前非对地支路导纳累加到I,J节点的自导纳上%
   D12(I)=D12(I)+BI;
   D11(J)=D11(J)+GJ;
   D12(J)=D12(J)+BJ;
   L=L+1;          %L是非零互导纳元素的计数单元%
   YZ1(L)=-GIJ;
   YZ2(L)=-BIJ;
   IY1(L)=I;
   IY2(L)=J;
   L=L+1;          %L是非零互导纳元素的计数单元%
   YZ1(L)=-GIJ;
   YZ2(L)=-BIJ;
   IY1(L)=J;
   IY2(L)=I;
end
end
J=0;                %J是有规则非零互导纳元素的计数单元，挨个累计%
K0=0;                %K0是有规则非零互导纳元素的计数单元，挨行累计%
for I=1:N             %I循环实现按行号由小到大将非零互导纳元素排列在Y1,Y2数组中%
J1=0;             %J1是当前行I非零互导纳元素的计数单元，开始置零%
for K=1:L       %K循环挨个检查不规则非零互导纳%
      if IY1(K)~=I
         continue
      end
      J3=IY2(K);    %IY1(K)如果等于I，则表明该非零互导纳YZ1(K),YZ2(K)是第I行的元素，将该非零互导纳的列号从IY2(K)单元取出，存于J3单元%
      Y=0;
      for K1=1:J1 %K1循环对当前行I已有规则排列在Y1,Y2数组中的互导纳元素进行挨个排查，是否有列号等于J3的元素%
         K2=K0+K1;
         if J3==IY(K2)
            Y=1;
            break
         end
      end
      if Y==0       %不存在列号等于J3,即非并联支路%
         J=J+1;
         Y1(J)=YZ1(K);
         Y2(J)=YZ2(K);
         IY(J)=J3;
         J1=J1+1;
         continue
      end
      if Y==1       %存在列号等于J3,即并联支路，将非零互导纳YZ1(K)，YZ2(K)累加到Y1(K2),Y2(K2)单元中%
         Y1(K2)=Y1(K2)+YZ1(K);
         Y2(K2)=Y2(K2)+YZ2(K);
      end
end
IN(I)=J1;          %将当前行I非零互导纳元素的个数J1存于IN数组，IN数组是用来存放正序导纳矩阵每行非零互导纳元素个数的%
K0=K0+J1;
end
LIP子程序（PV,PQ，平衡节点赋予初值）：

for I=1:N
U1(I)=UP;
U2(I)=0;
PD(I)=0;
QD(I)=0;
PF(I)=0;
QF(I)=0;
IVI(I)=0;
end
for I=1:IP                %PD,QD数组中将接有负荷的节点填上给定的负荷功率%
if ILP(I)==0
      continue
end
J=ILD(I);
PD(J)=WL1(I);
QD(J)=WL2(I);
end
for I=1:IQ                %在PF,QF数组中将接有发电机的节点填上给定的发电功率%
if IWGA(I)==0
      continue
end
J=IWG(I);
PF(J)=W1(I);
QF(J)=W2(I);
end
for I=1:N1                %给PV节点加标志1%
J=IPV(I);
U1(J)=PV(I);
IVI(J)=1;
end
U1(N0)=U0;

LJP子程序（牛顿拉夫逊法解雅可比矩阵）：

for IT=1:20                %IT是迭代次数储存单元%
AM=0;                   %AM是用来寄存节点功率误差的最大绝对值%
K0=1;                   %K0是导纳矩阵非零互导纳元素的指针，开始置1%
for I=1:N                %每次形成雅可比矩阵的一行元素及其相应的常数项，然后进行消去和规格化运算%
      A=D11(I)*U1(I);
      B=D12(I)*U1(I);
      R=D11(I)*U2(I);
      X=D12(I)*U2(I);
      A3=A-X;             %A3，B3单元寄存第I节点电流的实部和虚部%
      B3=R+B;
      J5=1;                %J5是雅可比矩阵第I行非零元素的计数单元，开始置1%
      for IG=1:IN(I)       %IG循环表示挨次取出导纳矩阵第I行的一个非零互导纳元素，形成雅可比矩阵相应的一个非零元素%
         J=IY(K0);
         A3=A3+Y1(K0)*U1(J)-Y2(K0)*U2(J);
         B3=B3+Y1(K0)*U2(J)+Y2(K0)*U1(J);
         if I~=N0&J~=N0
            J5=J5+1;
            JK(J5)=IY(K0);
            AK1(J5)=Y1(K0)*U2(I)-Y2(K0)*U1(I);
            AK3(J5)=Y1(K0)*U1(I)+Y2(K0)*U2(I);
            AK2(J5)=-AK3(J5);
            AK4(J5)=AK1(J5);
         end
         K0=K0+1;
      end
      if I==N0             %第I行如果是平衡节点%
         GQ(I)=A3*U2(I)-B3*U1(I);
         GP(I)=A3*U1(I)+B3*U2(I);
         CK1(I)=0;
         CK2(I)=0;
         JF(I)=0;
         continue
      end
      P2=A3*U1(I)+B3*U2(I)+PD(I);
      Q2=A3*U2(I)-B3*U1(I)+QD(I);
      GP(I)=P2;
      GQ(I)=Q2;
      P2=PF(I)-P2;
      Q2=QF(I)-Q2;
      P3=P2;
      Q3=Q2;
      if IVI(I)==1
         Q(I)=Q2;
         if Q(I)>0|IT-1<5
            Q3=0;
         end
      end
      AK3(1)=A3+A+X;
      AK4(1)=B3-B+R;
      JK(1)=I;
      CK2(I)=P2;
      Y=0;
      if IVI(I)~=1
         Y=1;
      end
      if Y==0
         Q(I)=Q2;
         if IT>=5
            if Q2<0
                  Y=1;
            end
         end
         if Y==0
            for J=2:J5
                  AK1(J)=0;
                  AK2(J)=0;
            end
            AK1(1)=2*U1(I);
            AK2(1)=2*U2(I);
            Q2=0;
            CK1(I)=Q2;
         end
      end
      if Y==1
         AK1(1)=R-B-B3;
         AK2(1)=A3-A-X;
         CK1(I)=Q2;
      end
      C5=abs(P3);             %判别节点功率误差的绝对值是否大于AM%
      D5=abs(Q3);
      if C5>AM
         I0=I;
         AM=C5;
      end
      if D5>AM
         I0=I;
         AM=D5;
      end
      K=1;
      for I1=1:I-1          %进行消去运算%
         X=0;
         for I3=2:J5
            if JK(I3)==I1
               X=1;
               break
            end
         end
         if X==0
            K=K+JF(I1);
            continue
         end
         for IG=1:JF(I1)
               Y=0;
               for I2=1:J5
                  if JK(I2)==IJ(K)
                     Y=1;
                     break
                  end
               end
               if Y==0
                  J5=J5+1;
                  AK1(J5)=-AK1(I3)*AJ1(K)-AK2(I3)*AJ3(K);
                  AK2(J5)=-AK1(I3)*AJ2(K)-AK2(I3)*AJ4(K);
                  AK3(J5)=-AK3(I3)*AJ1(K)-AK4(I3)*AJ3(K);
                  AK4(J5)=-AK3(I3)*AJ2(K)-AK4(I3)*AJ4(K);
                  JK(J5)=IJ(K);
               end
               if Y==1
               AK1(I2)=AK1(I2)-AK1(I3)*AJ1(K)-AK2(I3)*AJ3(K);
               AK2(I2)=AK2(I2)-AK1(I3)*AJ2(K)-AK2(I3)*AJ4(K);
               AK3(I2)=AK3(I2)-AK3(I3)*AJ1(K)-AK4(I3)*AJ3(K);
               AK4(I2)=AK4(I2)-AK3(I3)*AJ2(K)-AK4(I3)*AJ4(K);
               end
               K=K+1;
            end
            CK1(I)=CK1(I)-AK1(I3)*CK1(I1)-AK2(I3)*CK2(I1);
            CK2(I)=CK2(I)-AK3(I3)*CK1(I1)-AK4(I3)*CK2(I1);
         end
         K5=K;
         A=AK1(1)*AK4(1)-AK2(1)*AK3(1);
         if A==0
            CK1(I)=0;
            CK2(I)=0;
            JF(I)=0;
            continue
         end
         A3=AK4(1)/A;
         B3=-AK2(1)/A;
         C3=-AK3(1)/A;
         D3=AK1(1)/A;
      for J=2:J5
         if JK(J)>I
            IJ(K)=JK(J);
            AJ1(K)=A3*AK1(J)+B3*AK3(J);
            AJ2(K)=A3*AK2(J)+B3*AK4(J);
            AJ3(K)=C3*AK1(J)+D3*AK3(J);
            AJ4(K)=C3*AK2(J)+D3*AK4(J);
            K=K+1;
         end
      end
      A5=CK1(I);
      B5=CK2(I);
      CK1(I)=A3*A5+B3*B5;
      CK2(I)=C3*A5+D3*B5;
      JF(I)=K-K5;
end
      fprintf('%6.2',AM)
      fprintf('%6.2',IT)
      fprintf('%6.2',I0)
      if AM<=1.0E-4
         break
      end
      if IT>20
         break
      end
      for I=N-1:-1:1                   %I循环表示回代运算从n-1开始，倒推到第1行，每次计算一个节点电压的修正量%
         for IG=1:JF(I)
            K=K-1;
            J=IJ(K);
            CK1(I)=CK1(I)-AJ1(K)*CK1(J)-AJ2(K)*CK2(J);
            CK2(I)=CK2(I)-AJ3(K)*CK1(J)-AJ4(K)*CK2(J);
         end
      end
      for I=1:N                      %计算各节点的新电压%
         U1(I)=U1(I)+CK1(I);
         U2(I)=U2(I)+CK2(I);
   end
   for I=1:N1
         J=IPV(I);
         if Q(J)<0&IT>=5
            continue
         end
         C=PV(I)/sqrt(U1(J)*U1(J)+U2(J)*U2(J));
         U1(J)=C*U1(J);
         U2(J)=C*U2(J);
   end
  end

LRP子程序（输出节点信息和支路信息）：

GP(N0)=GP(N0)+PD(N0);
GQ(N0)=GQ(N0)+QD(N0);
PG=0;
QG=0;
PL=0;
QL=0;
PI=180/3.14159;
for I=1:N                   %I循环实现每次输出一个节点的信息，输出是按旧的节点号为顺序号%
I5=INB(I);
E=U1(I5);
F=U2(I5);
A=sqrt(E*E+F*F);
B=PI*atan(F/E);
PG=PG+GP(I5);
QG=QG+GQ(I5);
PL=PL+PD(I5);
QL=QL+QD(I5);
JD(I,1)=I;
JD(I,2)=A;
JD(I,3)=B;
JD(I,4)=GP(I5);
JD(I,5)=GQ(I5);
JD(I,6)=PD(I5);
JD(I,7)=QD(I5);
end
PG
QG
PL
QL
PLOSS=0;
QLOSS=0;
QB=0;
for K=1:M                   %K循环表示每次输出一条支路的信息%
IG=IZA(K);
I=IZ1(K);
J=IZ2(K);
R=Z1(K);
X=Z2(K);
B1=Z3(K);
if IG==0
      ZL(K,1)=K;
      I=INA(I);
      J=INA(J);
      ZL(K,2)=I;
      ZL(K,3)=I;
   continue
end
if R==0&X==0
      I=INA(I);
      J=INA(J);
      ZL(K,1)=K;
      ZL(K,2)=I
      ZL(K,3)=I;
      continue
end
E=U1(I);
F=U2(I);
if IG==4
      A=0;
      B=0;
end
if IG~=4
      A=U1(J);
      B=U2(J);
end
if IG==2|IG==3
      E=E/B1;
      F=F/B1;
      B1=0;
end
A3=R*R+X*X;
C3=E-A;
D3=F-B;
A5=(C3*R+D3*X)/A3;
B5=(D3*R-C3*X)/A3;
C3=(E*E+F*F)*B1/2;
D3=(A*A+B*B)*B1/2;
P1=E*A5+F*B5;
Q1=F*A5-E*B5-C3;
P2=-A*A5-B*B5;
Q2=-B*A5+A*B5-D3;
if IG~=4
      PLOSS=PLOSS+P1+P2;
      QLOSS=QLOSS+Q1+Q2;
      QB=QB-C3-D3;
end
I=INA(I);
if J~=0
      J=INA(J);
end
ZL(K,1)=K;
ZL(K,2)=I;
ZL(K,3)=J;
ZL(K,4)=P1;
ZL(K,5)=P2;
ZL(K,6)=Q1;
ZL(K,7)=Q2;
end
PLOSS
QLOSS
QB
xlswrite('outext.xls',JD);
xlswrite('outext1.xls',ZL);

gxuyang · 发表于 2009-11-28 08:50:54

我想看看，学习一下

zgk001 · 发表于 2009-11-29 10:02:05

那个excel文件读入是怎么处理的？对于节点数有没有要求？

lhw606 · 发表于 2009-11-29 16:33:06

不错，支持原创~~~~~~~~~~~~

klein · 发表于 2009-12-6 21:56:21

有些复杂，没弄明白excel是怎么调用的？

小草1987 · 发表于 2009-12-29 23:53:37

回复 5# klein

福州大学编写的《电力系统计算程序及其实现》书中的第二页给了个表格，你把它输入到excel文件中，matlab中xlsread就可读取

原随风 · 发表于 2010-1-5 11:45:42

这样的帖子要顶

fzhang1978 · 发表于 2010-1-5 18:31:32

学习一下，谢谢

rena1121 · 发表于 2010-3-24 13:40:56

好,力顶！！！

wz0430131 · 发表于 2010-3-28 18:54:24

看一下，呵呵

账号		自动登录	找回密码
密码			立即加入

[讨论] show下我自己用matlab编写的潮流程序

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

评分

相关帖子