5节点电力系统潮流计算matlab程序

合集下载

电力系统潮流计算matlab程序

电力系统潮流计算matlab程序电力系统潮流计算是电力系统运行和规划中的重要环节，它用于计算电力系统中各节点的电压、功率和电流等参数。

随着电力系统规模的不断扩大和复杂性的增加，传统的手工计算方法已经无法满足需求，因此，利用计算机编程进行潮流计算成为了一种必要的选择。

Matlab是一种功能强大的科学计算软件，它提供了丰富的数学函数和工具箱，可以方便地进行电力系统潮流计算。

下面我将介绍一下如何使用Matlab编写电力系统潮流计算程序。

首先，我们需要建立电力系统的节点模型。

节点模型是电力系统中各节点的电压、功率和电流等参数的数学表示。

在Matlab中，我们可以使用矩阵来表示节点模型。

假设电力系统有n个节点，我们可以定义一个n×n的复数矩阵Y来表示节点之间的导纳关系，其中Y(i,j)表示节点i和节点j之间的导纳。

同时，我们还需要定义一个n×1的复数向量V来表示各节点的电压，其中V(i)表示节点i的电压。

接下来，我们需要编写潮流计算的主程序。

主程序的主要功能是根据节点模型和潮流计算算法，计算出各节点的电压、功率和电流等参数。

在Matlab中，我们可以使用循环语句和矩阵运算来实现潮流计算。

具体的计算过程可以参考电力系统潮流计算的算法。

在编写主程序之前，我们还需要定义一些输入参数，如电力系统的节点数、发电机节点和负荷节点等。

这些参数可以通过用户输入或者读取文件的方式获取。

同时，我们还需要定义一些输出参数，如各节点的电压、功率和电流等。

这些参数可以通过矩阵运算和循环语句计算得到，并输出到文件或者显示在屏幕上。

最后，我们需要进行程序的测试和调试。

可以通过输入一些测试数据，运行程序并检查输出结果是否正确。

如果发现程序有错误或者结果不准确，可以通过调试工具和打印调试信息的方式进行调试。

总之，利用Matlab编写电力系统潮流计算程序可以提高计算效率和准确性，为电力系统的运行和规划提供有力的支持。

当然，编写一个完整的潮流计算程序需要考虑很多细节和特殊情况，这需要有一定的电力系统和编程知识。

潮流计算matlab程序

clear;%各节点参数：节点编号，类型，电压幅值，电压相位，注入有功，注入无功%类型：1＝PQ节点，2＝PV节点，3＝平衡节点%本程序中将最后一个节点设为平衡节点R_1=[1 1 1.0 0 0.2 0.2j;2 1 1.0 0 -0.45 -0.15j;3 1 1.0 0 -0.45 -0.05j;4 1 1.0 0 -0.6 -0.1j;5 3 1.0 0 0 0];%支路号首端节点末端节点支路导纳R_2=[1 5 2 1.25-3.75j;2 23 10.00-30.00j;3 34 1.25-3.75j;4 1 4 2.50-7.50j;5 1 5 5.00-15.00j;6 1 2 1.667-5.00j];n=5;L=6;%需要改变的到此为止i=0;j=0;a=0;precision=1;k=0;Y=zeros(n,n);u=zeros(1,n);delt=zeros(1,n);P=zeros(1,n);Q=zeros(1,n);G=[];B=[];PP=[];uu=[];U=[];dp=[];dq=[];for a=1:Li=R_2(a,2);j=R_2(a,3);Y(i,j)=-R_2(a,4);Y(j,i)=Y(i,j);endfor a=1:nfor b=1:nif a~=bY(a,a)=Y(a,a)+Y(a,b);endendendfor i=1:nfor j=1:nif i==jY(i,j)=-Y(i,j);endendendY %形成导纳矩阵for i=1:nfor j=1:nG(i,j)=real(Y(i,j));B(i,j)=imag(Y(i,j));endendfor a=1:nu(a)=R_1(a,3);P(a)=R_1(a,5);Q(a)=R_1(a,6);delt(a)=R_1(a,4);endwhile precision>0.0001 %判断是否满足精度要求for a=1:n-1for b=1:npt(b)=u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));qt(b)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b))); endpt,qtpi(a)=sum(pt);qi(a)=sum(qt); %计算PQ节点的注入功率dp(a)=P(a)-pi(a);dq(a)=Q(a)-qi(a); %计算PQ节点的功率不平衡量endfor a=1:n-1for b=1:n-1if a==bH(a,a)=-qi(a)-u(a)^2*B(a,a); N(a,a)=pi(a)+u(a)^2*G(a,a);J(a,a)=pi(a)-u(a)^2*G(a,a); L(a,a)=qi(a)-u(a)^2*B(a,a);JJ(2*a-1,2*a-1)=H(a,a); JJ(2*a-1,2*a)=N(a,a);JJ(2*a,2*a-1)=J(a,a); JJ(2*a,2*a)=L(a,a);elseH(a,b)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));J(a,b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));N(a,b)=-J(a,b);L(a,b)=H(a,b);JJ(2*a-1,2*b-1)=H(a,b);JJ(2*a-1,2*b)=N(a,b);JJ(2*a,2*b-1)=J(a,b); JJ(2*a,2*b)=L(a,b);endendend %计算jocbi各项,并放入统一矩阵JJ中，对JJ下标统一编号JJfor a=1:n-1PP(2*a-1)=dp(a);PP(2*a)=dq(a);end %按统一矩阵形成功率不平衡uu=inv(JJ)*PP';precision=max(abs(uu)); %判断是否收敛for b=1:n-1delt(b)=delt(b)+uu(2*b-1);u(b)=u(b)+uu(2*b)*u(b); %将结果分解为电压幅值和角度end %求解修正方程，得电压幅值变化量（标幺值）和角度变化量k=k+1;endfor a=1:nU(a)=u(a)*(cos(delt(a))+j*sin(delt(a)));endfor b=1:nI(b)=Y(n,b)*U(b);%求平衡节点的注入电流endS5=U(n)*sum(conj(I))%求平衡节点的注入功率for a=1:nfor b=1:nS(a,b)=U(a)*(conj(U(a))-conj(U(b)))*conj(-Y(a,b));endend %求节点i,j节点之间的功率，方向为由i指向j,S %显示支路功率。

Matlab实现潮流计算程序

程序代码如下:111111.%读入数据clcclearfilename='123.txt';a=textread(filename)n=a(1,1);pinghengjd=a(1,2);phjddianya=a(1,3);jingdu=a(1,4);b=zeros(1,9);j1=0;[m1,n1]=size(a);for i1=1:m1if a(i1,1)==0j1=j1+1;b(j1)=i1;endendb;%矩阵分块a1=a(b(1)+1:b(2)-b(1)+1,1:n1);a2=a(b(2)+1:b(3)-1,1:n1);a3=a(b(3)+1:b(4)-1,1:n1);a4=a(b(4)+1:b(5)-1,1:n1);a5=a(b(5)+1:b(6)-1,1:n1);%设置初值vcz=1;dcz=0;kmax=20;k1=0;%求节点导纳矩阵a11=zeros(4,6);for i0=1:3for j0=1:6a11(i0,j0)=a1(i0,j0);a11(4,j0)=a2(1,j0);endenda11;linei=a11(1:4,2);linej=a11(1:4,3);liner=a11(1:4,4);linex=a11(1:4,5);lineb=a11(1:4,6);branchi=0;branchj=0;branchb=0;G=zeros(4,4);B=zeros(4,4);for k=1:4i2=linei(k,1);j2=linej(k,1);r=liner(k,1);x=linex(k,1);b=0;GIJ=r/(r*r+x*x);BIJ=-x/(r*r+x*x);if k>=4 & lineb(k)~=0k0=lineb(k);G(i2,j2)=-GIJ/k0;G(j2,i2)=G(i2,j2);B(i2,j2)=-BIJ/k0;B(j2,i2)=B(i2,j2);G(i2,i2)=G(i2,i2)+GIJ/k0/k0; B(i2,i2)=B(i2,i2)+BIJ/k0/k0;elseG(j2,i2)=-GIJ;G(i2,j2)=G(j2,i2);B(j2,i2)=-BIJ;B(i2,j2)=B(j2,i2);G(i2,i2)=G(i2,i2)+GIJ;b=lineb(k);B(i2,i2)=B(i2,i2)+BIJ+b;endG(j2,j2)=G(j2,j2)+GIJ;B(j2,j2)=B(j2,j2)+BIJ+b;endG;B;B=B.*i;Yf=G+BY=abs(Yf);alf=angle(Yf);%赋Jacobian矩阵参数P=zeros(n,1);Q=zeros(n,1);Pd=zeros(1,n);Qd=zeros(1,n);dP=zeros(1,n);dQ=zeros(1,n);PG=a4(:,3);PD=a4(:,5);QG=a4(:,4);QD=a4(:,6);i8=a4(:,2);for j8=1:length(i8)P(i8(j8))=PG(i8(j8))-PD(i8(j8));Q(i8(j8))=QG(i8(j8))-QD(i8(j8));enddelt=zeros(n,1);V=ones(n,1);V(3)=1.10;V(4)=1.05;ddelt=zeros(n,1);dV=zeros(n,1);A=zeros(2*n,2*n);B=zeros(2*n,1);Jacobian=Jaco(V,delt,n,Y,alf)%求取矩阵功率for j5=1:kmaxdisp(['第' int2str(j5) '次计算结果'])if k>=kmaxbreakendfor i10=1:4Pd(i10)=0;Qd(i10)=0;for j10=1:nPd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(d elt(i10)-delt(j10)-alf(i10,j10));Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(d elt(i10)-delt(j10)-alf(i10,j10));endendfor i4=1:3dP(i4)=P(i4)-Pd(i4);endfor j4=1:2dQ(j4)=Q(j4)-Qd(j4);endA=Jaco(V,delt,n,Y,alf)for i14=1:nB(i14*2-1)=-dP(i14);B(i14*2)=-dQ(i14);endif max(abs(B))>jingduX=A\B;for i16=1:nddelt(i16)=X(2*i16-1);dV(i16)=X(2*i16)*V(i16);endV=V+dVdelt=delt+ddeltelsebreakenddisp('----------------')end%流氓算法% for ii=1:2% V(ii)=V(ii)+dV(ii);% end% V222222.function A=Jaco(V,delt,n,Y,alf)%计算Jacobian矩阵for i7=1:nHd1(i7)=0;Jd1(i7)=0;for j7=1:nHd1(i7)=Hd1(i7)+V(i7)*Y(i7,j7)*V(j7)*sin(delt(i7)-delt(j7)-alf(i7,j7));Jd1(i7)=Jd1(i7)+V(i7)*Y(i7,j7)*V(j7)*cos(delt(i7)-delt(j7)-alf(i7,j7));endendfor i6=1:nfor j6=1:nif i6~=j6H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-delt(j6)-alf(i6,j6));N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));J(i6,j6)=-N(i6,j6);L(i6,j6)=H(i6,j6);elseH(i6,i6)=Hd1(i6)-V(i6)*Y(i6,i6)*V(i6)*sin(delt(i6)-delt(j6)-alf(i6,j6));J(i6,j6)=-Jd1(i6)+V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));N(i6,j6)=-Jd1(i6)-V(i6)*Y(i6,i6)*V(i6)*cos(alf(i6,i6));L(i6,i6)=-Hd1(i6)+V(i6)*Y(i6,i6)*V(i6)*sin(alf(i6,i6));endendend%修正Jacobian矩阵for j9=3for i9=1:nN(i9,j9)=0;L(i9,j9)=0;J(j9,i9)=0;L(j9,i9)=0;endendL(j9,j9)=1;for j9=4for i9=1:nH(i9,j9)=0;N(i9,j9)=0;J(i9,j9)=0;L(i9,j9)=0;H(j9,i9)=0;N(j9,i9)=0;J(j9,i9)=0;L(j9,i9)=0;endendH(j9,j9)=1;L(j9,j9)=1;%Jaco=[H N;J L];%Jaco=zeros(2*n,2*n);for i11=1:nfor j11=1:nJaco(2*i11-1,2*j11-1)=H(i11,j11); Jaco(2*i11-1,2*j11)=N(i11,j11); Jaco(2*i11,2*j11-1)=J(i11,j11);Jaco(2*i11,2*j11)=L(i11,j11);endendA=Jaco;33333.数据:4 4 1.05 0.000011 12 0.1 0.40 0.015282 1 4 0.12 0.50 0.019203 24 0.08 0.40 0.014131 1 3 0 0.3 0.909090911 1 0 0 0.30 0.182 2 0 0 0.55 0.133 3 0.5 0 0 01 3 1.10 0 0。

基于Matlab的五节点潮流仿真

基于Matlab的五节点潮流仿真潮流计算是电力系统设计与运行的基本计算，可将网络中各节点电压和功率计算出来，继而进行损耗的计算。

通过潮流计算可以检测和预知电力系统中负荷和网络结构的变化，从而保障电力系统安全、稳定运行。

文章以电力系统五节点潮流算例为基础，应用牛顿潮流计算方法，构建基于Matlab的程序与仿真，并对程序与仿真结果进行对比。

标签：电力系统五节点；程序与仿真；牛顿潮流计算引言潮流计算是在电力系统稳态时的一种运算，由系统赋予的运行条件和整体结构来确定整个系统运行的状况。

同时它也是电力系统应用最基本和最广泛的电气计算。

较早的时候大多数的人用高斯——赛德尔迭代法求解潮流，因为它在运算的时候内部储存量比较小，但其收敛性却一般，所以在这以后人们就选用阻抗法，经过修改和运算方式的改进后，牛顿法就成为人们应用潮流计算的常用方法，其主要的特点就是将方程的非线性化变为方程的线性化。

文章重点以牛顿潮流计算作为潮流计算的方法，在Matlab中应用电力系统五节点算例进行仿真对比分析。

1 牛顿法潮流计算简言之，其核心就是把非线性方程式的求解过程变为相对应的线性求解的过程。

在牛顿潮流计算中形成雅克比矩阵之后，在此基础上建立相对应的修正方程，它的基本计算流程为：1.求解各个节点导纳，并组成相对应的矩阵YB；2、设置各个节点的初始值ei（0）、fi（0）；3、把各个节点的初始电压值带入式中求解修正方程式中的非平衡量ΔPi（0）、ΔQi（0）、Ui（0）2 ；4、求雅克比矩阵各元素Nij、Nij、Jij、Lij与Rij、Sij；5、求解修正量Δei（0）、Δfi（0）；6、再次计算各个节点的电压；7、在之前的基础上迭代；8、求解线路上的功率以及平衡节点的功率。

其中节点电压可采用直角坐标形式或者极坐标的形式。

在牛顿潮流计算的修正方程处理的方式上能联系修正方程式的求解，求解过程能采取形成、消元、储存同时进行的方式。

另外节点编号也能做有关的优化，其优化法通常有三种，它们为静态法、半动态法和动态法。

5节点电力系统潮流计算matlab程序汇总

%1.形成节点导纳矩阵, yb55=6.250-18.750j;yb51=-5.000+15.000j;yb52=-1.250+3.750j;yb53=0.000-0.000j;yb54=0.000-0.000j; yb15=-5.000+15.000j;yb11=10.834-32.500j;yb12=-1.667+5.000j;yb13=-1.667+5.000j;yb14=-2.500+7.500j; yb25=-1.250+3.750j;yb21=-1.667+5.000j;yb22=12.917-38.750j;yb23=-10.000+30.000j;yb24=0.000-0.000j; yb35=0.000-0.000j;yb31=-1.667+5.000j;yb32=-10.000+30.000j;yb33=12.917-38.750j;yb34=-1.250+3.750j; yb45=0.000-0.000j;yb41=-2.500+7.500j;yb42=0.000-0.000j;yb43=-1.250+3.750j;yb44=3.750-11.250j; YB=[yb11 yb12 yb13 yb14 yb15; yb21 yb22 yb23 yb24 yb25 ;yb31 yb32 yb33 yb34 yb35; yb41 yb42 yb43 yb44 yb45 ;yb51 yb52 yb53 yb54 yb55]; disp('节点导纳矩阵YB=');disp(YB); %计算各节点功率的不平衡量设U=E+jF ；Y=G+Bj；E(1)=1.00;E(2)=1.00;E(3)=1.00;E(4)=1.00; F(1)=0;F(2)=0;F(3)=0;F(4)=0;G=real(YB);B=imag(YB); %设S=P+Bj； S(1)=0.20+0.20i;S(2)=-0.45-0.15i;S(3)=-0.40-0.05i;S(4)=-0.60-0.10i; P=real(S);Q=imag(S); k=0;precision=1; N1=4; while precision > 0.00001 E(5)=1.06;F(5)=0; for m=1:N1 for n=1:N1+1 %计算Pi,Qi,设Pi=Pt;Qi=Qt Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n))+F(m)*(G(m,n)*F(n)+B(m,n)*E(n))); Qt(n)=(F(m)*(G(m,n)*E(n)-B(m,n)*F(n))-E(m)*(G(m,n)*F(n)+B(m,n)*E(n))); end %设P,Q的改变量为dP,dQ dP(m)=P(m)-sum(Pt); dQ(m)=Q(m)-sum(Qt); end for m=1:N1 for n=1:N1+1 %计算Hij Nij Jij Lij H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m);N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m);L(m,n)=G(m,n)*F(m)-B(m,n)*E(m); end end for m=1:N1 for n=1:N1+1Bi(n)=G(m,n)*F(n)+B(m,n)*E(n); Ai(n)=G(m,n)*E(n)-B(m,n)*F(n); end %sum(Ai)，sum(Bi)用于实现公式中的sigerma从j到n的求和; H(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))+2*G(m,m)*F(m); N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m))+2*G(m,m)*E(m); J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m)); L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))); end %设雅可比矩阵为JJ，以下语句用来实现雅可比矩阵中对角线上元素H N J L 的排列 for m=1:N1 JJ(2*m-1,2*m-1)=H(m,m); JJ(2*m-1,2*m)=N(m,m); JJ(2*m,2*m-1)=J(m,m); JJ(2*m,2*m)=L(m,m); end %以下语句用于实现雅可比矩阵非对角线上元素的排列 for m=1:N1 for n=1:N1 if m==n else H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m); N(m,n)=G(m,n)*E(m)+B(m,n)*F(m); J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m); L(m,n)=G(m,n)*F(m)-B(m,n)*E(m); JJ(2*m-1,2*n-1)=H(m,n); JJ(2*m-1,2*n)=N(m,n); JJ(2*m,2*n-1)=J(m,n); JJ(2*m,2*n)=L(m,n); end end end %设由P,Q的改变量组成的8×1矩阵为PQ，由E,F的改变量组成的8×1矩阵为dU for m=1:N1 PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m); end dU=inv(JJ)*PQ';precision=max(abs(dU)); for n=1:N1 F(n)=F(n)+dU(2*n-1); E(n)=E(n)+dU(2*n); end for n=1:N1+1 U(n)=E(n)+(F(n))*j; end k=k+1; k-1, dU=dU',PQ,U end。

matlab电力系统潮流计算程序

matlab电力系统潮流计算程序电力系统潮流计算是电力系统分析的关键步骤之一，用于确定电力系统各节点的电压和相角分布。

以下是一个简单的MATLAB电力系统潮流计算的基本步骤和代码示例：1.定义电力系统参数：-定义系统节点数量、支路数据、发电机数据、负荷数据等电力系统参数。

```matlab%电力系统参数busdata=[1,1.05,0,0,0,0,0,0;2,1.02,0,0,0,0,0,0;%...其他节点数据];linedata=[1,2,0.02,0.06,0.03;%...其他支路数据];gendata=[1,2,100,0,999,1.05,0.95;%...其他发电机数据];loaddata=[1,50,20;%...其他负荷数据];```2.构建潮流计算矩阵：-利用节点支路导纳、节点负荷和发电机功率等信息构建潮流计算的阻抗矩阵。

```matlabYbus=buildYbus(busdata,linedata);```3.迭代求解潮流方程：-利用迭代算法（如牛顿-拉夫森法）求解潮流方程，更新节点电压和相角。

```matlab[V,delta]=powerflow(Ybus,gendata,loaddata,busdata);```4.结果分析和可视化：-分析计算结果，可视化电压和相角分布。

```matlabplotVoltageProfile(busdata,V,delta);```这只是一个简单的潮流计算示例。

具体的程序实现可能涉及更复杂的算法和工程细节，取决于电力系统的复杂性和精确性要求。

您可能需要根据实际情况和数据格式进行调整和改进。

在实际工程中，也可以考虑使用专业的电力系统仿真软件。

电力系统潮流计算的MATLAB辅助程序设计-潮流计算程序

电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。

此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。

现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法.高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。

lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。

(读者可以到MATPOWER主页下载MATPOWER4.1，然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算）一、高斯—赛德尔法潮流计算使用的程序：高斯—赛德法的具体使用方法读者可参考后面的实例,这里仅介绍各程序的编写格式：lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。

程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。

根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。

对于PV节点，如发电机节点,要提供一个无功功率限定值。

当给定电压过高或过低时，无功功率可能超出功率限定值。

在几次迭代之后(高斯—塞德尔迭代为10次），需要检查一次发电机节点的无功出力，如果接近限定值,电压幅值进行上下5%的调整,使得无功保持在限定值内。

lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。

并将这些参数放在名为linedata的文件中。

这个程序将阻抗转换为导纳,并得到节点导纳矩阵.busout：该程序以表格形式输出结果,节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。

电力系统潮流分析matlab代码

functiontisco% 这是一个电力系统潮流计算的程序n=input( '\n 请输入节点数:n=') ;m =input( '\n 请输入支路数: m =') ;ph=input( '\n 请输入平衡母线的节点号:ph=') ;B1=input( '\n 请输入支路信息: B1=') ;% 它以矩阵形式存贮支路的情况, 每行存贮一条支路% 第一列存贮支路的一个端点% 第二列存贮支路的另一个端点% 第三列存贮支路的阻抗% 第四列存贮支路的对地导纳% 第五列存贮变压器的变比, 注意支路为1% 第六列存贮支路的序号B2=input( '\n 请输入节点信息: B2=') ;% 第一列为电源侧的功率% 第二列为负荷侧的功率% 第三列为该点的电压值% 第四列为该点的类型: 1 为PQ 节点,2 为PV 节点,3 为平衡节点A =input( '\n 请输入节点号及对地阻抗: A =') ;ip=input( '\n 请输入修正值: ip=') ;% ip为修正值Y =zeros( n) ;e=zeros( 1,n) ;f=zeros( 1,n) ;no=2*ph-1;fori=1:nif A( i,2) ~=0p=A( i,1) ;Y(p,p) =1./A( i,2) ;endendfori=1:mp=B1( i,1) ;q=B1( i,2) ;Y( p,p) =Y( p,p) +1./( B1( i,3) *B1( i,5) ^2+B1( i,4) ./2;Y( p,q) =Y( p,q) -1./( B1( i,3) *B1( i,5) ) ;Y(q,p) =Y( p,q) ;Y( q,q) =Y( q,q) +1./B1( i,3) +B1( i,4) ./2;endG =real( Y ) ;B=imag( Y ) ;fori=1:ne( i) =real( B2( i,3) ) ;f( i) =imag( B2( i,3) ) ;S(i) =B2( i,1) -B2( i,2) ;V(i) =B2( i,3) ;endP=real( S) ;Q =imag( S) ;[ C ,D ,D F ] =xxf( G ,B ,e,f,P ,Q ,n,B2,ph,V ,no) ;J=jacci( Y ,G ,B ,P ,Q ,e,f,V ,C ,D ,B2,n,ph,no) ; [ De,D f] =hxf( J,D F,ph,n,no) ;t=0;whilem ax( abs( D e) ) >ip& m ax( abs( D f) ) >ipt=t+1;e=e+D e;f=f+D f;[ C ,D ,D F] =xxf( G ,B ,e,f,P,Q ,n,B2,ph,V ,no) ;J=jacci( Y ,G ,B ,P,Q ,e,f,V ,C ,D ,B2,n,ph,no) ; [ De,D f] =hxf( J,D F ,ph,n,no) ;endv=e'+f'*j;fori=1:nhh( i) =conj( Y( ph,i) *v( i) ) ;endS(ph) =sum( hh) *v( ph) ;B2( ph,1) =S( ph) ;V =abs( v) ;jd=angle( v) *180/pi;result1=[ A( :,1) ,real( v) ,imag( v) ,V ,jd,real( S.') , imag( S.') ,real( B2( :,1) ) ,imag( B2( :,1) ) ,real ( B2( :,2) ) ,imag( B2( :,2) ) ] ;fori=1:ma( i) =conj( ( v( B1( i,1) ) /B1( i,5) -v( B1( i,2) ) ) /B1( i,3) ) ;b( i) =v( B1( i,1) ) *a( i) -j*B1( i,4) *v( B1( i,1) ) ^2/2;c( i) =-v( B1( i,2) ) *a( i) -j*B1( i,4) *v( B1( i,2) ) ^2/2;endresult2=[ B1( :,6) ,B1( :,1) ,B1( :,2) ,real( b.') ,imag ( b.') ,real( c.') ,imag( c.') ,real( b.'+c.') ,imag( b.'+c.') ] ;printout( t,result1,S,b,c,result2) ;typeresult.mfunction[ C,D ,D F] =xxf( G ,B,e,f,P,Q ,n,B2,ph,V ,no) % 该子程序是用来求取D Ffori=1:nifi~=phC(i) =0;D(i) =0;for j=1:nC( i) =C( i) +G( i,j) *e( j) -B( i,j) *f( j) ;D( i) =D( i) +G( i,j) *f( j) +B( i,j) *e( j) ;endP1=C(i) *e( i) +D( i) *f( i) ;Q 1=C(i) *f( i) -D( i) *e( i) ;V 1=e(i) ^2+f( i) ^2;if B2( i,4) ==2p=2*i-1;D F( p) =P( i) -P1;p=p+1;D F( p) =V( i) ^2-V 1^2;elsep=2*i-1;D F( p) =P( i) -P1;p=p+1;D F( p) =Q( i) -Q 1;endendendD F=D F';ifph~=nD F( no,:) =[ ] ;D F( no,:) =[ ] ;endfunction [ D e,D f] =hxf( J,D F,ph,n,no)% 该子函数是为求取D e D fD X =J\D F;D X 1=D X ;x1=length( D X 1) ;ifph~=nD X( no) =0;D X( no+1) =0;fori=( no+2) :( x1+2)D X(i) =D X 1( i-2) ;endelseD X =[ D X 1,0,0] ;endk=0;[ x,y] =size( D X ) ;fori=1:2:xk=k+1;D f( k) =D X( i) ;D e( k) =D X( i+1) ;endfunctionJ=jacci( Y ,G ,B,P,Q ,e,f,V ,C,D ,B2,n,ph,no) % 该子程序是用来求取jacci矩阵fori=1:nswitch B2( i,4)case 3continuecase 1for j=1:nif j- =i& j- =phX 1=G( i,j) *f( i) -B( i,j) *e( i) ;X 2=G( i,j) *e( i) +B( i,j) *f( i) ;X 3=-X 2;X 4=X 1;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;else if j==i& j~=phX 1=D( i) +G( i,i) *f( i) -B( i,i) *e( i) ;X 2=C( i) +G( i,i) *e( i) +B( i,i) *f( i) ;X 3=C( i) -G( i,i) *e( i) -B( i,i) *f( i) ;X 4=-D( i) +G( i,i) *f( i) -B( i,i) *e( i) ;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;endendendcase 2for j=1:nif j~=i& j~=phX 1=G( i,j) *f( i) -B( i,j) *e( i) ;X 2=G( i,j) *e( i) +B( i,j) *f( i) ;X 3=0;X 4=0;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;else if j==i& j~=phX 1=D( i) +G( i,i) *f( i) -B( i,i) *e( i) ; X 2=C( i) +G( i,i) *e( i) +B( i,i) *f( i) ; X 3=0;X 4=0;p=2*i-1;q=2*j-1;J(p,q) =X 1;m =p+1;J( m ,q) =X 3;q=q+1;J(p,q) =X 2;J( m ,q) =X 4;endendendendendifph~=nJ( no,:) =[ ] ;J( no,:) =[ ] ;J( :,no) =[ ] ;J( :,no) =[ ] ;end。

5节点电力系统牛顿-拉夫逊法潮流计算

（二〇一四年十二月课程论文学校代码： 10128 学号: 20141100304题目：五节点系统计算机潮流计算编程学生姓名：张佳羽学院：电力学院系别：电力系专业：电力系统及其自动化指导教师：郭力萍程序设计% 本程序的功能是用牛顿拉夫逊法进行潮流计算n=input('请输入节点数:n=’);nl=input(’请输入支路数：nl='）；isb=input('请输入平衡母线节点号:isb=’);pr=input('请输入误差精度:pr=’）；B1=input（’请输入由各支路参数形成的矩阵：B1=')；B2=input(’请输入各节点参数形成的矩阵：B2=’）；X=input（'请输入由节点号及其对地阻抗形成的矩阵:X=’）；Y=zeros(n）；e=zeros(1，n);f=zeros（1，n);V=zeros（1，n）；O=zeros(1，n）；S1=zeros(nl);for i=1:nlif B1(i，6)==0p=B1（i，1);q=B1(i,2)；else p=B1（i,2);q=B1(i，1)；endY（p，q)=Y(p，q）—1。

/(B1（i,3）*B1（i,5））;Y(q,p）=Y（p，q）；Y(q，q）=Y(q，q）+1。

/(B1(i,3)＊B1（i,5)^2）+B1（i,4)./2；Y(p,p)=Y（p，p)+1。

/B1（i,3）+B1(i，4)./2;end%求导纳矩阵disp('导纳矩阵Y=')；disp(Y）;G=real(Y）；B=imag(Y）;for i=1：ne(i）=real(B2（i,3))；f(i)=imag（B2（i，3));V(i)=B2(i,4）;endfor i=1:nS（i）=B2（i，1)—B2（i,2）；B(i,i)=B(i，i)+B2（i,5）；endP=real(S);Q=imag(S)；ICT1=0；IT2=1;N0=2*n;N=N0+1;a=0；while IT2~=0IT2=0；a=a+1;for i=1：nif i~=isbC（i）=0；D（i）=0；for j1=1：nC（i）= C(i）+G（i，j1）*e（j1)—B(i，j1)＊f(j1）;D(i)= D（i）+G(i,j1)*f(j1）+B(i,j1)*e(j1);endP1=C（i）*e(i)+f（i)＊D(i);Q1=f（i)＊C(i)—D(i)*e(i);V2=e(i）^2+f(i）^2;if B2（i,6）~=3DP=P（i)-P1;DQ=Q(i)-Q1;for j1=1：nif j1~=isb＆j1~=iX1=-G(i，j1）*e(i）—B（i,j1）*f(i）;X2=B(i,j1）*e(i)—G(i，j1）＊f(i);X3=X2;X4=-X1；p=2＊i-1；q=2*j1—1；J（p，q）=X3;J（p,N)=DQ;m=p+1;J(m,q)=X1；J（m,N）=DP;q=q+1;J(p,q)=X4；J（m，q）=X2;elseif j1==i&j1～=isbX1=—C（i）-G(i，i）*e（i)—B（i,i）*f(i);X2=-D(i)+B(i,i)*e(i)-G（i，i）＊f（i）；X3=D（i)+B（i，i)*e(i）-G（i,i）*f(i）；X4=-C(i)+G（i,i）*e（i）+B(i，i）*f（i）；p=2*i-1;q=2＊j1—1；J(p，q）=X3;J（p,N)=DQ;m=p+1;J(m，q)=X1;J（m，N)=DP；q=q+1;J(p，q)=X4;J(m，q）=X2;endendelseDP=P(i）—P1;DV=V(i）^2—V2；for j1=1：nif j1～=isb＆j1～=iX1=—G(i，j1)*e(i）-B(i，j1)＊f（i）；X2=B（i,j1)*e（i)—G（i,j1）*f（i）；X5=0;X6=0；p=2*i-1;q=2＊j1-1；J(p,q)=X5；J(p，N)=DV；m=p+1；J(m，q)=X1；J(m，N)=DP；q=q+1；J(p,q）=X6;J（m，q）=X2;elseif j1==i＆j1～=isbX1=—C(i）—G（i，i）＊e(i）—B（i，i）*f(i）；X2=-D(i）+B(i，i)＊e（i)-G(i，i)＊f（i)；X5=-2＊e(i）；X6=—2*f（i);p=2*i-1；q=2*j1—1；J（p，q)=X5;J（p,N）=DV;m=p+1；J(m，q)=X1；J(m，N)=DP；q=q+1；J(p,q）=X6；J(m，q）=X2;endendendendend％求雅可比矩阵for k=3：N0k1=k+1;N1=N；for k2=k1：N1J(k,k2)=J(k,k2)./J(k，k）；endJ(k,k)=1;if k~=3;k4=k—1；for k3=3:k4for k2=k1：N1J（k3,k2）= J(k3，k2)—J(k3,k)＊J(k，k2）；endJ(k3，k)=0；endif k==N0，break；endfor k3=k1:N0for k2=k1:N1J(k3，k2)=J(k3,k2)-J(k3，k）*J(k,k2);endJ（k3,k）=0；endelsefor k3=k1：N0for k2=k1:N1J（k3,k2）= J（k3,k2）—J（k3,k)*J（k，k2)；endJ（k3，k）=0;endendendfor k=3：2:N0—1L=(k+1)。

电力系统潮流计算的MATLAB辅助程序设计,潮流计算程序(精编文档).doc

【最新整理，下载后即可编辑】电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。

此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。

现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法。

高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。

lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。

(读者可以到MATPOWER主页下载MATPOWER4.1,然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算)一、高斯-赛德尔法潮流计算使用的程序：高斯-赛德法的具体使用方法读者可参考后面的实例，这里仅介绍各程序的编写格式:lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。

程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。

根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。

对于PV节点，如发电机节点，要提供一个无功功率限定值。

当给定电压过高或过低时，无功功率可能超出功率限定值。

在几次迭代之后(高斯-塞德尔迭代为10次)，需要检查一次发电机节点的无功出力，如果接近限定值，电压幅值进行上下5%的调整，使得无功保持在限定值内。

lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。

并将这些参数放在名为linedata的文件中。

这个程序将阻抗转换为导纳，并得到节点导纳矩阵。

busout：该程序以表格形式输出结果，节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。

matlab潮流计算程序

基于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 )))+imag(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(B2(i,3))*(real(Y(i,j))*imag(B2(j,3))+imag(Y(i,j))*real(B2(j,3)));endendendfor i=1:n %对PV节点的处理，注意这时不可再将h初始化为0if 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 )))+imag(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(B2(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初始化为0if 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 )))+imag(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(B2(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 )))+imag(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(B2(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】亲爱的朋友，上文已完，为感谢你的阅读，特加送另一篇范文，如果下文你不需要，可以下载后编辑删除，谢谢！矿井水灾事故专项应急演练方案1 应急演练目的、意义和目标1.1应急演练目的①评估我矿水灾事故的应急准备状态，发现并修改我矿水灾事故专项应急预案和执行程序中存在的缺陷和不足；②评估我矿在发生水灾事故时的应急能力，识别处理水灾事故的资源需求，澄清相关单位和人员的应急职责，改善水灾事故应急救援中的组织协调问题；③检验应急响应人员对水灾事故应急预案及执行程序的了解程度和实际操作技能；同时，通过调整演练难度，进一步培训和提高应急响应人员的业务素质和能力；④提高全员安全意识。

第六章MATLAB在电力系统潮流计算中的应用实例.

图4-1
2机5节点系统
1)字段baseMVA是一个标量，用来设置基准容量，如100MVA。
2)字段bus是一个矩阵，用来设置电网中各母线参数。
① bus_i用来设置母线编号(正整数)。
② type用来设置母线类型, 1为PQ节点母线, 2为PV节点母线, 3为
平衡(参考)节点母线，4为孤立节点母线。
④ Qmax和Qmin用来设置接入发电机(电源)的无功功率最大、最
小允许值。
⑤ Vg用来设置接入发电机(电源)的工作电压。
⑥ mBase用来设置接入发电机(电源)的功率基准,如果为默认值,就
是baseMVA变量的值。
⑦ status用来设置发电机(电源)工作状态, 1表示投入运行, 0表示退
出运行。
③ Pd和Qd用来设置母线注入负荷的有功功率和无功功率。
④ Gs、Bs用来设置与母线并联电导和电纳。
⑤ baseKV用来设置该母线基准电压。
⑥ Vm和Va用来设置母线电压的幅值、相位初值。
⑦ Vmax和Vmin用来设置工作时母线最高、最低电压幅值。
⑧ area和zone用来设置电网断面号和分区号，一般都设置为1，前
-0．06900
0
1．05000
0．435704
2
0．875129
-0．076619
1．0793
52
0．314147
1．037934
-0．07406
2
1．05000
0．383817
3
0．862445
-0．083228
1．0779
37
0．311611
1．036437
-0．07472
0
1．05000

基于matlab符号运算的电力系统潮流计算

基于matlab符号运算的电力系统潮流计算
电力系统潮流计算是一种用来分析电力系统的负载和电能输送的方法。

在Matlab中，可以使用符号运算的工具来进行电力系统潮流计算，下面是一个基于Matlab符号运算的电力系统潮流计算的步骤示例：
步骤1：定义系统节点参数和系统支路参数，包括节点电压和相角、支路阻抗等。

步骤2：根据系统节点参数和支路参数，建立节点电压方程组和功率方程组，利用功率方程组求解节点电压方程组。

步骤3：利用符号运算的工具，将节点电压方程组和功率方程组转换为矩阵形式，形成潮流计算的数学模型。

步骤4：使用Matlab的符号运算工具求解潮流计算的数学模型，得到节点电压和功率的解。

步骤5：根据节点电压和功率的解，计算系统负载和电能输送情况，进行电力系统的潮流分析和评估。

需要注意的是，电力系统潮流计算涉及大量的矩阵运算和复杂的方程组求解，因此在实际计算中可能需要对问题进行简化和近似处理，以提高计算的效率和精度。

基于Matlab的两机五节点网络潮流仿真计算—牛拉法项目计划书

基于Matlab的两机五节点网络潮流仿真计算—牛拉法计划书第一章电力系统潮流计算概述1.1 潮流计算简介电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。

在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性。

可靠性和经济性。

此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。

所以潮流计算是研究电力系统的一种很重要和基础的计算。

电力系统潮流计算也分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。

利用电子数字计算机进行电力系统潮流计算从50年代中期就已经开始。

在这20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。

对潮流计算的要求可以归纳为下面几点：（1）计算方法的可靠性或收敛性；（2）对计算机内存量的要求；（3）计算速度；（4）计算的方便性和灵活性。

电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题，其解法都离不开迭代。

因此，对潮流计算方法，首先要求它能可靠地收敛，并给出正确答案。

由于电力系统结构及参数的一些特点，并且随着电力系统不断扩大，潮流计算的方程式阶数也越来越高，对这样的方程式并不是任何数学方法都能保证给出正确答案的。

这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。

1.2 潮流计算的意义及其发展（1)在电网规划阶段,通过潮流计算,合理规划电源容量及接入点,合理规划网架,选择无功补偿方案,满足规划水平的大、小方式下潮流交换控制、调峰、调相、调压的要求。

(2)在编制年运行方式时,在预计负荷增长及新设备投运基础上,选择典型方式进行潮流计算,发现电网中薄弱环节,供调度员日常调度控制参考,并对规划、基建部门提出改进网架结构,加快基建进度的建议。

基于Matlab的两机五节点网络潮流仿真计算—牛拉法项目计划书

基于Matlab的两机五节点网络潮流仿真计算—牛拉法项目计划书基于Matlab的两机五节点网络潮流仿真计算—牛拉法计划书第一章电力系统潮流计算概述1.1 潮流计算简介电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。

在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性。

可靠性和经济性。

此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。

所以潮流计算是研究电力系统的一种很重要和基础的计算。

利用电子数字计算机进行电力系统潮流计算从50年代中期就已经开始。

在这20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。

电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题，其解法都离不开迭代。

因此，对潮流计算方法，首先要求它能可靠地收敛，并给出正确答案。

这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。

电力网潮流计算的MATLAB总程序(课程组自编)

电力网潮流计算的MATLAB总程序cleardisp('此程序适用可带有双绕组变压器的任意节点数目的通用网络的潮流计算')n=input('请输入节点数目n; n=');m=input('请输入PQ节点数目m; m=');disp('请将节点按如下规则排序,n个节点中前m个节点为PQ节点,第m+1至n-1个节点为PV节点,第n个节点为平衡节点.')t1=menu('请选择系统的给定参数形式:有名值选1,标幺值选2;','1','2')switch t1case 1disp('给定参数为有名值')t2=menu('选择所求系统是否含有变压器;是选1,否选2.','1','2')switch t2case 1disp('系统含有变压器')SB=input('请输入基准功率(MVA);SB=');disp('请输入系统共有几段电压等级Nv,单位:千伏.')Nv=input('Nv=');disp('为了输入快捷,在以下程序中要求输入"数组"或"矩阵"时请按如下方法输入')disp('输入的数据必须方在方括号内,一行中有多个元素时用空格或逗号隔开,行与行之间用分号或回车隔开.')disp('请为段排序后按序输入各段额定电压VN,单位:千伏;数组形式,如[220,110,10]') VN=input('VN=');[a,c1]=size(VN);while a~=1|c1~=Nvdisp('请为段排序后正确输入各段额定电压VN,单位:千伏;数组形式,如[220,110,10]') VN=input('VN=');[a,c1]=size(VN);enddisp('请按序对应输入各段电压的基准电压VB,数组形式,单位:千伏;如[230,115,10.5]') VB=input('VB=');[a,c2]=size(VB);while a~=1| c2~=Nvdisp('请正确输入基准电压VB,数组形式,单位:千伏;如[230,115,10.5]')VB=input('VB=');[a,c2]=size(VB);endb=1;for a=1:Nvdisp('请输入运行在第'),disp(a),disp('段中的节点号.数组形式,如:[2,3,5]') Nvn=input('请输入:'); %运行在第I,II,III,...段的节点编号[c1,c2]=size(Nvn);for b=1:c2NVN(a,b)=Nvn(b);endendNVNYT=zeros(n);for a=1:ndisp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,其形式如(冒号之内)"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');while (min(F)<=a|max(F)>n)&F~=0disp('请正确输入')disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a) disp('的节点号,数组形式如"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');endif F~=0[f1,f]=size(F); %求F中有几个元素for c=1:fb=F(c);[x1,x2]=find(a==NVN);[x3,x4]=find(b==NVN);if x1==x3disp('请输入节点'),disp(a),disp('与),disp(b),disp('之间的线路阻抗Zij单位:欧姆,电纳B,单位:S')disp('注:可以输入算式,如:(a+j*b)*c/d;复数虚部必须输入j*a,a 代表数值;电纳只输入数值即可')Zij(a,b)=input('Zij=');BB(a,b)=input('B=');while BB(a,b)^2<0disp('输入错误,此值不必输入虚数单位')BB(a,b)=input('B=');endZij(a,b)=Zij(a,b)*SB/VB(x1)^2; % 线路阻抗标幺值BB(a,b)=BB(a,b)*VB(x1)^2/SB; % 对地导纳标幺值Y(a,b)=-1/Zij(a,b);Y(b,a)=Y(a,b); %互导纳y0(a,b)=j*BB(a,b)/2; %对地导纳矩阵元素(虚数)elsedisp('节点'),disp(a),disp('与'),disp(b),disp('不在同一段电压等级,通过变压器相连;请输入它们之间的变压器参数')ST=input('额定容量,单位:兆伏安;S=');Vs100=input('短路电压百分数Vs%=');disp('请输入变压器两侧额定电压VTN;单位:千伏;数组形式,如:[10.5 121]')VTN=input('VTN=');while [1 2]~=size(VTN)disp('输入错误,请输入变压器两侧额定电压,单位:千伏;数组形式,如:[10.5 121]')VTN=input('VTN=');endif VN(x1)<VN(x3)if VTN(1)<VTN(2)ZT(a,b)=j*Vs100*VTN(2)^2*SB/(100*ST*VB(x3)^2); %变压器等效阻抗标幺值kT(b,a)=VTN(2)/VTN(1)/(VB(x3)/VB(x1)); %变压器变比标幺值elseZT(a,b)=j*Vs100*VTN(1)^2*SB/(100*ST*VB(x3)^2) ;%变压器等效阻抗标幺值kT(b,a)=VTN(1)/VTN(2)/(VB(x3)/VB(x1));%变压器变比标幺值 endkT(a,b)=kT(b,a)y0(b,a)=(1-kT(a,b))/ZT(a,b);y0(a,b)=kT(a,b)*(kT(a,b)-1)/ZT(a,b);elseif VTN(1)<VTN(2)ZT(a,b)=j*Vs100*VTN(2)^2*SB/(100*ST*VB(x1)^2);kT(a,b)=VTN(2)/VTN(1)/(VB(x1)/VB(x3));elseZT(a,b)=Vs100*VTN(1)^2*SB/(100*ST*VB(x1)^2);kT(a,b)=VTN(1)/VTN(2)/(VB(x1)/VB(x3));endy0(a,b)=(1-kT(a,b))/ZT(a,b);y0(b,a)=kT(a,b)*(kT(a,b)-1)/ZT(a,b);endY(a,b)=-kT(a,b)/ZT(a,b);Y(b,a)=Y(a,b);YT(a,b)=Y(a,b);YT(b,a)=YT(a,b);endendendY(a,a)=sum(y0(a,:))-sum(Y(a,:));endcase 2disp('系统不含变压器')SB=0;disp('为了输入快捷,在以下程序中要求输入"数组"或"矩阵"时请按如下方法输入')disp('输入的数据必须方在方括号内,一行中有多个元素时用空格或逗号隔开,行与行之间用分号或回车隔开.') ;for a=1:ndisp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,数组形式如"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');while (min(F)<=a|max(F)>n)&F~=0disp('请正确输入')disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a) disp('的节点号,数组形式如"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');endif F~=0[f1,f]=size(F); %求F中有几个元素for c=1:fb=F(c);disp('请输入节点'),disp(a),disp('与'),disp(b),disp('之间的线路阻抗Zij单位:欧姆,电纳B,单位:S')disp('注:可以输入算式,如:(r+j*x)*l;复数虚部必须输入j*a,a 代表数值;电纳只输入数值即可')Zij(a,b)=input('Zij=');BB(a,b)=input('B=');while BB(a,b)^2<0disp('输入错误,此值不必输入虚数单位')BB(a,b)=input('B=');endY(a,b)=-1/Zij(a,b);Y(b,a)=Y(a,b);y0(a,b)=j*BB(a,b)/2;endY(a,a)=sum(y0(a,:))-sum(Y(a,:))endendendcase 2disp('给定参数为标幺值')SB=0;t3=menu('请选择所求系统是否含有变压器,是选1,否选2', '1', '2');switch t3case 1 %对含变支路的输入和处理disp('系统含有变压器')nT=input('请输入含有变压器支路数目nT;nT=');YT=zeros(n); %产生一个n*n的零矩阵for a=1:nTdisp('请输入第'),disp(a),disp('条含有变压器支路的节点编号,其形式如"[x y]"或"[x,y]"(元素间用空格或逗号隔开)')disp('x为低压侧节点,y为高压侧节点')Tij=input('Tij=');x1=Tij(1);x2=Tij(2);disp('请输入节点'),disp(x1),disp('与'),disp(x2),disp('之间的阻抗Z(虚部j*a,a 代表数值)及变压器变比k')Zij(x1,x2)=input('Z=');kT=input('k=');t4=menu('请选择变压器阻抗Z是折算到高压侧的值还是低压侧,高压侧选1,低压侧选2','1','2');switch t4case 1Y(x1,x2)=-kT/Zij(x1,x2); %%此为含变支路互导纳Y(x2,x1)=Y(x1,x2); %对称性y0(x1,x2)=kT*(kT-1)/Zij(x1,x2);y0(x2,x1)=(1-kT)/Zij(x1,x2);case 2Y(x1,x2)=-1/(kT*Zij(x1,x2)); %%此为含变支路互导纳Y(x2,x1)=Y(x1,x2); %对称性y0(x1,x2)=(kT-1)/kT*Zij(x1,x2);y0(x2,x1)=(1-kT)/(kT^2*Zij(x1,x2));endendcase 2 %不含变disp('系统不含有变压器')YT=zeros(n); %产生一个n*n的零矩阵endfor a=1:n %导纳矩阵计算主程序disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,其形式如(冒号之内)"[3 5 6]"或"[3,5,6]",(变压器支路不必再输);若无则输入0')F=input('请输入:');while (min(F)<=a|max(F)>n)&F~=0disp('请正确输入')disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,其形式如(冒号之内)"[3 5 6]"或"[3,5,6]",(变压器支路不必再输);若无则输入0')F=input('请输入:');endif F~=0[f1,f]=size(F); %求F中有几个元素for c=1:fb=F(c);disp('请输入节点'),disp(a),disp('与'),disp(b),,disp('之间的线路阻抗Zij,线路对地支路导纳y0(虚部j*a,a 代表数值)')Zij(a,b)=input('Zij=');y0(a,b)=input('y0=');while y0(a,b)^2>0disp('输入错误,请重新输入;虚数须加虚数单位"j*"')y0(a,b)=input('y0=');endy0(b,a)=y0(a,b);Y(a,b)=-1/Zij(a,b); Y(b,a)=Y(a,b); %互导纳endendY(a,a)=sum(y0(a,:))-sum(Y(a,:));endendYG=real(Y);B=imag(Y);disp('以数组形式(形如"[50 1.5 75]")按节点顺序输入各节点注入的有功功率,无功功率和节点电压') disp('请按节点顺序输入PQ,PV节点的有功功率,可以输入算式形如[3+2,5+8]')Ps=input('Pis=');[x1,x2]=size(Ps);while x1~=1 | x2~=n-1disp('输入错误!请按节点顺序输入PQ,PV节点的有功功率,也可以输入算式形如[3+2,5+8]')Ps=input('Pis=');[x1,x2]=size(Ps);enddisp('请按节点顺序输入PQ节点的无功功率,也可以输入算式形如[3+2,5+8]')Qs=input('Qis=');[x1,x2]=size(Qs);while x1~=1 | x2~=mdisp('输入错误!请按节点顺序输入PQ节点的无功功率,也可以输入算式形如[3+2,5+8]')Qs=input('Qis=');[x1,x2]=size(Qs);enddisp('请按节点顺序输入PV节点及平衡节点电压幅值(从第m+1节点开始)')V=input('V=');[x1,x2]=size(V);while x1~=1 | x2~=n-mdisp('输入错误!请按节点顺序输入PV节点及平衡节点电压幅值')V=input('V=');[x1,x2]=size(V);enddisp('请输入平衡节点电压相角')r=input('r=');[x1,x2]=size(r);while x1~=1 | x2~=1disp('输入错误!请输入平衡节点电压相角')r=input('r=');[x1,x2]=size(r);endt10=menu('是否考虑PV节点无功越界?,是选1,否选2','1','2');switch t10case 1Qv=input('请按节点顺序输入PV节点可提供的最大无功功率,Qmax=');case 2Qv=inf;endswitch t1case 1switch t2case 1Ps=Ps/SB;Qs=Qs/SB;Vsd=ones(1,m);V=[Vsd,V];rsd=zeros(1,n-1);r=[rsd,r];for a=m+1:n-1[x1,x2]=find(a==NVN);V(a)=V(a)/VB(x1);endQv=Qv/SB;case 2[c1,c2]=size(V);Vsd=V(c2)*ones(1,m);V=[Vsd,V];rsd=zeros(1,n-1);r=[rsd,r];endcase 2Vsd=ones(1,m);V=[Vsd,V]rsd=zeros(1,n-1);r=[rsd,r]endk=0;jdgd=input('输入要求精度;jd=');jd=inf;while jd>jdgdfor a=1:n-1for b=1:nPn(b)=V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));Qn(b)=V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));endPi(a)=sum(Pn(:)); %PiQH(a)=sum(Qn(:));if a<=mQi(a)=sum(Qn(:)); %求QiendendDPi=Ps-Pi;DQi=Qs-Qi; %求得DPifor a=1:n-1for b=1:n-1if a~=bH(a,b)=-V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));if b<=mN(a,b)=-V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));endif a<=mK(a,b)=V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));if b<=mL(a,b)=-V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));endendelseH(a,a)=V(a)^2*B(a,a)+QH(a);if b<=mN(a,a)=-V(a)^2*G(a,a)-Pi(a);endif a<=mK(a,a)=V(a)^2*G(a,a)-Pi(a);if b<=mL(a,a)=V(a)^2*B(a,a)-QH(a);endendendendendfor a=1:n-1for b=1:n-1if a<=m & b<=mJJ(2*a-1,2*b-1)=H(a,b);JJ(2*a-1,2*b)=N(a,b);JJ(2*a,2*b-1)=K(a,b); JJ(2*a,2*b)=L(a,b);endif a<=m & b>mJJ(2*a-1,m+b)=H(a,b);JJ(2*a,m+b)=K(a,b);endif a>m & b<=mJJ(m+a,2*b-1)=H(a,b);JJ(m+a,2*b)=N(a,b);endif a>m & b>mJJ(m+a,m+b)=H(a,b);endendendJJfor a=1:n-1if a<=mDPQ(2*a-1)=DPi(a);DPQ(2*a)=DQi(a);elseDPQ(a+m)=DPi(a);endendDPQDrV=-inv(JJ)*(DPQ');for a=1:n-1if a<=mDr(a)=DrV(2*a-1);DVV(a)=DrV(2*a);elseDr(a)=DrV(a+m);endendDrj=Dr/pi*180for a=1:mVd2(a,a)=V(a);endDV=Vd2*(DVV')for a=1:n-1r(a)=r(a)+Dr(a);if a<=mV(a)=V(a)+DV(a);endendjd=max(abs(DPQ))switch t10case 1for a=m+1:n-1if Qv(a-m)<QH(a)disp(a),disp('节点无功越限')pvwgyj=1;m=m+1;for b=1:nif Y(a,b)~=0Yc(a,b)=Y(a,b);Y(a,b)=Y(m,b);Y(m,b)=Yc(a,b);Y(b,a)=Y(a,b);Y(b,m)=Y(m,b);endendQ(a)=Q(m);Qs(m)=Qv(a);Ps(m)=P(a);P(a)=P(m)Vc(a)=V(a);V(a)=V(m);V(m)=Vc(a);rc(a)=r(a);r(a)=r(m);r(m)=rc(a) elsepvwgyj=0endendendk=k+1;k=k+1;enddisp('迭代次数')k=k-1format short edisp('节点功率不平衡量')DPiDQiformat shortdisp('各节点电压相角; 单位:度;')rd=r/pi*180for a=1:nfor b=1:nPn(b)=V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));Qn(b)=V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));endPN=sum(Pn(:)); %P(n)QN=sum(Qn(:)); %Q(n)endPi(n)=PN;QH(n)=QN;Sn=Pi+QH*j;for a=1:nU(a)=V(a)*(cos(r(a))+j*sin(r(a)));endfor a=1:nfor b=1:nS(a,b)=V(a)^2*conj(y0(a,b))+U(a)*(conj(U(a))-conj(U(b)))*conj(-Y(a,b));endendif SB==0disp('各节点电压幅值')Vdisp('各节点的注入功率')Sndisp('各支路的功率分布')Selsefor a=1:n[x1,x2]=find(a==NVN);V(a)=V(a)*VB(x1);enddisp('各节点电压幅值,单位:千伏')Vdisp('各节点的注入功率,单位:兆伏安')Sn=Sn/SBdisp('各支路的功率分布,单位:兆伏安')S=S/SBend。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

N(m,n)=G(m,n)*E(m)+B(m,n)*F(m);
J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m);
L(m,n)=G(m,n)*F(m)-B(m,n)*E(m);
JJ(2*m-1,2*n-1)=H(m,n);
%设S=P+Bj；
S(1)=0.20+0.20i;S(2)=-0.45-0.15i;S(3)=-0.40-0.05i;S(4)=-0.60-0.10i;
P=real(S);Q=imag(S);
k=0;precision=1;
N1=4;
while precision > 0.00001
N(m,m)=sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m))+2*G(m,m)*E(m);
J(m,m)=-2*B(m,m)*F(m)+sum(Ai)-(G(m,m)*E(m)-B(m,m)*F(m));
L(m,m)=-2*B(m,m)*E(m)-(sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m)));
disp('节点导纳矩阵YB=');
disp(YB);
%计算各节点功率的不平衡量设U=E+jF ；Y=G+Bj；
E(1)=1.00;E(2)=1.00;E(3)=1.00;E(4)=1.00;
F(1)=0;F(2)=0;F(3)=0;F(4)=0;
G=real(YB);B=imag(YB);
yb25=-1.250+3.750j;yb21=-1.667+5.000j;yb22=12.917-38.750j;yb23=-10.000+30.000j;yb24=0.000-0.000j;
yb35=0.000-0.000j;yb31=-1.667+5.000j;yb32=-10.000+30.000j;yb33=12.917-38.750j;yb34=-1.250+3.750j;
%设由P,Q的改变量组成的8×1矩阵为PQ，由E,F的改变量组成的8×1矩阵为dU
for m=1:N1
PQ(2*m-1)=dP(m); PQ(2*m)=dQ(m);
end
dU=inv(JJ)*PQ';
precision=max(abs(dU));
yb45=0.000-0.000j;yb41=-2.500+7.500j;yb42=0.000-0.000j;yb43=-1.250+3.750j;yb44=3.750-11.250j;
YB=[yb11 yb12 yb13 yb14 yb15; yb21 yb22 yb23 yb24 yb25 ;yb31 yb32 yb33 yb34 yb35; yb41 yb42 yb43 yb44 yb45 ;yb51 yb52 yb53 yb54 yb55];
end
%设雅可比矩阵为JJ，以下语句用来实现雅可比矩阵中对角线上元素H N J L 的排列
for m=1:N1
JJ(2*m-1,2*m-1)=H(m,m);
JJ(2*m-1,2*m)=N(m,m);
JJ(2*m,2*m-1)=J(m,m);
end
for m=1:N1
for n=1:N1+1
%计算Hij Nij Jij Lij
H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m);
N(m,n)=G(m,n)*E(m)+B(m,n)*F(m);
Ai(n)=G(m,n)*E(n)-B(m,n)*F(n);
end
%sum(Ai)，sum(Bi)用于实现公式中的sigerma从j到n的求和;
H(m,m)=sum(Bi)-(B(m,m)*E(m)+G(m,m)*F(m))+2*G(m,m)*F(m);
E(5)=1.06;F(5)=0;
for m=1:N1
for n=1:N1+1
%计算Pi,Qi,设Pi=Pt;Qi=Qt
Pt(n)=(E(m)*(G(m,n)*E(n)-B(m,n)*F(n))+F(m)*(G(m,n)*F(n)+B(m,n)*E(n)));
J(m,n)=-B(m,n)*F(m)-G(m,n)*E(m);
L(m,n)=G(m,n)*F(m)-B(m,n)*E(m);
end
end
r m=1:N1
for n=1:N1+1
Bi(n)=G(m,n)*F(n)+B(m,n)*E(n);
JJ(2*m,2*m)=L(m,m);
end
%以下语句用于实现雅可比矩阵非对角线上元素的排列
for m=1:N1
for n=1:N1
if m==n
else
H(m,n)=-B(m,n)*E(m)+G(m,n)*F(m);
end
%1.形成节点导纳矩阵,
yb55=6.250-18.750j;yb51=-5.000+15.000j;yb52=-1.250+3.750j;yb53=0.000-0.000j;yb54=0.000-0.000j;
yb15=-5.000+15.000j;yb11=10.834-32.500j;yb12=-1.667+5.000j;yb13=-1.667+5.000j;yb14=-2.500+7.500j;
for n=1:N1
F(n)=F(n)+dU(2*n-1);
E(n)=E(n)+dU(2*n);
end
for n=1:N1+1
U(n)=E(n)+(F(n))*j;
end
k=k+1;
k-1, dU=dU',PQ,U
Qt(n)=(F(m)*(G(m,n)*E(n)-B(m,n)*F(n))-E(m)*(G(m,n)*F(n)+B(m,n)*E(n)));
end
%设P,Q的改变量为dP,dQ
dP(m)=P(m)-sum(Pt);
dQ(m)=Q(m)-sum(Qt);
JJ(2*m-1,2*n)=N(m,n);
JJ(2*m,2*n-1)=J(m,n);
JJ(2*m,2*n)=L(m,n);
end
end
end