谁能帮我看看这个潮流程序错在哪？
我自己编了个潮流的程序，现在只能在书上一个特定的算例上能够收敛，如果换成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);
=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);
=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;
=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);
=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=;
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=;
t
result1
S
b
c
result2
function =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 =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=;
end
k=0;
=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 头大啊，呵呵
页:
[1]
