潮流计算PQ分解法

%----------节点导纳矩阵程序lfybus.m----------
tic   % 计时开始
i=sqrt（-1）;
j=sqrt（-1）;
nl=linedata（:1）; 
nr=linedata（:2）; 
R=linedata（:3）;
X=linedata（:4）; 
Bc=j*linedata（:5）; 
a=linedata（: 6）;
nbr=length（linedata（:1））; 
nbus=max（max（nl） max（nr））;
Z=R + j*X; 
y=ones（nbr1）./Z; %支路导纳
for n=1:nbr
    if a（n）<=0  
       a（n）=1; 
    else end
    Ybus=zeros（nbusnbus）; % 将Ybus初始化为0
% ---非对角元素的形成---
    for k=1:nbr;
        Ybus（nl（k）nr（k））=Ybus（nl（k）nr（k））-y（k）/a（k）;
        Ybus（nr（k）nl（k））=Ybus（nl（k）nr（k））;
    end
end
%---对角元素的形成---
for  n=1:nbus
     for k=1:nbr
         if nl（k）==n
            Ybus（nn）=Ybus（nn）+y（k）/（a（k）^2） + Bc（k）;
         elseif nr（k）==n
                Ybus（nn）=Ybus（nn）+y（k） +Bc（k）;
         else end
     end
end


%----------PQ分解法潮流程序decouple.m----------
%1.快速解耦法电力系统潮流计算（PQ分解法）
%2.数据准备;用户需要在命令窗口编写busdata和linedata两个数据文件
basemva=100;  
accuracy=0.0001; 
maxiter=50;
ns=0; 
Vm=0; 
delta=0; 
yload=0; 
deltad=0;
nbus = length（busdata（:1））;
kb=[];
Vm=[]; 
delta=[]; 
Pd=[]; 
Qd=[]; 
Pg=[]; 
Qg=[]; 
Qmin=[]; 
Qmax=[];
Pk=[]; 
P=[]; 
Qk=[]; 
Q=[]; 
S=[]; 
V=[];
for k=1:nbus
    n=busdata（k1）;
    kb（n）=busdata（k2）; 
    Vm（n）=busdata（k3）; 
    delta（n）=busdata（k4）;
    Pd（n）=busdata（k5）; 
    Qd（n）=busdata（k6）; 
    Pg（n）=busdata（k7）; 
    Qg（n） = busdata（k8）;
    Qmin（n）=busdata（k 9）; 
    Qmax（n）=busdata（k10）;
    Qsh（n）=busdata（k11）;
    if Vm（n） <= 0  
        Vm（n） = 1.0; 
        V（n） = 1 + j*0;
    else delta（n） = pi/180*delta（n）;
         V（n） = Vm（n）*（cos（delta（n）） + j*sin（delta（n）））;
         P（n）=（Pg（n）-Pd（n））/basemva;
         Q（n）=（Qg（n）-Qd（n）+ Qsh（n））/basemva;
         S（n） = P（n） + j*Q（n）;
    end
 if kb（n） == 1
      ns = ns+1; 
 else end
nss（n） = ns;
end
Ym = abs（Ybus）; 
t = angle（Ybus）;
ii=0;
for ib=1:nbus
     if kb（ib） == 0 | kb（ib） == 2
     ii = ii+1;
     jj=0;
         for jb=1:nbus
             if kb（jb） == 0 | kb（jb） == 2
                jj = jj+1;
                B1（iijj）=imag（Ybus（ibjb））;
             elseend
         end
     else end
end

ii=0;
for ib=1:nbus
    if kb（ib） == 0
       ii = ii+1;
       jj=0;
       for jb=1:nbus
            if kb（jb） == 0
               jj = jj+1;
               B2（iijj）=imag（Ybus（ibjb））;
             elseend
       end
    else end
end
B1inv=inv（B1）; 
B2inv = inv（B2）;
maxerror = 1; 
converge = 1; 
iter = 0;
mline=ones（nbr1）;
for k=1:nbr
    for m=k+1:nbr
        if（（nl（k）==nl（m）） & （nr（k）==nr（m）））;
           mline（m）=2;
         elseif （（nl（k）==nr（m）） & （nr（k）==nl（m）））;
                 mline（m）=2;
         else end
    end
end

% decouple程序开始迭代
while maxerror >= accuracy & iter <= maxiter   %检验不平衡功率
      iter = iter+1;
      id=0; 
      iv=0;
for n=1:nbus
    nn=n-nss（n）;
    J11=0;  
    J33=0;
	for ii=1:nbr
		if mline（ii）==1   
     		if nl（ii） == n | nr（ii） == n
        	    if nl（ii） == n