PQ分解法算潮流

%一.读取参数
clear  
format short
f=load（‘F:\电力系统分析作业\123.txt‘）;%打开路径下存放的IEEE14节点的参数数据
zljz=zeros（208）;
zljz=f（1:201:8）;
jdjz=zeros（148）;
jdjz=f（21:341:8）;
gnd=zeros（13）; 
gnd=f（351:3）;
%对地支路数
ngnd=f（361）;
%节点个数 
npoint=f（371）;  
%线路支路数 
nzl=f（381）;
error=1e-12;%容许误差
Nmax=600;%最大迭代次数

%二.形成节点导纳矩阵
  %PV节点个数
pvjd=0;
for i=1:npoint          
       if（jdjz（i2）==2）
        pvjd=pvjd+1;
       end
end
   %PQ节点个数 
 pqjd=npoint-pvjd-1;
 %矩阵维数
 ndim=npoint-1+pqjd; 
%找出平衡节点序号
 for i=1:npoint            
       if（jdjz（i2）==1）
           nphjd=i;
       end
 end  
 %将平衡节点调整至最后一行
   jdjz1=zeros（18）;    
   jdjz1（1:）=jdjz（nphjd:）;
   jdjz（nphjd:）=jdjz（npoint:）;
   jdjz（npoint:）=jdjz1（1:）;
   for k=1:nzl             %更新线路
      i=1;
      for i=1:2
          if（zljz（ki）==nphjd）
          zljz（ki）=npoint;
          elseif（zljz（ki）==npoint）
              zljz（ki）=nphjd;
            end
        end
   end
%初始化各导纳矩阵，Y，YR（实部），YI（虚部）  
YR=zeros（npointnpoint）;
YI=zeros（npointnpoint）;
Y=zeros（npointnpoint）;
%1.不含变压器的线路导纳矩阵
k=1;
for k=1:nzl                
   i=zljz（k1）;
   j=zljz（k2）;
   mo=zljz（k3）^2+zljz（k4）^2;
    if（zljz（k7）==0）
   YR（ii）=YR（ii）+zljz（k3）/mo+0.5*zljz（k5）;
   YI（ii）=YI（ii）-zljz（k4）/mo+0.5*zljz（k6）;
   YR（ij）=YR（ij）-zljz（k3）/mo;
   YI（ij）=YI（ij）+zljz（k4）/mo;
   YR（ji）=YR（ji）-zljz（k3）/mo;
   YI（ji）=YI（ji）+zljz（k4）/mo;
   YR（jj）=YR（jj）+zljz（k3）/mo+0.5*zljz（k5）;
   YI（jj）=YI（jj）-zljz（k4）/mo+0.5*zljz（k6）;
   end
end
Y=YR+YI*1i;
%2.含变压支路的线路导纳矩阵
i=1;
for i=1:nzl  
   if（zljz（i7）==1）
       mo=zljz（i3）^2+zljz（i4）^2;
       YR（zljz（i1）zljz（i1））=（zljz（i8）^2）*zljz（i3）/mo+YR（zljz（i1）zljz（i1））;
       YI（zljz（i1）zljz（i1））=-（zljz（i8）^2）*zljz（i4）/mo+YI（zljz（i1）zljz（i1））;
       YR（zljz（i2）zljz（i2））=zljz（i3）/mo+YR（zljz（i2）zljz（i2））;
       YI（zljz（i2）zljz（i2））=-zljz（i4）/mo+YI（zljz（i2）zljz（i2））;
       YR（zljz（i1）zljz（i2））=-zljz（i8）*zljz（i3）/mo+YR（zljz（i1）zljz（i2））;
       YI（zljz（i1）zljz（i2））=-zljz（i8）*-zljz（i4）/mo+YI（zljz（i1）zljz（i2））;
       YR（zljz（i2）zljz（i1））=-zljz（i8）*zljz（i3）/mo+YR（zljz（i2）zljz（i1））;
       YI（zljz（i2）zljz（i1））=-zljz（i8）*-zljz（i4）/mo+YI（zljz（i2）zljz（i1））;
   end
end
  Y=YR+YI*1i; 
%3.对地并联支路的导纳矩阵修正
i=1;
  for i=1:ngnd               
   c=gnd（i1）;
    YR（cc）=YR（cc）+gnd（i2）;
    YI（cc）=YI（cc）-gnd（i3）;
end
   Y=YR+YI*1i;
disp（‘节点导纳矩阵为：‘） 
Y
%三.B’B’’矩阵的生成
%B’矩阵生成
B=imag（Y）; 
ph=find（jdjz（:2）==1）; 
BB=B; 
BB（:ph）=[]; BB（ph:）=[]; 
B1=BB;
%B’’矩阵的生成
phpv=find（jdjz（:2）~=3）; 
       BB=B;  
 BB（:phpv）=[]; 
 BB（phpv:）=[];      
  B2=BB;  
%四.计算不平衡量dPdQ
f=0;%迭代标志位
N=1;%迭代次数初始化
%计算各节点净有功
Pjd=zeros（1npoint）;
i=1;
for i=1:npoint
    Pjd（i）=jdjz（i3）-jdjz（i5）;
    end
%节点净无功功率
Qjd=zeros（1npoint）;
i=1;
for i=1:npoint
    Qjd（i）=jdjz（i4）-jdjz（i6）;
    end
%计算各支路功率有功和无功:统一潮流方程计算每个节点的功率大小
while（f==