matlab电力系统快速解耦法潮流计算及短路计算程序.doc

电力系统快速解耦法潮流分析及短路计算 一程序设计的基本思想：（1） 由于电力系统潮流分析中要利用到矩阵运算，复数运算，故采用matlab编程。采用文件输入，将系统的各个参数以文件的形式输入，便于程序的通用化。（2） 本程序共有两个输入文件，分别为线路参数的文件，和已知的节点状态文件（PQ）（3） 为了使程序不仅仅局限于计算9节点网络，在形成节点导纳的函数Yn（）中，利用循环，找出线路首节点中的最大编号数，自动确定节点导纳矩阵的维数。故对于任意n节点网络，均可以计算出节点导纳矩阵（4） 在（3）的前提下，为了使程序支持系统增加节点，增加负荷等造成的PQ参数改变，或者PQ表的加长。对程序做了如下优化。首先，程序执行的基础是PQ表中平衡节点在第一行，接下来是PV节点，最后是PQ节点，如果系统添加节点，或者删除节点，均在PQ表的末端操作，会造成PQ表的顺序不是平衡节点、PV节点、PQ节点的顺序。故引入了seqencing（）函数，其作用就是不论输入的PQ表是什么顺序，在程序读入后均按平衡-PV-PQ的顺序排列。其次，顺序打乱的PQ表必须与支路参数表对应，故在Yn（）函数中加入了两段循环体，使之对应（见相应函数体注释）（5） 在满足了上述4个条件后，程序便可以通用化了。当然，由于水平有限，且程序未能由大量数据测试，故缺陷在所难免，这里仅是做了通用化的尝试。在本文最末附加了该程序通用化的实例。二、潮流计算框图输入系统的支路参数和PQ参数计算节点导纳矩阵设置初始电压值U0设置收敛条件及最大迭代次数Kmax判断是否收敛求解修正方程得到U和对U0进行修正，U0=U0+UKKmax ?NY输出结果结束输出潮流不收敛信息YN开始三定义相应的函数1.形成节点导纳矩阵的函数Yn（）function Y=Yn(x,y) %定义一名为Yn的函数，其功能是自动识别输入表中节点的个数，形成相应的节点导纳矩阵fid,message=fopen(x,r) ; %从x文件中读入支路参数if fid=-1; %判断文件是否正确打开 error(message);end;HeadPoint,HeadNumber, EndPoint,EndNumber,R,X,B,k=textread(x,%s %d %s %d %f %f %f %f); %将读入的参数处理为以列为向量的数组fclose(fid);%关闭文件L=length(HeadNumber); %确定输入表的行数 Pointstyle,PointNumber,Ps,Qs,Uk,Ok=seqencing(y); %调用seqencing函数，引入y文件中的PQ参数A=PointNumber;for i=1:L; %通过以下两循环体，实现PQ参数与支路参数的编号对应 for j=1:L; if HeadNumber(i)=j; HeadNumber(i)=A(j); break; end; end;end;for i=1:L; for j=1:L; if EndNumber(i)=j; EndNumber(i)=A(j); break; end; end;end;Y=zeros(L,L); %根据txt文件中数据表的长度建立空的节点导纳矩阵for i=1:L m=HeadNumber(i);n=EndNumber(i); if k(i)=0; %判断是否何种元件，为输电线元件 if n=0; Y(m,m)=Y(m,m)+1j*B(i)+1/(R(i)+1j*X(i); Y(n,n)=Y(n,n)+1j*B(i)+1/(R(i)+1j*X(i); Y(m,n)=Y(m,n)-1/(R(i)+1j*X(i); Y(n,m)=Y(n,m)-1/(R(i)+1j*X(i); else Y(m,m)=Y(m,m)+R(i)+1j*X(i); end; else %为变压器元件 if n=0; Y(m,m)=Y(m,m)+1/(R(i)+1j*X(i); Y(m,n)=Y(m,n)-1/(k(i)*(R(i)+1j*X(i); Y(n,n)=Y(n,n)+1/(k(i)*k(i)*(R(i)+1j*X(i); Y(n,m)=Y(n,m)-1/(k(i)*(R(i)+1j*X(i); else Y(m,m)=Y(m,m)+R(i)+1j*X(i); end; end; end;maxm=HeadNumber(1);%通过下面两个循环体，确定输入表中节点编号的最大值，及为节点导纳矩阵的维数for i=1:L; if maxm=HeadNumber(i); maxm=HeadNumber(i); end;end;maxn=EndNumber(1);for i=1:L; if maxnPointstyle(j+1); t=Pointstyle(j+1); Pointstyle(j+1)=Pointstyle(j); Pointstyle(j)=t; t=PointNumber(j+1); PointNumber(j+1)=PointNumber(j); PointNumber(j)=t; t=Ps(j+1); Ps(j+1)=Ps(j); Ps(j)=t; t=Qs(j+1); Qs(j+1)=Qs(j); Qs(j)=t; t=Uk(j+1); Uk(j+1)=Uk(j); Uk(j)=t; end; end;end; 3、形成解耦算法B矩阵的函数 formB1（）function B1=formB1(x,y)%定义名为B1的函数形成解耦算法中的B矩阵，得到的B矩阵用B1表示fid,message=fopen(x,r) ; %从x文件中读入支路参数if fid=-1; %判断文件是否正确打开 error(message);end;HeadPoint,HeadNumber,EndPoint,EndNumber,R,X,B,k=textread(x,%s %d %s %d %f %f %f %f); %将读入的参数处理为以列为向量的数组fclose(fid);L=length(HeadNumber); Pointstyle,PointNumber,Ps,Qs,Uk,Ok=seqencing(y);%调用seqencing函数，引入y文件中的PQ参数A=PointNumber;%通过以下两循环体，实现PQ参数与支路参数的编号对应for i=1:L; for j=1:L; if HeadNumber(i)=j; HeadNumber(i)=A(j); break; end; end;end;for i=1:L; for j=1:L; if EndNumber(i)=j; EndNumber(i)=A(j); break; end; end;end;B1=zeros(L,L);for i=1:L %以行为单位，通过循环，用支路参数对B1进行修改，形成B矩阵 m=HeadNumber(i);n=EndNumber(i); B1(m,m)=B1(m,m)-1/X(i); B1(n,n)=B1(n,n)-1/X(i); B1(m,n)=B1(m,n)+1/X(i); B1(n,m)=B1(n,m)+1/X(i);endmaxm=HeadNumber(1);for i=1:L; if maxm=HeadNumber(i); maxm=HeadNumber(i); end;end;maxn=EndNumber(1);for i=1:L; if maxn=EndNumber(i); maxn=EndNumber(i); end;end;B1=B1(2:max(maxm,maxn),2:max(maxm,maxn); %形成B矩阵 4、形成解耦算法B矩阵的函数 formB11（）function B11=formB11(x,y)%定义名为B11的函数形成解耦算法中B矩阵，用B11表示从x文件中读入支路参数确定Y，从y文件中读入PQ参数确定B11的维数，即除去平衡节点和pv节点，此处要求PQ参数录入时，将平衡节点和PQ节点放在前排，这一要求在Yn函数中通过seqencing函数已经满足Y=Yn(x,y);B=imag(Y);Pointstyle,PointNumber,Ps,Qs,Uk,Ok=seqencing(y);i=1;j=1;while Pointstyle(i)=1; i=i+1; j=j+1;end;B11=B(j:end,j:end); %形成B矩阵5、计算正常情况下系统节点电压的函数 powerflow（）function U0,O0=powerflow(x,y)%定义名为powerflow的函数，利用快速解耦算法来计算正常情况下系统内各个节点的电压和相角Pointstyle,PointNumber,Ps,Qs,Uk,Ok=seqencing(y);% 调用seqencing函数对PQ参数表进行排序 Y=Yn(x,y); %形成节点导纳矩阵，Yn为n维 B1=formB1(x,y); %形成解耦算法中的B矩阵，B1为n-1维 B11=formB11(x,y); %形成解耦算法中的B矩阵，B为m维 G=real(Y); %取Y的实部 B=imag(Y); %取Y的虚部 U0=Uk; O0=Ok; L=length(PointNumber); P=zeros(L,1); Q=zeros(L,1); dP=zeros(L,1); dQ=zeros(L,1); number=1; i=1;k=1; while Pointstyle(i)=1; %通过k值确定系统中PQ节点的个数 i=i+1; k=k+1;end; while number100 %定义迭代次数上限为100次for i=2:L; sum1=0; for j=1:L; sum1=sum1+U0(j)*(G(i,j)*cos(O0(i)-O0(j)+B(i,j)*sin(O0(i)-O0(j); %潮流方程，n-1维 end; dP(i)=Ps(i)-U0(i)*sum1;endfor i=k:L; sum2=0; for j=1:L; sum2=sum2+U0(j)*(G(i,j)*sin(O0(i)-O0(j)-B(i,j)*cos(O0(i)-O0(j); %潮流方程，m维 end; dQ(i)=Qs(i)-U0(i)*sum2;end dP1=dP(2:L)./U0(2:L); dQ1=dQ(k:L)./U0(k:L); a=max(norm(dP1,inf); b=max(norm(dQ1,inf); if max(a,b)0.00001 %判断是否收敛 break;disp(迭代)disp(k);disp(次后收敛); else %如不收敛， dO=-inv(B1)*dP1; %dO为n-1维 dU=-inv(B11)*dQ1; %dU为m维 zero1=zeros(k-1,1); zero2=0; DU=zero1;dU; DO=zero2;dO; U0=U0+DU; O0=O0+DO; number=number+1; end; if number=100; disp(迭代100次后不收敛，迭代结束); end;end;四对相应的系统进行潮流分析和短路计算定义完上述函数之后，可直接调用函数形成导纳矩阵，计算正常情况下的节点电压，进行短路计算计算短路电流，短路后各个节点电压以及支路潮流分布。程序的输入表共有三个，为network.txt ，network2.txt ，PQ参数.txt。其中network.txt存放不含发电机节点导纳和负荷导纳的系统支路参数，network2.txt存放含发电机节点导纳和负荷导纳的系统支路参数，PQ参数.txt存放系统的PQ参数，现列表如下第一列：首节点名称 第二列：首节点编号 第三列：尾节点名称 第四列：尾节点编号第五列：支路参数R 第六列：支路参数X 第七列：支路参数B/2 第八列：变比K第一列：节点类型，0-平衡节点，1-PV节点，2-PQ节点 第二列：节点编号第三列：节点的有功P 第四列：节点的无功Q 第五列：节点电压，PQ节点置为0第六列：节点电压的相角，初始值为0程序代码如下：%正常情况下节点导纳矩阵Y=Yn(network.txt,PQ参数.txt);disp(节点导纳矩阵)disp(sparse(Y);%正常情况下潮流计算结果兵绘制潮流计算框图 U0,O0=powerflow(network.txt,PQ参数.txt) %各个节点电压为U0，相角为O0disp(正常情况下系统各个节点电压);disp(U0);disp(正常情况下系统各个节点相角);disp(O0);O0=O0*180/pidisp(正常情况下系统各个节点相角角度值)%计算平衡节点功率和PV节点无功功率：for m=1:9 p1=p1+U0(1)*U0(m)*(G(1,m)*cos(O0(1)-O0(m)+B(1,m)*sin(O0(1)-O0(m); q1=q1+U0(1)*U0(m)*(G(1,m)*sin(O0(1)-O0(m)-B(1,m)*cos(O0(1)-O0(m);q2=q2+U0(2)*U0(m)*(G(2,m)*sin(O0(2)-O0(m)-B(2,m)*cos(O0(2)-O0(m);q3=q3+U0(3)*U0(m)*(G(3,m)*sin(O0(3)-O0(m)-B(3,m)*cos(O0(3)-O0(m);end;disp(平衡节点有功、无功，两个电压节点有功分别为)disp(p1); %有功功率和无功功率disp(q1); %平衡节点无功功率disp(q2); %PV节点2无功功率disp(q3); %PV节点3无功功率%增加发电机导纳yi和负荷导纳yLDi后的发电机节点和负荷节点的自导纳 Y=Yn(network2.txt,PQ参数.txt) %network2文件中存放的是经修改的支路参数disp(加发电机导纳yi和负荷导纳yLDi后的发电机节点和负荷节点的自导纳);diag(Y)%精确法计算短路电流，短路后各节点电压以及网络中各支路的电流分布Y=Yn(network2.txt,PQ参数.txt); %形成包括发电机内阻抗和符合阻抗节点导纳矩阵 Z=inv(Y); %形成节点阻抗矩阵 %题中给出的是4节点短路Zf=Z(:,4);U0,O0=powerflow(network.txt,PQ参数.txt);U01=U0.*cos(O0)+1j*U0.*sin(O0);If=U01(4)/Z(4,4); %计算短路电流IfU=U01-Zf*If; %计算短路后的各个节点电压fid,message=fopen(network2.txt,r) ; %从x文件中读入支路参数if fid=-1; %判断文件是否正确打开 error(message);end;HeadPoint,HeadNumber,EndPoint,EndNumber,R,X,B,k=textread(network2.txt,%s %d %s %d %f %f %f %f); row=HeadNumber,EndNumber;for i=1:9 %计算系统中各个支路的潮流分布 I(row(i,1),row(i,2)=(U(row(i,1)-U(row(i,2)/Z(row(i,1),row(i,2);end;disp(阻抗矩阵中的第f列);disp(精确算法短路电流 If);disp(If模值);disp(abs(If);disp(相角为);disp(angle(If);disp(短路后各个节点的电压);disp(U);disp(模值为);disp(短路后各个支路的电流值 I) ;disp(模值为 I);disp(相角为（度）); %近似计算短路电流，短路后各节点电压以及网络中各支路的电流分布If2=1/Z(4,4); %计算短路电流IfU2=1-Zf*If; %计算短路后的各个节点电压disp(近似算法短路电流 If2);disp(If2模值);disp(abs(If);disp(相角为（度）);disp(angle(If2);disp(短路后各个节点的电压);disp(U2);disp(模值为);程序运行结果如下：节点导纳矩阵Y = (1,1) 0 -17.3611i (4,1) 0 +17.3611i (2,2) 0 -16.0000i (7,2) 0 +16.0000i (3,3) 0 -17.0648i (9,3) 0 +17.0648i (1,4) 0 +17.3611i (4,4) 3.3074 -39.3089i (5,4) -1.3652 +11.6041i (6,4) -1.9422 +10.5107i (4,5) -1.3652 +11.6041i (5,5) 2.5528 -17.3382i (7,5) -1.1876 + 5.9751i (4,6) -1.9422 +10.5107i (6,6) 3.2242 -15.8409i (9,6) -1.2820 + 5.5882i (2,7) 0 +16.0000i (5,7) -1.1876 + 5.9751i (7,7) 2.8047 -35.4456i (8,7) -1.6171 +13.6980i (7,8) -1.6171 +13.6980i (8,8) 2.7722 -23.3032i (9,8) -1.1551 + 9.7843i (3,9) 0 +17.0648i (6,9) -1.2820 + 5.5882i (8,9) -1.1551 + 9.7843i (9,9) 2.4371 -32.1539i迭代 k=9次后收敛正常情况下系统各个节点电压U0 = 1.0400 1.0250 1.0250 1.0258 0.9956 1.0127 1.0258 1.0159 1.0324正常情况下系统各个节点相角O0 = 0 0.1620 0.0814 -0.0387 -0.0696 -0.0644 0.0649 0.0127 0.0343正常情况下系统各个节点相角角度值O0 = 0 9.2800 4.6647 -2.2168 -3.9888 -3.6875 3.7197 0.7275 1.9667平衡节点有功、无功，两个电压节点有功分别为 p1=0.7164 q1= 0.2705q2= 0.0665 q3= -0.1086加发电机导纳yi和负荷导纳yLDi后的发电机节点和负荷节点的自导纳ans = 0 -20.6944i 0 -19.3333i 0 -20.3982i 3.3074 -39.3089i 1.2917 -16.8338i 2.3466 -15.5484i 2.8047 -35.4456i 1.8033 -22.9641i2.4371 -32.1539i 阻抗矩阵中的第f列Zf = -0.0630 + 0.1475i -0.0645 + 0.0860i -0.0634 + 0.0864i -0.0751 + 0.1758i -0.0892 + 0.1534i -0.0854 + 0.1532i -0.0779 + 0.1040i -0.0832 + 0.1032i -0.0758 + 0.1032i精确算法短路电流 If = -2.2963 - 4.8486iIf模值为 5.3649 相角为 -115.3422短路后各个节点的电压 U = 0.1801 + 0.0333i 0.4463 + 0.0502i 0.4573 - 0.0257i -0.0000 - 0.0000i 0.0446 - 0.1498i 0.0716 - 0.1271i 0.3406 - 0.0725i 0.3243 - 0.1533i 0.3572 - 0.0950i模值为0.1831 0.4491 0.4580 0.0000 0.1563 0.1459 0.3482 0.3587 0.3696短路后各个支路的电流值 I(1,4) -0.2501 - 1.1140i (4,5) 0.8560 - 0.2072i (4,6) 0.8320 + 0.0041i (2,7) 0.5086 - 0.9694i (5,7) 0.6818 + 1.8688i (7,8) 0.3965 - 0.3175i (3,9) 0.2141 - 0.7887i (6,9) 0.8241 + 1.7182i (8,9) -0.2355 + 0.3676i模值为 I(1,4) 1.1417 (4,5) 0.8807 (4,6) 0.8320 (2,7) 1.0947 (5,7) 1.9893 (7,8) 0.5079 (3,9) 0.8172 (6,9) 1.9056 (8,9) 0.4365相角为（度）(1,4) -102.6530 (4,5) -13.6077 (4,6) 0.2855 (2,7) -62.3148 (5,7) 69.9555 (7,8) -38.6844 (3,9) -74.8141 (6,9) 64.3752 (8,9) 122.6426近似算法短路电流 If2 = 0.2105 - 5.2357iIf2模值 5.2399 相角为（度） -87.6981短路后各个节点电压值U2 = 0.1611 + 0.0000i 0.4785 - 0.0351i 0.4773 - 0.0384i 0.0000 0.1155 - 0.0063i 0.1178 - 0.0098i 0.3698 - 0.0424i 0.3673 - 0.0436i 0.3752 - 0.0459i模值为 0.1611 0.4797 0.4788 0.0000 0.1157 0.1182 0.3722 0.3698 0.3780精确算法与近似算法的比较分析（1） 短路电流误差（If精-If近）/If精=2.33%精-近=27.6441（度）两电流相差进2.33%，可见在精度要求不是很高的前提下，近似算法是完全可取的。（2） 电压模值误差比较（U精-U近）/U精= 0.1205 -0.0681 -0.0455 0.0000 0.2597 0.1900 -0.0690 -0.0311 -0.0227可见电压模值的误差均在25%以内。至此，上述程序已完成题目的所有要求。为了实现程序的通用行，在编写程序的时候对程序进行了优化，使得上述程序不但能够满足本题要求，而且对于任意多节点的系统，只要将支路参数，PQ参数列表输入txt文件中，有上述程序均可完成正常情况下潮流计算和短路计算。而且，上述程序还支持支路参数的修改，添加或删除支路、节点。四附加： 程序通用化实例若在原题的基础上，在8母线和9母线之间增加一条母线10,，相应的支路参数如下表支路R（p.u.）X（p.u.）B/2(TK) （p.u.）1400.05761.02700.06251.03900.05861.0450.010.0850.088460.0170.0920.079570.0320.1610.153690.0390.170.179780.00850.0720.07458100.00590.05080.1045910 0.0050 0.0050 0.104510新增10节点为PQ节点 P=0，Q=0；由上述条件可以写出添加节点后的支路参数文件为PQ参数文件为以这两个表作为输入，执行上述程序，计算潮流，得到各个节点的电压和相角为U0 = 1.0400 1.0250 1.0250 1.0278 0.9988 1.0168 1.0304 1.0265 1.0397 1.0385O0 = 0 0.1601 0.0795 -0.0386 -0.0695 -0.0643 0.0635 0.0113 0.03270.0215即可得到十个节点的电压和相角，完成了程序的通用化。五结语学习是一个理论与实践相结合的过程，学习电力系统分析，对相应的电力系统进行简单的潮流计算和短路分析是我们必须掌握的基本知识。在用计算机编程是手段实现系统分析的功能时，我们的问题是只着眼于既定的题目和系统。所以，怎样将程序通用化是一个值得考虑的问题，由于电力系统的复杂性以及个人知识水平的限制，谨在此方面做了一些改进，真正的知识体系还亟待在今后的学习中得以建立和提高。20