show下我自己用matlab编写的潮流程序
本文的的案例是在福州大学编写的《电力系统计算程序及其实现》书中的第二页，所使用的方法也是本书中的第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); 我想看看，学习一下 那个excel文件读入是怎么处理的？对于节点数有没有要求？ 不错，支持原创~~~~~~~~~~~~ 有些复杂，没弄明白excel是怎么调用的？ 回复 5# klein
福州大学编写的《电力系统计算程序及其实现》书中的第二页给了个表格，你把它输入到excel文件中，matlab中xlsread就可读取 这样的帖子要顶 学习一下，谢谢 好,力顶！！！ 看一下，呵呵
