首先声明一下，这些是从网站上转载的，不是本人上编写的基于MATLAB的电力系统潮流计算%简单潮流计算的小程序，相关的原始数据数据数据输入格式如下：%B1是支路参数矩阵，第一列和第二列是节点编号。

节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点%编号，将变压器的串联阻抗置于低压侧处理。

%第三列为支路的串列阻抗参数。

%第四列为支路的对地导纳参数。

%第五烈为含变压器支路的变压器的变比%第六列为变压器是否是否含有变压器的参数，其中“1”为含有变压器，%“0”为不含有变压器。

%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点%负荷功率参数；第三列为节点电压参数；第六列为节点类型参数，其中%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

%X为节点号和对地参数矩阵。

其中第一列为节点编号，第二列为节点对地%参数。

n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');X=input('节点号和对地参数:X=');Y=zeros(n);Times=1; %置迭代次数为初始值%创建节点导纳矩阵for i=1:n1if B1(i,6)==0 %不含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else %含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);Y(q,q)=Y(q,q)+1/(B1(i,5)^2*B1(i,3));endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将OrgS、DetaS初始化%创建OrgS，用于存储初始功率参数h=0;j=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+im ag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B 2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+im ag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B 2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendOrgS%创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)==3t=t+1;PVU(t,1)=B2(i,3);endendPVU%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:n %对PV节点的处理，注意这时不可再将h初始化为0 if i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endendDetaS%创建I，用于存储节点电流参数i=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endendI%创建Jacbi(雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对PQ节点的处理if B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n %对PV节点的处理if B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));else %非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendJacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:n %对PQ节点处理if B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:n %对PV节点的处理if B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2%开始循环**********************************************************************while abs(max(DetaU))>prOrgS=zeros(2*n-2,1); %!!!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+im ag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B 2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:nif i~=isb&B2(i,6)==3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))+im ag(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3))*(real(Y(i,j))*real(B2(j,3))-imag(Y(i,j))*imag(B2(j,3)))-real(B 2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendOrgS%创建DetaSh=0;for i=1:nif i~=isb&B2(i,6)==2h=h+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2))-OrgS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,6)==3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2))-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1))^2+imag(PVU(t,1))^2-real(B2(i,3))^2-imag(B2(i,3))^2;endendDetaS%创建Ii=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1))/conj(B2(i,3));endendI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)==3h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3))+imag(I(h,1));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3))+real(I(h,1));Jacbi(2*h,2*k-1)=2*imag(B2(i,3));Jacbi(2*h,2*k)=2*real(B2(i,3));elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j))*real(B2(i,3))+real(Y(i,j))*imag(B2(i,3));Jacbi(2*h-1,2*k)=real(Y(i,j))*real(B2(i,3))+imag(Y(i,j))*imag(B2(i,3));Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbiDetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;for i=1:nif B2(i,6)==2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)==3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2Times=Times+1; %迭代次数加1endTimes一个原始数据的例子节点数5支路数5平衡节点编号5精度pr 0.000001B1(支路参数矩阵)[1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1]B2(节点参数矩阵)[0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1]X(节点号和对地参数)[1 0;2 0;3 0;4 0;5 0]程序2：%本程序采用牛顿-拉夫讯法对一个五节点电力系统进行潮流计算% n=input('请输入系统节点数目:n=');% nl=input('请输入系统支路数目:nl=');% ph=input('请输入平衡节点号:ph=');% jd=input('请输入误差精度：jd=');% B=input('请输入由支路参数形成的矩阵:B=');% A=input('请输入各节点参数形成的矩阵:A='); clear alln=5;nl=5;ph=5;jd=1E-6;%支路参数矩阵Br=[1 2 0.04+0.25*j 0.25j 1;1 3 0.1+0.35*j 0 1;2 3 0.08+0.30*j 0.25*j 1;2 4 0.015*j 0 1.05;3 5 0.03*j 0 1.05];%节点参数矩阵:2-PQ节点，3-PV节点，4-平衡节点A=[-1.6-0.8*j 1 0 2;-2-1*j 1 0 2;-3.7-1.3*j 1 0 2;5 1.05 1.05 3;0 1.05 0 4];%雅克比矩阵形成Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);for i=1:nlp=Br(i,1);q=Br(i,2);Y(p,q)=Y(p,q)-1./(Br(i,3)*Br(i,5));Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(Br(i,3))+Br(i,4);Y(p,p)=Y(p,p)+1./(Br(i,3)*Br(i,5)^2)+Br(i,4);enddisp('节点导纳矩阵为：Y=');disp(Y);G=real(Y);B=imag(Y);for i=1:ne(i)=real(A(i,2));f(i)=imag(A(i,2));V(i)=A(i,3);endfor i=1:nS(i)=A(i);endP=real(S);Q=imag(S);%雅克比矩阵求取Ci=0;a=1;NO=2*n;N=NO-1;while a~=0a=0;for i=1:nif i~=phC(i)=0;D(i)=0;for p=1:nC(i)=C(i)+G(i,p)*e(p)-B(i,p)*f(p);D(i)=D(i)+G(i,p)*f(p)+B(i,p)*e(p);endP1=e(i)*C(i)+f(i)*D(i);Q1=f(i)*C(i)-e(i)*D(i);V2=e(i)^2+f(i)^2;if A(i,4)~=3DP=P(i)-P1;DQ=Q(i)-Q1;for k=1:n%非平衡节点时，非对角线元素if k~=ph & k~=iX1=-G(i,k)*f(i)+B(i,k)*e(i);X2=-G(i,k)*e(i)-B(i,k)*f(i);X3=-X2;X4=X1;p=2*i-1;q=2*k-1;J(p,q)=X1;J(p,N)=DP;m=p+1;J(m,q)=X3;J(m,N)=DQ;q=q+1;J(p,q)=X2;J(m,q)=X4;%非平衡节点时对角线元素elseif k~=ph & k==iX1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);X2=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);X3=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);X4=D(i)+B(i,i)*e(i)-G(i,i)*f(i);p=2*i-1;q=2*k-1;J(p,q)=X1;J(p,N)=DP;m=p+1;J(m,q)=X3;J(m,N)=DQ;q=q+1;J(p,q)=X2;J(m,q)=X4;endendelseDP=P(i)-P1;DV=V(i)^2-V2;for s=1:nif s~=ph & s~=iX1=-G(i,s)*f(i)+B(i,s)*e(i);X2=-G(i,s)*e(i)-B(i,s)*f(i);X5=0;X6=0;p=2*i-1;q=2*s-1;J(p,q)=X1;J(p,N)=DP;m=p+1;J(m,q)=X5;J(m,N)=DV;q=q+1;J(p,q)=X2;J(m,q)=X6;elseif s~=ph & s==iX1=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);X2=-C(i)-G(i,s)*e(i)-B(i,s)*f(i);X5=-2*f(i);X6=-2*e(i);p=2*i-1;q=2*s-1;J(p,q)=X1;J(p,N)=DP;m=p+1;J(m,q)=X5;J(m,N)=DV;q=q+1;J(p,q)=X2;J(m,q)=X6;endendendendenddisp('系统雅克比矩阵为：J=');disp(J);%对于高压电力系统，最大元素出现在对角线位置，简单计算如下DJ=J(:,[1:NO-2]);DP=J([1:NO-2],N);DU=DJ\DP;disp('************************************');disp(' ')disp('第几次的修正值DU？');disp(Ci);disp('***********************************');disp(DU);for i=1:NO-2eps=abs(DP(i));if eps>=jd;a=a+1;endendCi=Ci+1;for i=1:n-1e(i)=e(i)-DU(2*i);f(i)=f(i)-DU(2*i-1);endE=e+j*f;disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');disp(' ');disp('第几次迭代后的近似计算值');disp(Ci);disp('&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&');for i1=1:nV(i1)=sqrt(e(i1)^2+f(i1)^2);ANGLE(i1)=atan(f(i1)./e(i1))*(180./pi);end%下列输出结果均按节点号从小到大排列disp('各个节点的电压实际值为：E=');disp(E);disp('各个节点电压幅值为：V=');disp(V);disp('各个节点的电压相交为：ANGLE=');disp(ANGLE);disp('各个节点的功率误差变化');disp(J([1:8],9)');MDU(Ci)=max(abs(J([1:8],9)'));for p=1:nSm(p)=0;for q=1:nSm(p)=Sm(p)+conj(Y(p,q))*conj(E(q));endS(i)=E(p)*Sm(p);endenddisp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @'); disp(' ');disp('计算的最终功率近似值如下各项值');disp(' ');disp('@ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @'); disp('各个节点功率为:S=');disp(S);disp('各条支路的首端功率Si为:');for i=1:nlp=Br(i,1);q=Br(i,2);A(i)=E(p)*[conj(E(p))-conj(E(q))]*conj(1./(Br(i,3)));SBS(i)=[abs(E(p))]^2*conj(Br(i,4))+A(i);enddisp(SBS);disp('各条支路的末端功率Sj为:');for i=1:nlp=Br(i,1);q=Br(i,2);A2=E(q)*[conj(E(q))-conj(E(p))]*conj(1./(Br(i,3)));SBE(i)=[abs(E(q))]^2*conj(Br(i,4))+A2;enddisp(SBE);disp('各条支路的功率损耗为:');for i=1:nlp=Br(i,1);q=Br(i,2);SBL(i)=SBS(i)+SBE(i);enddisp(SBL);%绘制功率误差曲线，先用最小二乘法进行曲线拟合for i=1:CiCSH(i)=i;endP=polyfit(CSH,MDU,3);plot(CSH,MDU,'-O')。