潮流计算上机报告 华电
潮流上机课程设计-华电
课程设计报告( 2011—2012年度第一学期)名称：电力系统潮流上机院系：电气与电子工程学院班级：学号：学生：指导教师：设计周数：两周成绩：日期： 2011年12月19日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细容见附录）1．手算2．计算机计算3．思考题三、课程设计总结或结论潮流计算是研究电力系统稳态运行的一种基本计算，最初求解电力系统潮流时大多使用手算，但随着电力系统结构的日趋复杂，计算量也越来越大。
复杂电力系统潮流计算中，由于节点数量巨大，所形成的修正方程已经无法通过手算方式解决，尤其是需要迭代次数较多时，手算所需要的时间太长，计算机潮流计算无疑为解决这一问题提供了极大的便利。
计算机潮流计算可以迅速解决复杂网络的潮流计算问题，这是由于无论系统的复杂程度如何，其节点与支路的类型是固定的，所以只需要输入节点与支路的数据，就可以解决任何一个复杂网络的潮流计算问题。
即只需要一次编程，就可以基本上解决所有复杂网络的计算。
需要注意的是，在使用牛顿—拉弗逊发计算潮流时，对于初值要选择比较接近它们的精确解，否则迭代过程可能不收敛。
潮流计算Ｃ语言程序编程过程中需要注意的是，Ｃ语言无法实现复数运算，需要将得到的值的实部与虚部分开储存并计算。
这个过程复杂并且容易出错，编写程序是需要注意。
另外需要注意的一点是：Ｃ语言数组的编号是从零开始的，在程序编写过程中应注意下标的对应。
通过这一次的电力系统潮流计算编程，我不仅对Ｃ语言的编程有了更深刻的理解，也对《电力系统分析》这门课程进行了查漏补缺和巩固，对电力系统的运行也有了更加深入的了解，受益匪浅。
四、参考文献1.《电力系统计算：电子数字计算机的应用》，交通大学等合编。
：水利电力；2.《现代电力系统分析》，王锡凡主编，科学；3.《电力系统稳态分析》，珩，中国电力，2007年，第三版；附录（包括：1.手算过程及结果；2.计算机计算流程图、表格、数据；3.思考题答案）附录一、手算过程及其结果附录二、计算机计算流程图、表格、数据编写潮流计算程序简单系统如下图所示，支路数据如下： 支路14，27，39为变压器支路，参数为100.1,058.0114==K X ，050.1,063.0227==K X 100.1,059.0339==K X其余支路为线路支路，参数为075.02/,072.0019.07878=+=B j Z , 105.02/,101.0012.08989=+=B j Z 153.02/,161.0032.05757=+=B j Z 179.02/,170.0039.06969=+=B j Z 088.02/,085.0010.04545=+=B j Z 079.02/,092.0017.04646=+=B j Z节点数据如下：o U 004.11∠=•025.1,63.122==U P ，025.1,85.033==U P5.025.15j S --= ,3.09.06j S --= ,35.00.18j S --=节点数据线路数据节点导纳矩阵Y B迭代次数：Count_Num=3各节点电压：平衡节点功率：节点编号1: 0.72887-j0.10579 PV节点无功功率：节点编号2: j0.30332PV节点无功功率：节点编号3: j0.17199线路功率及损耗：线路总损耗：ΔS= 0.05887-j0.78048附录三、思考题1、潮流计算的方法有哪些？各有何特点？答：潮流计算分为简单电力网络的手算和复杂电力网络的机算两大类，其中机算又有高斯-赛德尔迭代法，牛顿-拉夫逊迭代法和P-Q分解法。
潮流计算实验
电力系统分析实验报告实验一：潮流计算的计算机算法>> clear;n=10;nl=10;isb=1;pr=0.00001;B1=[120.03512+0.08306i0.13455i10;230.0068+0.18375i0 1.023811;140.05620+0.13289i0.05382i10;450.00811+0.24549i0 1.023811;160.05620+0.13289i0.05382i10;460.04215+0.09967i0.04037i10;670.0068+0.18375i0 1.023811;680.02810+0.06645i0.10764i10;8100.00811+0.24549i011;890.03512+0.08306i0.13455i10] B2=[00 1.1 1.101;001002;00.343+0.21256i1002;001002;00.204+0.12638i1002;001002;00.306+0.18962i1002;001002;0.50 1.1 1.103;00.343+0.21256i1002]Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);elsep=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;enddisp('导纳矩阵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);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej) endP1=C(i)*e(i)+f(i)*D(i);Q1=C(i)*f(i)-e(i)*D(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;q=q+1;J(p,q)=X4;J(m,N)=DP;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;endendendendendfor 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~=3k4=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==N0break;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)./2;e(L)=e(L)-J(k,N);k1=k+1;f(L)=f(L)-J(k1,N);endfor k=3:N0DET=abs(J(k,N));if DET>=prIT2=IT2+1;endendICT2(a)=IT2;ICT1=ICT1+1;enddisp('迭代次数：');disp(ICT1);disp('没有达到精度要求的个数：');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;enddisp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率S为(节点号从小到大排列)：');disp(S);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);elseSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));Siz(i)=Si(p,q);enddisp(Si(p,q));SSi(p,q)=Si(p,q);ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);elseSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));Sjy(i)=Sj(q,p);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p);ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i);ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);disp('-----------------------------------------------------');endfigure(1);subplot(1,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(1,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;figure(2);subplot(2,2,1);P=real(S);Q=imag(S);bar(P);xlabel('节点号');ylabel('节点注入有功');grid on;subplot(2,2,2);bar(Q);xlabel('节点号');ylabel('节点注入无功');grid on;subplot(2,2,3);P1=real(Siz);Q1=imag(Siz);bar(P1);xlabel('支路号');ylabel('支路首端注入有功');grid on;subplot(2,2,4);bar(Q1);xlabel('支路号');ylabel('支路首端注入无功');grid on;B1 =1.00002.0000 0.0351 + 0.0831i 0 + 0.1346i 1.0000 02.00003.0000 0.0068 + 0.1838i 0 1.0238 1.00001.0000 4.0000 0.0562 + 0.1329i 0 + 0.0538i 1.0000 04.00005.0000 0.0081 + 0.2455i 0 1.0238 1.00001.0000 6.0000 0.0562 + 0.1329i 0 + 0.0538i 1.0000 04.0000 6.0000 0.0422 + 0.0997i 0 + 0.0404i 1.0000 06.00007.0000 0.0068 + 0.1838i 0 1.0238 1.00006.0000 8.0000 0.0281 + 0.0664i 0 + 0.1076i 1.0000 08.0000 10.0000 0.0081 + 0.2455i 0 1.0000 1.00008.0000 9.0000 0.0351 + 0.0831i 0 + 0.1346i 1.0000 0 B2 =0 0 1.1000 1.1000 0 1.00000 0 1.0000 0 0 2.00000 0.3430 + 0.2126i 1.0000 0 0 2.00000 0 1.0000 0 0 2.00000 0.2040 + 0.1264i 1.0000 0 0 2.00000 0 1.0000 0 0 2.00000 0.3060 + 0.1896i 1.0000 0 0 2.00000 0 1.0000 0 0 2.0000 0.5000 0 1.1000 1.1000 0 3.00000 0.3430 + 0.2126i 1.0000 0 0 2.0000导纳矩阵Y=Columns 1 through 69.7177 -22.8591i -4.3185 +10.2135i 0 -2.6996 + 6.3834i 0 -2.6996 + 6.3834i-4.3185 +10.2135i 4.5104 -15.3311i -0.1964 + 5.3083i 0 0 00 -0.1964 + 5.3083i 0.2011 - 5.4347i 0 0 0-2.6996 + 6.3834i 0 0 6.4271 -18.7292i -0.1313 + 3.9744i -3.5993 + 8.5110i0 0 0 -0.1313 + 3.9744i 0.1344 - 4.0690i 0-2.6996 + 6.3834i 0 0 -3.5993 + 8.5110i 0 11.8891 -32.7444i0 0 0 0 0 -0.1964 + 5.3083i0 0 0 0 0 -5.3984 +12.7660i0 0 0 0 0 00 0 0 0 0 0Columns 7 through 100 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0-0.1964 + 5.3083i -5.3984 +12.7660i 0 00.2011 - 5.4347i 0 0 00 9.8514 -26.9275i -4.3185 +10.2135i -0.1344 + 4.0690i0 -4.3185 +10.2135i 4.3185 -10.1462i 00 -0.1344 + 4.0690i 0 0.1344 - 4.0690i迭代次数：4没有达到精度要求的个数：17 18 17 0各节点的实际电压标幺值E为(节点号从小到大排列)：Columns 1 through 61.1000 1.0757 - 0.0207i 1.0050 - 0.0780i 1.0772 - 0.0175i 1.0171 - 0.0631i 1.0762 - 0.0152iColumns 7 through 101.0112 - 0.0666i 1.0778 - 0.0051i 1.0996 + 0.0304i 1.0177 - 0.0814i-----------------------------------------------------各节点的电压大小V为(节点号从小到大排列)：1.1000 1.0759 1.0080 1.0773 1.0191 1.0763 1.0134 1.0778 1.1000 1.0209-----------------------------------------------------各节点的电压相角sida为(节点号从小到大排列)：0 -1.1046 -4.4373 -0.9283 -3.5503 -0.8106 -3.7665 -0.27181.5822 -4.5707各节点的功率S为(节点号从小到大排列)：Columns 1 through 60.7165 + 0.2587i 0.0000 + 0.0000i -0.3430 - 0.2126i -0.0000 + 0.0000i -0.2040 - 0.1264i -0.0000 + 0.0000iColumns 7 through 10-0.3060 - 0.1896i -0.0000 + 0.0000i 0.5000 + 0.0089i -0.3430 - 0.2126i-----------------------------------------------------各条支路的首端功率Si为(顺序同输入B1时一致)：0.3485 + 0.0932iS(1,2)=0.3485+0.093157i-----------------------------------------------------0.3441 + 0.2420iS(2,3)=0.34409+0.24201i-----------------------------------------------------0.1904 + 0.0760iS(1,4)=0.19038+0.07599i-----------------------------------------------------0.2044 + 0.1400iS(4,5)=0.20445+0.13999i-----------------------------------------------------0.1777 + 0.0895iS(1,6)=0.17767+0.089525i------------------------------------------------------0.0163 - 0.0055iS(4,6)=-0.016305-0.0054856i-----------------------------------------------------0.3069 + 0.2128iS(6,7)=0.30686+0.21281i------------------------------------------------------0.1477 - 0.0234iS(6,8)=-0.14767-0.02338i-----------------------------------------------------0.3443 + 0.2509iS(8,10)=0.34427+0.25091i------------------------------------------------------0.4925 - 0.1508iS(8,9)=-0.49251-0.15077i-----------------------------------------------------各条支路的末端功率Sj为(顺序同输入B1时一致)：-0.3441 - 0.2420iS(2,1)=-0.34409-0.24201i------------------------------------------------------0.3430 - 0.2126iS(3,2)=-0.343-0.21256i------------------------------------------------------0.1881 - 0.1345iS(4,1)=-0.18815-0.13451i------------------------------------------------------0.2040 - 0.1264iS(5,4)=-0.204-0.12638i------------------------------------------------------0.1755 - 0.1482iS(6,1)=-0.17551-0.14815i-----------------------------------------------------0.0163 - 0.0413iS(6,4)=0.016326-0.041272i------------------------------------------------------0.3060 - 0.1896iS(7,6)=-0.306-0.18962i-----------------------------------------------------0.1482 - 0.1001iS(8,6)=0.14824-0.10014i------------------------------------------------------0.3430 - 0.2126iS(10,8)=-0.343-0.21256i-----------------------------------------------------0.5000 + 0.0089iS(9,8)=0.5+0.0089402i-----------------------------------------------------各条支路的功率损耗DS为(顺序同输入B1时一致)：0.0044 - 0.1488iDS(1,2)=0.0044095-0.14885i-----------------------------------------------------0.0011 + 0.0294iDS(2,3)=0.0010897+0.029445i-----------------------------------------------------0.0022 - 0.0585iDS(1,4)=0.0022306-0.058518i-----------------------------------------------------0.0004 + 0.0136iDS(4,5)=0.00044972+0.013613i-----------------------------------------------------0.0022 - 0.0586iDS(1,6)=0.0021584-0.058629i-----------------------------------------------------0.0000 - 0.0468iDS(4,6)=2.1344e-005-0.046758i-----------------------------------------------------0.0009 + 0.0232iDS(6,7)=0.00085804+0.023186i-----------------------------------------------------0.0006 - 0.1235iDS(6,8)=0.00056584-0.12352i-----------------------------------------------------0.0013 + 0.0384iDS(8,10)=0.001267+0.038353i-----------------------------------------------------0.0075 - 0.1418iDS(8,9)=0.0074931-0.14183i----------------------------------------------------- >>。
电力系统潮流计算实验报告
电力系统潮流上机计算实验报告１1．手算过程已知：节点1：PQ 节点，节点， s(1)= s(1)= -0.5000-j0.3500 节点2：PV 节点，节点， p(2)=0.4000 v(2)=1.0500 p(2)=0.4000 v(2)=1.0500 节点3：平衡节点，：平衡节点，U(3)=1.0000U(3)=1.0000U(3)=1.0000∠∠0.0000 网络的连接图：0.0500+j0.2000 1 0.0500+j0.2000231）计算节点导纳矩阵由2000.00500.012j Z +=Þ71.418.112j y -=; 2000.00500.013j Z +=Þ71.418.113j y -=;\导纳矩阵中的各元素：42.936.271.418.171.418.1131211j j j y y Y -=-+-=+=;71.418.11212jy Y +-=-=; 71.418.11313j y Y +-=-=; =21Y71.418.11212j y Y +-=-=; 71.418.12122j y Y -==; 002323j y Y +=-=;=31Y 71.418.11313j y Y +-=-=; =32Y 002323j y Y +=-=; 71.418.13133j y Y -==;\形成导纳矩阵B Y ：úúúûùêêêëé-++-+-+-+-+--=71.418.10071.418.10071.418.171.418.171.418.171.418.142.936.2j j j j j j j j j Y B 2）计算各PQ PQ、、PV 节点功率的不平衡量，及PV 节点电压的不平衡量：取:000.0000.1)0(1)0(1)0(1j jf e U +=+=000.0000.1)0(2)0(2)0(2j jf e U +=+=节点3是平衡节点，保持000.0000.1333j jf e U +=+=为定值。
华电潮流上机计算实验报告全解
院系：电气与电子工程学院班级：电气1205学号：1121181223学生姓名：王城指导教师：孙英云设计周数：两周成绩：日期：2015年7月7日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1．手算2．计算机计算3．思考题三、课程设计总结或结论（详细内容见附录）四、参考文献1.《电力系统计算：电子数字计算机的应用》，西安交通大学等合编。
北京：水利电力出版社；2.《现代电力系统分析》，王锡凡主编，科学出版社；3.《电力系统稳态分析》，陈珩，中国电力出版社，1995年，第三版；附录（设计流程图、程序、表格、数据等）4.机算潮流程序及结果// dierti.cpp : 定义控制台应用程序的入口点。
//#include "stdafx.h"struct Line //线路结构体{int Num,NumI,NumJ; //线路号左节点名右节点名float R,X,B,K; //电阻电抗电纳变比（K等于1为普通支路，不等于1为变压器支路的变比） };struct Bus //节点结构体{int Num ;float Volt,Phase,GenP,GenQ,LoadP,LoadQ;int Type;};#include"stdio.h"#include"string.h"#include"math.h"#include"stdlib.h"#define NBUS 4#define NLINE 4/* Global variables */int nL,nB,nVA,nSH;float X[NBUS];int L;double def[2*NBUS];double mn[50];void Gauss(double a[50][50],double b[50], int n) /*定义高斯法 */{int JS[50];int i,j,k;float d,t,x[50];FILE *fp;int L=1;for(i=0;i<50;i++) JS[i]=0;for(k=0;k<n;k++){d=0.0;for(j=k;j<n;j++)if(fabs(a[k][j])>d){ /*在一行中找到一个最大值赋值d,并用JS[K]记住这个最大值所在的列号*/ d=fabs(a[k][j]);JS[k]=j;}if(fabs(d)<0.000001) /*如果d的数值太小，做为被除数将带来很大的误差 */L=0;else {if(JS[k]!=k)for(i=0;i<n;i++){t=a[i][k];a[i][k]=a[i][JS[k]]; /*进行列交换，让最大值始终在对角元上*/a[i][JS[k]]=t;}}if(L==0)break;for(j=k+1;j<n;j++)a[k][j]=a[k][j]/a[k][k]; /*对角元上的元素消为1*/b[k]=b[k]/a[k][k];for(i=k+1;i<n;i++){for(j=k+1;j<n;j++)a[i][j]=a[i][j]-a[i][k]*a[k][j]; /*使下三角阵的元素为0*/b[i]=b[i]-a[i][k]*b[k];}}if(fabs(a[n-1][n-1])>0.00001){ /*用追赶法，解方程组，求未知数x*/ x[n-1]=b[n-1];for(i=n-2;i>=0;i--){t=0.0;for(j=i+1;j<n;j++)t=t+a[i][j]*x[j];x[i]=(b[i]-t);}}if((fp=fopen("gauss.txt","w"))==NULL) /*将结果写到TXT文件中*/{printf("err");exit(0);}for(i=0;i<n;i++){fprintf(fp,"%lf",x[i]);mn[i]=x[i];fprintf(fp,"\n");}fclose(fp);if(fp!=NULL) fclose(fp);}int _tmain(int argc, _TCHAR* argv[]){FILE *fp;FILE *fpout;int i,j,k,l,h,n,v;int i1,i2,i3,kp,kq;float d1,d2,d3,d4,d5,d6,r,x,g,b,tt,LL,e,ps,qs,shsh,m;struct Line sL[NLINE];struct Bus sB[NBUS];float YG[NBUS+1][NBUS+1],YB[NBUS+1][NBUS+1];double u[50][2];i1=i2=i3=0;d1=d2=d3=d4=d5=d6=ps=qs=0.0;for(i=0;i<NBUS;i++)if((fp=fopen("in.txt","r"))==NULL){ printf("Can not open the file named 'in.txt' \n");exit(0);}fscanf(fp,"%d,%d,%d",&nB,&nL,&nSH);for(i=0;i<nB;i++){sB[i].Num=sB[i].Type=0;sB[i].Volt=1.0;sB[i].Phase=sB[i].GenP=sB[i].GenQ=sB[i].LoadP=sB[i].LoadQ=0.0;fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&i1,&d1,&d2,&d3,&d4,&d5,&d6,&i2);sB[i].Num=i1;sB[i].Volt=d1;sB[i].Phase=d2;sB[i].GenP=d3;sB[i].GenQ=d4;sB[i].LoadP=d5,sB[i].LoadQ=d6;sB[i].T ype=i2;};for(i=0;i<nL;i++){sL[i].Num=sL[i].NumI=sL[i].NumJ=0;sL[i].R=sL[i].X=sL[i].B=0.0;sL[i].K=1.0;fscanf(fp,"%2d %3d %3d %f %f %f %f",&i1,&i2,&i3,&d1,&d2,&d3,&d4);sL[i].Num=i1;sL[i].NumI=i2;sL[i].NumJ=i3;sL[i].R=d1;sL[i].X=d2;sL[i].B=d3;sL[i].K=d4;}if(fp!=NULL) fclose(fp);/*Make Y Matrix*/for(i=1;i<nB+1;i++)for(j=1;j<nB+1;j++){YG[i][j]=0.0;YB[i][j]=0.0;};for(l=0; l<nL; l++){i=sL[l].NumI;j=sL[l].NumJ;r=sL[l].R;x=sL[l].X;d1=r*r+x*x;g=r/d1;b=-x/d1;m=sL[l].K;if(fabs(sL[l].K-1.0)<0.000001) //普通支路 {YG[i][i]=YG[i][i]+g;YG[j][j]=YG[j][j]+g;YB[i][i]=YB[i][i]+b+sL[l].B;YB[j][j]=YB[j][j]+b+sL[l].B;YG[i][j]=YG[i][j]-g;YG[j][i]=YG[j][i]-g;YB[i][j]=YB[i][j]-b;YB[j][i]=YB[j][i]-b;}else //变压器支路{YG[i][i]=YG[i][i]+g/m+g*(m-1)/m;YG[j][j]=YG[j][j]+g/m+g*(1-m)/m/m;YB[i][i]=YB[i][i]+b/m+b*(m-1)/m;YB[j][j]=YB[j][j]+b/m+b*(1-m)/m/m;YG[i][j]=YG[i][j]-g/m;YG[j][i]=YG[j][i]-g/m;YB[i][j]=YB[i][j]-b/m;YB[j][i]=YB[j][i]-b/m; }}/* Check the Y matrix */if((fp=fopen("GGBB.txt","w"))==NULL){printf("Can not open the file named 'GGBB.txt' \n");exit(0);}fprintf(fp,"---Y Matrix---\n");for(i=1;i<nB+1;i++)for(j=1;j<nB+1;j++)if(fabs(YB[i][j]-0.0)>0.000001) fprintf(fp,"Y(%3d,%-3d)=(%10.5f,%10.5f)\n",i,j,YG[i][j],YB[i][j]);if(fp!=NULL) fclose(fp);/* 节点电压附初值 */for(i=1;i<nB+1;i++){if(sB[i-1].Type==0){u[i][0]=0.0;u[i][1]=1.0;}else if(sB[i-1].Type==1){u[i][1]=sB[i-1].Volt;u[i][0]=0.0;}else if(sB[i-1].Type==2){u[i][1]=sB[i-1].Volt;u[i][0]= sB[i-1].Phase;}}for(v=1;;v++)/* 迭代次数可以无限大 */{/* 节点电压附初值 */printf("迭代第%d次赋予的电压初值为e+jf:\n",v); for(i=1;i<nB+1;i++)printf("%lf,%lf\n",u[i][1],u[i][0]);printf("\n");printf("\n");/* 求偏移量 */double P_P[10];double P_Q[10];double P_UU[10];for(i=1;i<nB+1;i++){if(sB[i-1].Type==2){P_P[i]=0.0;P_Q[i]=0.0;P_UU[i]=1.05;}if(sB[i-1].Type==0){double tempP=0.0;double tempQ=0.0;for(j=1;j<nB+1;j++){tempP+=YG[i][j]*u[j][1]-YB[i][j]*u[j][0];tempQ+=YG[i][j]*u[j][0]+YB[i][j]*u[j][1];}P_P[i]=(sB[i-1].GenP-sB[i-1].LoadP)-tempP*u[i][1]-tempQ*u[i][0]; P_Q[i]=(sB[i-1].GenQ-sB[i-1].LoadQ)-tempP*u[i][0]+tempQ*u[i][1]; P_UU[i]=0.0;}if(sB[i-1].Type==1){double tempP=0.0;double tempQ=0.0;for(j=1;j<nB+1;j++){tempP+=YG[i][j]*u[j][1]-YB[i][j]*u[j][0];tempQ+=YG[i][j]*u[j][0]+YB[i][j]*u[j][1];P_P[i]=(sB[i-1].GenP-sB[i-1].LoadP)-tempP*u[i][1]-tempQ*u[i][0]; }P_UU[i]=sB[i-1].Volt*sB[i-1].Volt-u[i][1]*u[i][1]-u[i][0]*u[i][0]; P_Q[i]=0.0;}}/* 偏移量阵 */double P_PQ[6];int a=0;for(i=1;i<3;i++){P_PQ[a]=P_P[i];a=a+2;}a=1;for(i=1;i<3;i++){P_PQ[a]=P_Q[i];a=a+2;P_PQ[4]=P_P[3];P_PQ[5]=P_UU[3];printf("迭代第%d次的偏移量为：\n",v);for(i=0;i<6;i++){printf("%f",P_PQ[i]);printf("\n");}printf("\n");printf("\n");/* 雅可比矩阵 */double H[6][6],N[6][6],J[6][6],L[6][6],R[6][6],S[6][6],aa[6],bb[6]; for(i=1;i<5;i++){ if(fabs(sB[i-1].Type-2.0)<0.000001)continue;else{for(j=1;j<5;j++)if(i!=j){H[i][j]=-YB[i][j]*u[i][1]+YG[i][j]*u[i][0];N[i][j]=YG[i][j]*u[i][1]+YB[i][j]*u[i][0];J[i][j]=-N[i][j];L[i][j]=H[i][j];R[i][j]=0;S[i][j]=0;}else{aa[i]=bb[i]=0.0;aa[i]+=YG[i][n]*u[n][1]-YB[i][n]*u[n][0];bb[i]+=YG[i][n]*u[n][0]+YB[i][n]*u[n][1];}H[i][i]=-YB[i][i]*u[i][1]+YG[i][i]*u[i][0]+bb[i]; N[i][i]=YG[i][i]*u[i][1]+YB[i][i]*u[i][0]+aa[i]; J[i][i]=-YG[i][i]*u[i][1]-YB[i][i]*u[i][0]+aa[i]; L[i][i]=YG[i][i]*u[i][0]-YB[i][i]*u[i][1]-bb[i]; R[i][i]=2*u[i][0];S[i][i]=2*u[i][1];}}}double ss[50][50];for(i=0;i<6;i++)for(j=0;j<6;j++)ss[i][j]=0.0;for(i=1;i<3;i++)for(j=1;j<4;j++){ss[2*i-2][2*j-2]=H[i][j];ss[2*i-2][2*j-1]=N[i][j];ss[2*i-1][2*j-2]=J[i][j];ss[2*i-1][2*j-1]=L[i][j];}i=3;for(j=1;j<4;j++){ss[2*i-2][2*j-2]=H[i][j];ss[2*i-2][2*j-1]=N[i][j];ss[2*i-1][2*j-2]=R[i][j];ss[2*i-1][2*j-1]=S[i][j];}printf("迭代第%d次的雅可比矩阵为：\n",v);for(i=0;i<6;i++){for(j=0;j<6;j++)printf("%10f",ss[i][j]);printf("\n");}printf("\n");printf("\n");Gauss(ss,P_PQ,6);for(i=1;i<nB;i++){u[i][0]=u[i][0]+mn[2*(i-1)];u[i][1]=u[i][1]+mn[2*i-1];}double max;max=fabs(P_PQ[0]);for(i=0;i<=5;i++)if (max<fabs(P_PQ[i]))max=fabs(P_PQ[i]);if(fabs(max)<0.0001){printf("满足精度要求，迭代终止，迭代次数为%d\n",v); printf("\n");}/* 叠代循环的括号 */printf("最终求得的节点电压值为e+jf:\n");for(i=1;i<nB+1;i++)printf("%lf,%lf\n",u[i][1],u[i][0]);printf("\n");printf("\n");double uu[5],Phase[5];for(i=1;i<nB+1;i++){uu[i]=sqrt(u[i][1]*u[i][1]+u[i][0]*u[i][0]); Phase[i]=atan(u[i][0]/u[i][1]);}for(i=1;i<nB+1;i++)printf("%lf,%lf\n",uu[i],Phase[i]);*计算线路功率和平衡节点 PV节点功率*/double P[5],Q[5];double tempP=0.0;double tempQ=0.0;for(i=1;i<nB+1;i++){for(j=1;j<nB+1;j++){tempP+=YG[i][j]*u[j][1]-YB[i][j]*u[j][0]; tempQ+=YG[i][j]*u[j][0]+YB[i][j]*u[j][1];}P[i]=tempP*u[i][1]+tempQ*u[i][0];Q[i]=tempP*u[i][0]-tempQ*u[i][1];}for(i=1;i<nB+1;i++)printf("节点%d注入功率为%lf,%lf\n",i,P[i],Q[i]);/* 支路功率 */double V[4][2];for(i=1;i<5;i++)for(j=0;j<3;j++)V[i][j]=u[i][j];double sP[5][5],sQ[5][5];double dsq,dsp,dp,sumgen;for(i=1;i<NBUS+1;i++){for(j=1;j<NBUS+1;j++){sP[i][j]=0.0;sQ[i][j]=0.0;}}for(l=0; l<nL; l++){i=sL[l].NumI;j=sL[l].NumJ;r=sL[l].R;x=sL[l].X;d1=r*r+x*x;if(fabs(sL[l].K-1.0)<0.000001){/*Normal lines or transformers*/sP[i][j]=V[i][1]*V[i][1]*g-V[i][1]*V[j][1]*(g*cos(V[i][0]-V[j][0])+b*sin(V[i][0]-V[j][0]));sQ[i][j]=-(V[i][1]*V[i][1]*sL[l].B+V[i][1]*V[i][1]*b+V[i][1]*V[j][1]*(g*sin(V[i][0]-V[j][0])-b*cos(V[i ][0]-V[j][0])));sP[j][i]=V[j][1]*V[j][1]*g-V[i][1]*V[j][1]*(g*cos(V[j][0]-V[i][0])+b*sin(V[j][0]-V[i][0]));sQ[j][i]=-(V[j][1]*V[j][1]*sL[l].B+V[j][1]*V[j][1]*b+V[i][1]*V[j][1]*(g*sin(V[j][0]-V[i][0])-b*cos(V[j ][0]-V[i][0])));}else{/*abnormal transformer ratio*/sP[i][j]=V[i][1]*V[i][1]*g/sL[l].B/sL[l].B-V[i][1]*V[j][1]*(g*cos(V[i][0]-V[j][0])/sL[l].B+b*sin(V[i][ 0]-V[j][0])/sL[l].B);sQ[i][j]=-(V[i][1]*V[i][1]*b/sL[l].B/sL[l].B+V[i][1]*V[j][1]*(g*sin(V[i][0]-V[j][0])/sL[l].B-b*cos(V[i ][0]-V[j][0])/sL[l].B));sP[j][i]=V[j][1]*V[j][1]*g-V[i][1]*V[j][1]*(g*cos(V[j][0]-V[i][0])/sL[l].B+b*sin(V[j][0]-V[i][0])/sL[l ].B);sQ[j][i]=-(V[i][1]*V[i][1]*b+V[i][1]*V[j][1]*(g*sin(V[j][0]-V[i][0])/sL[l].B-b*cos(V[j][0]-V[i][0])/sL [l].B));}}/* 输电效率 */dsp=P[4];sumgen=P[4];for(i=0;i<NBUS;i++){dsp+=sB[i].GenP-sB[i].LoadP;dsq+=sB[i].GenQ-sB[i].LoadQ;sumgen+=sB[i].GenP;}dp=dsp/sumgen*100;/* 输出功率情况 */if((fp=fopen("功率情况.txt","w"))==NULL){printf("Can not open the file named '功率情况.txt' \n");exit(0);}fprintf(fp,"---功率情况---\n");fprintf(fp,"平衡节点功率S=%10.5f+ j%10.5f\n",P[4],Q[4]);for(i=1;i<NBUS+1;i++)for(j=1;j<NBUS+1;j++)if(fabs(sP[i][j]-0.0)>0.000001)fprintf(fp,"S(%3d,%-3d)=(%10.5f,j%10.5f)\n",i,j,sP[i][j],sQ[i][j]); fprintf(fp,"网损为%10.5f+j%10.3f,输电效率为%10.3f\n",dsp,dsq,100-dp);if(fp!=NULL) fclose(fp);return 0;}结果：1．导纳阵Y( 1,1 )=( 1.01534, -8.19201) Y( 1,2 )=( -0.56148, 2.30208) Y( 1,3 )=( 0.00000, 3.66667) Y( 1,4 )=( -0.45386, 1.89107) Y( 2,1 )=( -0.56148, 2.30208) Y( 2,2 )=( 1.04225, -4.67651) Y( 2,4 )=( -0.48077, 2.40385) Y( 3,1 )=( 0.00000, 3.66667) Y( 3,3 )=( 0.00000, -3.33333) Y( 4,1 )=( -0.45386, 1.89107) Y( 4,2 )=( -0.48077, 2.40385) Y( 4,4 )=( 0.93463, -4.26159)2．设定电压初值01.1;01;01)0(3)0(3)0(2)0(2)0(1)0(1j jf e j jf e j jf e +=++=++=+ 3．计算功率和电压偏移；27731.0])()([41)0(11)0(1)0(141)0(11)0(1)0(11)0(11)0(1-=++--=-=∆∑∑==j j j j j j jj s s e B f G ff B eG e P P P P0.05097])()([41)0(11)0(1)0(141)0(11)0(1)0(11)0(11)0(1-=----=-=∆∑∑==jj j j jj jj s s e B f G ef B e G f Q Q Q Q同理可算出52596.0)0(22)0(2-=-=∆P P P s ,0196.0)0(22)0(2=-=∆Q Q Q s 5.0)0(33)0(3=-=∆P P P s ,0.02)0(3232)0(3=-=∆U U U s 4．根据求的第一次迭代时雅可比矩阵各元素的公式计算雅可比矩阵各个元素的具体值：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤-----------⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡20000.200066667.300003334.40052691.406629.130208.256148.00001821.182612.456148.030208.266667.3030208.256148.006298.803803.1066667.356148.030208.299265.032104.85．求高斯计算后的修正量：⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆∆∆∆∆=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∆∆∆∆∆∆-0.0000000.1276520.023566-0.108546-0.006511-0.007919-2)0(3)0(3)0(2)0(2)0(1)0(11)0()0(3)0(3)0(2)0(2)0(1)0(1U P Q P Q P e f e f e f J 6．计算各节点电压的一次近似值：12765.010855.000792.010000.197643.099349.0)0(3)0(3)1(3)0(2)0(2)1(2)0(1)0(1)1(1)0(3)0(3)1(3)0(2)0(2)1(2)0(1)0(1)1(1=∆+=-=∆+=-=∆+==∆+==∆+==∆+=f f f f f f fffe e e e e e e e e返回第三步重新迭代，并校验收敛与否，令410-=ε。
华北电力大学潮流上机课程设计报告
课程设计报告( 2014—2015年度第一学期) 名称：电力系统潮流上机院系：电气与电子工程学院班级：电气1215 学号：学生姓名：指导教师：王莉丽设计周数：两周成绩：日期： 2015年1月4日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1．手算2．计算机计算3．思考题三、课程设计总结或结论四、参考文献1.《电力系统分析基础》，李庚银，机械工业出版社，2011年，第一版；2.《电力系统稳态分析》，陈珩，中国电力出版社，2007年，第三版；附录（包括：1.手算过程及结果；2.计算机计算流程图、表格、数据；3.思考题答案）附录1.手算过程及结果2.计算机计算流程图3.完整的潮流程序头文件（定义变量）：#define Bus_Num 9#define Line_Num 9#define Precision 1e-5struct Bus{int No ;float Voltage,Phase,GenP,GenQ,LoadP,LoadQ;int Type;}gBus[Bus_Num];struct Line{int No,No_I,No_J;float R,X,B,k;}gLine[Line_Num];float gY_G[Bus_Num][Bus_Num],gY_B[Bus_Num][Bus_Num];float gDelta_P[Bus_Num-1],gDelta_Q[Bus_Num-1],gDelta_PQ[2*(Bus_Num-1)]; float gJaccobi[2*(Bus_Num-1)][2*(Bus_Num-1)];float gDelta_f[Bus_Num-1],gDelta_e[Bus_Num-1],gDelta_fe[2*(Bus_Num-1)]; float gf[Bus_Num],ge[Bus_Num];主程序：// flow.cpp: 主项目文件#include"stdafx.h"#include"NEquation.h"#include"math.h"#include"stdio.h"#include"config.h"using namespace System;void test(){NEquation ob1;ob1.SetSize(2);ob1.Data(0,0)=1;ob1.Data(0,1)=2;ob1.Data(1,0)=2;ob1.Data(1,1)=1;ob1.Value(0)=4;ob1.Value(1)=6;ob1.Run();printf("x1=%f\n",ob1.Value(0));printf("x2=%f\n",ob1.Value(1));}void GetData()//Read the data{FILE *fp;int i;fp=fopen("E:\\1121960726\\flow\\data\\data.txt","r");if(fp==NULL){printf("Can not open the file named 'data.txt' \n");return;}for(i=0;i<=Bus_Num-1;i++){fscanf(fp,"%d,%f,%f,%f,%f,%f,%f,%d",&gBus[i].No,&gBus[i].Voltage,&gBus[i].Phase, &gBus[i].GenP,&gBus[i].GenQ,&gBus[i].LoadP,&gBus[i].LoadQ,&gBus[i].Type);}for(i=0;i<=Line_Num-1;i++){fscanf(fp,"%d,%d,%d,%f,%f,%f,%f",&gLine[i].No,&gLine[i].No_I,&gLine[i].No_J, &gLine[i].R,&gLine[i].X,&gLine[i].B,&gLine[i].k);}fclose(fp);}void GetYMatrix(){int i,j,bus1,bus2;float r,x,d,g,b,g1,b1,g2,b2,g3,b3;FILE *fp;for(i=0;i<=Bus_Num-1;i++){for(j=0;j<=Bus_Num-1;j++){gY_G[i][j]=0;gY_B[i][j]=0;}}for(i=0; i<=Line_Num-1; i++){if(gLine[i].k==0){bus1=gLine[i].No_I-1;bus2=gLine[i].No_J-1;r=gLine[i].R;x=gLine[i].X;d=r*r+x*x;g=r/d;b=-x/d;gY_G[bus1][bus1]=gY_G[bus1][bus1]+g;gY_G[bus2][bus2]=gY_G[bus2][bus2]+g;gY_G[bus1][bus2]=gY_G[bus1][bus2]-g;gY_G[bus2][bus1]=gY_G[bus2][bus1]-g;gY_B[bus1][bus1]=gY_B[bus1][bus1]+b+gLine[i].B;gY_B[bus2][bus2]=gY_B[bus2][bus2]+b+gLine[i].B;gY_B[bus1][bus2]=gY_B[bus1][bus2]-b;gY_B[bus2][bus1]=gY_B[bus2][bus1]-b; }else{bus1=gLine[i].No_I-1;bus2=gLine[i].No_J-1;x=gLine[i].X;d=r*r+x*x;g=r/d;b=-x/d;g1=g/gLine[i].k;b1=b/gLine[i].k;g2=g*(1-gLine[i].k)/(gLine[i].k*gLine[i].k);b2=b*(1-gLine[i].k)/(gLine[i].k*gLine[i].k);g3=g*(gLine[i].k-1)/gLine[i].k;b3=b*(gLine[i].k-1)/gLine[i].k;gY_G[bus1][bus1]=gY_G[bus1][bus1]+g1+g2;gY_G[bus2][bus2]=gY_G[bus2][bus2]+g1+g3;gY_G[bus1][bus2]=gY_G[bus1][bus2]-g1;gY_G[bus2][bus1]=gY_G[bus2][bus1]-g1;gY_B[bus1][bus1]=gY_B[bus1][bus1]+b1+b2;gY_B[bus2][bus2]=gY_B[bus2][bus2]+b1+b3;gY_B[bus1][bus2]=gY_B[bus1][bus2]-b1;gY_B[bus2][bus1]=gY_B[bus2][bus1]-b1;}}// output the Y matrixfp=fopen("E:\\1121960726\\flow\\data\\ymatrix.txt","w");if(fp==NULL){printf("Can not open the file named 'ymatrix.txt' \n");return ;}fprintf(fp,"---Y Matrix---\n");for(i=0;i<=Bus_Num-1;i++){for(j=0;j<=Bus_Num-1;j++){fprintf(fp,"Y(%d,%d)=(%10.5f,%10.5f)\n",i+1,j+1,gY_G[i][j],gY_B[i][j]);}}fclose(fp);}void SetInitial(){int i;for(i=0;i<=Bus_Num-1;i++){if(gBus[i].Type==3){gf[i]=gBus[i].Voltage*sin(gBus[i].Phase);ge[i]=gBus[i].Voltage*cos(gBus[i].Phase);}else{gf[i]=0;ge[i]=1;}}}void GetUnbalance(){int i,j;FILE *fp;for(i=0;i<=Bus_Num-2;i++){gDelta_P[i]=gBus[i+1].GenP-gBus[i+1].LoadP;if(gBus[i+1].Type==2) //PV节¨2点ì?gDelta_Q[i]=gBus[i+1].Voltage*gBus[i+1].Voltage-(ge[i+1]*ge[i+1]+gf[i+1]*gf[i+1]);elsegDelta_Q[i]=gBus[i+1].GenQ-gBus[i+1].LoadQ;for(j=0;j<=Bus_Num-1;j++){gDelta_P[i]=gDelta_P[i]-ge[i+1]*(gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j])-gf[i+1]*(gY_G[i+1][j]*gf[j]+gY_B[i+1][j]*ge[j]);if(gBus[i+1].Type==1) //PQ节¨2点ì?gDelta_Q[i]=gDelta_Q[i]-gf[i+1]*(gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j])+ge[i+1]*(gY_G[i+1][j]*gf[j]+gY_B[i+1][j]*ge[j]);}}for(i=0;i<=Bus_Num-2;i++) //合?并?é{gDelta_PQ[2*i]=gDelta_P[i];gDelta_PQ[2*i+1]=gDelta_Q[i];}fp=fopen("E:\\1121960726\\flow\\data\\unbalance.txt","w");if(fp==NULL){printf("无法打开文件：”'unbalance.txt' \n");return ;}fprintf(fp,"---Unbalance---\n");for(i=0;i<=2*Bus_Num-3;i++){fprintf(fp,"Unbalance[%d]=%10.5f\n",i+1,gDelta_PQ[i]);}fclose(fp);}void GetJaccobi(){int i,j;float ga[Bus_Num-1],gb[Bus_Num-1];FILE *fp;for(i=0;i<=Bus_Num-2;i++) //计算注入电流{ga[i]=0;gb[i]=0;for(j=0;j<=Bus_Num-1;j++){ga[i]=ga[i]+gY_G[i+1][j]*ge[j]-gY_B[i+1][j]*gf[j];gb[i]=gb[i]+gY_G[i+1][j]*gf[j]+gY_B[i+1][j]*ge[j];}}for(i=0;i<=Bus_Num-2;i++){for(j=0;j<=Bus_Num-2;j++){if(i!=j){gJaccobi[2*i][2*j]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1];gJaccobi[2*i][2*j+1]=gY_G[i+1][j+1]*ge[i+1]+gY_B[i+1][j+1]*gf[i+1];if(gBus[i+1].Type==2) //PV节¨2点ì?{gJaccobi[2*i+1][2*j]=0;gJaccobi[2*i+1][2*j+1]=0;}else//PQ{gJaccobi[2*i+1][2*j]=-gJaccobi[2*i][2*j+1];gJaccobi[2*i+1][2*j+1]=gJaccobi[2*i][2*j];}}else{gJaccobi[2*i][2*j]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1]+gb[i];gJaccobi[2*i][2*j+1]=gY_G[i+1][j+1]*ge[i+1]+gY_B[i+1][j+1]*gf[i+1]+ga[i];if(gBus[i+1].Type==2) //PV节¨2点ì?{gJaccobi[2*i+1][2*j]=2*gf[i+1];gJaccobi[2*i+1][2*j+1]=2*ge[i+1];}else//PQ节点{gJaccobi[2*i+1][2*j]=-gY_G[i+1][j+1]*ge[i+1]-gY_B[i+1][j+1]*gf[i+1]+ga[i];gJaccobi[2*i+1][2*j+1]=-gY_B[i+1][j+1]*ge[i+1]+gY_G[i+1][j+1]*gf[i+1]-gb[i];}}}}fp=fopen(,"w");if(fp==NULL){printf("无法打开文件：” 'jaccobi.txt' \n");return ;}fprintf(fp,"---Jaccobi Matrix---\n");for(i=0;i<=2*Bus_Num-3;i++){for(j=0;j<=2*Bus_Num-3;j++){fprintf(fp,"jaccobi(%d,%d)=%10.5f\n",i+1,j+1,gJaccobi[i][j]);}}fclose(fp);}void GetRevised(){int i,j;FILE *fp;NEquation ob1; //解矩阵方程ob1.SetSize(2*(Bus_Num-1));for(i=0;i<=2*Bus_Num-3;i++)for(j=0;j<=2*Bus_Num-3;j++)ob1.Data(i,j)=gJaccobi[i][j];for(i=0;i<=2*Bus_Num-3;i++)ob1.Value(i)=gDelta_PQ[i];ob1.Run();for(i=0;i<=Bus_Num-2;i++){gDelta_f[i]=ob1.Value(2*i);gDelta_e[i]=ob1.Value(2*i+1);gDelta_fe[2*i]=gDelta_f[i];gDelta_fe[2*i+1]=gDelta_e[i];}fp=fopen("E:\\1121960726\\flow\\data\\revised.txt","w");if(fp==NULL){printf("无法打开文件：” 'revised.txt' \n");return ;}fprintf(fp,"---Revised---\n");for(i=0;i<=2*Bus_Num-3;i++){fprintf(fp,"revised[%d]=%10.5f\n",i+1,gDelta_fe[i]);}fclose(fp);}void GetNewValue(){int i;FILE *fp;for(i=0;i<=Bus_Num-2;i++){gf[i+1]=gf[i+1]+gDelta_f[i];ge[i+1]=ge[i+1]+gDelta_e[i];}fp=fopen("E:\\1121960726\\flow\\data\\newvalue.txt","w");if(fp==NULL){printf("无法打开文件：”'newvalue.txt' \n");return ;}fprintf(fp,"---New Value---\n");for(i=0;i<=Bus_Num-2;i++){fprintf(fp,"f(%d)=%10.5f,e(%d)=%10.5f\n",i+1,gf[i+1],i+1,ge[i+1]); }fclose(fp);}int main(array<System::String ^> ^args){int i,Count_Num;float maxValue;//test();GetData();GetYMatrix();SetInitial();for(Count_Num=0;Count_Num<=100;Count_Num++) {GetUnbalance();GetJaccobi();GetRevised();GetNewValue();maxValue=fabs(gDelta_fe[0]);for(i=1;i<=2*(Bus_Num-1)-1;i++){if(maxValue<fabs(gDelta_fe[i])){maxValue=fabs(gDelta_fe[i]);}}if(maxValue<Precision){break;}}printf("%d\n",Count_Num);for(i=0;i<=Bus_Num-1;i++){printf("%10.5f\n",sqrt(ge[i]*ge[i]+gf[i]*gf[i]));}while(true){}return 0;}计算结果：（1）原始数据：（2）节点导纳矩阵：(3)Jacobi矩阵：（4）NewVoltageValue:(5)不平衡量：（6）修正量：（6）结果：3.回答思考题1.潮流计算的方法有哪些？各有何特点？2. 如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）3. 设计中遇到的问题和解决的办法。
华中科技大学现代电力系统分析潮流计算作业
现代电力系统分析作业------基于Matpower的电力系统潮流计算专业：班级：姓名：学号：目录基于Matpower的电力系统潮流计算 (1)1.本次潮流计算的目的及意义 (1)2.电力系统潮流计算及其意义 (1)3.电力系统潮流计算常规方法 (1)3.1 牛顿-拉夫逊法 (1)3.2 节点电压用直角坐标表示时的牛顿—拉夫逊潮流计算 (3)3.3 牛顿—拉夫逊法潮流计算程序框图 (4)4.选用的潮流计算的系统 (5)5.利用软件matpower计算潮流 (7)5.1 matpower简介 (7)5.2 基态潮流计算 (7)5.21 基态潮流计算条件说明及数据输入 (7)5.22 潮流计算结果 (9)5.23 基态潮流结果计算分析 (10)5.3 最优潮流计算 (11)5.31 最优潮流计算条件说明及数据输入 (12)5.32 最优潮流的理论结果 (12)5.33 最优潮流仿真计算结果 (12)5.23 最优潮流结果计算分析 (13)6.感想与小结 (14)参考文献 (15)附录一、基态潮流计算文件 (15)附录二、最优潮流计算文件 (17)基于Matpower的电力系统潮流计算1.本次潮流计算的目的及意义本次潮流计算的目的及意义主要是了解电力系统潮流计算及其意义，在此基础上，了解电力系统潮流计算的模型以及常规的潮流计算的方法，掌握并熟练使用电力系统潮流计算软件Matpower。
最后，利用Matpower计算美国西部电网WSCC 三机九节点系统的静态潮流及最优潮流并给出分析。
2.电力系统潮流计算及其意义电力系统潮流计算是研究电力系统稳态运行情况的基本电气计算，电力系统潮流计算的任务是根据给定的网络结构及运行条件，求出电网的运行状态，其中包括各母线的电压、各支路的功率分布以及功率损耗等。
潮流计算分为离线计算和在线计算两大类。
离线计算可以用于电力系统调度，并确定系统的运行方式；离线计算的结果还可以用于电力系统规划方案的分析以及优化系统的运行状态；此外离线潮流计算可以作为初值，用于配合系统的故障分析以及稳定性分析。
华北电力大学 电力系统基础 第五章 简单电力系统潮流计算
dU 3ZIZ 3Z ( 3SU22 )
(R jX )( P2 jQ2 ) U2
P2R Q2 X j P2X Q2R U jU
U2
U2
U1 (U2 U)2 U2 arctg U
U2 U
近似计算：
(21.10
j 21.35)
0.0745
j0.0754MVA
S
'
2
(5.62
j3.34) (0.0745
j0.0754)
5.55
j3.26MVA
S
"
1
S'Fra bibliotek2 S1
(5.55
j3.26) (20.15
j13.96)
14.60
j10.70MVA
2 N
1000SN2
QZT
S22
U
2 2
*
U
k
%U
2 N
100SN
PyT
U12
*
P0 1000U
2 N
QyT
U12
*
I0 %SN
100U
2 N
发电厂
PZT
PkU
2 N
S1
2
1000U12 S N2
QZT
U
k
%U
2 N
S1
2
100U
2 2
S
N
PyT
电
电
所
厂
的
的
变
变
压
压
华电潮流上机课程设计报告程序..doc
院系：电气与电子工程学院班级：学号：学生姓名：指导教师：刘宝柱设计周数：成绩：日期：2012年1月5日课程课程设计报告一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1．手算2．计算机计算3．思考题3.1潮流计算的方法有哪些？各有何特点？答：潮流计算的方法主要有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法和P-Q分解法。
它们各自的特点如下：（1）高斯-赛德尔迭代法分为以节点导纳矩阵为基础的高斯-赛德尔迭代法和以以节点阻抗矩阵为基础的高斯-赛德尔迭代法的原理比较简单，要求的数字计算机的内存量也比较小，但收敛性差，当系统规模变大时，迭代次数急剧上升，往往出现迭代不收敛的情况；而阻抗法改善了电力系统潮流计算导纳法德收敛性问题，在当时获得了广泛的应用，但是，阻抗法的主要缺点是占用计算机的内存很大，每次迭代的计算量很大。
当系统不断扩大时，这些缺点就更加突出。
（2）牛顿-拉夫逊法是数学中求解非线性方程式的典型方法，有较好的收敛性。
只要在迭代过程中尽可能保持方程式系数矩阵的稀疏性，就可以大大提高牛顿潮流计算程序的计算效率，牛顿法在收敛性、内存要求、计算速度方面都超过了阻抗法，成为知道目前仍被广泛采用的方法。
（3）P-Q分解法潮流计算派生于以极坐标表示时的牛顿-拉夫逊法，它根据电力系统的特点，抓住主要矛盾，对纯数学的牛顿法进行了改造。
与牛顿法相比，P-Q分解法的修正方程的系数矩阵B’和B”分别是（n-1）和（m-1）的方阵，替代了原有的（n+m-2）阶系数矩阵J；B’、B”在迭代过程中可以保持不变且为对称的系数矩阵，提高了计算速度，降低了对存储容量的要求。
P-Q分解法在计算速度方面有显著地提高，迅速得到了推广。
3.2如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作？（收集哪些数据，如何整理，计算结果如何分析）答：(1)在进行北京城市电网的潮流计算之前需要了解北京城市电网中所有的节点支路的相关数据，并对节点和支路分类。
华北电力大学潮流上机课程设计报告
课程设计报告( 2011-- 2012年度第一学期)名称：电力系统潮流上机院系：电气与电子工程学院班级：电管0902学号：1091140202 学生：指导教师：麻秀设计周数：两周成绩：日期：2011年1月6日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细容见附录）1、手算在此我们同样采取计算机计算这种方便快捷的方法。
具体的方法及结果见附录。
2、计算机计算通过计算机编程，我们可以很方便的求解出相关量，其具体计算方法及计算结果见附录。
3、思考题答案见附录。
三、课程设计总结或结论潮流计算是研究电力系统稳态运行的一种基本运算方法，通过对于潮流的计算，实现对潮流的控制，维护电力系统稳定，防止出现超过约束条件的情况，从而使电力系统免遭崩溃的危险。
最初求解电力系统潮流时大多是用手算，这种计算方法不但耗费时间长，精度差，而且对于大容量的电力系统无法做到很好的控制。
随着电力系统结构日趋复杂，计算量也越来越大，仅靠手算难以求解一些复杂网络的潮流。
计算机作为一种处理数据的工具，其计算速度快，准确率高，存储量大，因此它广泛的应用于大规模复杂化的电力系统网络中。
为了能使计算机能进行潮流计算，就必须编制相应的程序，使其能按照人们在程序中设定的计算流程一步一步的执行，直到计算结果达到要求的精度，得到最终想要的结果。
C语言是一种简单实用的语言，利用它进行程序设计可以大大提高电力系统潮流控制的自动化程度。
本次的潮流上机，就是要求我们利用C语言编程，进行计算机潮流计算。
这既是对我们电力系统潮流计算知识的考验，也是对我们C语言编程能力的考验。
是一次各种知识的交汇综合的运用。
通过这次潮流上机，我对电力系统潮流计算有了更深刻的认识，同时，我的C语言编程能力也得到了提高。
我明白了学习不能仅限于课本，还要做到“学以致用”的道理。
在这次潮流设计中，我遇到了由于多问题，其中既包括由于对潮流知识掌握不牢而导致的问题，也包括由于C语言编程能力有限而导致的问题。
电力系统潮流上机计算报告书写要求
电力系统潮流上机计算报告
系别：
班级：
姓名：
学号：
2012年×月×日
一、程序说明
包括：程序设计思想、程序流程图、程序使用说明。
二、给定题目的手算过程（迭代两次）
包括：原题目、节点导纳矩阵、雅克比矩阵、第一次和第二次迭代
结果。
三、给定题目的程序计算结果
包括：原题目、节点导纳矩阵、雅克比矩阵、程序输入和输出文件（误差0.0001）。
四、编程特色与创新
包括：程序能够完成的基本功能；程序能够完成的高级功能（如：是否包括非标准变比变压器支路，是否采用了稀疏矩阵技术，
是否增加了人机对话界面，程序的通用性和实用性如何）。
五、总结
包括：手算结果与程序计算结果的分析比较；本次上机体会，如：独立编程体会、跟踪调试技能的掌握情况、C语言中结构体、
指针、文件输入输出的掌握情况等。
报告要求：
1．报告中除上面的第“三”项外，其他部分必须手写（最好使用黑色水笔）。
2．报告统一采用A4打印纸书写（留出页边距: 1.5~2厘米）。
不使用实验报告纸。
3．封面按上述格式书写。
4．装订统一在左侧1厘米，二个钉。
5．上述五部分内容必须齐全，各部分内容可以扩充。
6．报告书写要求字迹清楚，不得潦草，流程图中的框要用尺子画。
7．报告必须与本人提交程序吻合，否则取消成绩。
8．报告不得有雷同，否则全部取消成绩。
潮流计算实验报告分析
一、实验背景与目的电力系统潮流计算是电力系统分析中的一个重要环节，它通过对电力系统网络中功率和电压的分布进行计算，以评估系统的运行状态。
本实验旨在通过实际操作，加深对电力系统潮流计算原理和方法的理解，并掌握使用PSASP、ETAP等软件进行潮流计算的基本技能。
二、实验原理与方法1. 基本原理潮流计算主要基于基尔霍夫电流定律和基尔霍夫电压定律，通过求解电力系统网络中的功率和电压分布，得到各节点电压、线路电流和设备功率等参数。
2. 计算方法常用的潮流计算方法包括牛顿-拉夫逊法、快速分解法、迭代法等。
本实验采用牛顿-拉夫逊法进行潮流计算。
3. 实验步骤（1）建立电力系统网络模型，包括节点、线路、变压器等元件；（2）设置各节点电压初始值和负荷功率；（3）计算网络中各支路功率和节点电压，判断是否满足功率平衡和电压平衡；（4）根据功率平衡和电压平衡条件，修正节点电压，重复步骤（3）直至满足收敛条件。
三、实验过程与结果分析1. 实验数据本实验采用某实际电力系统网络进行计算，网络包括10个节点、15条线路和3个变压器。
2. 实验步骤（1）根据实验数据，建立电力系统网络模型；（2）设置各节点电压初始值和负荷功率；（3）使用PSASP软件进行潮流计算；（4）分析计算结果，包括节点电压、线路电流和设备功率等。
3. 结果分析（1）节点电压分布合理，各节点电压满足运行要求；（2）线路电流分布均匀，线路负载率在合理范围内；（3）设备功率分配合理，满足电力系统运行需求。
四、实验总结与讨论1. 实验总结本实验通过实际操作，加深了对电力系统潮流计算原理和方法的理解，掌握了使用PSASP软件进行潮流计算的基本技能。
2. 讨论（1）实验中，节点电压初始值设置对计算结果有较大影响，需要根据实际情况进行设置；（2）潮流计算结果受网络拓扑结构、元件参数和负荷分布等因素的影响，需要综合考虑；（3）在实际工程应用中，应根据具体情况选择合适的潮流计算方法，以保证计算结果的准确性和可靠性。
华北电力大学 电力系统稳态潮流上机计算程序结果
unsigned int i,j,k;
//消元过程
for(i=0;i<Dimension;i++)
{
//规格化过程:(每次消元前,先规格化,以便以下各行消元时,消元系数直接取待消
//列元素的值即可,也便于回代过程,而运算量并不增加)
for( j = i+1; j < Dimension; j++ )
D[y][y].diandao+=lineg;//尾节点自导纳中电导
D[y][x].dianna=D[x][y].dianna;//首节点与尾节点之间导纳与尾节点与首节点之间导纳相等
D[y][x].diandao=D[x][y].diandao;
}
if(gLineData[i].Type==2)//对线路为变压器的类型 对变压器进行π型等值 计算方法与导线型相同
{
if( FactorMatrix[j][i] != 0 ) //如果第j行第i列元素本就是0,则不需本列对应的消元过程
{
for( k = i + 1; k < Dimension; k++ ) //当FactorMatrix[i][k]=0,a[j][k]值不变,可省去运算
{
D[x][y].dianna+=(-1)*lineb/gLineData[i].Conductance;
D[x][y].diandao+=(-1)*lineg/gLineData[i].Conductance;
D[x][x].diandao+=lineg/gLineData[i].Conductance+(1-gLineData[i].Conductance)*lineg/(gLineData[i].Conductance*gLineData[i].Conductance);
电力网潮流计算实训报告
一、实训目的本次电力网潮流计算实训旨在使学生掌握电力系统潮流计算的基本原理和方法，提高学生对电力系统运行状态的分析能力，培养学生运用所学知识解决实际问题的能力。
通过实训，使学生能够熟练运用MATLAB等工具进行电力系统潮流计算，并能够根据计算结果对电力系统运行状态进行分析和评价。
二、实训内容1. 电力系统潮流计算的基本原理电力系统潮流计算是一种分析电力系统稳态运行状态的方法，其主要任务是确定电力系统中各节点的电压、电流和功率分布。
电力系统潮流计算的基本原理如下：（1）节点电压方程：根据节点电压和注入功率之间的关系，建立节点电压方程。
（2）支路功率方程：根据支路两端电压和电流之间的关系，建立支路功率方程。
（3）潮流计算方法：采用迭代方法求解节点电压方程和支路功率方程，得到电力系统中各节点的电压、电流和功率分布。
2. 电力系统潮流计算实训本次实训以MATLAB为工具，进行电力系统潮流计算。
实训步骤如下：（1）建立电力系统模型：根据实训要求，建立电力系统模型，包括节点、支路、变压器等。
（2）设置电力系统参数：根据电力系统模型，设置各节点的注入功率、支路参数、变压器参数等。
（3）进行潮流计算：运用MATLAB中的电力系统分析工具箱（Power System Toolbox）进行潮流计算。
（4）分析计算结果：对潮流计算结果进行分析，包括节点电压、电流、功率分布、网络损耗等。
三、实训结果与分析1. 潮流计算结果通过潮流计算，得到以下结果：（1）节点电压：各节点电压的幅值和相角。
（2）电流：各支路电流的幅值和相角。
（3）功率分布：各支路的有功功率和无功功率。
（4）网络损耗：电力系统中的网络损耗。
2. 结果分析（1）节点电压分析：分析各节点电压是否满足要求，是否存在过电压或欠电压现象。
（2）电流分析：分析各支路电流是否满足要求，是否存在过载现象。
（3）功率分布分析：分析电力系统中功率分布是否合理，是否存在不平衡现象。
第一章华北电力大学 电力系统潮流计算1new
U i
k 1
(1-17) 上式是该算法最基本的迭代计算公式。 其迭代收敛的判据是 maxU k 1 U k
i i i
i 1 n 1 Pi jQi Yi1 U1s ( Yij U i( k 1) Yij U i( k ) i 2,3,, n Yii ( k ) j 2 j i 1 Ui
第三节 潮流计算的几种基本方法
一 高斯-塞德尔法 以导纳矩阵为基础，并应用高斯-塞 德尔迭代的算法是电力系统应用最早的 潮流计算方法。
三．潮流计算的几种基本方法
高斯—塞德尔迭代法原理
已知方程组 (1) 1
( 2)
x1 0 3 0.3333
( x21) 0 2 0.6667 3
第二节 潮流计算问题的数学模型
对这样的线性网络一般采用节点电压 法进行分析。节点电压与节点注入电流 之间的关系为:
或
YU I
U Z I
第二节 潮流计算问题的数学模型
式中：
I1 I I 2 , In . U1 . U U 2 U. n
j 1
Ui
n
或
U i Z ij
j 1
Pj jQ j
i 1,2,, n
Uj
(1-7)
第二节 潮流计算问题的数学模型
上两式是潮流计算问题的基本方程式， 是一个以节点电压为变量的非线性代数 方程组。而采用节点功率作为节点注入 量是造成方程组呈非线性的根本原因。 由于方程组为非线性的，因此必须采用 迭代方法进行数值求解。 根据对方程组的不同处理方式，形成 了不同的潮流算法。
潮流计算实验报告
潮流计算实验报告潮流计算实验报告潮流计算是电力系统运行中的重要工具，用于分析电力系统中各节点的电压、功率等参数，以确保电力系统的稳定运行。
本次实验旨在通过潮流计算方法，对一个简化的电力系统进行分析，探讨电力系统的稳定性和可靠性。
1. 实验背景电力系统是一个复杂的网络，由发电厂、输电线路、变电站和用户组成。
在电力系统中，电流和电压的分布是非常重要的，因为它们直接影响到电力系统的稳定性和可靠性。
潮流计算是一种基于电力系统的拓扑结构和电气参数，通过求解节点电压和功率的方程组，来分析电力系统中各节点的电压、功率等参数的方法。
2. 实验目的本次实验的目的是通过潮流计算方法，对一个简化的电力系统进行分析，了解电力系统的稳定性和可靠性。
具体目标包括：- 分析电力系统中各节点的电压、功率等参数；- 研究电力系统中负荷变化对电压和功率的影响；- 探讨电力系统中的潮流分布情况。
3. 实验过程本次实验采用Matlab软件进行潮流计算。
首先，根据给定的电力系统拓扑结构和电气参数，建立电力系统的节点电压和功率方程组。
然后，通过求解该方程组，得到电力系统中各节点的电压和功率等参数。
最后，根据求解结果，分析电力系统中的潮流分布情况。
4. 实验结果通过潮流计算，得到了电力系统中各节点的电压和功率等参数。
根据实验结果，可以得出以下结论：- 在电力系统中，电压和功率的分布是不均匀的，不同节点的电压和功率存在差异；- 负荷变化会对电力系统中的电压和功率产生影响，负荷增加会导致电压下降，功率增加；- 电力系统中存在潮流集中的现象，即部分节点的潮流较大，而其他节点的潮流较小。
5. 实验分析通过对实验结果的分析，可以得出以下结论：- 电力系统中的电压和功率分布不均匀，这是由于电力系统中各节点的拓扑结构和电气参数的差异所导致的；- 负荷变化对电力系统的稳定性和可靠性具有重要影响，负荷增加会导致电力系统中的电压下降，功率增加，从而可能引发电力系统的故障；- 电力系统中的潮流集中现象可能会导致部分节点的负荷过载，从而影响电力系统的稳定运行。
潮流计算报告
潮流计算报告一、系统结构图：二、网络参数：1、支路参数：网络类型支路编号电压等级（kV）装机容量(MW)导线的技术参数1r（Ω）1x（Ω）1b（）环网1-222010013.6 125.5 67.85 1-3 8.321 130.5 52.24 3-5 10.2 128.8 74.99 2-3 8.5 105.4 28.36 1-4 7.579 129.6 51.45 4-5 13.84 125.31 2.78辐射网1-2————2 4——2-3 3 63-4 4 82-5 1 25-6 4 42、节点参数：节点类型节点编号发电功率（MW）负荷视在功率环网1 0 未知(平衡节点)2 100 0（PV节点）3 0 15+9.4i4 0 27+6i5 0 35.5+25.5i辐射网1 0 未知(平衡节点)2 0 4+2i3 0 6+3.2i4 0 3+1.44i5 0 4+3.2i6 0 2+1.1i三、潮流计算流程图:Un 2 1.1];FT=[%首端末端4 33 26 55 22 1];RX=[% R X4 83 64 41 22 4];NN=size(PQ,1); %节点数NB=size(FT,1); %支路数数V=PQ(:,1); %V初始电压相量maxd=1k=1while maxd>0.0001k=k+1;PQ2=PQ; %每一次迭代各节点的注入有功和无功相同PL=0.0;for i=1:NBkf=FT(i,1); %前推始节点号kt=FT(i,2); %前推终节点号x=(PQ2(kf,2)^2+PQ2(kf,3)^2)/V(kf)/V(kf);%计算沿线电流 /平方Alosss(i,1)=RX(i,1)*x; %计算线路有功损耗 /MWlosss(i,2)=RX(i,2)*x; %计算线路无功损耗/MWPQ1(i,1)=PQ2(kf,2)+RX(i,1)*x; %计算支路首端有功/MW RX(i,1)*RPQ1(i,2)=PQ2(kf,3)+RX(i,2)*x; %计算沿支路的无功/MW RX(i,2)*XPQ2(kt,2)= PQ2(kt,2)+PQ1(i,1); %用PQ1去修正支路末端节点的有功P 单位MWPQ2(kt,3)= PQ2(kt,3)+PQ1(i,2); %用PQ1去修正支路末端节点的有功Q 单位Mvarendangle(1)=0.0;for i=NB:-1:1kf=FT(i,2); %回代始节点号kt=FT(i,1); %回代终节点号dv1=(PQ1(i,1)*RX(i,1)+PQ1(i,2)*RX(i,2))/V(kf); %计算支路电压损耗的纵分量dv1dv2=(PQ1(i,1)*RX(i,2)-PQ1(i,2)*RX(i,1))/V(kf); %计算支路电压损耗的横分量dv2V2(kt)=sqrt((V(kf)-dv1)^2+dv2^2); %计算支路末端电压/kV angle(kt)=angle(kf)+atand(dv2/(V(kf)-dv1)); %计算支路endmaxd=abs(V2(2)-V(2));V2(1)=V(1);for i=3:1:NNif abs(V2(i)-V(i))>maxd;maxd=abs(V2(i)-V(i));endendfullloss(1,1)=0;%计算线路总损耗fullloss(1,2)=0;finalPQ=max(PQ1);for i=1:NBfullloss(1,1)=fullloss(1,1)+losss(i,1);fullloss(1,2)=fullloss(1,2)+losss(i,2);enddisp('辐射网迭代次数：')kdisp('辐射网系统电压差精度：')maxddisp('辐射网系统末端节点有功和无功：')finalPQ %潮流分布即支路首端潮流MVAdisp('辐射网系统总功率损耗：')fullloss %线路总损耗MVAdisp('辐射网系统各支路功率损耗：')losss %各支路损耗MVAdisp('辐射网系统各节点电压幅值：')V=V2 %节点电压模计算结果kVdisp('辐射网系统各节点电压相角：')angle %节点电压角度计算结果单位度endclcdisp('辐射网迭代次数：')kdisp('辐射网系统电压差精度：')maxddisp('辐射网系统末端节点有功和无功/MVA：')FinPQ=finalPQ(1,1)+finalPQ(1,2)*j %潮流分布即支路首端潮流MVA disp('辐射网系统总功率损耗/MVA：')Fulloss=fullloss (1,1)+fullloss(1,2)*j %线路总损耗MVA disp('辐射网系统各支路功率损耗/MVA：')for(a=1:5)LOSS=losss (a,1)+losss(a,2)*j %各支路损耗MVAenddisp('辐射网系统各节点电压幅值/KV：')V=V2 %节点电压模计算结果kVdisp('辐射网系统各节点电压相角：')angle %节点电压角度计算结果单位度n=5; %input('节点数');nl=6; %input('支路数');isb=1; %input('平衡母线节点号');pr=0.000001; %input('误差精度：pr=');B1=[1,2,13.6+125.5i,0.00006785i,1,0;1,3,8.321+130.5i,0.00005224i,1,0;3,5,10.2+128.8i,0.00007499i,1,0;2,3,8.5+105.4i,0.00002836i,1,0;1,4,7.579+129.6i,0.00005145i,1,0;4,5,13.84+125.31i,0.0000278i,1,0]; %input('由支路参数形成的矩阵');B2=[-FinPQ,0,Un,0,0,1;100,0,Un,Un,0,3;0,15+9.4i,Un,0,0,2;0,27+6i,Un,0,0,2;0,35.5+25.5i,Un,0,0,2]; %input('各节点参数形成的矩阵');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1= zeros(nl);%对各矩阵置零%－－－－－－－修改部分－－－－－－－－－－－－ym=1;SB=100;UB=Un; %定义视在功率和电压基值if ym~=0 %若不是标幺值YB=SB./UB./UB; %定义导纳标幺值BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB; %切换为阻抗标幺值B1(i,4)=B1(i,4)./YB; %切换为导纳标幺值enddisp('支路矩阵B1=');sparseB1=sparse(B1);disp(sparseB1) %输出标幺值稀疏矩阵B1disp('-----------------------------------------------------');for i=1:nB2(i,1)=B2(i,1)./SB; %切换为视在功率标幺值B2(i,2)=B2(i,2)./SB; %切换为视在功率标幺值B2(i,3)=B2(i,3)./UB; %切换为电压标幺值B2(i,4)=B2(i,4)./UB; %切换为电压标幺值B2(i,5)=B2(i,5)./SB; %切换为视在功率标幺值enddisp('节点矩阵B2=');sparseB2=sparse(B2);disp(sparseB2) %输出标幺值稀疏矩阵B2enddisp('-----------------------------------------------------');% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于1侧p=B1(i,1);q=B1(i,2);elsep=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; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵disp('导纳矩阵 Y=');sparseY=sparse(Y);disp(sparseY) %输出导纳稀疏矩阵disp('-----------------------------------------------------');%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部for i=1:ne(i)=real(B2(i,3)); %给定i节点初始电压的实部f(i)=imag(B2(i,3)); %给定i节点初始电压的虚部V(i)=B2(i,4); %PV节点电压给定模值endfor i=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%=================================================================== P=real(S);Q=imag(S); %定义有功功率和无功功率ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %定义迭代次数ICT1和不满足精度要求的节点个数IT2while IT2~=0 %仍有不满足精度要求的节点IT2=0;a=a+1; %IT2置零for i=1:nif i~=isb %非平衡节点C(i)=0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fi Σ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-ei Σ(Gij*fj+Bij*ej)%求P',Q'V2=e(i)^2+f(i)^2; %电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素 =========if B2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算 =================%================= 求取Jacobi矩阵 ===================for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %X1=N(i,j1)=dDP(i)/de(j1)X2=B(i,j1)*e(i)-G(i,j1)*f(i); %X2=H(i,j1)=dDP(i)/df(j1)X3=X2; %X2=H(i,j1)=dDP(i)/df(j1)=X3=M(i,j1)=dDQ(i)/de(j1)X4=-X1; %X1=N(i,j1)=dDP(i)/de(j1)=-X4=-L(i,j1)=-dDQ(i)/df(j1)p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;%扩展列△QJ(m,q)=X1;J(m,N)=DP;q=q+1;%扩展列△PJ(p,q)=X4;J(m,q)=X2; %对Jacobi矩阵赋值elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%X1=N(i,i)=dDP(i)/de(i)X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%X2=H(i,i)=dDP(i)/df(i)X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); %X3=M(i,i)=dDQ(i)/de(i)X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%X4=L(i,i)=dDQ(i)/df(i)p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q m=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2; %对Jacobi矩阵赋值endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========DP=P(i)-P1; % PV节点有功误差DV=V(i)^2-V2; % PV节点电压误差for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %X1=N(i,j1)=dDP(i)/de(j1)X2=B(i,j1)*e(i)-G(i,j1)*f(i); %X2=H(i,j1)=dDP(i)/df(j1)X5=0;X6=0; %X5=R(i,j1)=X6=S(i,j1)=0p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;%扩展列△V m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;%扩展列△PJ(m,q)=X2; %对Jacobi矩阵赋值elseif j1==i&j1~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%X1=N(i,i)=dDP(i)/de(i)X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%X2=H(i,i)=dDP(i)/df(i)X5=-2*e(i); % X5=R(i,i)=-2e(i) X6=-2*f(i); % X6=F(i,i)=-2f(i) p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;%扩展列△V m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;%扩展列△PJ(m,q)=X2; %对Jacobi矩阵赋值endendendendend%========= 以上为求雅可比矩阵的各个元素 =====================for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即 N=2*n+1扩展列△P、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%========================================================== ==k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行前的k4行消去for k2=k1:N1 % k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for 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;endendend%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵(利用线性代数求解电压实部与虚部)=====for k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N); %修改节点电压实部k1=k+1;f(L)=f(L)-J(k1,N); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解"w=-J*V"disp('迭代次数');disp(ICT1);disp('没有达到精度要求的个数');disp(ICT2);disp('-----------------------------------------------------'); for k=1:nV(k)=sqrt(e(k)^2+f(k)^2); %计算实际电压大小sida(k)=atan(f(k)./e(k))*180./pi;%计算实际电压相角E(k)=e(k)+f(k)*j; %计算实际电压相量end%=============== 计算各输出量 ===========================disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');sparseE=sparse(E);disp(sparseE);EE=E*UB;disp('-----------------------------------------------------'); disp('各节点的实际电压EE为(节点号从小到大排列)：');disp(sparseEE);disp('-----------------------------------------------------'); disp('各节点的电压标幺值幅值V为(节点号从小到大排列)：');sparseV=sparse(V);disp(sparseV);disp('-----------------------------------------------------'); VV=V*UB;disp('各节点的电压幅值VV为(节点号从小到大排列)：');sparseVV=sparse(VV);disp(sparseVV);disp('-----------------------------------------------------'); disp('各节点的电压相角为(节点号从小到大排列)：');sparsesida=sparse(sida);disp(sparsesida)disp('-----------------------------------------------------'); for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));%计算电流的共轭endS(p)=E(p)*C(p);%计算节点的视在功率enddisp('各节点的功率标幺值S为(节点号从小到大排列)：');disp(sparseS);disp('-----------------------------------------------------');disp('各节点的功率实际值SS为(节点号从小到大排列)：');SS=S*SB;sparseSS=sparse(SS);disp(sparseSS);disp('-----------------------------------------------------');disp('各条支路的功率损耗S标幺值和实际值SS为(顺序支路参数矩阵顺序一致)：'); HDDS=0;for i=1:nlp=B1(i,1);q=B1(i,2);Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(E(p)-E(q))*(conj(E(p))-con j(E(q)))*conj(1./(B1(i,3))))+E(q)*(conj(E(q))*conj(B1(i,4)./2));ZF1=['S(',num2str(p),',',num2str(q),')=',num2str(Si(p,q))];disp(ZF1);SSi(p,q)=Si(p,q)*SB;%计算各条支路的消耗功率实际值SSiHDDS=HDDS+SSi(p,q);ZF=['SS(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);enddisp('环网总网损为；');ZH=['HDDS=',num2str(HDDS)];disp(ZH);%环网总网耗ZSS=HDDS+Fulloss;disp('总网损为：');ZSS五、额定电压不同时对系统参数的分析：（1）额定电压为240V：辐射网：环网：（2）额定电压为220V：辐射网：环网：（3）额定电压为200V：辐射网：环网：结果分析：240V时总损耗为0.8856KVA 220V时为0.95KVA 200V时为0.9828KVA电压等级越高，损耗越小。
潮流上机课程设计报告华电
《电力系统潮流上机》课程设计报告院系班级：学号：学生姓名：指导教师：设计周数成绩：日期：年月日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录）1.手算： 要求应用牛顿-拉夫逊法或P-Q 分解法手算求解，要求精度为0.001MW 。
节点1为平衡节点，电压︒∠=00.11U ，节点2为PQ 节点，负荷功率6.08.0~2j S +=，节点3是PV 节点，1.1,4.033==U P ，两条支路分别为04.001.013j Z +=，2.005.012j Z +=，对地支路33.030j y =。
2.计算机计算：编写潮流计算程序，要求如下：2.1据给定的潮流计算任务书整理潮流计算的基础数据：节点的分类，线路模型，等值变压器模型，电压等级的归算，标幺值的计算；2.2基础数据的计算机存储：节点数据，支路数据（包括变压器）； 2.3用牛顿-拉夫逊法计算；2.4根据所选潮流计算方法画流程图，划分出功能模块，有数据输入模块，导纳阵形成模块，解线性方程组模块，计算不平衡功率模块，形成雅可比矩阵模块，解修正方程模块，计算线路潮流，网损，PV 节点无功功率和平衡节点功率，数据输出模块； 2.5据上述模块编制程序并上机调试程序，得出潮流计算结果； 2.6源程序及其程序中的符号说明集、程序流图简单系统如下图所示，支路数据如下：41.01.012j z +=，3.013j z =，5.012.014j z +=，40.008.024j z +=01528.01,202,10j y y ==，0192.01,404,10j y y ==，01413.02,404,20j y y ==1.1=k节点数据如下：18.030.01j S --= ，13.055.02j S --= ， 5.03=S ，10.13=U ，o U 005.14∠= 324,10y 1,40y1)节点导纳阵#include <stdio.h> #include <math.h> #include <fstream.h> #include "LF.h"//form node conductance matrixintMakeY( intnB, intnL, Line* sL, double** YG, double** YB ) {inti,j,l;double r,x,d1,g,b,t;for(i=0;i<nB;i++)for(j=0;j<nB;j++){ YG[i][j]=0.0;YB[i][j]=0.0;}for(i=0;i<nL;i++){r=sL[i].R;x=sL[i].X;g=r/(r*r+x*x);b=-x/(r*r+x*x);switch(sL[i].Type){case 1://Linebreak;case 2://Transformerg*=1/sL[i].K;b*=1/sL[i].K;break;}YG[sL[i].NumI][sL[i].NumI]+=g;YG[sL[i].NumJ][sL[i].NumJ]+=g;YG[sL[i].NumI][sL[i].NumJ]-=g;YG[sL[i].NumJ][sL[i].NumI]-=g;YB[sL[i].NumI][sL[i].NumI]+=b+sL[i].B;YB[sL[i].NumJ][sL[i].NumJ]+=b+sL[i].B;YB[sL[i].NumI][sL[i].NumJ]-=b;YB[sL[i].NumJ][sL[i].NumI]-=b; }printf("实部：\n");for(i=0;i<nB;i++){for(j=0;j<nB;j++)printf("%lf\t",YG[i][j]);printf("\n");}printf("虚部：\n");for(i=0;i<nB;i++){for(j=0;j<nB;j++)printf("%lf\t",YB[i][j]);printf("\n");}/* Check the Y matrix */ofstreamfout("out.txt");fout<< "--------------Y Matrix--------------------" <<endl;for(i=0;i<nB;i++){for(j=0;j<nB;j++)fout<< YG[i][j] << "+j" << YB[i][j] << "\t";fout<<endl;}fout.close();return 0;}2）计算功率不平衡量#include <stdio.h>#include <math.h>#include <fstream.h>#include "LF.h"//form delta p and delta qintCalDeltaPQ( intnpv, intnpq, Bus* bus, double** YG, double** YB, int* p_Jtobus, double* deltaf ){intk,i,j;for(k=0;k<npv+npq*2;k++){ i=p_Jtobus[k];if(k<npv){ deltaf[k]=bus[i].GenP-bus[i].LoadP;for(j=0;j<npv+npq+1;j++){deltaf[k]-=bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j ]*sin(bus[i].Phase-bus[j].Phase));}printf("PV节点%d的有功功率是%lf\n",i,deltaf[k]);}if(k<npq+npv&&k>=npv){ deltaf[k]=bus[i].GenP-bus[i].LoadP;for(j=0;j<npv+npq+1;j++){deltaf[k]-=bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j ]*sin(bus[i].Phase-bus[j].Phase));}printf("PQ节点%d的有功功率是%lf\n",i,deltaf[k]);}if(k<npq*2+npv&&k>=npv+npq){deltaf[k]=bus[i].GenQ-bus[i].LoadQ;for(j=0;j<npv+npq+1;j++){deltaf[k]-=bus[i].Volt*bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j ]*cos(bus[i].Phase-bus[j].Phase));}printf("PQ节点%d的无功功率是 %lf\n",i,deltaf[k]);}}return 0;}3）雅各比矩阵的计算/*Purpose: for undergraduate courseTask: Load FlowCopyright @ NCEPU, Liu Chongru*/#include <stdio.h>#include <math.h>#include <fstream.h>#include "LF.h"//form Jacobian matrixintFormJacobian( intnpv, intnpq, Bus* bus, double** YG, double** YB, int* p_Jtobus, double** Jac ){ intnp = npv+npq,j,k,i,m;//TODOdouble a[14],q[14];for(k=0;k<npv+npq*2;k++){i=p_Jtobus[k];a[i]=0;q[i]=0;if(k<np)//H N{for(j=0;j<np+1;j++)if(j!=i){a[i]+=bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j]*cos(bus[i].Phase-bu s[j].Phase));q[i]+=bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*sin(bus[i].Phase-bu s[j].Phase));}for(m=0;m<npv+npq*2;m++){j=p_Jtobus[m];if(j!=i){if(m<np)Jac[k][m]=bus[i].Volt*bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j]*cos (bus[i].Phase-bus[j].Phase));//Form HelseJac[k][m]=bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*sin (bus[i].Phase-bus[j].Phase));//Form N}else if(j==i){if(m<np)Jac[k][m]=-bus[i].Volt*a[i];//Form HelseJac[k][m]=bus[i].Volt*q[i]+2*bus[i].Volt*bus[i].Volt*YG[i][j];//Form N}}}else{for(j=0;j<np+1;j++)if(j!=i){a[i]+=bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j]*cos(bus[i].Phase-bu s[j].Phase));q[i]+=bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*sin(bus[i].Phase-bu s[j].Phase));}for(m=0;m<npv+npq*2;m++){j=p_Jtobus[m];if(j!=i){if(m<np)Jac[k][m]=-bus[i].Volt*bus[j].Volt*(YG[i][j]*cos(bus[i].Phase-bus[j].Phase)+YB[i][j]*si n(bus[i].Phase-bus[j].Phase)); //Form JelseJac[k][m]=bus[i].Volt*bus[j].Volt*(YG[i][j]*sin(bus[i].Phase-bus[j].Phase)-YB[i][j] *cos(bus[i].Phase-bus[j].Phase)); //Form L }else if(j==i){ if(m<np)Jac[k][m]=bus[i].Volt*q[i];elseJac[k][m]=bus[i].Volt*a[i]-2*bus[i].Volt*bus[i].Volt*YB[i][j];} }}}for(i=0;i<np+npq;i++){for(int j=0;j<np+npq;j++){printf("%d %d %f ",i,j,Jac[i][j]);}printf("\n");}//Output the matrix to check the Jacobian matrixofstreamfout("out.txt",ios::app);fout<< "--------------Jacobian Matrix--------------------" <<endl;for(i=0; i<np+npq;i++ ){for(j=0; j<np+npq; j++ )fout<<Jac[i][j] << "\t";fout<<endl;}fout.close();return 0;}4）线路损耗//8.calculate the power flowdouble* p_Pij, *p_Qij, *p_Pji, *p_Qji;p_Pij = new double[nL];p_Qij = new double[nL];p_Pji = new double[nL];p_Qji = new double[nL];int x1,x2;for( i=0; i<nL; i++ ){ x1=line[i].NumI;x2=line[i].NumJ;if(line[i].Type==1){p_Pij[i]=bus[x1].Volt*bus[x1].Volt*(-YG[x1][x2])-bus[x1].Volt*bus[x2].Volt*((-YG[x1][x2 ])*cos(bus[x1].Phase-bus[x2].Phase)+(-YB[x1][x2])*sin(bus[x1].Phase-bus[x2].Phase));p_Qij[i]=-bus[x1].Volt*bus[x1].Volt*(line[i].B+(-YB[x1][x2]))-bus[x1].Volt*bus[x2]. Volt*((-YG[x1][x2])*sin(bus[x1].Phase-bus[x2].Phase)-(-YB[x1][x2])*cos(bus[x1].Phase-bu s[x2].Phase));p_Pji[i]=bus[x2].Volt*bus[x2].Volt*(-YG[x2][x1])-bus[x2].Volt*bus[x1].Volt*((-YG[x2][x1 ])*cos(bus[x2].Phase-bus[x1].Phase)+(-YB[x2][x1])*sin(bus[x2].Phase-bus[x1].Phase));p_Qji[i]=-bus[x2].Volt*bus[x2].Volt*(line[i].B+(-YB[x2][x1]))-bus[x2].Volt*bus[x1]. Volt*((-YG[x2][x1])*sin(bus[x2].Phase-bus[x1].Phase)-(-YB[x2][x1])*cos(bus[x2].Phase-bu s[x1].Phase));}else{p_Pij[i]=bus[x1].Volt*bus[x1].Volt*(-YG[x1][x2])/line[i].K-bus[x1].Volt*bus[x2].Vol t*((-YG[x1][x2])*cos(bus[x1].Phase-bus[x2].Phase)+(-YB[x1][x2])*sin(bus[x1].Phase-bus[x 2].Phase));p_Qij[i]=-bus[x1].Volt*bus[x1].Volt*((-YB[x1][x2])/line[i].K+line[i].B)-bus[x1].Volt*bu s[x2].Volt*((-YG[x1][x2])*sin(bus[x1].Phase-bus[x2].Phase)-(-YB[x1][x2])*cos(bus[x1].Ph ase-bus[x2].Phase));p_Pji[i]=bus[x2].Volt*bus[x2].Volt*(-YG[x2][x1]*line[i].K)-bus[x2].Volt*bus[x1].Volt*(( -YG[x2][x1])*cos(bus[x2].Phase-bus[x1].Phase)+(-YB[x2][x1])*sin(bus[x2].Phase-bus[x1].P hase));p_Qji[i]=-bus[x2].Volt*bus[x2].Volt*((-YB[x2][x1])*line[i].K+line[i].B)-bus[x2].Volt*bu s[x1].Volt*((-YG[x2][x1])*sin(bus[x2].Phase-bus[x1].Phase)-(-YB[x2][x1])*cos(bus[x2].Ph ase-bus[x1].Phase));}}//p and q of PH bus and PV busint s=0;double p[9],q[9],Ps[9],Qs[9],PS=0,QS=0;for( i=0; i<nB; i++ ){p[i]=0;q[i]=0;for(int j=0; j<nB; j++ ){p[i]+=(bus[j].Volt*(YG[i][j])*cos(bus[j].Phase)-bus[j].Volt*(YB[i][j])*sin(bus[j].Phas e));q[i]-=(bus[j].Volt*(YG[i][j])*sin(bus[j].Phase)+bus[j].Volt*(YB[i][j])*cos(bus[j].Phase ));}Ps[i]=bus[i].Volt*cos(bus[i].Phase)*p[i]-bus[i].Volt*sin(bus[i].Phase)*q[i]; Qs[i]=bus[i].Volt*cos(bus[i].Phase)*q[i]+bus[i].Volt*sin(bus[i].Phase)*p[i];} for(i=0;i<nB;i++){ PS+=Ps[i];QS+=Qs[i];}printf("PS=%7.7f,QS=%7.7f\n",PS,QS);}//lossdoublePloss=0, Qloss=0;for( i=0; i<nB; i++ ){Ploss+=p_Pij[i]+p_Pji[i];Qloss+=p_Qij[i]+p_Qji[i];}5）程序流图如下6--------------Y Matrix--------------------0+j-17.3611 0+j0 0+j0 0+j17.3611 0+j0 0+j0 0+j0 0+j0 0+j00+j0 0+j-16 0+j0 0+j0 0+j0 0+j0 0+j16 0+j0 0+j00+j0 0+j0 0+j-17.0648 0+j0 0+j0 0+j0 0+j0 0+j0 0+j17.06480+j17.3611 0+j0 0+j0 3.30738+j-39.3089 -1.36519+j11.6041 -1.94219+j10.5107 0+j0 0+j0 0+j00+j0 0+j0 0+j0 -1.36519+j11.6041 2.55279+j-17.3382 0+j0 -1.1876+j5.97513 0+j0 0+j00+j0 0+j0 0+j0 -1.94219+j10.5107 0+j0 3.2242+j-15.8409 0+j0 0+j0 -1.28201+j5.588240+j0 0+j16 0+j0 0+j0 -1.1876+j5.97513 0+j0 2.80473+j-35.4456 -1.61712+j13.698 0+j00+j0 0+j0 0+j0 0+j0 0+j0 0+j0 -1.61712+j13.698 2.77221+j-23.3032 -1.15509+j9.784270+j0 0+j0 0+j17.0648 0+j0 0+j0 -1.28201+j5.58824 0+j0 -1.15509+j9.78427 2.4371+j-32.1539--------------Jacobian Matrix--------------------16.4 0 0 0 0 -16.4 0 0 0 0 0 0 0 00 17.4915 0 0 0 0 0 -17.4915 0 0 0 0 0 00 0 40.1703 -11.6041 -10.5107 0 0 0 3.30738 -1.36519 -1.94219 0 0 00 0 -11.6041 17.5792 0 -5.97513 0 0 -1.36519 2.55279 0 -1.1876 0 00 0 -10.5107 0 16.0989 0 0 -5.58824 -1.94219 0 3.2242 0 0 -1.28201-16.4 0 0 -5.97513 0 36.0731 -13.698 0 0 -1.1876 0 2.80473 -1.61712 00 0 0 0 0 -13.698 23.4822 -9.78427 0 0 0 -1.61712 2.77221 -1.155090 -17.4915 0 0 -5.58824 0 -9.78427 32.864 0 0 -1.28201 0 -1.15509 2.43710 0 -3.30738 1.36519 1.94219 0 0 0 38.4474 -11.6041 -10.5107 0 0 00 0 1.36519 -2.55279 0 1.1876 0 0 -11.6041 17.0972 0 -5.97513 0 00 0 1.94219 0 -3.2242 0 0 1.28201 -10.5107 0 15.5829 0 0 -5.588240 0 0 1.1876 0 -2.80473 1.61712 0 0 -5.97513 0 34.8181 -13.698 00 0 0 0 0 1.61712 -2.77221 1.15509 0 0 0 -13.698 23.1242 -9.784270 0 0 0 1.28201 0 1.15509 -2.4371 0 0 -5.58824 0 -9.78427 31.4437--------------Jacobian Matrix--------------------16.9269 0 0 0 0 -16.9269 0 0 0 0 0 1.68793 0 00 18.1691 0 0 0 0 0 -18.1691 0 0 0 0 0 0.8836270 0 41.9297 -12.1301 -11.1536 0 0 0 3.54272 -1.0628 -1.76646 0 0 00 0 -12.0455 18.0609 0 -6.01539 0 0 -1.78138 1.30819 0 -2.10262 0 00 0 -11.0484 0 16.8144 0 0 -5.76607 -2.33608 0 2.42598 0 0 -1.97778-16.9269 0 0 -6.36224 0 37.9476 -14.6585 0 0 -0.357534 0 3.05959 -0.930027 00 0 0 0 0 -14.4721 24.8873 -10.4152 0 0 0 -2.509 1.86088 -1.473890 -18.1691 0 0 -6.05157 0 -10.4721 34.6928 0 0 -0.733327 0 -0.9919732.66270 0 -3.52149 1.0628 1.76646 0 0 0 42.0299 -12.1301 -11.1536 0 0 00 0 1.78138 -3.884 0 2.10262 0 0 -12.0455 17.2037 0 -6.01539 0 00 0 2.33608 0 -4.31386 0 0 1.97778 -11.0484 0 16.2993 0 0 -5.766071.68793 0 0 0.357534 0 -2.97549 0.930027 0 0 -6.36224 0 38.3226 -14.6585 00 0 0 0 0 2.509 -3.98289 1.47389 0 0 0 -14.4721 24.2355 -10.41520 0.883627 0 0 0.733327 0 0.991973 -2.60893 0 0 -6.05157 0 -10.4721 34.8585 --------------Jacobian Matrix--------------------16.7457 0 0 0 0 -16.7457 0 0 0 0 0 1.63043 0 00 18.0388 0 0 0 0 0 -18.0388 0 0 0 0 0 0.8501960 0 41.3695 -11.8919 -10.9686 0 0 0 3.48069 -1.02775 -1.73712 0 0 00 0 -11.8057 17.6918 0 -5.8861 0 0 -1.7602 1.28091 0 -2.0217 0 00 0 -10.8651 0 16.5476 0 0 -5.68251 -2.29737 0 2.40655 0 0 -1.91027-16.7457 0 0 -6.21183 0 37.3041 -14.3465 0 0 -0.382862 0 2.95313 -0.937485 00 0 0 0 0 -14.1704 24.4052 -10.2348 0 0 0 -2.42909 1.86079 -1.433530 -18.0388 0 0 -5.94693 0 -10.2872 34.273 0 0 -0.757656 0 -0.9892082.598470 0 -3.48089 1.02775 1.73712 0 0 0 41.3703 -11.8919 -10.9686 0 0 00 0 1.7602 -3.78189 0 2.0217 0 0 -11.8057 16.6941 0 -5.8861 0 00 0 2.29737 0 -4.20764 0 0 1.91027 -10.8651 0 15.9488 0 0 -5.682511.63043 0 0 0.382862 0 -2.95077 0.937485 0 0 -6.21183 0 37.3083 -14.3465 00 0 0 0 0 2.42909 -3.86262 1.43353 0 0 0 -14.1704 23.7059 -10.23480 0.850196 0 0 0.757656 0 0.989208 -2.59706 0 0 -5.94693 0 -10.2872 34.2743 --------------Jacobian Matrix--------------------16.7435 0 0 0 0 -16.7435 0 0 0 0 0 1.63 0 00 18.0374 0 0 0 0 0 -18.0374 0 0 0 0 0 0.850 0 41.3625 -11.8888 -10.9664 0 0 0 3.48016 -1.02713 -1.73662 0 0 00 0 -11.8026 17.6871 0 -5.8845 0 0 -1.76008 1.28053 0 -2.02045 0 00 0 -10.8628 0 16.5444 0 0 -5.68158 -2.29703 0 2.40632 0 0 -1.90929-16.7435 0 0 -6.20987 0 37.296 -14.3426 0 0 -0.383399 0 2.95114 -0.937742 00 0 0 0 0 -14.1667 24.3994 -10.2326 0 0 0 -2.42794 1.86097 -1.433020 -18.0374 0 0 -5.94567 0 -10.285 34.2681 0 0 -0.758139 0 -0.9892032.597340 0 -3.48016 1.02713 1.73662 0 0 0 41.3625 -11.8888 -10.9664 0 0 00 0 1.76008 -3.78053 0 2.02045 0 0 -11.8026 16.6871 0 -5.8845 0 00 0 2.29703 0 -4.20632 0 0 1.90929 -10.8628 0 15.9444 0 0 -5.681581.63 0 0 0.383399 0 -2.95114 0.937742 0 0 -6.20987 0 37.296 -14.3426 00 0 0 0 0 2.42794 -3.86097 1.43302 0 0 0 -14.1667 23.6994 -10.23260 0.85 0 0 0.758139 0 0.989203 -2.59734 0 0 -5.94567 0 -10.285 34.2681--------------iteration---------------iteration = 4--------------voltage magnitude and angle--------------------1.04 0 01.025 0.161967 9.280011.025 0.0814153 4.664761.02579 -0.0386902 -2.216790.995631 -0.0696178 -3.988811.01265 -0.0643572 -3.68741.02577 0.064921 3.71971.01588 0.0126979 0.7275371.03235 0.0343257 1.96672-------------bus P and Q------------------1 0.71641 0.2704592 1.63 0.06653663 0.85 -0.1085974 0 05 -1.25 -0.56 -0.9 -0.37 0 08 -1 -0.359 0 0----------------line flow--------------------NUM------i-------j-----------begin------------end--------------1 4 1 -0.71641+j-0.239232 0.71641+j0.270462 7 2 -1.63+j0.0917816 1.63+j0.06653653 9 3 -0.85+j0.149553 0.85+j-0.1085974 7 8 0.763799+j-0.00797398 -0.759046+j-0.1070415 9 8 0.241834+j0.0311946 -0.240954+j-0.2429586 7 5 0.866201+j-0.0838079 -0.843202+j-0.1131287 9 6 0.608166+j-0.180748 -0.594627+j-0.1345668 5 4 -0.406798+j-0.386872 0.409373+j0.2289319 6 4 -0.305372+j-0.165433 0.307037+j0.0102993--------------------Ploss and Qloss-----------------Ploss = 0.0471901 Qloss = -0.9574833.思考题3.1潮流计算的方法有哪些？各有何特点？答：潮流计算分为手算和机算两大类，其中机算又有高斯-赛德尔迭代法、牛顿-拉夫逊迭代法、P-Q 分解法等算法。
潮流上机课程设计华电
潮流上机课程设计－华电————————————————————————————————作者：————————————————————————————————日期：ﻩ课程设计报告( 2０１1—２0１２年度第一学期)名称：电力系统潮流上机院系：电气与电子工程学院班级:学号:学生姓名:指导教师：设计周数:两周成绩：日期: 20１1年12月19日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力,掌握计算机潮流计算的相关知识二、设计正文（详细内容见附录)1．手算2．计算机计算3．思考题三、课程设计总结或结论潮流计算是研究电力系统稳态运行的一种基本计算,最初求解电力系统潮流时大多使用手算,但随着电力系统结构的日趋复杂，计算量也越来越大。
复杂电力系统潮流计算中，由于节点数量巨大,所形成的修正方程已经无法通过手算方式解决，尤其是需要迭代次数较多时，手算所需要的时间太长,计算机潮流计算无疑为解决这一问题提供了极大的便利。
计算机潮流计算可以迅速解决复杂网络的潮流计算问题，这是由于无论系统的复杂程度如何,其节点与支路的类型是固定的,所以只需要输入节点与支路的数据，就可以解决任何一个复杂网络的潮流计算问题。
即只需要一次编程，就可以基本上解决所有复杂网络的计算。
需要注意的是，在使用牛顿—拉弗逊发计算潮流时，对于初值要选择比较接近它们的精确解,否则迭代过程可能不收敛。
潮流计算C语言程序编程过程中需要注意的是,C语言无法实现复数运算，需要将得到的值的实部与虚部分开储存并计算。
这个过程复杂并且容易出错，编写程序是需要注意。
另外需要注意的一点是：C语言数组的编号是从零开始的，在程序编写过程中应注意下标的对应。
通过这一次的电力系统潮流计算编程，我不仅对Ｃ语言的编程有了更深刻的理解,也对《电力系统分析》这门课程进行了查漏补缺和巩固,对电力系统的运行也有了更加深入的了解，受益匪浅。
四、参考文献1.《电力系统计算：电子数字计算机的应用》,西安交通大学等合编。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《电力系统潮流上机》课程设计报告院系：电气与电子工程学院班级：学号：学生姓名：指导教师：刘宝柱设计周数：两周成绩：日期：2012年1月5日一、课程设计的目的与要求培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识 二、设计正文（详细内容见附录）1． 手算： 要求应用牛顿-拉夫逊法或P-Q 分解法手算求解，要求精度为。
节点1为平衡节点，电压︒∠=00.11U &，节点2为PQ 节点，负荷功率6.08.0~2j S +=，节点3是PV 节点，1.1,4.033==U P ，两条支路分别为04.001.013j Z +=，2.005.012j Z +=，对地支路33.030j y =。
2． 计算机计算 1．导纳阵Y( 1,1 )=( , Y( 1,2 )=( , Y( 1,3 )=( , Y( 1,4 )=( , Y( 2,1 )=( , Y( 2,2 )=( , Y( 2,4 )=( , Y( 3,1 )=( , Y( 3,3 )=( , Y( 4,1 )=( , Y( 4,2 )=( , Y( 4,4 )=( , 2．设定电压初值3．计算功率和电压偏移；同理可算出52596.0)0(22)0(2-=-=∆P P P s ,0196.0)0(22)0(2=-=∆Q Q Q s5.0)0(33)0(3=-=∆P P P s ,0.02)0(3232)0(3=-=∆U U U s4．根据求的第一次迭代时雅可比矩阵各元素的公式计算雅可比矩阵各个元素的具体值：5．求高斯计算后的修正量： 6．计算各节点电压的一次近似值： 返回第三步重新迭代，并校验收敛与否，令410-=ε。
经过 四 次迭代后，收敛条件满足，停止迭代，求出的电压终值：7．计算出平衡节点4的注入功率。
8．各节点间功率9．网损为:04487.001787.0~j S-= 10．网损效率为：% 3． 思考题潮流计算的方法有哪些各有何特点答：潮流计算方法主要包括：高斯-赛德尔迭代法、牛顿-拉夫逊迭代法、P-Q 分解法等算法。
各方法特点如下所示：如果交给你一个任务，请你用已有的潮流计算软件计算北京城市电网的潮流，你应该做哪些工作（收集哪些数据，如何整理，计算结果如何分析）答：1.所需要收集的数据：A.电网中所有节点的数据：a．各节点的类型,包括平衡节点、PV 节点、PQ 节点b. 对于平衡节点要了解节点的电压大小相位、及节点所能提供的最大最小有功无功功率c. PV节点要知道节点电压大小注入有功功率及节点所能提供的最大和最小无功功.率d. PQ节点要知道节点的注入有功和无功功率B．电网中所有支路的数据：a.各支路类型，即是否含有变压器b.各支路的电阻、电感、电纳c.各变压器的变比。
2．数据整理：将上述数据资料进行分类整理，并为每个节点及支路编上编号。
将整理的结果写成本实验中所要求的格式（原始数据的 txt 文档），再用本实验所编制的程序进行求解，得到各节点电压、相位，各线路传输功率、损耗，平衡节点注入功率等数值。
3.计算结果分析：可以大致了解整个电网的运行状况并按照电力系统知识依次进行所需的调整。
如调节各节点电压的大小以防止超出规定限度，调整各线路的有功无功分布状况分析整个系统的网损是否达到标准等。
设计中遇到的问题和解决的办法。
1.在设计最开始遇到的问题就是由于以前C语言所学习的知识遗忘较多，对于程序的使用不太熟练，通过老师讲解、看书、询问同学逐渐扫除了这个障碍。
2.设计的时候对于电力系统整体的知识脉络不是很清晰，编写程序时所定义的变量容易弄混淆，编写的程序产生了较多的错误，经过仔细学习和对每一排程序的排查，找出了产生错误的地方。
同时也对潮流计算的过程有了更清楚的认识。
三、课程设计总结或结论该门课程的主要目的是培养学生的电力系统潮流计算机编程能力，掌握计算机潮流计算的相关知识。
潮流计算是电力系统非常重要的分析计算，是电力系统最基本的计算，主要用以研究系统规划和运行中提出的各种问题。
现代电力系统潮流就算过程比较复杂，手算有很大的困难，而计算机作为一种处理数据的工具，具有计算速度快，准确率高的优点，因此现在电力系统的潮流大多也都采用计算机的算法求解。
对于写程序而言，C 语言是一种简单但很实用的语言，利用 C 语言编程实现潮流计算不仅简单易行，而且还能加深对电力系统潮流这一物理感念的理解。
通过这次潮流对 c++的面向对象编程及其过程有了更深入的了解。
对于潮流计算，我更深刻的掌握了牛拉法解潮流方程的方法以及对数学高次多元方程线性化在迭代求解的方法有了更深入的认识。
通过这次课程设计，自己动手编写潮流计算的程序，仔细分析了潮流计算各个步骤，理论与实践相结合,对潮流计算的思想和方法以及过程有了更深入的理解。
在这两周的学习里，我不仅对C语言的学习有了更深层次的理解，还对自己在《电力系统分析》这门课程里所学到的知识进行了查缺补漏，对电力系统的运行分析等有了更加深入的了解，受益匪浅。
四、参考文献1.《电力系统计算：电子数字计算机的应用》，西安交通大学等合编。
北京：水利电力出版社；2.《现代电力系统分析》，王锡凡主编，科学出版社；3.《电力系统稳态分析》，陈珩，中国电力出版社，1995年，第二版；附录（设计流程图、程序、表格、数据等）1.2.程序设计:#include""#include""#include""#include""#define NBUS 5#define NLINE 7/* Global variables */int nL,nSH,nBUS,nVA;double d,t;int L;struct Line{int Num,NumI,NumJ;double R;double X;double B;double K;};struct Bus{int Num ;double Volt,Phase,GenP,GenQ,LoadP,LoadQ; int Type;};struct Shunt{int Num,NumI;double G,B;};void main(){FILE *fp;FILE *fpout;int i,j,k,l,h,c;int i1,i2,i3,kp,kq,LH[100];double d,d1,d2,d3,d4,d5,d6,t,r,x,g,b,tt,LL,e,ps,qs,shsh,sum,w;struct Line sL[NLINE];struct Bus sBus[NBUS];struct Shunt sSH[NBUS];double YG[NBUS][NBUS],YB[NBUS][NBUS];doubleV[100][2],U[100],dPQ[100],PQ[100][2],JJ[100][100]={0},H[100][100],J[100][1 00],N[100][100],L[100][100],sP[NBUS][NBUS]={0},sQ[NBUS][NBUS]={0},dsp,dsq, sumgen,dp;i1=i2=i3=0;d1=d2=d3=d4=d5=d6=ps=qs=;e=;for(i=0;i<NBUS;i++)/*Read the input data*/if((fp=fopen("","r"))==NULL){ printf("Can not open the file named '' \n");exit(0);}fscanf(fp,"%d,%d,%d",&nBUS,&nL,&nSH);for(i=0;i<nBUS;i++){sBus[i].Num=sBus[i].Type=0;sBus[i].Volt=;sBus[i].Phase=sBus[i].GenP=sBus[i].GenQ=sBus[i].LoadP=sBus[i].LoadQ=;fscanf(fp,"%d,%lf,%lf,%lf,%lf,%lf,%lf,%d",&i1,&d1,&d2,&d3,&d4,&d5,&d6,&i2) ;sBus[i].Num=i1;sBus[i].Volt=d1;sBus[i].Phase=d2;sBus[i].GenP=d3;sBus[i].Ge nQ=d4;sBus[i].LoadP=d5,sBus[i].LoadQ=d6;sBus[i].Type=i2;};for(i=0;i<nL;i++){sL[i].Num=sL[i].NumI=sL[i].NumJ=0;sL[i].R=sL[i].X=sL[i].B=;sL[i].K=;fscanf(fp,"%2d %3d %3d %lf %lf %lf %lf",&i1,&i2,&i3,&d1,&d2,&d3,&d4);sL[i].Num=i1;sL[i].NumI=i2;sL[i].NumJ=i3;sL[i].R=d1;sL[i].X=d2;sL[i].B=d3; sL[i].K=d4;}for(i=0;i<nSH;i++){sSH[i].Num=sSH[i].NumI=0;sSH[i].G=sSH[i].B=; fscanf(fp,"%2d %3d %lf",&i1,&i2,&d1);sSH[i].Num=i1;sSH[i].NumI=i2;sSH[i].B=d1; }if(fp!=NULL) fclose(fp);/*Make Y Matrix*/for(i=1;i<nBUS+1;i++)for(j=1;j<nBUS+1;j++){ YG[i][j]=;YB[i][j]=;};for(l=0; l<nL; l++){i=sL[l].NumI;j=sL[l].NumJ;r=sL[l].R;x=sL[l].X;d1=r*r+x*x;g=r/d1;b=-x/d1;if(fabs(sL[l].<{/*Normal lines or transformers*/YG[i][i]=YG[i][i]+g;YG[j][j]=YG[j][j]+g;YB[i][i]=YB[i][i]+b+sL[l].B;YB[j][j]=YB[j][j]+b+sL[l].B;YG[i][j]=YG[i][j]-g;YG[j][i]=YG[j][i]-g;YB[i][j]=YB[i][j]-b;YB[j][i]=YB[j][i]-b;}else{/*abnormal transformer ratio*/if(fabs(sL[l].B)>{YG[i][i]=YG[i][i]+g/sL[l].B/sL[l].B;YB[i][i]=YB[i][i]+b/sL[l].B/sL[l].B;YG[j][j]=YG[j][j]+g;YB[j][j]=YB[j][j]+b;YG[i][j]=YG[i][j]-g/sL[l].B;YG[j][i]=YG[j][i]-g/sL[l].B;YB[i][j]=YB[i][j]-b/sL[l].B;YB[j][i]=YB[j][i]-b/sL[l].B;}elseprintf("节点%2d,%-2d之间变压器变比为0 \n",i,j); }}ype==0){V[i][0]=;V[i][1]=;}for(i=1;i<nBUS+1;i++)if(sBus[i-1].Type==1){V[i][1]=sBus[i-1].Volt;V[i][0]=;}for(i=1;i<nBUS+1;i++)if(sBus[i-1].Type==2){V[i][1]=sBus[i-1].Volt;V[i][0]=sBus[i-1].Phase;}/* 输出电压初值if((fp=fopen("d:\\lx\\电压初值.txt","w"))==NULL){printf("Can not open the file named '电压初值.txt' \n");exit(0);}fprintf(fp,"---电压初值---\n");for(i=1;i<nBUS+1;i++)for(j=1;j<2;j++)fprintf(fp,"Y(%2d )=(%,%\n",i,V[i][0],V[i][1]);if(fp!=NULL) fclose(fp);*/for(c=1;;c++){/* 计算偏移量 */for(i=1;i<nBUS+1;i++){if(sBus[i-1].Type!=2){PQ[i][0]=0;PQ[i][1]=0;for(j=1;j<nBUS+1;j++){PQ[i][0]-=V[i][1]*V[j][1]*(YG[i][j]*cos(V[i][0]-V[j][0])+YB[i][j]*sin( V[i][0]-V[j][0]));PQ[i][1]-=V[i][1]*V[j][1]*(YG[i][j]*sin(V[i][0]-V[j][0])-YB[i][j]*cos( V[i][0]-V[j][0]));}}}h=1;for(i=1;i<nBUS+1;i++){if(sBus[i-1].Type==0){dPQ[h]=PQ[i][0]+sBus[i-1].GenP-sBus[i-1].LoadP;h++;dPQ[h]=PQ[i][1]+sBus[i-1].GenQ-sBus[i-1].LoadQ;h++;}if(sBus[i-1].Type==1){dPQ[h]=PQ[i][0]+sBus[i-1].GenP-sBus[i-1].LoadP;h++;}}/* 输出偏移量if((fp=fopen("d:\\lx\\偏移量 .txt","w"))==NULL){printf("Can not open the file named '偏移量.txt' \n"); exit(0);}fprintf(fp,"--- 偏移量 ---\n");for(i=1;i<h;i++)fprintf(fp,"dPQ(%2d )=(%\n",i,dPQ[i]);if(fp!=NULL) fclose(fp);*//* 计算雅克比矩阵 */for(i=1;i<nBUS+1;i++){for(j=1;j<nBUS+1;j++){if(i==j){H[i][j]=PQ[i][1]-V[i][1]*V[i][1]*YB[i][j];J[i][j]=-PQ[i][0]-V[i][1]*V[i][1]*YG[i][j];N[i][j]=-PQ[i][0]+V[i][1]*V[i][1]*YG[i][j];L[i][j]=-PQ[i][1]-V[i][1]*V[i][1]*YB[i][j];}else{H[i][j]=V[i][1]*V[j][1]*(YG[i][j]*sin(V[i][0]-V[j][0])-YB[i][j]*cos(V[ i][0]-V[j][0]));J[i][j]=-V[i][1]*V[j][1]*(YG[i][j]*cos(V[i][0]-V[j][0])+YB[i][j]*sin(V [i][0]-V[j][0]));N[i][j]=V[i][1]*V[j][1]*(YG[i][j]*cos(V[i][0]-V[j][0])+YB[i][j]*sin(V[ i][0]-V[j][0]));L[i][j]=V[i][1]*V[j][1]*(YG[i][j]*sin(V[i][0]-V[j][0])-YB[i][j]*cos(V[ i][0]-V[j][0]));}}}h=1;for(i=1;i<nBUS+1;i++){l=1;if(sBus[i-1].Type==0){for(j=1;j<nBUS+1;j++){if(sBus[j-1].Type==0){JJ[h][l]=H[i][j];l++;JJ[h][l]=N[i][j];l++;}if(sBus[j-1].Type==1){JJ[h][l]=H[i][j];l++;}}h++;l=1;for(j=1;j<nBUS+1;j++) {if(sBus[j-1].Type==0){JJ[h][l]=J[i][j];l++;JJ[h][l]=L[i][j];l++;}if(sBus[j-1].Type==1){JJ[h][l]=J[i][j];l++;}}h++;}if(sBus[i-1].Type==1){for(j=1;j<nBUS+1;j++){if(sBus[j-1].Type==0){JJ[h][l]=H[i][j];l++;JJ[h][l]=N[i][j];l++;}if(sBus[j-1].Type==1){JJ[h][l]=H[i][j];l++;}}h++;}}/* 输出雅克比矩阵if((fp=fopen("d:\\lx\\雅克比矩阵 .txt","w"))==NULL){printf("Can not open the file named '雅克比矩阵.txt' \n"); exit(0);}fprintf(fp,"--- 雅克比矩阵 ---\n");fprintf(fp," ");for(i=1;i<h;i++){fprintf(fp," %d\t\t",i);}fprintf(fp,"\n");*/for(i=1;i<h;i++){fprintf(fp,"%2d ",i);for(j=1;j<h;j++){fprintf(fp,"%\t",JJ[i][j]);}fprintf(fp,"\n");}if(fp!=NULL) fclose(fp);/*高斯法求解方程组*/l=1;for(i=1;i<h;i++){LH[i]=0;}for(k=1;k<h;k++){d=;for(j=k;j<h;j++){if(fabs(JJ[k][j])>d){d=fabs(JJ[k][j]); /*在一行中找到一个最大值赋值d,并用JS[K]记住这个最大值所在的列号*/LH[k]=j;}}if(fabs(d)< /*如果d的数值太小，做为被除数将带来很大的误差*/{l=0;}else if(LH[k]!=k){for(i=1;i<h;i++){t=JJ[i][k];JJ[i][k]=JJ[i][LH[k]]; /*进行列交换，让最大值始终在对角元上*/JJ[i][LH[k]]=t;}}if(l==0){break;}for(j=k+1;j<h;j++){JJ[k][j]=JJ[k][j]/JJ[k][k]; /*对角元上的元素消为1*/ }dPQ[k]=dPQ[k]/JJ[k][k];for(i=k+1;i<h;i++){for(j=k+1;j<h;j++){JJ[i][j]=JJ[i][j]-JJ[i][k]*JJ[k][j]; /*使下三角阵的元素为0*/}dPQ[i]=dPQ[i]-JJ[i][k]*dPQ[k];}}if(fabs(JJ[h-1][h-1])>{ /*用追赶法，解方程组，求未知数x*/ U[h-1]=dPQ[h-1];for(i=h-2;i>=0;i--){t=;for(j=i+1;j<h;j++){t=t+JJ[i][j]*U[j];}U[i]=(dPQ[i]-t);}}/*输出高斯结果if((fp=fopen("d:\\lx\\高斯.txt","w"))==NULL) {printf("err");exit(0);}for(i=1;i<h;i++){fprintf(fp,"%f",U[i]);fprintf(fp,"\n");}fclose(fp);*//*得到电压值*/h=1;for(i=1;i<nBUS+1;i++){if(sBus[i-1].Type==0){V[i][0]+=U[h];h++;V[i][1]+=U[h];h++;}if(sBus[i-1].Type==1){V[i][0]+=U[h];h++;}}/* 输出电压值if((fp=fopen("d:\\lx\\电压值.txt","w"))==NULL){printf("Can not open the file named '电压值.txt' \n"); exit(0);}fprintf(fp,"---电压初值---\n");for(i=1;i<nBUS+1;i++)for(j=1;j<2;j++)fprintf(fp,"Y(%2d )=(%,%\n",i,V[i][0],V[i][1]);if(fp!=NULL) fclose(fp);*//*求最大变化值*/w=0;for(i=1;i<h;i++){if(U[i]>w){w=U[i];}}if(w<{break;}}/* 输出电压终值if((fp=fopen("d:\\lx\\电压终值.txt","w"))==NULL){printf("Can not open the file named '电压终值.txt' \n"); exit(0);}fprintf(fp,"---电压终值---\n");fprintf(fp,"循环%d次\n",c);for(i=1;i<nBUS+1;i++)for(j=1;j<2;j++)fprintf(fp,"Y(%2d )=(%,%\n",i,V[i][0]*180/,V[i][1]);if(fp!=NULL) fclose(fp);*/ps=0;qs=0;for(i=1;i<nBUS+1;i++)if(sBus[i-1].Type==2){for(j=1;j<nBUS+1;j++){ps+=V[i][1]*V[j][1]*(YG[i][j]*cos(V[i][0]-V[j][0])+YB[i][j]*sin(V[i][0 ]-V[j][0]));qs+=V[i][1]*V[j][1]*(YG[i][j]*sin(V[i][0]-V[j][0])-YB[i][j]*cos(V[i][0 ]-V[j][0]));}}}for(i=1;i<nBUS+1;i++){for(j=1;j<nBUS+1;j++){sP[i][j]=0;sQ[i][j]=0;}}for(l=0; l<nL; l++)i=sL[l].NumI;j=sL[l].NumJ;r=sL[l].R;x=sL[l].X;d1=r*r+x*x;g=r/d1;b=-x/d1;if(fabs(sL[l].<{/*Normal lines or transformers*/sP[i][j]=V[i][1]*V[i][1]*g-V[i][1]*V[j][1]*(g*cos(V[i][0]-V[j][0])+b*s in(V[i][0]-V[j][0]));sQ[i][j]=-(V[i][1]*V[i][1]*sL[l].B+V[i][1]*V[i][1]*b+V[i][1]*V[j][1]*( g*sin(V[i][0]-V[j][0])-b*cos(V[i][0]-V[j][0])));sP[j][i]=V[j][1]*V[j][1]*g-V[i][1]*V[j][1]*(g*cos(V[j][0]-V[i][0])+b*s in(V[j][0]-V[i][0]));sQ[j][i]=-(V[j][1]*V[j][1]*sL[l].B+V[j][1]*V[j][1]*b+V[i][1]*V[j][1]*( g*sin(V[j][0]-V[i][0])-b*cos(V[j][0]-V[i][0])));}else{/*abnormal transformer ratio*/sP[i][j]=V[i][1]*V[i][1]*g/sL[l].B/sL[l].B-V[i][1]*V[j][1]*(g*cos(V[i] [0]-V[j][0])/sL[l].B+b*sin(V[i][0]-V[j][0])/sL[l].B);sQ[i][j]=-(V[i][1]*V[i][1]*b/sL[l].B/sL[l].B+V[i][1]*V[j][1]*(g*sin(V[ i][0]-V[j][0])/sL[l].B-b*cos(V[i][0]-V[j][0])/sL[l].B));sP[j][i]=V[j][1]*V[j][1]*g-V[i][1]*V[j][1]*(g*cos(V[j][0]-V[i][0])/sL[ l].B+b*sin(V[j][0]-V[i][0])/sL[l].B);sQ[j][i]=-(V[i][1]*V[i][1]*b+V[i][1]*V[j][1]*(g*sin(V[j][0]-V[i][0])/s L[l].B-b*cos(V[j][0]-V[i][0])/sL[l].B));}}dsp=ps;dsq=qs;sumgen=ps;for(i=0;i<nBUS;i++){dsp+=sBus[i].GenP-sBus[i].LoadP;dsq+=sBus[i].GenQ-sBus[i].LoadQ;sumgen+=sBus[i].GenP;}dp=dsp/sumgen*100;/* 输出功率情况 */if((fp=fopen("功率情况.txt","w"))==NULL){printf("Can not open the file named '功率情况.txt' \n");exit(0);}fprintf(fp,"---功率情况---\n");fprintf(fp,"平衡节点功率S=%+ j%\n",ps,qs);for(i=1;i<nBUS+1;i++)for(j=1;j<nBUS+1;j++)if(fabs(sP[i][j]>fprintf(fp,"S(%3d,%-3d)=(%,j%\n",i,j,sP[i][j],sQ[i][j]); fprintf(fp,"网损为%+j%,输电效率为%\n",dsp,dsq,100-dp);if(fp!=NULL) fclose(fp);/* 输出最终结果 */if((fp=fopen("最终结果.txt","w"))==NULL){printf("Can not open the file named '最终结果.txt' \n");exit(0);}fprintf(fp,"---最终结果---\n");fprintf(fp,"\n\n\n");fprintf(fp,"一共循环迭代计算%d次\n",c);fprintf(fp,"\n\n\n");fprintf(fp,"---电压终值---\n");for(i=1;i<nBUS+1;i++)for(j=1;j<2;j++)fprintf(fp,"U(%2d )=(%,%\n",i,V[i][0]*180/,V[i][1]);fprintf(fp,"\n\n\n");fprintf(fp,"平衡节点功率S=%+j%\n",ps,qs);fprintf(fp,"\n\n\n");fprintf(fp,"---各节点间功率---\n");for(i=1;i<nBUS+1;i++)for(j=1;j<nBUS+1;j++)if(fabs(sP[i][j]>fprintf(fp,"S(%2d,%-2d)=(%,j%\n",i,j,sP[i][j],sQ[i][j]); fprintf(fp,"\n\n\n");fprintf(fp,"网损为：%+j%。