中国石油大学电力系统分析大作业

1、 电力系统分析大作业 编程实现N-R迭代和P-Q迭代潮流计算 院 系： 信息与控制工程学院 班 级： 电气12-04班 学 号： 120534 姓 名： 日期：2015年 06 月 17 日要求：利用NR迭代和P-Q迭代计算所给系统的潮流，设发电机G1的端电压为1 p.u.，发出的有功、无功可调；发电机G2的端电压为1 p.u.，按指定的有功P=0.5 p.u.发电，取=10xe-4。1:牛顿-拉夫逊法潮流算法1.1牛顿-拉夫逊法潮流算法的基本原理 电力系统潮流计算是电力系统运行和规划中最基本和最经常的计算，其任务是在已知某些运行参数的情况下，计算出系统中全部的运行参数，一般来说，各个母线所供

2、负荷的功率是已知的，各个节点电压是未知的（平衡节点除外），可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵和网络拓扑结构列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初

3、值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。1.2节点导纳矩阵的算法 根据定义直接求取节点导纳矩阵时，注意以下几点： 1) 节点导纳矩阵是方阵，其阶数就等于网络中除去参考节点外的节点数。参考节点一般取大地，编号为零。 2) 节点导纳矩阵是稀疏矩阵，其各行非零非对角元素就等于与该行相对应节点所连接的不接地支路数。 3) 节点导纳矩阵的对角元素就等于该节点所连接导纳的总和。因此，与没有接地支路的节点对应的行或列中，对角元素为非对角元素之和的负值。 4) 节点导纳矩阵的非对角元素等于连接节点i，j支路导纳的负值。因此，在

4、一般情况下，节点导纳矩阵的对角元素往往大于非对角元素的负值。 5) 节点导纳矩阵一般是对称矩阵，这是网络的互易特性所决定的。从而，一般只要求取这个矩阵的上三角或下三角部分。 6) 网络中的变压器。1.3 牛顿拉夫逊迭代法的步骤及算法流程图：1）形成各节点导纳矩阵Y。2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。3）计算各个节点的功率不平衡量。4）根据收敛条件判断是否满足，若不满足则向下进行。5）计算雅可比矩阵中的各元素。6）修正方程式个节点电压7）利用新值自第3）步开始进入下一次迭代，直至达到精度退出循环。8）计算平衡节点输出功率和各线路功率流程图如下：形成节点导纳矩阵设非平

5、衡节点电压初值、e输入原始数据启动对PQ节点计算(对PV节点计算) 计算雅可比矩阵各元素Hij(k)、Nij(k)、Jij(k)、Lij(k)解修正方程，由及雅可比矩阵用牛顿-拉夫逊法求各节点的ei（k）、Ui(k)/Ui计算节点的新电压增加迭代次数k=k+1令迭代次数k=0计算平衡节点的功率及线路功率输出是否1.4程序代码及详细注释：clear;disp('【输出结果如下】：')m=5;m1=5;isb=5;pr=0.0001;A1=1 2 0.025+0.08i 0.07i 1 0; 1 3 0.03+0.1i 0.09i 1 0; 2 3 0.02+0.06i 0.05i

6、 1 0; 4 2 0.1905i 0 1.0522 1; 5 3 0.1905i 0 1.0522 1;%A1是支路参数矩阵，第一列和第二列是节点编号.%第三列为之路串联阻抗参数.%第四列为含支路的对地导纳矩阵.%第五列为含变压器之路的变压器变比.%第六列表式是否含有变压器的参数，1表示含有变压器，0表示不含有.%变压器侧为1，否则为0.A2=0 -0.8055-0.5320i 1 0 0 2; 0 -0.18-0.12i 1 0 0 2; 0 0 1 0 0 2; 0 0.5 1 1.0522 0 3; 0 0 1 1.0522 0 1;% A2为节点参数矩阵.%第二列为节点负荷功率参数.

7、%第三列为节点电压参数.%第六列为节点类型参数其中1为平衡节点，2为PQ节点，3为PV节点.disp('-【 已知量 ：】-')disp('【支路参数矩阵】')disp(A1)disp('【节点参数矩阵】')disp(A2)%=计算导纳矩阵=X=1 0;2 0;3 0;4 0;5 0;Y=zeros(m);U=zeros(1,m);cta=zeros(1,m);V=zeros(1,m);O=zeros(1,m);S1=zeros(m1);for i=1:m if X(i,2)=0; p=X(i,1); Y(p,p)=X(i,2); endendf

8、or i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end Y(p,q)=Y(p,q)-1./(A1(i,3)*A1(i,5); Y(q,p)=Y(p,q); Y(q,q)=Y(q,q)+1./A1(i,3)+A1(i,4)./2; Y(p,p)=Y(p,p)+1./(A1(i,3)* A1(i,5)2) +A1(i,4)./2;enddisp('【导纳矩阵】 Y=');disp(Y)%输出导纳矩阵%=给定各节点初始电压及给定各节点注入功率=G=real(Y);B=imag(Y);%分解出导

9、纳阵的实部和虚部for i=1:m cta(i)=angle(A2(i,3); U(i)=abs(A2(i,3); V(i)=A2(i,4);%PV、平衡节点及PQ节点电压模值endfor i=1:m %给定各节点注入功率 S(i)=A2(i,1)-A2(i,2);%i节点注入功率SG-SL B(i,i)=B(i,i)+A2(i,5);%i节点无功补偿量(电纳值)end%=用牛顿-拉夫逊法迭代求解非线性代数方程（功率方程）= P=real(S);Q=imag(S);%分解出各节点注入的有功和无功功率ICT1=0;IT2=1;%迭代次数ICT1；不满足收敛要求的节点数IT2while IT2=0

10、 IT2=0;t1=1;t2=1; for i=1:m if i=isb C(i)=0; D(i)=0; for j1=1:m %第i行共m列(m个节点间互导纳及节点电压相乘即电流) C(i)=C(i)+U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1); D(i)=D(i)+U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1); end DP(t1)=P(i)-C(i); t1=t1+1; if A2(i,6)=2 DQ(t2)=Q(i

11、)-D(i); t2=t2+1; %计算失配功率 end endenddisp('失配功率: ') %输出失配功率 disp(DP); disp(DQ); t1=t1-1;t2=t2-1; DPQ=DP'DQ' %求DP,DQ %= 以下为形成完整的雅克比矩阵 = for i=1:t1+t2 if abs(DPQ(i)>pr IT2=IT2+1; end end H=zeros(t1,t1);N=zeros(t1,t2);K=zeros(t2,t1);L=zeros(t2,t2); for i=1:t1 for j1=1:t1 if j1=isb&

12、j1=i H(i,j1)=0-U(i)*U(j1)*(G(i,j1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1); elseif j1=isb&j1=i H(i,j1)=U(i)2*B(i,j1)+D(i); end end end %计算H矩阵 for i=1:t1 for j1=1:t2 if j1=isb&j1=i N(i,j1)=0-U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1); elseif j1=isb&j1=i N(i,j1

13、)=0-U(i)2*G(i,j1)-C(i); end end end %计算N矩阵 for i=1:t2 for j1=1:t1 if j1=isb&j1=i K(i,j1)= U(i)*U(j1)*(G(i,j1)*cos(cta(i)-cta(j1)+B(i,j1)*sin(cta(i)-cta(j1); elseif j1=isb&j1=i K(i,j1)=U(i)2*G(i,j1)-C(i); end end end %计算K矩阵 for i=1:t2 for j1=1:t2 if j1=isb&j1=i L(i,j1)=0-U(i)*U(j1)*(G(i,j

14、1)*sin(cta(i)-cta(j1)-B(i,j1)*cos(cta(i)-cta(j1); elseif j1=isb&j1=i L(i,j1)=U(i)2*B(i,j1)-D(i); end end end %计算L矩阵 J=H,N;K,L;%求雅可比矩阵 disp('【 雅可比矩阵 ：】'); disp(J); modify=-JDPQ; Dcta=modify(1:t1,:); t3=U(:,1:t2); DU=diag(t3,0)*modify(t1+1:t1+t2,:); t4=1; for i=1:t1 if A2(i,6)=1 cta(1,i)=c

15、ta(1,i)+Dcta(t4,1); t4=t4+1; end end t5=1; for i=1:t2 if A2(i,6)=2 U(1,i)=U(1,i)+DU(t5,1); t5=t5+1; end end ICT1=ICT1+1;end %修正原值for i=1:m UU(i)=U(i)*cos(cta(i)+1i*U(i)*sin(cta(i);endfor p=1:m c(p)=0; for q=1:m c(p)=c(p)+conj(Y(p,q)*conj(UU(q); end s(p)=UU(p)*c(p);enddisp('*');disp('【 各节

16、点电压U为（节点从小到大排列）:】');disp(UU);disp('*');disp('【 各节点电压相角为（节点从小到大排列） ：】');disp(180*angle(UU)/pi);disp('*');disp('【 按公式计算全部线路功率，结果如下:】');for i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end Si(p,q)=UU(p)*(conj(UU(p)*conj(A1(i,4)./2)+(conj(UU(p)*A

17、1(i,5)-conj(UU(q)*conj(1./(A1(i,3)*A1(i,5);%各条支路首端功率Si f=p,q,Si(p,q); disp(f);endfor i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end Sj(q,p)=UU(q)*(conj(UU(q)*conj(A1(i,4)./2)+(conj(UU(q)./A1(i,5)-conj(UU(p)*conj(1./(A1(i,3)*A1(i,5);%各条支路末端功率Sj f=q,p,Sj(q,p); disp(f);enddisp(&

18、#39;*');disp('【 各条支路的功率损耗DS为(顺序同输入A1时一样):】');for i=1:m1 if A1(i,6)=0 p=A1(i,1);q=A1(i,2); else p=A1(i,2);q=A1(i,1); end DS(i)=Si(p,q)+Sj(q,p);%各条支路功率损耗DS disp(DS(i);endSp=0;for i=1:m Sp=Sp+UU(isb)*conj(Y(isb,i)*conj(UU(i);enddisp('【 平衡节点的功率：】');disp(Sp); 1.6输出结果：>> dlxt【输出结

19、果如下】：-【 已知量 ：】-【支路参数矩阵】1.0000 2.0000 0.0250 + 0.0800i 0 + 0.0700i 1.0000 0 1.0000 3.0000 0.0300 + 0.1000i 0 + 0.0900i 1.0000 0 2.0000 3.0000 0.0200 + 0.0600i 0 + 0.0500i 1.0000 0 4.0000 2.0000 0 + 0.1905i 0 1.0522 1.0000 5.0000 3.0000 0 + 0.1905i 0 1.0522 1.0000 【节点参数矩阵】0 -0.8055 - 0.5320i 1.0000 0

20、0 2.0000 0 -0.1800 - 0.1200i 1.0000 0 0 2.0000 0 0 1.0000 0 0 2.0000 0 0.5000 1.0000 1.0522 0 3.0000 0 0 1.0000 1.0522 0 1.0000 【导纳矩阵】 Y=6.3110 -20.4822i -3.5587 +11.3879i -2.7523 + 9.1743i 0 0 -3.5587 +11.3879i 8.5587 -31.0693i -5.0000 +15.0000i 0 + 4.9889i 0 -2.7523 + 9.1743i -5.0000 +15.0000i 7.7

21、523 -28.8457i 0 0 + 4.9889i 0 0 + 4.9889i 0 0 - 5.2493i 0 0 0 0 + 4.9889i 0 0 - 5.2493i失配功率: 0.8055 0.1800 0 -0.5000 0.6120 0.4275 0.3175【 雅可比矩阵 ：】 -20.5622 11.3879 9.1743 0 -6.3110 3.5587 2.7523 11.3879 -31.3768 15.0000 4.9889 3.5587 -8.5587 5.0000 9.1743 15.0000 -29.1632 0 2.7523 5.0000 -7.7523 0

22、4.9889 0 -4.9889 0 0 0 6.3110 -3.5587 -2.7523 0 -20.4022 11.3879 9.1743 -3.5587 8.5587 -5.0000 0 11.3879 -30.7618 15.0000 -2.7523 -5.0000 7.7523 0 9.1743 15.0000 -28.5282失配功率: -0.2695 0.0231 0.0748 0.0799 -0.1030 -0.1277 -0.1137【 雅可比矩阵 ：】 -28.8103 15.9814 12.8289 0 -10.1478 4.5599 3.4378 15.7341 -41

23、.6855 20.1846 5.7668 5.3512 -11.7083 6.6232 12.6028 20.1216 -38.4668 0 4.1917 6.8122 -10.2937 0 5.7668 0 -5.7668 0 0.5799 0 7.9977 -4.5599 -3.4378 0 -30.0803 15.9814 12.8289 -5.3512 11.3945 -6.6232 0.5799 15.7341 -42.1809 20.1846 -4.1917 -6.8122 10.4433 0 12.6028 20.1216 -38.6942失配功率: -0.0107 0.0024

24、 0.0044 0.0015 -0.0018 -0.0041 -0.0042【 雅可比矩阵 ：】 -27.3240 15.1545 12.1695 0 -9.3998 4.4171 3.3502 14.9730 -39.8614 19.2510 5.6374 4.9977 -11.1924 6.3412 12.0040 19.2055 -36.8266 0 3.9019 6.4777 -9.8938 0 5.6374 0 -5.6374 0 0.5015 0 7.7673 -4.4171 -3.3502 0 -28.3915 15.1545 12.1695 -4.9977 10.8373 -6

25、.3412 0.5015 14.9730 -40.1097 19.2510 -3.9019 -6.4777 9.9027 0 12.0040 19.2055 -36.8350失配功率: 1.0e-004 * -0.1524 0.0448 0.0746 0.0085 1.0e-005 * -0.0652 -0.4685 -0.5571【 雅可比矩阵 ：】 -27.2778 15.1288 12.1490 0 -9.3743 4.4143 3.3490 14.9504 -39.8057 19.2220 5.6333 4.9855 -11.1784 6.3329 11.9862 19.1774 -3

26、6.7767 0 3.8915 6.4669 -9.8837 0 5.6333 0 -5.6333 0 0.5000 0 7.7633 -4.4143 -3.3490 0 -28.3418 15.1288 12.1490 -4.9855 10.8184 -6.3329 0.5000 14.9504 -40.0457 19.2220 -3.8915 -6.4669 9.8838 0 11.9862 19.1774 -36.7767*【 各节点电压U为（节点从小到大排列）:】 1.1586 + 0.1243i 1.1292 + 0.0995i 1.1251 + 0.0951i 1.0000 - 0

27、.0007i 1.0000 *【 各节点电压相角为（节点从小到大排列） ：】 6.1212 5.0334 4.8335 -0.0388 0*【 按公式计算全部线路功率，结果如下:】 1.0000 2.0000 0.4176 + 0.2856i 1.0000 3.0000 0.3879 + 0.2464i 2.0000 3.0000 0.0923 + 0.0216i 2.0000 4.0000 0.5000 + 1.1124i 3.0000 5.0000 0.4747 + 1.0795i 2.0000 1.0000 -0.4123 - 0.3613i 3.0000 1.0000 -0.3825

28、- 0.3468i 3.0000 2.0000 -0.0921 - 0.0851i 4.0000 2.0000 -0.5000 - 0.8919i 5.0000 3.0000 -0.4747 - 0.8717i*【 各条支路的功率损耗DS为(顺序同输入A1时一样):】 0.0053 - 0.0757i 0.0054 - 0.1004i 0.0002 - 0.0635i 0 + 0.2205i -0.0000 + 0.2078i【 平衡节点的功率：】 -0.4747 - 0.3638i>> 7、结果截图：2:P-Q解耦迭代法2.1 程序代码cleardisp('【电力系统混合

29、坐标下的PQ解耦迭代潮流计算】:');n=5;m=3;%= 参数初始化= Y=zeros(n,n); %导纳矩阵U=ones(n,1); %电压矢量a=zeros(n,1); %相角矢量Ps=zeros(n,1);Qs=zeros(n,1);P=zeros(n,1);Q=zeros(n,1);p=zeros(n-1,1);q=zeros(m,1);aa=zeros(n-1,1);u=zeros(m,1);k=0;%=节点导纳矩阵的生成=I=sqrt(-1);%由于后面用到ij作为循环变量，故重定义虚数单位Y=6.311-I*20.3722 -2.7523+I*9.1743 -3.558

30、7+I*11.3879 0 0; -3.5587+I*11.3879 8.5587-I*30.9993 -5+I*15 I*4.9889 0; -2.7523+I*9.1743 -5+I*15 7.7523-I*28.7557 0 I*4.9889; 0 I*4.9889 0 -I*5.2493 0; 0 0 I*4.9889 0 -I*5.2493;%=初值的设定=U=1;1;1;1;1;a=0;0;0;0;0;Ps=-0.8055;-0.18;0;0.5;0;Qs=-0.532;-0.12;0;0;0;%求失配功率P(k)和Q(k)for i=1:n-1 s=0; for j=1:n s=

31、s+U(j)*(real(Y(i,j)*cos(a(i)-a(j)+imag(Y(i,j)*sin(a(i)-a(j); end P(i)=U(i)*s;endfor i=1:n-1 s=0; for j=1:n s=s+U(j)*(real(Y(i,j)*sin(a(i)-a(j)-imag(Y(i,j)*cos(a(i)-a(j); end Q(i)=U(i)*s;endfor i=1:n-1 detp(i)=Ps(i)-P(i);endfor i=1:m detq(i)=Qs(i)-Q(i);end %生成B'和B"矩阵 Bp=zeros(n-1,n-1); Bpp=z

32、eros(m,m); for i=1:n-1; for j=1:n-1; Bp(i,j)=imag(Y(i,j); end; end; for i=1:m; for j=1:m; Bpp(i,j)=imag(Y(i,j); end; end; k=0;while (max(detp)>=0.0001 | max(detq)>=0.0001) %设定精确度 aa=(-1)*inv(Bp)*detp(1)/U(1);detp(2)/U(2);detp(3)/U(3);detp(4)/U(4); u=(-1)*inv(Bpp)*detq(1)/U(1);detq(2)/U(2);detq(3)/U(3); %求解修正方程，得到修正量 aa,u for i=1:n-1 %计算 a(电压相角),U(电压幅值) a(i)=a(i)+aa(i); end for i=1:m U(i)=U(i)+u(i); end% detp'detq'迭代过程中失配功率变化情况 k=k+1;%(记录迭代次数) for i=1:n-1 %计算功率误