有关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

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

页: [1]

电力研学论坛's Archiver

有关MATLABdPQ分解法的潮流程序