潮流计算(matlab)实例计算

合集下载

潮流计算(matlab)实例计算

潮流计算(matlab)实例计算

潮流例题:根据给定的参数或工程具体要求(如图),收集和查阅资料;学习相关软件(软件自选:本设计选择Matlab进行设计)。

2.在给定的电力网络上画出等值电路图。

3.运用计算机进行潮流计算。

4.编写设计说明书。

一、设计原理1.牛顿-拉夫逊原理牛顿迭代法是取x0 之后,在这个基础上,找到比x0 更接近的方程的跟,一步一步迭代,从而找到更接近方程根的近似跟。

牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程f(x) = 0 的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根。

电力系统潮流计算,一般来说,各个母线所供负荷的功率是已知的,各个节点电压是未知的(平衡节点外)可以根据网络结构形成节点导纳矩阵,然后由节点导纳矩阵列写功率方程,由于功率方程里功率是已知的,电压的幅值和相角是未知的,这样潮流计算的问题就转化为求解非线性方程组的问题了。

为了便于用迭代法解方程组,需要将上述功率方程改写成功率平衡方程,并对功率平衡方程求偏导,得出对应的雅可比矩阵,给未知节点赋电压初值,一般为额定电压,将初值带入功率平衡方程,得到功率不平衡量,这样由功率不平衡量、雅可比矩阵、节点电压不平衡量(未知的)构成了误差方程,解误差方程,得到节点电压不平衡量,节点电压加上节点电压不平衡量构成新的节点电压初值,将新的初值带入原来的功率平衡方程,并重新形成雅可比矩阵,然后计算新的电压不平衡量,这样不断迭代,不断修正,一般迭代三到五次就能收敛。

牛顿—拉夫逊迭代法的一般步骤:(1)形成各节点导纳矩阵Y。

(2)设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。

(3)计算各个节点的功率不平衡量。

(4)根据收敛条件判断是否满足,若不满足则向下进行。

(5)计算雅可比矩阵中的各元素。

(6)修正方程式个节点电压(7)利用新值自第(3)步开始进入下一次迭代,直至达到精度退出循环。

(8)计算平衡节点输出功率和各线路功率2.网络节点的优化1)静态地按最少出线支路数编号这种方法由称为静态优化法。

matlab潮流计算仿真方法

matlab潮流计算仿真方法

matlab潮流计算仿真方法
MATLAB 是一种强大的编程语言和环境,可用于执行各种仿真和计算任务,包括电力系统潮流计算。

以下是一个简单的 MATLAB 潮流计算仿真方法的
示例:
1. 定义系统参数:首先,你需要定义电力系统的参数,如发电机、负荷、变压器等。

这些参数通常包括额定电压、额定功率、电抗、电阻等。

2. 建立系统模型:使用这些参数,你可以在 MATLAB 中建立电力系统的模型。

这通常涉及到定义节点和支路,以及为它们分配相应的参数。

3. 编写潮流计算函数:接下来,你需要编写一个用于执行潮流计算的函数。

这个函数应该能够接收系统的模型和参数,并返回计算出的潮流结果,如电压、电流、功率等。

4. 运行仿真:最后,你可以运行仿真并调用你编写的潮流计算函数。

这将返回计算出的潮流结果,你可以使用这些结果进行进一步的分析或可视化。

这只是一个简单的示例,实际上在编写 MATLAB 潮流计算仿真方法时可能
需要考虑更多因素,例如系统的约束条件、初始条件、迭代算法的收敛性等。

如果你需要具体的 MATLAB 代码示例或更详细的指导,我建议你查阅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实现潮流计算程序

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 极坐标求解潮流计算

电力系统分析大作业——潮流计算班级:XXX姓名:XXX学号:XXX程序清单:%直角坐标法求解潮流计算clear;close;clc;Branch=[1 2 0.1+0.4j 0.01528j 1 03 1 0.3j inf 1.1 11 4 0.12+0.50j 0.01920j 1 02 4 0.08+0.40j 0.01413j 1 0]%Branch矩阵:1、支路首端号;2、支路末端号;3、支路阻抗;4、支路对地导纳;5、支路的变化;6、支路首端处于K侧为1,1侧为0Y=zeros(4); %节点导纳矩阵for i=1:4if Branch(i,6)==0 %不含变压器的支路j=Branch(i,1);k=Branch(i,2);Y(j,k)=Y(j,k)-1/Branch(i,3);Y(k,j)=Y(j,k);Y(j,j)=Y(j,j)+1/Branch(i,3)+Branch(i,4);Y(k,k)=Y(k,k)+1/Branch(i,3)+Branch(i,4);elsej=Branch(i,1);k=Branch(i,2);Y(j,k)=Y(j,k)-Branch(i,5)/Branch(i,3);Y(k,j)=Y(j,k);Y(j,j)=Y(j,j)+1/Branch(i,3);Y(k,k)=Y(k,k)+Branch(i,5)^2/Branch(i,3);endenddisp('节点导纳矩阵');YG=real(Y);B=imag(Y);V=[1;1;1.1;1.05];%给定V的初始计算值e=real(V);f=imag(V);disp('节点电压的实部:')edisp('节点电压的虚部:')fdisp('节点注入有功功率:')Ps=[-0.3;-0.55;0.5;0]disp('节点注入无功功率:')Qs=[-0.18;-0.13;0;0]%由各节点电压向量(状态变量)可得各节点注入功率:P=e.*(G*e-B*f)+f.*(G*f+B*e);Q=f.*(G*e-B*f)-e.*(G*f+B*e);del_W=[-0.30-P(1);-0.18-Q(1);-0.55-P(2);-0.13-Q(2);0.5-P(3);1.1^2-e(3)^2-f(3)^2];n=0;%辅助循环计数变量while (any(del_W>1e-5))&(n<5)n=n+1;disp('迭代次数:');disp(n)%--------------------------------------------------------%修正方程式:% [-0.30-P(1) ] [del_e(1)]% [-0.18-Q(1) ] [del_f(1)]% [-0.55-P(2) ] [del_e(2)]% [-0.13-Q(2) ] + J * [del_f(2)] = 0 % [0.5-P(3) ] [del_e(3)]% [1.1^2-e(3)^2-f(3)^2] [del_f(3)]%雅克比矩阵为:J=[P1_e1 P1_f1 P1_e2 P1_f2 P1_e3 P1_e3;% Q1_e1 Q1_f1 Q1_e2 Q1_f2 Q1_e3 Q1_f3% P2_e1 P2_f1 P2_e2 P2_f2 P2_e3 P2_f3% Q2_e1 Q2_f1 Q2_e2 Q2_f2 Q2_e3 Q2_f3% P3_e1 P3_f1 P3_e2 P3_f2 P3_e3 P3_f3% V3_e1 V3_f1 V3_e2 V3_f2 V3_e3 V3_f3]%--------------------------------------------------------%-------------------求雅可比矩阵的参数-----------------------%GeBf=G*e-B*f;%辅助计算函数GfBe=G*f+B*e;%辅助计算函数for i=1:2for j=1:3if i==jP_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);elseP_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);Q_f(i,j)=-P_e(i,j);P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);Q_e(i,j)=P_f(i,j);endendendfor j=1:2P_e(3,j)=-G(3,j)*e(3)-B(3,j)*f(3);P_f(3,j)=B(3,j)*e(3)-G(3,j)*f(3);V_e(3,j)=0;V_f(3,j)=0;endP_e(3,3)=-GeBf(3)-G(3,3)*e(3)-B(3,3)*f(3);P_f(3,3)=-GfBe(3)-G(3,3)*f(3)+B(3,3)*e(3);V_e(3,3)=-2*e(3);V_f(3,3)=-2*f(3);%----------------------------------------------------%%--------------------求出雅克比矩阵-------------------% J=zeros(6,6);for i=1:3for j=1:3J(2*i-1,2*j-1)=P_e(i,j);J(2*i-1,2*j)=P_f(i,j);endendfor i=1:2for j=1:3J(2*i,2*j-1)=Q_e(i,j);J(2*i,2*j)=Q_f(i,j);endendfor j=1:3J(6,2*j-1)=V_e(3,j);J(6,2*j)=V_f(3,j);enddisp('雅可比矩阵为:')J%-----------------------------------------------------%%------------解修正方程并得出修正后的电压向量-----------% del_V=-J\del_W;clear i;for k=1:3V(k)=V(k)+del_V(2*k-1)+i*del_V(2*k);enddisp('节点电压为:')disp(V)%-----------------------------------------------------%%------------------计算节点不平衡量--------------------%e=real(V);f=imag(V);P=e.*(G*e-B*f)+f.*(G*f+B*e);Q=f.*(G*e-B*f)-e.*(G*f+B*e);del_W=[-0.30-P(1);-0.18-Q(1);-0.55-P(2);-0.13-Q(2);0.5-P(3);1.1^2-e(3)^2-f(3)^2];disp('节点不平衡量为:')disp(del_W)end%------------------------最终结论-------------------------%disp('最终各节点电压幅值为:')disp(abs(V))disp('最终各节点电压相角(度)为:')disp(180*angle(V)/pi)disp('最终各节点注入功率:')S=P+i*Qdisp('最终各节点注入有功功率为:')Pdisp('最终各节点注入无功功率为:')Q运行结果:Branch =1.00002.0000 0.1000 + 0.4000i 0 + 0.0153i 1.0000 03.0000 1.0000 0 + 0.3000i Inf 1.1000 1.00001.0000 4.0000 0.1200 + 0.5000i 0 + 0.0192i1.0000 02.0000 4.0000 0.0800 + 0.4000i 0 + 0.0141i 1.0000 0节点导纳矩阵Y =1.0421 - 8.2429i -0.5882 +2.3529i 0 +3.6667i -0.4539 + 1.8911i-0.5882 + 2.3529i 1.0690 - 4.7274i 0 -0.4808 + 2.4038i0 + 3.6667i 0 0 - 3.3333i 0-0.4539 + 1.8911i -0.4808 + 2.4038i 0 0.9346 - 4.2616i节点电压的实部:e =1.00001.00001.10001.0500节点电压的虚部:f =节点注入有功功率:Ps =-0.3000-0.55000.5000节点注入无功功率:Qs =-0.1800-0.1300迭代次数:1雅可比矩阵为:J =-1.0194 -8.3719 0.5882 2.3529 0 3.6667 -8.1138 1.0648 2.3529 -0.5882 3.6667 00.5882 2.3529 -1.0450 -4.8770 0 02.3529 -0.5882 -4.5778 1.0930 0 00 4.0333 0 0 0 -3.66670 0 0 0 -2.2000 0节点电压为:0.9935 - 0.0088i0.9763 - 0.1078i1.1000 + 0.1267i1.0500节点不平衡量为:-0.0013-0.0028-0.0135-0.05470.0030-0.0160迭代次数:2雅可比矩阵为:J =-0.8091 -8.3613 0.6052 2.3325 0.0324 3.6429 -7.9992 1.4071 2.3325 -0.6052 3.6429 -0.03240.8280 2.2338 -1.0190 -4.6364 0 02.2338 -0.8280 -4.3641 2.0878 0 0-0.4644 4.0333 0 0 -0.0324 -3.64290 0 0 0 -2.2000 -0.2533节点电压为:0.9847 - 0.0086i0.9590 - 0.1084i1.0924 + 0.1289i1.0500节点不平衡量为:-0.0000-0.0000-0.0003-0.00110.0001-0.0001迭代次数:3雅可比矩阵为:J =-0.7940 -8.2936 0.5995 2.3120 0.0315 3.6107 -7.9228 1.4000 2.3120 -0.5995 3.6107 -0.03150.8191 2.1927 -0.9865 -4.6144 0 02.1927 -0.8191 -4.2210 2.0885 0 0-0.4728 4.0056 0 0 -0.0315 -3.61070 0 0 0 -2.1849 -0.2579。

matlab电力系统潮流计算程序

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

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

现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和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 粗略程序

%========================================================================== %========================================================================== %========================================================================== %潮流计算MATLAB 粗略程序 C.zhou 2009.3.27%========================================================================== %========================================================================== %========================================================================== %creat a new_datat=0;s=0;r=0;w=0;number=input('How many node are there=');% Convert Pq to a new arrayfor ii=1:numberif data(ii,4)==1t=t+1;for jj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;s=s+1; %record the number of the PQ node end;end;%Convert pv to a new arrayfor ii=1:numberif data(ii,4)==2t=t+1;for jj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;r=r+1; %record the number of the PV node end;end;%Convert set_v to a new arrayfor ii=1:numberif data(ii,4)==3t=t+1;for jj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;w=w+1;end;end;%creat a new_data2[x,y]=size(data2)for ii=1:xfor jj=1:2for mm=1:numberif data2(ii,jj)==a(1,mm)new_data2(ii,jj)=mm;end;end;end;end;for ii=1:xfor jj=3:14new_data2(ii,jj)=data2(ii,jj);end;end;%creat a YY=zeros(number,number);YY=zeros(number,number);yy=zeros(number,number);for ii=1:x% for jj=1:14iii=new_data2(ii,1);jjj=new_data2(ii,2);if new_data2(ii,5)==2sub=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6 ))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;Y(iii,jjj)=-sub./new_data2(ii,14);YY(iii,jjj)=sub./new_data2(ii,14);Y(jjj,iii)=-sub/new_data2(ii,14);YY(jjj,iii)=sub./new_data2(ii,14);yy(iii,jjj)=(1.-new_data2(ii,14))./(new_data2(ii,14).*new_data2(ii,14)).*sub;yy(jjj,iii)=(new_data2(ii,14)-1)./(new_data2(ii,14)).*sub;elseY(iii,jjj)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2 (ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;YY(iii,jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data 2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;Y(jjj,iii)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2 (ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;YY(jjj,iii)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;yy(iii,jjj)=new_data2(ii,8)./2.*i;yy(jjj,iii)=new_data2(ii,8)./2.*i;end;%end;end;for iii=1:numberY(iii,iii)=0;end;%for ii=1:x% for jj=1:14for iii=1:numberfor jj=1:number% if iii~=jjY(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);% end;end;end;%creat B, Gfor ii=1:numberfor jj=1:numberG(ii,jj)= real(Y(ii,jj));B(ii,jj)= imag(Y(ii,jj));end;end;%creat Initial_P Initial_Q Initial_Vfor ii=1:(s+r)set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7))./100;end;for ii=1:s;set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8))./100;end;for ii=1:rset_V(ii,1)=new_data1(ii+s,12).*new_data1(ii+s,12);%try to modify for sike of correcting end;Initial_p_q_v=[set_P;set_Q;set_V];disp(Initial_p_q_v);%creat Initial_e,Initial_ffor ii=1:number-1e(ii,1)=1;f(ii,1)=0.0;%change f to test used to be 1.0end;e(number,1)=new_data1(number,12);f(number,1)=0;% e(64,1)=0.88;%test 118ieee% f(64,1)=0.39395826829394;% f(14,1)=0;% e(10,1)=1.045;%e(11,1)=1.01;%e(12,1)=1.07;%e(13,1)=1.09;%//////////////////////////////////////////////////////////////////////////%/////////////////////////////////////////////////////////////////////////%//////////////////////////////////////////////////////////////////////////%//////////////////////////////////////////////////////////////////////////% Start NEWTOWN CALULATIONfor try_time=1:25%Creat every node consume P Q and Un=s;m=r;for ii=1:(n+m)sum1=0;for jj=1:(n+m+1)sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1))+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));end;p(ii,1)=sum1;end;for ii=1:nsum2=0;for jj=1:(n+m+1)sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1))-e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));end;q(ii,1)=sum2;end;disp('q=');disp(q);u=zeros((n+m),1);for ii=(n+1):(n+m)u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;for ii=n+1:(n+m)extra_u((ii-n),1)=u(ii,1);end;disp('extra_u=');disp(extra_u);sum=[p;q;extra_u];disp(sum)disp(s);disp(p);%creat Jacobiandisp(n);disp(m);for ii=1:(n+m)for jj=1:(n+m)if (ii~=jj)PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);elsess=0;qq=0;for num=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1PE(ii,jj)=-qq-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1end;end;end;copy=3.14159;disp('================copy================')for ii=1:nfor jj=1:m+nif (ii~=jj)QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1elsess=0;qq=0;for num=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1end;end;end;%disp('QF');%disp(QF);%disp('QE');%disp(QE);UE=zeros((n+m),(n+m));UF=zeros((n+m),(n+m));for ii=n+1:n+mfor jj=1:(n+m)if (ii~=jj)UE(ii,jj)=0;UF(ii,jj)=0;elsess=0;qq=0;for num=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;UF(ii,jj)=-2.*f(ii,1);UE(ii,jj)=-2.*e(ii,1);end;end;end;for ii=(n+1):(n+m)for jj=1:(n+m)extra_UE((ii-n),jj)=UE(ii,jj);extra_UF((ii-n),jj)=UF(ii,jj);end;end;%disp('extra_UE');%disp(extra_UE);%disp('extra_Uf');%disp(extra_UF);Jacobian=[PF,PE;QF,QE;extra_UF,extra_UE];%disp('Jacobian=');%disp(Jacobian);%creat substract resultsubstract_result=Initial_p_q_v-sum;%disp('substract_result');%disp(substract_result);%calculate delta_f_edelta_f_e=-inv(Jacobian)*substract_result; %disp(delta_f_e);for ii=1:number-1;f(ii,1)=f(ii,1)+delta_f_e(ii,1);e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1); end;if max(substract_result)<1e-4break;end ;end;%disp('substract_result');%disp(substract_result);%disp('e=');%disp(e);%disp('f=');%disp(f);for ii=1:numberuuu(ii,1)= e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);U_RESULT(ii,1)=sqrt(uuu(ii,1));end;for ii=1:numberfor jj=1:numberif ii==a(1,jj)Old_Uresult(ii,1)=U_RESULT(jj,1)end;end;end;for ii=1:numberOld_Uresult(ii,2)=ii;end;%disp('U_result');%disp(U_RESULT);disp('====================================='); disp('The last result is :')disp('===========U===================BUS-NO.'); disp('U=')disp(Old_Uresult);%calculate the anglePI=3.141592for ii=1:numberAngle(ii,1)=atan(f(ii,1)./e(ii,1))./PI*180;end;for ii=1:numberfor jj=1:numberif ii==a(1,jj)Old_Angle(ii,1)=Angle(jj,1);Old_Angle(ii,2)=ii;end;end;end;disp('=======Angle===================BUS-NO.'); disp('Angle=');disp(Old_Angle);disp('=====Try-times=======================') disp('Try-times=')disp(try_time);%disp('p====================');%disp(p);% for jj=1:number% if a(1,jj)==118% man=jj% end;%end;%disp('man=========');%disp(man)sum4=0;for jj=1:numberY_conj(number,jj)=conj(Y(number,jj));sum4=sum4+Y_conj(number,jj).*(e(jj,1)-f(jj,1)*i);end;%sum4=sum4*e(number,1);disp('===============Balance P Q=========BUS-NO');%disp(sum4);Blance_Q(1,1)=imag(sum4)*100;Blance_Q(1,2)=a(1,number);Blance_P(1,1)=real(sum4)*100;Blance_P(1,2)=a(1,number);disp('Q Of the Balance node= ');disp(Blance_Q);disp('P Of the Balance node= ')disp(Blance_P);disp('=================================BUS-NO');%calculate the Q of the P-V nodeQ_PV_node=zeros(number,2);Y_conj=conj(Y);for ii=(s+1):(s+r)for jj=1:numberQ_PV_node(ii,1)=Q_PV_node(ii,1)+(e(ii,1)+f(ii,1)*i).*(Y_conj(ii,jj).*(e(jj,1)-f(jj,1)*i));end;end;for ii=(s+1):(s+r);Q_PV_node(ii,1)=Q_PV_node(ii,1).*100+new_data1(ii,8)*i;end;disp('This program is from /breadwinner') ;for ii=1:numberold_number=a(1,ii);Q_PV_node_old(old_number,1)=Q_PV_node(ii,1);end;for ii=1:numberQ_PV_node_old(ii,1)=imag(Q_PV_node_old(ii,1));end;for ii=1:numberQ_PV_node_old(ii,2)=ii;end;disp('Q gen=');disp(Q_PV_node_old);。

matlab潮流计算程序

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实验电力系统潮流计算

matlab实验电力系统潮流计算

matlab实验电力系统潮流计算电力系统潮流计算是电力系统运行分析的基础,它通过计算电力系统中各节点的电压幅值和相角,以及各支路的功率和电流,来研究电力系统的稳态工作状态。

本文将介绍潮流计算的原理及其在电力系统中的应用。

潮流计算的基本原理是基于电力系统节点间的功率平衡方程,即节点注入功率等于节点负荷消耗功率加上节点发电机注入功率和节点之间传输输电功率的代数和。

潮流计算通常分为直流潮流计算和交流潮流计算两种方法。

直流潮流计算是指忽略电流相位差的计算方法。

在直流潮流计算中,电力系统的节点电压幅值和相角可以用复数来表示,节点注入功率和节点负荷消耗功率也采用复功率的形式。

直流潮流计算的基本方程为:P+iQ = V*(Gcosθ+Bsinθ)其中,P和Q分别表示节点注入有功功率和无功功率,V表示节点电压幅值,θ表示节点电压相角,G和B分别表示节点导纳矩阵的实部和虚部。

交流潮流计算考虑了电流相位差的影响,是更为准确的潮流计算方法。

交流潮流计算通常采用牛顿-拉夫逊法(Newton-Raphson method)进行迭代求解。

该方法以功率不平衡最小为目标,通过迭代计算更新节点电压幅值和相角,直到收敛为止。

潮流计算在电力系统运行和规划中具有重要的应用价值。

首先,潮流计算可以用来评估电力系统的稳态工作状态,包括节点电压幅值和相角、支路功率和电流等信息。

通过分析潮流计算结果,可以发现电力系统中潜在的潮流瓶颈和潮流分布情况,为电网调度和运行提供参考依据。

其次,潮流计算可以用来优化电力系统的运行和规划。

通过分析潮流计算结果,可以确定潮流分布不均衡的节点和支路,进而优化电力系统的输电和变电容量配置,提高电力系统的可靠性、经济性和稳定性。

此外,潮流计算还可用于电力系统的故障分析和重构,对于故障点的电压幅值、相角以及故障后的支路功率和电流进行分析,有助于电力系统故障的定位和恢复。

总的来说,电力系统潮流计算是电力系统运行分析的重要工具,通过计算电力系统中各节点的电压幅值和相角,以及各支路的功率和电流,可以评估电力系统的稳态工作状态,优化电力系统的运行和规划,实现电力系统的安全、稳定和高效运行。

前推回代法潮流计算matlab

前推回代法潮流计算matlab

前推回代法潮流计算matlab随着计算机技术的不断发展,代数计算在科学计算和工程领域中越来越重要。

代数计算能够快速精确地求解方程组、积分、微分等数学问题,提高了计算效率,也为科研工作提供了有力的支持。

在代数计算中,代数推导是一种重要的方法,它可以通过代数等式的变换推导出新的等式,从而简化复杂的计算过程。

代数推导方法有很多,其中推回代法是一种常用的方法。

推回代法是指将已知的一些等式反过来使用,以便求解未知的变量。

在MATLAB中,使用推回代法可以快速地求解一些复杂的方程组、积分和微分等问题。

以求解方程组为例,假设有以下方程组:x + y = 52x – y = 1使用推回代法求解该方程组,可以先将第一个方程式化为y = 5 – x,然后将其代入第二个方程式中,得到:2x – (5 – x) = 1化简得到3x = 6,即x = 2。

将x = 2代入第一个方程式中,得到y = 3。

因此,该方程组的解为x = 2,y = 3。

在MATLAB中,可以使用syms命令定义符号变量,使用solve命令求解方程组。

例如,求解上述方程组的代码如下:syms x ysolve(x + y == 5, 2*x - y == 1)运行代码后,MATLAB会输出方程组的解:x = 2,y = 3。

除了求解方程组,推回代法在积分、微分等计算中也有广泛的应用。

例如,对于以下积分:∫(x2 + 2x + 1)dx可以使用推回代法将其化简为∫(x + 1)2dx,然后使用二次积分公式求解。

推回代法是一种常用的代数推导方法,在MATLAB中也有广泛的应用。

使用推回代法可以简化复杂的计算过程,提高计算效率,为科研工作提供有力的支持。

基于matlab的潮流计算

基于matlab的潮流计算

课程设计任务书第学期学院:信息商务学院专业:电气工程及其自动化学生姓名:学号课程设计题目:基于MATLAB的潮流计算起迄日期:月日~月日课程设计地点:电气工程系工程实验室指导教师:系主任:课程设计任务书课程设计任务书1 引言潮流计算是电力系统最基本最常用的计算。

根据系统给定的运行条件,网络接线及元件参数,通过潮流计算可以确定各母线的电压(幅值和相角),各支路流过的功率,整个系统的功率损耗。

潮流计算是实现电力系统安全经济发供电的必要手段和重要工作环节。

因此,潮流计算在电力系统的规划计算,生产运行,调度管理及科学计算中都有着广泛的应用。

潮流计算在数学上是多元非线性方程组的求解问题,牛顿—拉夫逊Newton-Raphson法是数学上解非线性方程组的有效方法,有较好的收敛性。

运用电子计算机计算一般要完成以下几个步骤:建立数学模型,确定解算方法,制订计算流程,编制计算程序。

1.1 设计任务与要求利用MATLAB/SIMULINK仿真软件,搭建多电压等级电力网络潮流计算模型。

掌握各种潮流计算方法,掌握MATLAB/SPS工具箱使用方法及其在潮流计算中的应用。

MATLAB作为先进的仿真软件,在电力系统潮流计算方面,尤其是矩阵运算上具有突出的优秀性。

本设计要求利用MATLAB编写程序,输入数据并完成潮流分布计算。

1.2实用价值与理论意义在电力系统的正常运行中,随着用电负荷的变化和系统运行方式的改变,网络中的损耗也将发生变化。

要严格保证所有的用户在任何时刻都有额定的电压是不可能的,因此系统运行中个节点出现电压的偏移是不可避免的。

为了保证电力系统的稳定运行,要进行潮流调节。

随着电力系统及在线应用的发展,计算机网络已经形成,为电力系统的潮流计算提供了物质基础。

电力系统潮流计算是电力系统分析计算中最基本的内容,也是电力系统运行及设计中必不可少的工具。

根据系统给定的运行条件、网络接线及元件参数,通过潮流计算可以确定各母线电压的幅值及相角、各元件中流过的功率、整个系统的功率损耗等。

matlab潮流计算

matlab潮流计算

jds=IEEE(1,1);%节点数phjdh=IEEE(1,2);%平衡节点号phv=IEEE(1,3);%平衡节点电压jd=IEEE(1,4);%精度IN=size(IEEE);aa=0;ln=1;for li=3:IN(1)%读取线路参数if IEEE(li,1)>0Line(ln,:)=IEEE(li,:); ln=ln+1;endif IEEE(li,1)==0breakendendLI=li;%读取变压器参数cn=1;for ci=LI+1:IN(1)if IEEE(LI+1,1)==0trans=0;endif IEEE(ci,1)>0trans(cn,:)=IEEE(ci,:); cn=cn+1;endif IEEE(ci,1)==0breakendendif ci==LI+1CI=ci+1;else CI=ci;end%接地支路数jdn=1;for jdi=CI+1:IN(1)if IEEE(CI+1,1)==0jdzls=0;endif IEEE(jdi,1)>0jdzls(jdn,:)=IEEE(jdi,:);jdn=jdn+1;endif IEEE(jdi,1)==0breakendendif jdi==CI+1JDI=jdi+1;else JDI=jdi;end%读取节点功率数据ppn=1;for ppi=JDI+1:IN(1)if IEEE(ppi,1)>0jdgl(ppn,:)=IEEE(ppi,:);ppn=ppn+1;endif IEEE(ppi,1)==0breakendendPPI=ppi;%读取PV节点数据pvn=1;for pvi=PPI+1:IN(1)if IEEE(pvi,1)>0pvsj(pvn,:)=IEEE(pvi,:);pvn=pvn+1;endif IEEE(pvi,1)==0breakendendnline=size(Line,1);%读取线路个数ntrans=size(trans,1);%变压器个数npq=size(jdgl);Y=zeros(jds);if nline>=1for k=1:nlinet1=Line(k,2); t2=Line(k,3); b2=Line(k,6); Yl=1/(Line(k,4)+j*Line(k,5));Y(t1,t1)=Y(t1,t1)+Yl+j*b2;Y(t1,t2)=Y(t1,t2)-Yl;Y(t2,t1)=Y(t2,t1)-Yl;Y(t2,t2)=Y(t2,t2)+Yl+j*b2;endendif ntrans>=1for k=1:ntranst1=trans(k,2);t2=trans(k,3);t3=trans(k,6); Yt=1/(trans(k,4)+j*trans(k,5));Yt1=Yt/t3;Yt2=Yt*(1-t3)/(t3*t3);Yt3=Yt*(t3-1)/t3;Y(t1,t1)=Y(t1,t1)+Yt1+Yt2;Y(t1,t2)=Y(t1,t2)-Yt1;Y(t2,t1)=Y(t2,t1)-Yt1;Y(t2,t2)=Y(t2,t2)+Yt1+Yt3;endendG=real(Y);B=imag(Y);%设各个节点电压for ii=1:jdsf(ii)=0;u(ii)=1;endfor ii=1:pvsju(pvsj(ii,2))=pvsj(ii,3);endu(phjdh)=real(phv);f(phjdh)=imag(phv);npv=size(pvsj);npq=size(jdgl)-npv(1);for ii=1:size(jdgl)p0(jdgl(ii,2))=jdgl(ii,3)-jdgl(ii,5);q0(jdgl(ii,2))=jdgl(ii,4)-jdgl(ii,6);enddmax=1;%误差上限kn=0;while dmax>IEEE(1,4)for i1=1:jdsp(i1)=0;q(i1)=0;for j1=1:jdsif i1==phjdhbreakendp1(j1)=u(i1)*(G(i1,j1)*u(j1)-B(i1,j1)*f(j1))+f(i1)*(G(i1,j1)*f(j1 )+B(i1,j1)*u(j1));p(i1)=p1(j1)+ p(i1);q1(j1)=f(i1)*(G(i1,j1)*u(j1)-B(i1,j1)*f(j1))-u(i1)*(G(i1,j1)*f(j1 )+B(i1,j1)*u(j1));q(i1)=q1(j1)+ q(i1);endend%计算PQ节点有功,无功功率不平衡量for i1=1:jdsif i1==phjdhii1=i1;breakenddp(i1)=p0(i1)-p(i1);dq(i1)=q0(i1)-q(i1);U1(i1)=u(i1)-f(i1)*i;I(i1)=(p(i1)-q(i1)*i)./U1(i1);a(i1)=real(I(i1));b(i1)= imag(I(i1));endfor i1=ii1+1:jdsdp(i1-1)=p0(i1)-p(i1);dq(i1-1)=q0(i1)-q(i1);U1(i1-1)=u(i1)-f(i1)*i;I(i1-1)=(p(i1)-q(i1)*i)./U1(i1-1);a(i1-1)=real(I(i1-1));b(i1-1)= imag(I(i1-1));end%pq对角元素for i1=1:jds-1if i1<=npq(1)H(i1,i1)=-B(i1,i1)*u(i1)+G(i1,i1)*f(i1)+b(i1);N(i1,i1)=G(i1,i1)*u(i1)+B(i1,i1)*f(i1)+a(i1);J(i1,i1)=-G(i1,i1)*u(i1)-B(i1,i1)*f(i1)+a(i1);L(i1,i1)=-B(i1,i1)*u(i1)+G(i1,i1)*f(i1)-b(i1);YKB(i1*2-1,i1*2-1)= H(i1,i1);YKB(i1*2-1,i1*2)= N(i1,i1);YKB(i1*2,i1*2-1)= J(i1,i1);YKB(i1*2,i1*2)=L(i1,i1);endif i1>npq(1)H(i1,i1)=-B(i1,i1)*u(i1)+G(i1,i1)*f(i1)+b(i1); N(i1,i1)=G(i1,i1)*u(i1)+B(i1,i1)*f(i1)+a(i1); R(i1,i1)=2*f(i1);S(i1,i1)=2*u(i1);YKB(i1*2-1,i1*2-1)= H(i1,i1);YKB(i1*2-1,i1*2)= N(i1,i1);YKB(i1*2,i1*2-1)= R(i1,i1);YKB(i1*2,i1*2)=S(i1,i1);endend%pq非对角元素for i1=1:jds-1for j1=1:jds-1if i1~=j1if i1<=npqH(i1,j1)=-B(i1,j1)*u(i1)+G(i1,j1)*f(i1);N(i1,j1)=G(i1,j1)*u(i1)+B(i1,j1)*f(i1);J(i1,j1)=-B(i1,j1)*f(i1)-G(i1,j1)*u(i1);L(i1,j1)=G(i1,j1)*f(i1)-B(i1,j1)*u(i1);YKB(i1*2-1,j1*2-1)= H(i1,j1);YKB(i1*2-1,j1*2)= N(i1,j1);YKB(i1*2,j1*2-1)= J(i1,j1);YKB(i1*2,j1*2)=L(i1,j1);endif i1>npq(1)H(i1,j1)=-B(i1,j1)*u(i1)+G(i1,j1)*f(i1);N(i1,j1)=G(i1,j1)*u(i1)+B(i1,j1)*f(i1);R(i1,j1)=0;S(i1,j1)=0;YKB(i1*2-1,j1*2-1)= H(i1,j1);YKB(i1*2-1,j1*2)= N(i1,j1);YKB(i1*2,j1*2-1)= R(i1,j1);YKB(i1*2,j1*2)=S(i1,j1);endendendendYKB%求修正向量for i1=1:jds-1if i1<=npqdpq(2*i1-1,1)=dp(i1);dpq(2*i1,1)=dq(i1);elsedpq(2*i1-1)=dp(i1);dpq(2*i1,1)=0;endenddfu=YKB\dpqdmax=max(abs(dfu));nn=1;%nn避免与平衡节点重叠for i1=1:jdsif i1==phjdhii1=i1+1;breakendif i1~=phjdhdf(nn)=dfu(2*i1-1,1);du(nn)=dfu(2*i1,1);endnn=nn+1;endfor i1=ii1:jdsdf(i1-1)=dfu(2*i1-1,1); du(i1)=dfu(2*i1-1,1);endfor i1=1:jdsf(i1)=f(i1)+df(i1);u(i1)=u(i1)+du(i1);endkn=kn+1end%while的%求平衡节点功率for i1=1:jdsU(i1)=u(i1)+f(i1)*i;endS11=0;for i1=1:jdsss=conj(Y(1,i1)).*conj(U(i1));S11=S11+ss;endS1=u(1)*S11;%平衡节点功率nn=1;for i1=1:jdsif i1~=phjdhp0(nn)=p0(i1);q0(nn)=q0(i1);endif i1==phjdhp0(nn)=real(S1);q0(nn)=imag(S1);endnn=nn+1;endP(i1)=p0(i1)+q0(i1)*i;%求网络总损耗SP=0;for i1=1:jdsSP=P(i1)+SP;enddS=S1+SP;%网络总损耗%求输电效率a=0;b=0;for i1=1:jdsif p0(i1)>0a=a+p0(i1);elseb=b+p0(i1);endendn=abs(b)/(real(S1)+a)*100;%输电效率%求线路功率for i1=1:jdsfor j1=1:jdsS(i1,j1)=U(i1)*(conj(U(i1))-conj(U(j1)))*conj(-Y(i1,j1));%线路功率endendS。

基于MATLAB的潮流计算

基于MATLAB的潮流计算

目录1.摘要 (3)2.题目原始资料 (4)3.题目分析 (6)4.题目求解 (7)1)根据题意要求画出等值电路 (7)2)读程序画出拉夫逊法的流程图 (8)3)变电所负荷为题目所给数据进行求解 (8)4)编写程序并运行 (10)5)具体调压调损耗过程 (10)1.改变变压器变比调压 (10)2.改变发电机机端电压调压 (12)3.负荷按照一定比例变化的潮流计算分析 (15)4.轮流断开支路双回线中的一条的潮流计算 (19)5.仿真并比较 (28)6.设计心得 (30)7.参考文献 (31)摘要本文运用MATLAB软件进行潮流计算,对给定题目进行分析计算,再应用DDRTS软件,构建系统图进行仿真,最终得到合理的系统潮流。

潮流计算是电力系统最基本最常用的计算。

根据系统给定的运行条件,网络接线及元件参数,通过潮流计算可以确定各母线的电压幅值和相角,各元件流过的功率,整个系统的功率损耗。

潮流计算是实现电力系统安全经济发供电的必要手段和重要工作环节。

因此,潮流计算在电力系统的规划计算,生产运行,调度管理及科学计算中都有着广泛的应用。

首先,画出系统的等效电路图,在计算出各元件参数的基础上,应用牛顿—拉夫逊Newton-Raphson法以及MATLAB软件进行计算对给定系统图进行了四种不同负荷下的潮流计算,经过调节均得到符合电压限制及功率限制的潮流分布。

其次,牛顿—拉夫逊Newton-Raphson法具有较好的收敛性,上述计算过程经过四到五次迭代后均能收敛。

根据运算结果,分析各支路损耗和系统总损耗。

最后,应用DDRTS软件,构建系统图,对给定负荷重新进行分析,潮流计算后的结果也能满足相应的参数要求。

关键词:牛顿-拉夫逊法 MATLAB DDRTS 潮流计算一、 题目原始资料:1.系统图:两个发电厂分别通过变压器和输电线路与四个变电所相连。

2、发电厂资料:母线1和2为发电厂发高压母线,发电厂一总装机容量为(300MW ),母线3为机压母线上装机容量为(100MW ),最大负荷和最小负荷分别为50MW 和20MW ,发电厂二总装机容量为(200MW )3、变电所资料:(一)变电所1、2、3、4低压母线的电压等级分别为:35KV 10KV 35KV 10KV (二)变电所的负荷分别为: 50MW 50MW 40MW 70MW变电所1 变电所2母线电厂一 电厂二(三)每个变电所的功率因数均为cos φ=0.85;(四)变电所1和变电所3分别配有两台容量为75MVA 的变压器,短路损耗414KW ,短路电压(%)=16.7;变电所2和变电所4分别配有两台容量为63MVA 的变压器,短路损耗为245KW ,短路电压(%)=10.5;4、输电线路资料:发电厂和变电所之间的输电线路的电压等级及长度标于图中,单位长度的电阻为Ω17.0,单位长度的电抗为Ω0.402,单位长度的电纳为S -610*2.78。

基于MATLAB电力系统潮流计算和分析

基于MATLAB电力系统潮流计算和分析

一、实验目的了解计算机潮流分析的基本原理、主要步骤;熟悉Matlab运行环境;了解MATLAB潮流分析的步骤;对给定网络的运行方式做潮流分析,并初步分析计算结果。

二、实验原理:实验原理如下图:图1 系统原理图三、实验仪器、设备:一台装有MATLAB R2010a的个人计算机三相同步发电机模型,变压器模型,负荷模型,线路元件模型四、实验步骤:(1)熟悉原始资料:根据计算要求,整理数据,包括:计算网络中线路、变压器的参数(以上数据均采用有名值计算)(2)上机调试:熟悉Matlab的运行环境,准确输入原始数据、节点编号、节点注入功率等信息;(3)整理计算结果:根据计算结果作电网潮流分布图。

五、实验网络接线图及原始数据如图所示,3为平衡节点,1、2为P、Q节点,电压等级为110kV,节点处功率已将各线路充电功率考虑在内,3节点电压为115kV,角度为0。

原始数据各参数是以其自身额定功率和额定电压为基准的标幺值。

发电机参数 Pn=100MV·A,Un=10.5KV ,fn=50Hz, 变压器参数采用Y-Y 连接方式 T1的参数 Pn=100MV·A,fn=50Hz,一次绕组:U1=10.5KV ,R1=0.002,L1=0,二次绕组:U2=121KV ,R2=0.002,L2=0.015, Rm=5000,Lm=5000 T2的参数 Pn=100MV·A,fn=50Hz,一次绕组:U1=10.5KV ,R1=0.002,L1=0, 二次绕组:U2=121KV ,R2=0.002,L2=0.03,25+j18MV A3225+j18MVA32Rm=5000,Lm=5000线路参数L23:R*=0.08,X*=0.30,Y*=0.5L31:R*=0.10,X*=0.35,Y*=0L12:R*=0.04,X*=0.25,Y*=0.5六、实验数据记录及处理:从各个Scope中可以看到输电线π型等值电路两端的有功与无功功率的波形,具体操作方法是从Workspace中读出记录的数据(如图三、图四),数据有多组,取其平均值,分别得到各输电线π型等值电路两端的有功和无功功率。

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

电力网潮流计算的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。

matlab潮流计算

matlab潮流计算

matlab潮流计算%线路数据line=[1-节点编号2-线路首端节点号3-线路末端节点号4-支路电阻5-支路电抗6-支路电纳(注意:此处取的是b/2)]line=[1120.040.250.252130.10.3503230.080.300.25];%变压器数据transform=[1-支路编号2-支路首节点编号3-支路末节点编号4-支路电阻(p.u.)5-支路电抗(p.u.)6-变压器变比(p.u.)]transform=[12400.0151.0523500.031.05];%数据预处理nbus=5;%节点数nline=size(line,1);%线路个数ntrans=size(transform,1);%变压器个数slack=5;%均衡节点号npq=4;%pq节点的个数ss=0;%计算节点导纳矩阵y=zeros(nbus);ifnline>=1%推论与否存有线路fork=1:nline%以下处理线路t1=line(k,2);t2=line(k,3);b2=line(k,6);%分别抽出线路的首端节点编号t1、末端节点编号t2和对地电纳b2yl=1/(line(k,4)+j*line(k,5));%计算线路的支路电导yly(t1,t1)=y(t1,t1)+yl+j*b2;%修正第k条线路首端节点的自导纳y(t1,t2)=y(t1,t2)-yl;%修正第k条线路首端节点与末端节点之间的互磁滞y(t2,t1)=y(t2,t1)-yl;%修正第k条线路末端节点与首端节点之间的互磁滞y(t2,t2)=y(t2,t2)+yl+j*b2;%修正第k条线路末端节点的自磁滞endendifntrans>=1%判断是否存在变压器fork=1:ntrans%以下处理变压器t1=transform(k,2);t2=transform(k,3);t3=transform(k,6);%分别抽出变压器的首端节点编号t1、末端节点编号t2和变比t3yt=1/(transform(k,4)+j*transform(k,5));yt1=yt/t3;yt2=yt*(1-t3)/(t3*t3);yt3=yt*(t3-1)/t3;y(t1,t1)=y(t1,t1)+yt1+yt2;y(t1,t2)=y(t1,t2)-yt1;y(t2,t1)=y(t2,t1)-yt1;y(t2,t2)=y(t2,t2)+yt1+yt3;endendg=real(y);b=imag(y);%区分节点磁滞矩阵的实部和虚部gb%赋初值delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;u(1)=1;u(2)=1;u(3)=1;u(4)=1;p(1)=-0.30;q(1)=-0.18;p(2)=-0.55;q(2)=-0.13;p(3)=0;q(3)=0;p(4)=0.5;q(4)=1.10;p(5)=0.8;q(5)=0.50;k=0;precision=1;npq=4;%npq分别就是网络中的pq节点数%[unbalance]=-[jacobi][correction]whileprecision>0.00001%设定误差上限,判断是否继续迭代u(5)=1.06;delt(5)=0;%设定平衡节点电压相角与幅值k;u;delt;form=1:npqforn=1:nbuspt(n)=u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));%由节点电压求出的pq节点转化成军功功率qt(n)=u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));%由节点电压求得的pq节点注入无功功率endunbalance(2*m-1)=p(m)-sum(pt);%排序pq节点军功功率不平来衡量unbalance(2*m)=q(m)-sum(qt);%排序pq节点无功功率不平来衡量end%[unbalance]是节点不平衡量矩阵form=1:npqforn=1:nbush0(n)=u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));n0(n)=-u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));j0(n)=-u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));l0(n)=-u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));endh(m,m)=sum(h0)-u(m)^2*(g(m,m)*sin(delt(m)-delt(m))-b(m,m)*cos(delt(m)-delt(m)));n(m,m)=sum(n0)+u(m)^2*(g(m,m)*cos(delt(m)-delt(m))+b(m,m)*sin(delt(m)-delt(m)))-2*u(m)^2*g(m,m);j(m,m)=sum(j0)+u(m)^2*(g(m,m)*cos(delt(m)-delt(m))+b(m,m)*sin(delt(m)-delt(m)));l(m,m)=sum(l0)+u(m)^2*(g(m,m)*sin(delt(m)-delt(m))-b(m,m)*cos(delt(m)-delt(m)))+2*u(m)^2*b(m,m);jacobi(2*m-1,2*m-1)=h(m,m);jacobi(2*m-1,2*m)=n(m,m);jacobi(2*m,2*m-1)=j(m,m);jacobi(2*m,2*m)=l(m,m);end%计算m=n情况下的jacobi矩阵中的子矩阵元素form=1:npqforn=1:npqifm==nelseh(m,n)=-u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));j(m,n)=u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));n(m,n)=-j(m,n);l(m,n)=h(m,n);jacobi(2*m-1,2*n-1)=h(m,n);jacobi(2*m-1,2*n)=n(m,n);jacobi(2*m,2*n-1)=j(m,n);jacobi(2*m,2*n)=l(m,n);endendend%排序m≠n情况下的jacobi矩阵中的子矩阵元素correction=-jacobi\\(unbalance');%计算电压相角和幅值的修正量precision=max(abs(correction));%取误差最大值form=1:npqdelt(m)=delt(m)+correction(2*m-1);%修正pq节点电压相角u(m)=u(m)+correction(2*m);%修正pq节点电压幅值endk=k+1;%运算轮数+1endk,u,delt,jacobi,precisionform=1:nbusu(m)=u(m)*(cos(delt(m))+j*sin(delt(m)));%采用直角坐标系表示电压i(m)=y(nbus,m)*u(m);%计算注入平衡节点的电流endsslack=u(nbus)*sum(conj(i))%排序转化成均衡节点的功率form=1:nbusforn=1:nbuss(m,n)=u(m)*(conj(u(m))-conj(u(n)))*conj(-y(m,n));%排序线路功率ploss(m,n)=u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));endends%表明线路功率矩阵,s(m,n)则表示假设功率由节点m流向节点n%若数值为+则说明实际功率流向与假设方向相同,若数值为-则说明实际功率流向与假设方向相反。

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

潮流例题:根据给定的参数或工程具体要求(如图),收集和查阅资料;学习相关软件(软件自选:本设计选择Matlab进行设计)。

2.在给定的电力网络上画出等值电路图。

3.运用计算机进行潮流计算。

4.编写设计说明书。

一、设计原理
1.牛顿-拉夫逊原理
牛顿迭代法是取x0 之后,在这个基础上,找到比x0 更接近的方程的跟,一步一步迭代,从而找到更接近方程根的近似跟。

牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程f(x) = 0 的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根。

电力系统潮流计算,一般来说,各个母线所供负荷的功率是已知的,各个节点电压是未知的(平衡节点外)可以根据网络结构形成节点导纳矩阵,然后由节点导纳矩阵列写功率方程,由于功率方程里功率是已知的,电压的幅值和相角是未知的,这样潮流计算的问题就转化为求解非线性方程组的问题了。

为了便于用迭代法解方程组,需要将上述功率方程改写成功率平衡方程,并对功率平衡方程求偏导,得出对应的雅可比矩阵,给未知节点赋电压初值,一般为额定电压,将初值带入功率平衡方程,得到功率不平衡量,这样由功率不平衡量、雅可比矩阵、节点电压不
平衡量(未知的)构成了误差方程,解误差方程,得到节点电压不平衡量,节点电压加上节点电压不平衡量构成新的节点电压初值,将新的初值带入原来的功率平衡方程,并重新形成雅可比矩阵,然后计算新的电压不平衡量,这样不断迭代,不断修正,一般迭代三到五次就能收敛。

牛顿—拉夫逊迭代法的一般步骤:
(1)形成各节点导纳矩阵Y。

(2)设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。

(3)计算各个节点的功率不平衡量。

(4)根据收敛条件判断是否满足,若不满足则向下进行。

(5)计算雅可比矩阵中的各元素。

(6)修正方程式个节点电压
(7)利用新值自第(3)步开始进入下一次迭代,直至达到精度退出循环。

(8)计算平衡节点输出功率和各线路功率
2.网络节点的优化
1)静态地按最少出线支路数编号
这种方法由称为静态优化法。

在编号以前。

首先统计电力网络个节点的出线支路数,然后,按出线支路数有少到多的节点顺序编号。

当由n 个节点的出线支路相同时,则可以按任意次序对这n 个节点进行编号。

这种编号方法的根据是导纳矩阵中,出线支路数最少的节点所对应的行中非零元素也2)动态地按增加出线支路数最少编号在上述的方法中,各节点的出线支路数是按原始网络统计出来的,在编号过程中认为固定不变的,事实上,在节点消去过程中,每消去一个节点以后,与该节点相连的各节点的出线支路数将发生变化(增加,减少或保持不变)。

因此,如果每消去一个节点后,立即修正尚未编号节点的出线支路数,然后选其中支路数最少的一个节点进行编号,就可以预期得到更好的效果,动态按最少出线支路数编号方法的特点就是按出线最少原则编号时考虑了消去过程中各节点出线支路数目的变动情况。

3.MATLAB编程应用
Matlab 是“Matrix Laboratory”的缩写,主要包括:一般数值分析,矩阵运算、数字信号处理、建模、系统控制、优化和图形显示等应用程序。

由于使用Matlab 编程运算与人进行科学计算的思路和表达方式完全一致,所以不像学习高级语言那样难于掌握,而且编程效率和计算效率极高,还可在计算机上直接输出结果和精美的图形拷贝,所以它的确为一高效的科研助手。

二、设计内容
1.设计流程图
2.程序
clear;clc
%重新编号,把原题中的节点1,2,3,4,5重新依次编号为5,1,2,3,4,其中1-4号为PQ节点,5号为平衡节点
y=0;
y (1,2)=1/(0.06+0.18i); y (1,3)=1/(0.06+0.18i); y (1,4)=1/(0.04+0.12i);
y(1,5)=1/(0.02+0.06i);
y(2,3)=1/(0.01+0.03i);y(2,5)=1/(0.08+0.24i);
y(3,4)=1/(0.08+0.24i);
y(4,5)=0;
for i=1:5
for j=i:5
y(j,i)=y(i,j);
end
end
Y=0;
%求互导纳
for i=1:5
for j=1:5
if i~=j
Y(i,j)=-y(i,j);
end
end
end
%求自导纳
for i=1:5
Y(i,i)=sum(y(i,:));
end
Y %Y 为导纳矩阵
G=real(Y);
B=imag(Y);
%原始节点功率
S(1)=0.2+0.2i;
S(2)=-0.45-0.15i;
S(3)=-0.4-0.05i;
S(4)=-0.6-0.1i;
S(5)=0;
P=real(S);
Q=imag(S);
%赋初值
U=ones(1,5);U(5)=1.06;
e=zeros(1,5);
ox=ones(8,1);fx=ones(8,1);
count=0 %计算迭代次数
while max(fx)>1e-5
for i=1:4
for j=1:4
H(i,j)=0;N(i,j)=0;M(i,j)=0;L(i,j)=0;oP(i)=0;oQ(i)=0;
end
end
for j=1:5
oP(i)=oP(i)-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));
oQ(i)=oQ(i)-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)));
end
oP(i)=oP(i)+P(i); oQ(i)=oQ(i)+Q(i);
end
fx=[oP,oQ]';
%求雅克比矩阵
%当i~=j时候求H,N,M,L 如下:
for i=1:4
for j=1:4
if i~=j H(i,j)=-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)));
N(i,j)=-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));
L(i,j)=H(i,j);
M(i,j)=-N(i,j);
end
end
end
H,N,M,L
%当i=j 时H,N,M,L如下:
for i=1:4
for j=1:5
if i~=j
H(i,i)=H(i,i)+U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i, j)*cos (e(i)-e(j))); N(i,i)=N(i,i)-U(i)*U(j)*(G(i, j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));
M(i,i)=M(i,i)-U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)));
L(i,i)=L(i,i)-U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)));
end
end
N(i,i)=N(i,i)-2*(U(i))^2*G(i,i);
L(i,i)=L(i,i)+2*(U(i))^2*B(i,i);
end
J=[H,N;M,L] %J 为雅克比矩阵
ox=-((inv(J))*fx);
for i=1:4
oe(i)=ox(i); oU(i)=ox(i+4)*U(i);
end
for i=1:4
e(i)=e(i)+oe(i); U(i)=U(i)+oU(i);
end
count=count+1;
end
ox,U,e,count
%求节点注入的净功率
i=5;
P(i)=U(i)*U(j)*(G(i,j)*cos(e(i)-e(j))+B(i,j)*sin(e(i)-e(j)))+P(i);
Q(i)=U(i)*U(j)*(G(i,j)*sin(e(i)-e(j))-B(i,j)*cos(e(i)-e(j)))+Q(i);
end
S(5)=P(5)+Q(5)*sqrt(-1);
S
%求节点注入电流
I=Y*U'
3.运行结果
Y值:
迭代过程:
电压值:
平衡节点注入功率及电流:。

相关文档
最新文档