潮流计算(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 潮流计算仿真方法的
示例:
1. 定义系统参数:首先,你需要定义电力系统的参数,如发电机、负荷、变压器等。
这些参数通常包括额定电压、额定功率、电抗、电阻等。
2. 建立系统模型:使用这些参数,你可以在 MATLAB 中建立电力系统的模型。
这通常涉及到定义节点和支路,以及为它们分配相应的参数。
3. 编写潮流计算函数:接下来,你需要编写一个用于执行潮流计算的函数。
这个函数应该能够接收系统的模型和参数,并返回计算出的潮流结果,如电压、电流、功率等。
4. 运行仿真:最后,你可以运行仿真并调用你编写的潮流计算函数。
这将返回计算出的潮流结果,你可以使用这些结果进行进一步的分析或可视化。
这只是一个简单的示例,实际上在编写 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实现潮流计算程序

程序代码如下: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 极坐标求解潮流计算

电力系统分析大作业——潮流计算班级: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电力系统潮流计算的基本步骤和代码示例: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 粗略程序

%========================================================================== %========================================================================== %========================================================================== %潮流计算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);。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 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值:
迭代过程:
电压值:
平衡节点注入功率及电流:。