潮流计算的MATLAB源程序

1、%简单潮流计算的牛顿拉夫逊程序，相关的原始数据数据数据输入格式如下： %B1是支路参数矩阵，第一列和第二列是节点编号。节点编号由小到大编写 %对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点 %编号，将变压器的串联阻抗置于低压侧处理。%第三列为支路的串列阻抗参数。%第四列为支路的对地导纳参数。%第五列为含变压器支路的变压器的变比 %第六列为变压器是否含有变压器的参数，其中“ 1”为含有变压器， %“0”为不含有变压器。%B2为节点参数矩阵，其中第一列为节点注入发电功率参数；第二列为节点负荷 功率参数；第三列为节点电压参数；第六列为节点类型参数，其中“ 1”为平衡 节点，“2”为

2、PQ节点，“3”为PV节点参数。%X为节点号和对地参数矩阵。其中第一列为节点编号，第二列为节点对地 %参数。n=input(' 请输入节点数 :n='); n1=input(' 请输入支路数 :n1='); isb=input(' 请输入平衡节点号 :isb='); pr=input(' 请输入误差精度 :pr=');B1=input(' 请输入支路参数 :B1=');B2=input(' 请输入节点参数 :B2=');X=input(' 节点号和对地参数 :X=');Y=zeros

3、(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(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);

4、丫(q,q)=Y(q,q)+1/(B1(i,5F2*B1(i,3);endendYOrgS=zeros(2*n-2,1);DetaS=zeros(2*n-2,1); %将 OrgS、DetaS 初始化%创建OrgS,用于存储初始功率参数h=0;j=0;for i=1:n%对PQ节点的处理if 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)*r eal(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B2(i,3)*(real(Y(i,j)*

5、 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);endendh 初始化为 0endfor i=1:n %对PV节点的处理，注意这时不可再将if i=isb&B2(i,6)=3h=h+1;for j=1:nOrgS(2*h-1,1)=OrgS(2*h-1,

6、1)+real(B2(i,3)*(real(Y(i,j)*r eal(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);endendendOrgS%

7、创建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);endend PVU %创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:n %对PQ节点的处理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

8、初始化为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)F 2-imag(B2(i,3)F2;endendDetaS%创建 I ，用于存储节点电流参数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

9、( 雅可比矩阵 )Jacbi=zeros(2*n-2);h=0;k=0;for i=1:n % 对PQ节点的处理if B2(i,6)=2h=h+1;for j=1:n if j=isb k=k+1;if i=j % 对角元素的处理 Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(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)=-Jac

10、bi(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)+re al(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

11、);endif k=(n-1) % 将用于内循环的指针置于初始值，以确保雅可 比矩阵换行k=0;endendendend end k=0;for i=1:n % 对PV节点的处理if B2(i,6)=3h=h+1;for j=1:n if j=isbk=k+1;if i=j %对角元素的处理Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(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)+rea

12、l(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)+re al(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;endif k=(n-1) % 将用于内循环的指针置于初始值，以确保 雅可比

13、矩阵换行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=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);enden

14、dB2%开始循环*while abs(max(DetaU)>prOrgS=zeros(2*n-2,1); %!初始功率参数在迭代过程中是不累加的，以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i=isb&B2(i,6)=2 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,

15、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)*(real(Y(i,j)*r eal(B2(j,3)-imag(Y(i,j)*imag(B2(j

16、,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);endendendOrgS%创建 DetaSh=0;for i=1:nif i=isb&B2(i,6)=2h=h+1;DetaS(2*h-1,1)=

17、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)F2+imag(PVU(t,1)F2-real(B2(i,3)F2-imag(B2(i,3)F2;endendDetaS%创建 Ii=zeros(n-1,1);h=0;for i=1:nif i=isbh=h+1; I(h,1)

18、=(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)=2 h=h+1; for j=1:n if j=isb k=k+1; if i=j Jacbi(2*h-1,2*k-1)=-imag(Y(i,j)*real(B2(i,3)+re al(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

19、(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)+re al(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(

20、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)=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)+re al(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)*real(B2(i,3)+re al(Y(i,j)*imag(B2(i,3);Jacbi(2*h-1,2*k)=real(Y(i,j)*real(B2(i,3