潮流计算的MATLAB源程序Word下载.docx

1、%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（i,1）;q=B1（i,2）;Y（p,q）=Y（p,

2、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 %含有变压器的支路Y（p,q）=Y（p,q）-1/（B1（i,3）*B1（i,5）;Y（p,p）=Y（p,p）+1/B1（i,3）;Y（ q,q）二 Y（ q,q）+1/（B1（i,5F2*B1（i,3）;endYOrgS=zeros（2*n-2,1）;DetaS二zeros（2*n-2,1）; %将 OrgS、DetaS初始化%创建OrgS,用于存储初始功率参数h=0;j=0;if i=

3、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）*（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）*im

4、ag（B2（j,3）+imag（Y（i,j）*real（B2（j,3）;n %对 PV 节点的处理，注意这时不可再将 h 初始化为 0B2（i,6）=3OrgS%创建 PVU 用于存储 PV 节点的初始电压PVU=zeros（n-h-1,1）;t=0;if B2（i,6）=3t=t+1;PVU（t,1）=B2（i,3）;PVU%创建DetaS,用于存储有功功率、无功功率和电压幅值的不平衡量n %对 PQ 节点的处理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）;n %对 PV 节点

5、的处理,注意这时不可再将 h 初始化为 0DetaS（2*h,1）=real（PVU（t,1）A2+imag（PVU（t,1）F2-real（B2（i,3）A2-imag（B2（i,3）F2;DetaS%创建 I ，用于存储节点电流参数i=zeros（n-1,1）;if i=isbI（h,1）=（OrgS（2*h-1,1）-OrgS（2*h,1）*sqrt（-1）/conj（B2（i,3）;%创建Jacbi雅可比矩阵）Jacbi=zeros（2*n-2）;k=0;if B2（i,6）=2if j=isbk=k+1;if i=j %对角元素的处理Jacbi（2*h-1,2*k-1）=-imag（

6、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,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）;Jacbi（2*h-1,2*k-1）=-imag（Y（i,j）*real（B2（i,3）+real（Y（i,j）*imag（B2（i,3）

7、; 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）;if k=（n-1） % 将用于内循环的指针置于初始值，以确保雅可比矩阵换行 k=0;n %对 PV 节点的处理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）=real（Y（i,j

8、）*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;if k=（n-1） %将用于内循环的指针置于初始值，以确保雅可比矩阵换行Jacbi%求解修正方程，获取节点电压的不平衡量DetaU=zeros（2*n-2,1）;DetaU=inv（Jacbi）*DetaS;DetaU%修正节点电压 j=0;j=j+1;B2（i,3）=B2（i,3）+DetaU（2*j,1）+DetaU（2*j-1,1）*sqrt（-1）;B2%开始循环 * while abs（max（DetaU）pr %!初始功率参数在迭

9、代过程中是不累加的，所以在这里必须将其初始化为零矩阵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）;end for i=1:%创建 DetaSDetaS（2*h,1）=imag（B2（i,2）-OrgS（2*h,1）;%创建 I i=zeros（n-1,1）;I%创建 Jacbiif i=j（I（h,1）;Jacbi（2*h-1,2

10、*k-1）=-imag（Y（i,j）*real（B2（i,3）+real（Y（i,j）*imag（B2（i,3）+imagelseend if k=（n-1）if k=（n-1）%修正节点电压Times=Times+1; %迭代次数加 1Times一个原始数据的例子节点数 5支路数 5平衡节点编号 5精度 pr 0.000001B1（支路参数矩阵）1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 30.03i 0 1.05 1B2（节点参数矩阵） 0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1 X（ 节点号和对地参数 ）1 0;2 0;3 0;4 0;5 0