《潮流计算的MATLAB源程序》由会员分享，可在线阅读，更多相关《潮流计算的MATLAB源程序（7页珍藏版）》请在金锄头文库上搜索。
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(n);Times=1; % 置迭代次数为初始值%创建节点导纳矩阵for i=1:n1if B1(i,6)=0%不含变压器的支路p=B1
3、(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);丫(q,q)=Y(q,q)+1/(B1(i,5f2*B1(i,3);endendYOrgS=zeros(2*n-2,1);DetaS=zer
4、os(2*n-2,1);%将 OrgS、DetaS 初始化%创建 O rg S ，用于存储初始功率参数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)*real(B2(j,3)-imag(Y(i,j)*imag(B2(j,3)+imag(B 2(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
5、(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:n%对PV节点的处理，注意这时不可再将h初始化为0if 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(B 2(i
6、,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%创建 PVU 用于存储 PV 节点的初始电压PVU=zeros(n-h-1,1);t=0;for i=1:nif B2(i,6)=3t=t+1;PVU(t
7、,1)=B2(i,3);endendPVU%创建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初始化为0if 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);De
8、taS(2*h,1)=real(PVU(t,1)A2+imag(PVU(t,1)A2-real(B2(i,3)f2-imag(B2(i,3)A2; 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(雅可比矩阵)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=
9、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,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
10、(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)=-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
11、=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,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)=
12、-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;endif k=(n-1)%将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0;endendendendendJacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros(2*n-2,1);DetaU=inv(Jacbi)*DetaS;DetaU%修正节点电压j=0;fo
13、r 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);endendB2% 开始循环 *while abs(max(DetaU)prOrgS=zeros(2*n-2,1); %! 初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵 h=0;j=0;for i=1:ni
14、f 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(B 2(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);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)+
