谁能帮我看看这个潮流程序错在哪？

zhouyou529 · 发表于 2009-5-23 14:36:49

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

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

×

我自己编了个潮流的程序，现在只能在书上一个特定的算例上能够收敛，如果换成IEEE算例就不收敛了，不知道怎么回事情。
刚才看别人一个程序，它基本和我一样，只不过我用左除，他编写了一个高斯的程序。都没有考虑任何PQ。PV节点互换之类的问题。但是他的程序却能在IEEE算例收敛，不知道为什么？
我的程序如下
function newton
clear;
n=input('\n请输入节点数：n=');
m=input('\n请输入支路数：m=');
ph=input('\n请输入平衡节点的节点号:ph=');
B1=input('\n请输入支路信息：B1=');
% 他以矩阵的形式存贮支路的情况，每行贮存一条支路
% 第一列存贮支路的一个端点
% 第二列存贮支路的另一个端点
% 第三列存贮支路的阻抗
% 第四列存贮支路的对地导纳
% 第五列存贮支路的序号
B2=input('\n请输入节点信息：B2=');
% 第一列为电源侧的功率
% 第二列为负荷侧的功率
% 第三列为该点的电压值
% 第四列为该点的类型：1为PQ节点，2为PV节点，3为平衡节点
A=input('\n请输入节点号及对地阻抗：A=');
ip=input('\n请输入修正精度:ip=');
Y=zeros(n);
e=zeros(1,n);
f=zeros(1,n);
no=2*ph-1;
for i=1:n
if A(i,2)~=0;
   p=A(i,1);
   Y(p,p)=1./A(i,2);
end
end
for i=1:m
p=B1(i,1);
q=B1(i,2);
Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; % 自导纳形成
Y(p,q)=Y(p,q)-1./B1(i,3);
Y(q,p)=Y(p,q); % 互导纳形成
Y(q,q)=Y(q,q)+1./B1(i,3)+B1(i,4)./2;
end
G=real(Y);
B=imag(Y);
for i=1:n;
e(i)=real(B2(i,3));
f(i)=imag(B2(i,3));
S(i)=B2(i,1)-B2(i,2);
V(i)=B2(i,3);
end
P=real(S);
Q=imag(S);
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);
[De,Df]=hxf(J,DF,ph,n,no);
t=0;
while (max(abs(De))>ip|max(abs(Df))>ip)
t=t+1;
e=e+De;
f=f+Df;
[C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no);
J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no);
[De,Df]=hxf(J,DF,ph,n,no);
end
v=e'+f'*j;
for i=1:n
hh(i)=conj(Y(ph,i)*v(i));
end
S(ph)=sum(hh)*v(ph);
B2(ph,1)=S(ph);
V=abs(v);
jd=angle(v)*180/pi;
result1=[A(:,1),real(v),imag(v),V,jd,real(S.'),imag(S.'),real(B2(:,1)),imag(B2(:,1)),real(B2(:,2)),imag(B2(:,2))];
for i=1:m
a(i)=conj((v(B1(i,1))-v(B1(i,2)))/B1(i,3));
b(i)=v(B1(i,1))*a(i)-j*B1(i,4)*v(B1(i,1))^2/2;
c(i)=-v(B1(i,2))*a(i)-j*B1(i,4)*v(B1(i,2))^2/2;
end
result2=[B1(:,5),B1(:,1),B1(:,2),real(b.'),imag(b.'),real(c.'),imag(c.'),real(b.'+c.'),imag(b.'+c.')];
t
result1
S
b
c
result2

function [C,D,DF]=xxf(G,B,e,f,P,Q,n,B2,ph,V,no)  % 该子程序是用来求取DF
for i=1:n
if i~=ph
      C(i)=0;
      D(i)=0;
      for j=1:n
         C(i)=C(i)+G(i,j)*e(j)-B(i,j)*f(j);  % aii
         D(i)=D(i)+G(i,j)*f(j)+B(i,j)*e(j);  % bii
      end
      P1=C(i)*e(i)+D(i)*f(i);
      Q1=C(i)*f(i)-D(i)*e(i);
      V1=e(i)^2+f(i)^2;
      if B2(i,4)==2
         p=2*i-1;
         DF(p)=P(i)-P1;
         p=p+1;
         DF(p)=V(i)^2-V1^2;
      else
         p=2*i-1;
         DF(p)=P(i)-P1;
         p=p+1;
         DF(p)=Q(i)-Q1;
      end
end
end
DF=DF';
if ph~=n
DF(no,:)=[];
DF(no,:)=[];
end

function [De,Df]=hxf(J,DF,ph,n,no)  %  该子函数是为求取De Df
DX=J\DF;
DX1=DX;
x1=length(DX1);
if ph~=n
DX(no)=0;
DX(no+1)=0;
for i=(no+2):(x1+2)
      DX(i)=DX1(i-2);
end
else
DX=[DX1,0,0];
end
k=0;
[x,y]=size(DX);
for i=1:2:x
k=k+1;
Df(k)=DX(i);
De(k)=DX(i+1);
end

function J=jacci(Y,G,B,P,Q,e,f,V,C,D,B2,n,ph,no)  %  该子程序是用来求取jacci矩阵
for i=1:n
switch B2(i,4)
      case 3
         continue
      case 1
         for j=1:n
            if (j~=i&j~=ph)
                  X1=G(i,j)*f(i)-B(i,j)*e(i);
                  X2=G(i,j)*e(i)+B(i,j)*f(i);
                  X3=-X2;
                  X4=X1;
                  p=2*i-1;
                  q=2*j-1;
                  J(p,q)=X1;
                  m=p+1;
                  J(m,q)=X3;
                  q=q+1;
                  J(p,q)=X2;
                  J(m,q)=X4;
            else if (j==i&j~=ph)
                     X1=D(i)+G(i,i)*f(i)-B(i,i)*e(i);
                     X2=C(i)+G(i,i)*e(i)+B(i,i)*f(i);
                     X3=C(i)-G(i,i)*e(i)-B(i,i)*f(i);
                     X4=-D(i)+G(i,i)*f(i)-B(i,i)*e(i);
                     q=2*j-1;
                     p=2*i-1;
                     J(p,q)=X1;
                     m=p+1;
                     J(m,q)=X3;
                     q=q+1;
                     J(p,q)=X2;
                     J(m,q)=X4;
                  end
            end
         end
      case 2
         for j=1:n
         if (j~=i&j~=ph)
            X1=G(i,j)*f(i)-B(i,j)*e(i);
            X2=G(i,j)*e(i)+B(i,j)*f(i);
            X3=0;
            X4=0;
            p=2*i-1;
            q=2*j-1;
            J(p,q)=X1;
            m=p+1;
            J(m,q)=X3;
            q=q+1;
            J(p,q)=X2;
            J(m,q)=X4;
         else if j==i&j~=ph
                  X1=D(i)+G(i,i)*f(i)-B(i,i)*e(i);
                  X2=C(i)+G(i,i)*e(i)+B(i,i)*f(i);
                  X3=2*f(i);
                  X4=2*e(i);
                  p=2*i-1;
                  q=2*j-1;
                  J(p,q)=X1;
                  m=p+1;
                  J(m,q)=X3;
                  q=q+1;
                  J(p,q)=X2;
                  J(m,q)=X4;
            end
         end
end
end
if ph~=n
J(no,:)=[];
J(no,:)=[];
J(:,no)=[];
J(:,no)=[];  % 删除不需要的空间
end

fat.taotao · 发表于 2009-5-23 15:45:20

头大啊，呵呵

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

谁能帮我看看这个潮流程序错在哪？

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

相关帖子