pq.m

基于MATLAB的PQ分解法潮流计算程序

clear
yiobus=[0.01494i;0.00225i;0;0.01805i];
vbus=[1.0*exp(i*0);1.0*exp(i*0);1.1*exp(i*0);1.05*exp(i*0)];
ybus=[1.042093-8.242876i,-0.588235+2.352941i,3.666667i,-0.453858+1.891074i;-0.588235+2.352941i,1.069005-4.727377i,0,-0.480769+2.403846i;3.666667i,0,-3.333333i,0;-0.453858+1.891074i,-0.480769+2.403846i,0,0.934627-4.261590i];
pxs=[-0.30;-0.55;0.5];
qxs=[-0.18;-0.13;0];
sizey=size(ybus);
ybusi=sizey(1,1);
ybusj=sizey(1,2);
ny=ybusi;
my=2;
vbusi=ybusi;


flag=1;
k=0;


while	flag==1
nbus=zeros(3,1);
mbus=zeros(3,1);
k=k+1;

kp=1;
kq=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for x=1:4
   for y=1:4
       gbus(x,y)=real(ybus(x,y));
   end
end
for x=1:4
   for y=1:4
       bbus(x,y)=imag(ybus(x,y));
   end
end
for x=1:4
   for y=1:4
       smybus(x,y)=angle(ybus(x,y));
   end
end
for x=1:4
       smvbus(x,1)=angle(vbus(x,1));
   end
 for x=1:3
    for y=1:4
         n=abs(vbus(y,1))*(gbus(x,y)*cos(smvbus(x,1)-smvbus(y,1))+bbus(x,y)*sin(smvbus(x,1)-smvbus(y,1)));
         nbus(x,1)=nbus(x,1)+n;
     end
 end
for x=1:3
 pbus(x,1)=abs(vbus(x,1))*nbus(x,1);
end
for x=1:3
 dp(x,1)=pxs(x,1)-pbus(x,1);
end

for x=1:3
           adp(x,1)=abs(dp(x,1));
   end
rn=max(adp)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if rn>=0.00001

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
for x=1:4
      
         absvbus(x,1)=abs(vbus(x,1));
end
for x=1:3
       for y=1:3
               if x==y
         absv1(x,y)=abs(vbus(x,1));
 else
          absv1(x,y)=0;
  end
  end
end
absv=absv1(1:3,:)
for x=1:3
        dpv(x,1)=dp(x,1)/absvbus(x,1);
end
for x=1:3
        for y=1:3
               bibus(x,y)=bbus(x,y);
       end
end
for x=1:2
        for y=1:2
               biibus(x,y)=bbus(x,y);
       end
end


vds=inv(-bibus)*dpv;
ds=inv(absv)*vds;
for x=1:3
        smvbus(x,1)=smvbus(x,1)+ds(x,1);
end
kq=1;

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

else
kp=0;
if kq==0
flag=0;
else
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for x=1:2
   for y=1:4
        m=abs(vbus(y,1))*(gbus(x,y)*sin(smvbus(x,1)-smvbus(y,1))-bbus(x,y)*cos(smvbus(x,1)-smvbus(y,1)));
        mbus(x,1)=mbus(x,1)+m;
   end
end
for x=1:2
 qbus(x,1)=abs(vbus(x,1))*mbus(x,1);
end
for x=1:2
 dq(x,1)=qxs(x,1)-qbus(x,1);
end
 for x=1:2
           adq(x,1)=abs(dq(x,1));
   end
   rqn=max(adq)
   
%%%%%%%%%%%%%%%%%%%%%%%


if  rqn>=0.00001

   %%%%%%%%%%%%%%%%%%%%%%%

for x=1:2
        dqv(x,1)=dq(x,1)/absvbus(x,1);
end

dv=inv(-biibus)*dqv
for x=1:2
        absvbus(x,1)=absvbus(x,1)+dv(x,1);
end

for x=1:4
       vbus(x,1)=absvbus(x,1)*exp(i*smvbus(x,1));
end
kp=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
kq=0;
if kp==0
flag=0;
else

end
end

end
p4jq4=0;
for x=1:4
    p4jq4=p4jq4+vbus(4,1)*conj(ybus(4,x))*conj(vbus(x,1)); 
end

   p4jq4
   save p4jq4.dat p4jq4
   
   save vbusl.dat vbus -ascii
   save absvbus.dat absvbus -ascii