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

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

2、);isb=input( 请输入平衡节点号 :isb=);pr=input( 请输入误差精度 :pr=);B1=input( 请输入支路参数 :B1=);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/B1(i,3)+0.5*B1(i,4);Y(q,q)=Y(q,q)+1/B1(

3、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);Y(q,q)=Y(q,q)+1/(B1(i,5)2*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节点的处理if i=isb&B2(i,6)=2 h=h+1;

4、for j=1:n OrgS(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);end end end%三：对 PV节点的处理，注意这时不可再将 h 初始化为 0 for i=1:nif i=isb&B2(i,6)=3 h=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)*(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)=3 t=t+1; PVU(t,1)=B2(i,3);endendPVU;%五：创建 DetaS, 用于存储有功功率、无功功率和电压幅值的不平衡量 h=0;for i=1:n % 对 PQ节点的处理

7、 if i=isb&B2(i,6)=2 h=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:n%六：对 PV节点的处理，注意这时不可再将 h 初始化为 0 if i=isb&B2(i,6)=3 h=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)2+imag(PVU(t,1)2-real(B2(i, 3)2-imag( B2(i

8、,3)2;endend DetaS;%七：创建 I ，用于存储节点电流参数 i=zeros(n-1,1); h=0;for i=1:n if i=isb h=h+1;I(h,1)=(OrgS(2*h-1,1)-OrgS(2*h,1)*sqrt(-1)/conj(B2(i,3); endend I;%八：创建 Jacbi( 雅可比矩阵 )Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n %对 PQ节点的处理if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j % 对角元素的处理Jacbi(2*h-1,2*k-1)=-im

9、ag(Y(i,j)*real(B2(i,3)+real(Y(i,j)*imag( B2(i,3)+imag(I(h,1);Jacbi(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)+r

10、eal(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)=-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:n %对 PV节点的处理if B2(i,6)=3h=h+1;for j=1:nif j=isbk=k+1;if i=j % 对角

11、元素的处理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,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)*i

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

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

14、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)*(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)*

15、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)*(real(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

16、(B2(j,3)-real(B2(i,3)*(real(Y(i,j)*imag(B2( j,3)+imag(Y(i,j)*real(B2(j,3);endendendOrgS;%创建 DetaSh=0;for i=1:nif i=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-

17、1,1); DetaS(2*h,1)=real(PVU(t,1)2+imag(PVU(t,1)2-real(B2(i,3)2-im ag(B2(i,3)2;end end DetaS; %创建 I i=zeros(n-1,1); h=0; for i=1:n if i=isb h=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:n if B2(i,6)=2 h=h+1; for j=1:n if j=isb k=k+1

18、; 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);Jacbi(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); elseJacbi(2*h-1,2*k-1)=-imag(Y

19、(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)=-Jacbi(2*h-1,2*k); Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1); end if k=(n-1)k=0; end end end end end k=0; for i=1:n if B2(i,6)=3 h=h+1; for j=1:n if j=isb k=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);Jacbi(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)=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)*r