有关MATLABdPQ分解法的潮流程序

xiaotiejiang523 · 发表于 2010-12-25 19:45:49

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

您需要登录才可以下载或查看，没有账号？立即加入

×

如题，首先明确一点，我不太懂。我照原程序抄的程序，运行后发现迭代次数不一样，结果大多数相同，有个别的第四位小数不同，这是为什么呢？谢谢。

facat · 发表于 2010-12-26 01:40:05

你要把代码拿上来给大家看啊

jcxx01 · 发表于 2010-12-26 09:04:24

看看会不会是设置的迭代精度不一样我觉得如果代码没有问题的话有可能的就是你的matlab设置的计算精度和原来的不一样这个我觉得不是什么问题呵呵可以的重要的是理解算法本身~加油~

xiaotiejiang523 · 发表于 2010-12-29 12:02:06

本帖最后由 xiaotiejiang523 于 2010-12-29 12:03 编辑

以下是程序代码，其中从Np(K)=ICT2;这条语句开始到if与end之间没看懂，我觉得应该是LU分解，来解方程，但好像不是那样。求高人解答啊。在此先谢谢了！

%P-Q分解法进行潮流计算
data=zeros(1,4)
data=load('d:\MATLAB7\work\data.txt')
disp('节点数：')
n=data(1)
disp('支路数：')
nl=data(2)
disp('平衡节点编号：')
isb=data(3)
disp('误差精度：')
pr=data(4)
disp('PQ节点数：');
na=data(5)
disp('由支路参数形成的矩阵：')
B1=dlmread('d:\MATLAB7\work\B1data.txt')
disp('各节点参数形成的矩阵：')
B2=dlmread('d:\MATLAB7\work\B2data.txt')
Y=zeros(n);YI=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);th=zeros(1,n);
for i=1:nl
if B1(i,6)==0;
      p=B1(i,1);q=B1(i,2);
   else p=B1(i,2);q=B1(i,1);
   end
   Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));
   YI(p,q)=YI(p,q)-1./B1(i,3);
   Y(q,p)=Y(p,q);
   YI(q,p)=YI(p,q);
   Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;
   YI(q,q)=YI(q,q)+1./B1(1,3);
   Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;
   YI(p,p)=YI(p,p)+1./B1(i,3);
end
%求导纳矩阵
G=real(Y);B=imag(YI);BI=imag(Y);
for i=1:n
S(i)=B2(i,1)-B2(i,2);
BI(i,i)=BI(i,i)+B2(i,5);
end
P=real(S);Q=imag(S);
for i=1:n
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
V(i)=B2(i,4);
end
for i=1:n
if B2(i,6)==2
      V(i)=sqrt(e(i)^2+f(i)^2);
      th(i)=atan(f(i)./e(i));
end
end
for i=2:n
if i==n
      B(i,i)=1./B(i,i);
else IC1=i+1;
      for j1=IC1:n
         B(i,j1)=B(i,j1)./B(i,i);
      end
      B(i,i)=1./B(i,i);
      for k=i+1:n
         for j1=i+1:n
            B(k,j1)=B(k,j1)-B(k,i)*B(i,j1);
         end
      end
end
end
p=0;q=0;
for i=1:n
if B2(i,6)==2
      p=p+1;k=0;
      for j1=1:n
         if B2(j1,6)==2
            k=k+1;
            A(p,k)=BI(i,j1);
         end
      end
end
end
for i=1:na
if i==na
      A(i,i)=1./A(i,i);
else k=i+1;
      for j1=k:na
      A(i,j1)=A(i,j1)./A(i,i);
      end
      A(i,i)=1./A(i,i);
      for k=i+1:na
         for j1=i+1:na
            A(k,j1)=A(k,j1)-A(k,i)*A(i,j1);
         end
      end
end
end
ICT2=1;ICT1=0;kp=1;kq=1;K=1;DET=0;ICT3=1;
while ICT2~=0|ICT3~=0;
ICT2=0;ICT3=0;
for i=1:n
      if i~=isb
         C(i)=0;
         for k=1:n
            C(i)=C(i)+V(k)*(G(i,k)*cos(th(i)-th(k))+BI(i,k)*sin(th(i)-th(k)));
         end
         DP1(i)=P(i)-V(i)*C(i);
         DP(i)=DP1(i)./V(i);
         DET=abs(DP1(i));
         if DET>=pr
            ICT2=ICT2+1;
         end
      end
end
Np(K)=ICT2;
if ICT2~=0
      for i=2:n
         DP(i)=B(i,i)*DP(i);
         if i~=n
            IC1=i+1;
            for k=IC1:n
                  DP(k)=DP(k)-B(k,i)*DP(i);
            end
         else
            for LZ=3:i
                  L=i+3-LZ;
                  IC4=L-1;
                  for MZ=2:IC4
                     I=IC4+2-MZ;
                     DP(I)=DP(I)-B(I,L)*DP(L);
                  end
            end
         end
      end
      for i=2:n
         th(i)=th(i)-DP(i);
      end
      kq=1;L=0;
      for i=1:n
         if B2(i,6)==2
            C(i)=0;L=L+1;
            for k=1:n
                  C(i)=C(i)+V(k)*(G(i,k)*sin(th(i)-th(k))-BI(i,k)*cos(th(i)-th(k)));
            end
            DQ1(i)=Q(i)-V(i)*C(i);
            DQ(L)=DQ1(i)./V(i);
            DET=abs(DQ1(i));
            if DET>=pr
                  ICT3=ICT3+1;
            end
         end
      end
else kp=0;
      if kq~=0
         L=0;
         for i=1:n
            if B2(i,6)==2
                  C(i)=0;L=L+1;
                  for k=1:n
                     C(i)=C(i)+V(k)*(G(i,k)*sin(th(i)-th(k))-BI(i,k)*cos(th(i)-th(k)));
                  end
                  DQ1(i)=Q(i)-V(i)*C(i);
                  DQ(L)=DQ1(i)./V(i);
                  DET=abs(DQ1(i));
                  if DET>=pr
                     ICT3=ICT3+1;
                  end
            end
         end
      end
end
Nq(K)=ICT3;
if ICT3~=0;
      L=0;
      for i=1:na
         DQ(i)=A(i,i)*DQ(i);
         if i==na
            for LZ=2:i
                  L=i+2-LZ;
                  IC4=L-1;
                  for MZ=1:IC4
                     I=IC4+1-MZ;
                     DQ(I)=DQ(I)-A(I,L)*DQ(L);
                  end
            end
         else
            IC1=i+1;
            for k=IC1:na
                  DQ(k)=DQ(k)-A(k,i)*DQ(i);
            end
         end
      end
      L=0;
      for i=1:n
         if B2(i,6)==2
            L=L+1;
            V(i)=V(i)-DQ(L);
         end
      end
      kp=1;
      K=K+1;
else
      kq=0;
      if kp~=0
         K=K+1;
      end
end
end
disp('迭代次数');
disp(K);
disp('每次没有达到精度要求的有功功率个数为');
disp(Np);
disp('每次没有达到精度要求的无功功率个数为');
disp(Nq);
for k=1:n
E(k)=V(k)*cos(th(k))+V(k)*sin(th(k))*j;
th(k)=th(k)*180./pi;
end
disp('各节点的实际电压标幺值E为（节点号从小到大排列）：');
disp(E);
disp('各节点的电压大小V为（节点号从小到大排列）：');
disp(V);
disp('各节点的电压相角th为（节点号从小到大排列）：');
disp(th);
for p=1:n
C(p)=0;
for q=1:n
      C(p)=C(p)+conj(Y(p,q))*conj(E(q));
end
S(p)=E(p)*C(p);
end
disp('各节点的功率S为（节点号从小到大排列):');
disp(S);
disp('各条支路的首端功率Si为（顺序同您输入B1时一样）：');
for i=1:nl
if B1(i,6)==0
      p=B1(i,1);q=B1(i,2);
else p=B1(i,2);q=B1(i,1);
end
Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
disp(Si(p,q));
end
disp('各条支路的末端功率Sj为（顺序同您输入B1时一样）：');
for i=1:nl
if B1(i,6)==0
      p=B1(i,1);q=B1(i,2);
else p=B1(i,2);q=B1(i,1);
end
Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
disp(Sj(q,p));
end
disp('各条支路的功率损耗DS为（顺序同您输入B1时一样）：');
for i=1:nl
if B1(i,6)==0
      p=B1(i,1);q=B1(i,2);
else p=B1(i,2);q=B1(i,1);
end
DS(i)=Si(p,q)+Sj(q,p);
disp(DS(i));
end

gxuyang · 发表于 2011-1-7 10:22:52

en~~~~第四位小数点不同而已，这个已经是对的了，同一个数据，就算大家算法一样，都不可能让结果一模一样的。

zhangkefu · 发表于 2011-1-11 10:18:09

回复 4# xiaotiejiang523

你好，可以把你的数据'd:\MATLAB7\work\data.txt'发给我吗？
我最近也在看潮流计算。可以一起交流交流
我的邮箱：zhangkefulove@126.com
谢谢了

小虎 · 发表于 2011-5-30 02:58:38

这里真多牛人，太佩服你们了。

账号		自动登录	找回密码
密码			立即加入

有关MATLABdPQ分解法的潮流程序

马上加入，结交更多好友，共享更多资料，让你轻松玩转电力研学社区！

相关帖子