电力系统分析潮流计算的计算机算法

1、潮流计算的计算机算法实验报告:学号：班级：一、实验目的掌握潮流计算的计算机算法。熟悉MATLAB,并掌握MATLAB程序的基本调试方法。二、实验准备根据课程容，熟悉MATLAB软件的使用方法，自行学习MATLAB程序的基 础语法，并根据所学知识编写潮流计算牛顿拉夫逊法（或PQ分解法）的计算 程序，用相应的算例在MATLAB ±进行计算、调试和验证。三、实验要求每人一组，在实验课时，用MATLAB调试和修改运行程序，用算例计算输 出潮流结果。四、实验程序clear; %清空存n=input（1iW输入节点数:心）；n I=InputCiW 输入支路数:nl 二'）；isb=in

2、putfiW输入平衡节点号Iisb=*）;pr=input（,iW输入误差精度:Pu）;Bl=input（,iW输入支路参数:BU）;B2=input（,输入节点参数:B2二'）；X=iput（节点号和对地参数:X*Y=ZerOs();Times=I;%-：创建节点导纳矩阵for i=l:nlifBl(i,6)=0%不含变压器的支路P=BI (i,l);q=Bl(L2);Y (p,q)=Y(p,q)-lBl(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+lBl (i,3)+0.5*Bl (i,4);Y (q,q)=Y(q,q)+lBl(i,3)+0.5*Bl(L4)

3、;else%含有变压器的支路P=BI (i,l);q=Bl(L2);Y (p,q)=Y(p,q)-l/(Bl(L3)*B1(L5);Y (q,p)=Y(p,q);Y(p,p)=Y(p,p)+lBl(i,3);Y (q,q)=Y(q,q)+l/(Bl(L5)2*B1 (i,3);endendY;%将C)rgS、DetaS初始化OrgS=ZerOS(2*n-2,l);DetaS=ZerOS(2*n-2,l);%二：创建OrgS,用于存储初始功率参数h=0;j 二 0;for i=l:n%对PQ节点的处理if i 二 isb&B2(iz6)=2h=h+l;for j=l:nOrgS(2*h-

4、lJ)=OrgS(2*h-l,l)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-im ag(Y(izj)*lmag(B2(jz3)+imag(B2(iz3)*(real(Y(i,j)*imag(B2(jz3)+ima g(Y(Lj)*real(B2(L3);OrgS(2*h,l)=OrgS(2*hJ)+imag(B2(L3)*(real(Y(ij)*real(B2(j,3)-imag (Y(Li)*imag(B2(jz3)-real(B2(i,3)*(real(Y(ij)*imag(B2(L3)+imag(Y(i,j )*real(B2(j,3);endend

5、end%三：对PV节点的处理，注慧这时不可再将h初始化为0for i=l:nif i-=isb&B2(i,6)=3h 二 h+1;OrgS(2FJI)=OrgS(23lthjl )+reol(B2(L3)广(reQ(Y(ij)FeQ(B2(j,3)-im g(Y(U)*img(B2(j,3)+img(B2(i,3)*(rel(Y(i,j)*img(B2(j,3)+im g(Y(U)*real(B2(i,3);OrgS(2*h,l)=OrgS(2*hJ)+imag(B2(L3)*(real(Y(ij)*real(B2(j,3)-imag (Y(Li)*imag(B2(L3)-real(B

6、2(L3)*(real(Y(Lj)*imag(B2(L3)+imag(Y(Lj )*real(B2(j,3);endendendOrgS;%四：创建PVU用于存储PV节点的初始电压PVU=ZerOS(-h-l,l);t 二 0;for i=l:nif B2(iz6)=3M+l;PV(tJ)=B2(i,3);endendPV;%五：创建DetaSZ用于存储有功功率、无功功率和电压幅值的不平衡量if i 二 isb&B2(i,6)=2h=h+l;DetaS(2*h-l, 1 )=real(B2(i,2)-OrgS(2*h-lJ);DetaS (2*hz 1) =imag (B2(i,2)

7、-OrgS (2*h, 1):endendt=0;for i=l:n%六：对PV节点的处理，注意这时不可再将h初始化为Oif i 二 isb&B2(i,6)=3h=h+l;t=t+l;DetaS(2*h-l J )=real(B2(i,2)-OrgS(2*h-U);DetaS(2*M)=real(PV(t,l)2+imag(PV(tzl)2-real(B2(i,3)2-imag(B2(L3)2;endendDetaS;%七：创建I,用于存储节点电流参数if i =isbh=h+l; l(hJ)=(OrgS(2*h-lJ)-OrgS(2*hJ)*sqrt(-l)conj(B2(L3);e

8、ndendI；%八：创建JaCbi傩可比矩阵)JaCbi=ZerOS(2*n-2);h=0;k 二 0;for i=l:n%对PQ节点的处理if B2(L6)=2h=h+l;for j=l:nif j=isbk 二 k+1;ifi=j%对角元素的处理JaCbi(2*h-L2*k-l)=-imag(Y(ij)*real(B2(L3)+real(Y(i,j)*imag(B2(i,3)+imag(l(h,l);JaCbi(2*h-l,2*k)=l(Y(i,j)*rel(B2(L3)+img(Y(iJ)*img(B2(i,3)+real(l(h,l);JaCbi(2*h,2*k-l )=-Jocbi(

9、2*h-l ,2*k)+2*reol(l(h,l);JaCbi(2*h,2*k)=Jacbi(2*h-L2*k-l)-2*imag(l(hJ);else%非对角元素的处理JaCbi(2*h-L2*k-l)=-imag(Y(Lj)*real(B2(l3)+real(Y(ij)*imag(B2(L3)/JaCbi(2*h-b2*k)=real(Y(Lj)*real(B2(l3)+imag(Y(ij)*imag(B2(L3);JaCbi(2*h,2*k-l )=-Jacbi(2*h-l ,2*k);JaCbi(2*h,2*k)=Jacbi(2*h-lr2*k-l);endif k二=(Fl-I) %

10、将用于循环的指针責于初始值，以确保雅可比矩阵换行k 二 0;endendendk 二 O;for i=l:n%对PV节点的处理if B2(L6)=3h=h+l;for j=l :nif j =isbk=k+l;ifi=j %对角元素的处理JaCbi(2*h-lz2*k-l)=-img(Y(izj)*rel(B2(i,3)+rel(Y(i,j)*img(B2(L3)+imag(l(h,l);JaCbi(2*h-l,2*k)=al(Y(i,j)*real(B2(L3)+imag(Y(iJ)*imag(B2(i,3)+real(l(hzl);JaCbi(2*h,2*k-l)=2*imag(B2(i,

11、3);JaCbi (2*hz2*k) =2*real (B2(i,3):else%非对角元素的处理JaCbi(2*h-L2*k-l)=-imag(Y(iJ)*real(B2(L3)+real(Y(ij)*imag(B2(L3)JaCbi(2*h-L2*k)=rel(Y(ij)*rel(B2(L3)+img(Y(izj)*img(B2(ir3);JaCbi(2*h,2*k-l )=0;JaCbi(2*h,2*k)=0;endif k=(n-l)%将用于循环的指针責于初始值，以确保雅可比矩阵换行k 二 0;endendendendendJacbi;%九：求解修正方程，获取节点电压的不平衡量Deta

12、U=ZerOS(2*n-2z 1);DetaU=inv(Jacbi)*DetaS;Detail;%修正节点电压j 二 0;for i=l:n%对PQ节点处理if B2(L6)=2i=i÷i;B2(i,3)=B2(iz3)+Det U (2*jz 1)+Deta U (2*j- IJ) *sqrt (-1);endfor i=l:n%对PV节点的处理if B2(i,6)=3i=i+i;B2(i,3)=B2(iz3)+Deta U (2*jf 1)+DetaU (2*j- LI) *sqrt (-1);endendB2;%十：开始循环*While abs(ma刈DetoU)>prO

13、rgS=ZerOS(2*n-2zl);%!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵 h 二 0;j=0;for i=l:nif i 二 isb&B2(L6)=2h=h+l;for j=l:nOrgS(2*h-l ,1 )=OrgS(2*h-l,l )+rel(B2(L3)*(rel(Y(i,j)*rel(B2(jz3)-im ag(Y(Li)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+ima g(Y(Lj)*real(B2(j,3);OrgS(2*h,l)=OrgS(2*hzl)+imag(B2(

14、L3)*(real(Y(ij)*real(B2(jz3)-imag (Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imog(B2(j,3)+imag(Y(i,j )*real(B2(j,3);endendendfor i=l:nif i-=isb&B2(i,6)=3h=h+l;for j=l: nOrgS(2*h-l ,1 )=OrgS(2*h-l,l )+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-im ag(Y(ij)*iag(B2(j,3)+imag(B2(i,3)*(real(Y(Li)*imag

15、(B2(jz3)+ima g(Y(Lj)*real(B2(j,3);OrgS(2*h,l)=OrgS(2*hJ)+imag(B2(L3)*(real(Y(ij)*real(B2(j,3)-imag (Y(Lj)*imag(B2(L3)-real(B2(L3)*(real(Y(ij)*imag(B2(L3)+imag(Y(Lj )*real(B2(L3);endOrgS;%创建DetaSh=0;for i=l:nif i=isb&B2(i,6)=2h=h+l;DetaS (2*h-l, 1 )=real (B2(L2)-OrgS (2*h-l,l);DetaS (2*h, 1) =ima

16、g (B2(i,2) -OrgS (2*h, 1):endendt=0;for i=l :nif i-=isb 8<B2(L6)=3h 二 h+1;t=t+l;DetaS(2*h-l, 1 )=real(B2(L2)-OrgS(2*h-lJ);DetaS(2*h,l)=real(PV(tzl)2+imag(PV(t,l)2-real(B2(iz3)2-ima g(B2(i,3)2;endDetS;%创建IRZerOs (n-l,l);h=0;for i=l:nif i 二 isbh 二 h+1; l(hJ)=(OrgS(2*h-U)-OrgS(2*hJ)*sqrt(-l)conj(B2(

17、L3);endendI；%创建JaCbiJaCbi=ZerOS(2*-2);h=0;k 二 O;for i=l:nif B2(i,6)=2h 二 h+1;for j=l:nif j 二 isbk=k+l;if i=jJaCbi(2*h-l,2*k-l)=-img(Y(ij)*rel(B2(L3)+rel(Y(ij)*img(B2(L3)+imag(l(h,l);JaCbi(2*h-lz2*k)=real(Y(i,j)*real(B2(L3)+imag(Y(ij)*imag(B2(i,3)+real(l(hzl);JaCbi(2*hz2*k-l )=-Jacbi(2*h-l ,2*k)+2*re

18、al(l(hJ);JaCbi(2*h,2*k)=Jacbi(2*h-L2*k-l)-2*imag(l(hJ);elseJaCbi(2*h-b2*k-l)=-imag(Y(ij)*real(B2(i,3)+real(Y(Lj)*imag(B2(L3)JaCbi(2*h-b2*k)=real(Y(ij)*real(B2(iz3)+imag(Y(iJ)*imag(B2(i,3);JaCbi(2*h,2*k-l )=-Jacbi(2*h-l ,2*k);JaCbi(2*hz2*k)=Jacbi(2*h-L2*k-1);endif k=(-l)k 二 O;endendendendk 二 O;for i=

19、l:nif B2(L6)=3h=h+l;for j=l:nif j 二 isbk=k+l;if i=jJaCbi(2*h-L2*k-l)=-img(Y(Lj)*rel(B2(L3)+rel(Y(ij)*img(B2(iz3)+imag(l(h,l);JaCbi(2*h-L2*k)=real(Y(i,j)*real(B2(L3)+imag(Y(ij)*imag(B2(i,3)+r eal(l(hzl);JaCbi (2*h,2*k-1 )=2*imag (B2(i,3):JaCbi(2*h,2*k)=2*real(B2(L3);elseJaCbi(2*h-L2*k-l)=-img(Y(Lj)*r

20、el(B2(L3)+rel(Y(ij)*img(B2(iz3)JaCbi(2*h-L2*k)=real(Y(Lj)*real(B2(L3)+imag(Y(ij)*imag(B2(i,3);JaCbi(2*h,2*k-l )=0;JaCbi(2*h,2*k)=0;endif k=(-l)k 二 0;endendendendendJacbi;DetaU=ZerOS(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaiI;%修正节点电压j 二 O;for i=l:nif B2(L6)=2i=i÷i;B2(i3)二 B2(i3)+DetoU(2Sl)+DetoU(2牛

21、1广Sqrt 卜1);endendfor i=l:nif B2(i,6)二=3i=i÷i;B2(i,3)=B2(iz3)+Deta U (2*jf 1)+DetaU (2*j- LI) *sqrt (-1);endendB2;TimeS=TimeS+1;% 迭代次数加 1endTimes;五、实验流程六、实验结果参数输入:命令行SB话辎入节点g$：n=5i5 入支 83?i：nl-5谛输入平衢节点号:ib=5iS入逞差桶度:爼=0OOoOOI请笹入支SS:Bl=CI 2 0. 04*0. 25i 0. 6i 1 0;1 3 0. 1*0.36i 0 1 0:2 3 0. 08-0.

22、30i O. 5i 1 0;4 2 O. OlSi O 1.05 1；5 3 0.03i O 1.05 1请S入节点參数:B2=S -1.6-0. 81 1 O O 2:0 -2-Ii IoO 2:0 -3. 7-1.31 1 0 0 2;O SfOi 1.05 1.05 O 3.0 O 1.05 1.05 O 1节点号和对地琴数:X=U 0;2 0.3 0.4 0.5 OJ运行结果：l_f DetaS x!_J DetaU X l_J OrgSR-； 8x1 doublet=H 8x1 double 圧 8x1 dou1111234-7.0314.1-1.4654.11.6000-2.663

23、9.2-8.3690.2 -0.80001.3407e.31.9428e 3-2.0000-1.1136.4-6.71014-1.00005-3.6230.5-3.0339.56-3.70006 9612.6-7.3178.-1.300073.8685e.73.91 79e.785.00008-2.1913.8-T68291.81 314x1 COmPleXdout1-1.7721 + 1.0793i2.0505 + 03142i-3.4664 + 1.5174i5.0625 + 0.1690i1234511.3787-6.2917i-0.6240 + 3.9002i-0.7547 +2.6415i0.0000 + 0.0OOOi0.0000 + 0.0OOOi20.6240 十 3.9002i1.4539 669808i-0.8299 +3.1120i0.0000 + 63.4921 i0.0000 + 0.0OOOi3-0.7547 t 2.6415i-0.829