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

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

2、#39;）;isb=input（'请输入平衡节点号：isb='）;pr=input（'请输入误差精度:pr='）;B仁input('请输入支路参数:B仁');B2=input('请输入节点参数：B2=');X=input('节点号和对地参数:X=');Y=zeros(n);Times=1;%一：创建节点导纳矩阵for i=1:n1if B1(i,6)=0%不含变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/B1(i,3);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1

3、/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(i,3)+0.5*B1(i,4);else%含有变压器的支路p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)*B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/B1(i,3);丫(q,q)=Y(q,q)+1/(B1(i,5F2*B1(i,3);endendY;%将OrgS、DetaS初始化OrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1);%二：创建OrgS,用于存储初始功率参数h=0;j=0;for i=1:n%对PQ

4、节点的处理for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*I mag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+im

5、ag(Y(i,j)*real(B2(j,3);endendend%三：对 PV 节点的处理，注意这时不可再将 h 初始化为 0for i=1:nif i=isb&B2(i,6)=3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*im ag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,

6、3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendOrgS;%四：创建 PVU 用于存储 PV 节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)=3t=t+1;PVU(t,1)=B2(i,3);endendPVU;%五：创建DetaS用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n%对PQ节点的处理DetaS(2*h-1

7、,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);endendt=0;for i=1:n%六：对 PV 节点的处理，注意这时不可再将 h 初始化为 0if i=isb&B2(i,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1)F2+imag(PVU(t,1)F2-real(B2(i,3)F2-imag( B2(i,3)F2;endendDetaS;%七：创建 I ，用于存储节点电流

8、参数i=zeros(n-1,1);h=0;for i=1:nif i=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3);endendI;%八：创建Jacbi雅可比矩阵)Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n%对PQ节点的处理if B2(i,6)=2h=h+1;for j=1:nif j=isbk=k+1;if i=j%对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I

9、(h,1)JJacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1);else%非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B

10、2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);保雅可比矩阵换行k=0;endendendendendk=0;for i=1:n%对 PV 节点的处理if B2(i,6)=3h=h+1;for j=1:nif j=isbk=k+1;if i=j%对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1)Jacbi(2*h-1,

11、2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=2*imag(B2(i,3);Jacbi(2*h,2*k)=2*real(B2(i,3);else%非对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=0;Ja

12、cbi(2*h,2*k)=0;endif k=(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendJacbi;%九：求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU;%修正节点电压j=0;for i=1:n%对 PQ 节点处理if B2(i,6)=2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);for i=1:n%对 PV 节点的处理if B2(i,6)=3j=j+1;B2(i,3)=B2(i,3

13、)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2;%十：开始循环*while abs(max(DetaU)>prOrgS=zeros(2*n-2,1);%!初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i=isb&B2(i,6)=2h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(re

14、al(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendendfor i=1:nif i=isb&B2(i,6)=3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,1)+real(B2(i,3)*(rea

15、l(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);OrgS(2*h,1)=OrgS(2*h,1)+imag(B2(i,3)*(real(Y(i,j)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2(j,3)+imag(Y(i,j)*real(B2(j,3);endendOrgS;%创建 DetaSh=0;for i=1:nif i

16、=isb&B2(i,6)=2h=h+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=imag(B2(i,2)-OrgS(2*h,1);endendt=0;for i=1:nif i=isb&B2(i,6)=3h=h+1;t=t+1;DetaS(2*h-1,1)=real(B2(i,2)-OrgS(2*h-1,1);DetaS(2*h,1)=real(PVU(t,1)A2+imag(PVU(t,1)F2-real(B2(i,3)A2-imag(B2(i,3)F2;endendDetaS;%创建 Ii=zeros(n

17、-1,1);h=0;for i=1:nif i=isbh=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3);endendI;%创建 JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,6)=2h=h+1;for j=1:nif j=isbk=k+1;if i=jJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1)JJacbi(2*h-1,2*k)=real(Y(i,j)*

18、real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+real(I(h,1);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1);elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*

19、k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k=(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,6)=3h=h+1;for j=1:nif j=isbk=k+1;if i=jJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3)+imag(I(h,1)JJacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3)+re

20、al(I(h,1);Jacbi(2*h,2*k-1)=2*imag(B2(i,3);Jacbi(2*h,2*k)=2*real(B2(i,3);elseJacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3)+imag(Y(i,j)*imag(B2(i,3);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;end if k=(n-1)k=0;endendendendendJacbi;DetaU=zeros(

21、2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU;%修正节点电压j=0;for i=1:nif B2(i,6)=2j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendfor i=1:nif B2(i,6)=3j=j+1;B2(i,3)=B2(i,3)+DetaU(2*j,1)+DetaU(2*j-1,1)*sqrt(-1);endendB2;Times=Times+1;%迭代次数加 1endTimes;五、实验流程六、实验结果参数输入：静令行 口诘输入站辭mi好诰和节点号：i呦咄齢1八唉差瀚

22、度；时-匚CKJOCOl诒输A支踏垂散:肌=二1 2 0. 040. 251 0. 5i I 0; 1 3 0. 1>. 36i 0 I 0;2 3 O.MM).30i 0. Si 1 C;4 3 0,015i 0 l. 05 l；5 3 0. 031 0 1”。$ 1)请输人节点參戲:时-k6-0. 8i 1 0 0 2：C -'li 1 0 0 2. D -3- 7-1, 31 1002.Q i+Qi 1. Qa 1. Qo 0 2 0 0 1, Do 1,临 0 1韦点岸和对城夕齣:X-L 0 2 0,3 Q + 0;T 0运行结果：nJatbi:-j 3*8 double

23、12345671-0J357-3JDC0-0,0) GZ-2.2153-0.630 *a02-3.4D854.22720.GI62-3.30600.83812.2153005-420780.643669.517S-22.6941-3.4(5720.1770-&5.1433209B244064S6-4.207015«1-0.1770-3 +67-20324-65.14335-2.6716-O.9S44-3.1521-1.096533.3305O.936S00EO9S44-2P67161 09B5-7.fl69635 29570n700-61S30S?.3O460065.1433-?OQS?G0000000.78131.9492LJffiDeta8x1 doubleu田DetaU8x1 doubleOrgSI ； 8x1 doubku1 X L4x1 rfimnlmr rioub11111-L7721 + 1.0793i1| -7.0314. J1-1.4654.1| -1,60002-2.6639.22-o.sooo31.3407e.1.942fie.3-2,00002-2.0505 + 0.31424-1+1136,4-6.7101 “A1 .anno33.