Respected Sir, I have tried atleast manytimes....to run NR loadflow with two UPFCs but getting error again and agian. Main problem I got inf while inversing jaccobian.

1 Ansicht (letzte 30 Tage)
dhiman banerjee
dhiman banerjee am 12 Aug. 2025
Bearbeitet: Torsten am 21 Aug. 2025
Ran in:
l1=1;l2=1;m1=1;m2=1;n1=1;n2=1;o1=1;o2=1;r1=1;r2=1;s1=1;s2=1;t1=1;t2=1;u1=1;u2=1;w1=1;w2=1;x1=1;x2=1;y1=1;y2=1;z1=1;z2=1;
L2=1;M2=1;N2=1;O2=1;R2=1;T2=1;X2=1;Z2=1;SS2=1;
for i=2:busno
l2=1;m2=1;n2=1;o2=1;r2=1;t2=1;u2=1;w2=1;s2=1;x2=1;y2=1;z2=1;
S2=1;
for j=2:busno
if j==2
if options.upfcno==2
if (i==options.sbus(1)||i==options.ebus(1))
O2=1;
M2=1;
N2=1;
L2=1;
elseif (i==options.sbus(2)||i==options.ebus(2))
O2=2;
M2=2;
N2=2;
L2=2;
end
end
end
if j==i
A=0;
for k=1:busno
if k~=i
A=A+V(i)*V(k)*Y(i,k)*sin(theta(i,k)-delta(i)+delta(k));
end
end
if options.ds(i)==1
if i==options.sbus(L2)
j1(l1,l2)=A-V(i)*Vp(options.sbus(L2))*Yp(options.sbus(L2))*sin(thetap(options.sbus(L2))-delta(i)+deltap(options.sbus(L2)))+V(i)*Vs(options.sbus(L2))*Ys(options.sbus(L2))*sin(thetas(options.sbus(L2))-delta(i)+deltas(options.sbus(L2)))-V(i)*V(options.ebus(L2))*Ys(options.sbus(L2))*sin(thetas(options.sbus(L2))-delta(i)+delta(options.ebus(L2)));
l2=l2+1;
elseif i==options.ebus(L2)
j1(l1,l2)=A-V(i)*V(options.sbus(L2))*Ys(options.sbus(L2))*sin(thetas(options.sbus(L2))-delta(i)+delta(options.sbus(L2)))-V(i)*Vs(options.sbus(L2))*Ys(options.sbus(L2))*sin(thetas(options.sbus(L2))-delta(i)+deltas(options.sbus(L2)));
l2=l2+1;
end
else
j1(l1,l2)=A;
l2=l2+1;
end
A=2*V(i)*Y(i,i)*cos(theta(i,i));
for k=1:busno
if k~=i
A=A+V(k)*Y(i,k)*cos(theta(i,k)-delta(i)+delta(k));
end
end
if options.ds(i)==1 & i==options.ebus(M2)
j2(m1,m2)=A+2*V(i)*Gs(options.sbus(M2))-V(options.sbus(M2))*Ys(options.sbus(M2))*cos(thetas(options.sbus(M2))-delta(i)+delta(options.sbus(M2)))-Vs(options.sbus(M2))*Ys(options.sbus(M2))*cos(thetas(options.sbus(M2))-delta(i)+deltas(options.sbus(M2)));
m2=m2+1;
elseif options.dp(i)==1
j2(m1,m2)=-V(i)*Yp(options.sbus(M2))*cos(thetap(options.sbus(M2))-delta(i)+deltap(options.sbus(M2)));
m2=m2+1;
elseif d(i)==2
j2(m1,m2)=A;
m2=m2+1;
end
if d(i)==2
A=0;
for k=1:busno
if k~=i
A=A+V(i)*V(k)*Y(i,k)*cos(theta(i,k)-delta(i)+delta(k));
end
end
if options.ds(i)==1
if i==options.sbus(N2)
j3(n1,n2)=A-V(i)*Vp(options.sbus(N2))*Yp(options.sbus(N2))*cos(thetap(options.sbus(N2))-delta(i)+deltap(options.sbus(N2)))+V(i)*Vs(options.sbus(N2))*Ys(options.sbus(N2))*cos(thetas(options.sbus(N2))-delta(i)+deltas(options.sbus(N2)))-V(i)*V(options.ebus(N2))*Ys(options.sbus(N2))*cos(thetas(options.sbus(N2))-delta(i)+delta(options.ebus(N2)));
n2=n2+1;
elseif i==options.ebus(N2)
j3(n1,n2)=A-V(i)*V(options.sbus(N2))*Ys(options.sbus(N2))*cos(thetas(options.sbus(N2))-delta(i)+delta(options.sbus(N2)))-V(i)*Vs(options.sbus(N2))*Ys(options.sbus(N2))*cos(thetas(options.sbus(N2))-delta(i)+deltas(options.sbus(N2)));
n2=n2+1;
end
else
j3(n1,n2)=A;
n2=n2+1;
end
end
if d(i)==2
A=-2*V(i)*Y(i,i)*sin(theta(i,i));
for k=1:busno
if k~=i
A=A-V(k)*Y(i,k)*sin(theta(i,k)-delta(i)+delta(k));
end
end
if options.ds(i)==1& i==options.ebus(O2)
j4(o1,o2)=A-2*V(i)*Bs(options.sbus(O2))+V(options.sbus(O2))*Ys(options.sbus(O2))*sin(thetas(options.sbus(O2))-delta(i)+delta(options.sbus(O2)))+Vs(options.sbus(O2))*Ys(options.sbus(O2))*sin(thetas(options.sbus(O2))-delta(i)+deltas(options.sbus(O2)));
o2=o2+1;
elseif options.dp(i)==1
j4(o1,o2)=V(i)*Yp(options.sbus(O2))*sin(thetap(options.sbus(O2))-delta(i)+deltap(options.sbus(O2)));
o2=o2+1;
elseif d(i)==2
j4(o1,o2)=A;
o2=o2+1;
end
end
if options.dp(i)==1
js5(r1,r2)=+V(i)*Vp(options.sbus(R2))*Yp(options.sbus(R2))*sin(thetap(options.sbus(R2))-delta(i)+deltap(options.sbus(R2)));
r2=r2+1;
if d(i)==2
js8(t1,t2)=V(i)*Vp(options.sbus(R2))*Yp(options.sbus(R2))*cos(thetap(options.sbus(R2))-delta(i)+deltap(options.sbus(R2)));
t2=t2+1;
end
%T2=T2+1;
js12(x1,x2)=V(i)*Yp(options.sbus(R2))*cos(thetap(options.sbus(R2))-deltap(options.sbus(R2))+delta(i))-2*Vp(options.sbus(R2))*Gp(options.sbus(R2));
js17(x1,x2)=-V(i)*Yp(options.sbus(R2))*cos(thetap(options.sbus(R2))-delta(i)+deltap(options.sbus(R2)));
js22(x1,x2)=V(i)*Yp(options.sbus(R2))*sin(thetap(options.sbus(R2))-delta(i)+deltap(options.sbus(R2)));
x2=x2+1;
js13(z1,z2)=Vp(options.sbus(Z2))*V(i)*Yp(options.sbus(Z2))*sin(thetap(options.sbus(Z2))-deltap(options.sbus(Z2))+delta(i));
js18(z1,z2)=V(i)*Vp(options.sbus(Z2))*Yp(options.sbus(Z2))*sin(thetap(options.sbus(Z2))-delta(i)+deltap(options.sbus(Z2)));
js23(z1,z2)=V(i)*Vp(options.sbus(Z2))*Yp(options.sbus(Z2))*cos(thetap(options.sbus(Z2))-delta(i)+deltap(options.sbus(Z2)));
z2=z2+1;
end
if (options.dp(i)==1 && Z2<options.upfcno) && (R2<options.upfcno)
R2=R2+1;
Z2=Z2+1;
end
if options.ds(i)==1 && i==options.sbus(S2)
js6(s1,s2)=-V(i)*Vs(options.sbus(S2))*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js7(s1,s2)=V(i)*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
s2=s2+1;
if d(i)==2
js9(u1,u2)=-V(i)*Vs(options.sbus(S2))*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js10(u1,u2)=-V(i)*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
u2=u2+1;
end
js14(y1,y2)=-Vs(options.sbus(S2))*V(i)*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-deltas(options.sbus(S2))+delta(i))+Vs(options.sbus(S2))*V(options.ebus(S2))*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-deltas(options.sbus(S2))+delta(options.ebus(S2)));
js15(y1,y2)=-V(i)*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-deltas(options.sbus(S2))+delta(i))+V(options.ebus(S2))*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-deltas(options.sbus(S2))+delta(options.ebus(S2)))-2*Vs(options.sbus(S2))*Gs(options.sbus(S2));
js19(y1,y2)=-V(i)*Vs(options.sbus(S2))*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js20(y1,y2)=V(i)*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js24(y1,y2)=-V(i)*Vs(options.sbus(S2))*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js25(y1,y2)=-V(i)*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
y2=y2+1;
end
if options.ds(i)==1 && i==options.sbus(SS2)
js11(w1,w2)=Vs(options.sbus(SS2))*V(i)*Ys(options.sbus(SS2))*sin(thetas(options.sbus(SS2))-deltas(options.sbus(SS2))+delta(i))-Vp(options.sbus(SS2))*V(i)*Yp(options.sbus(SS2))*sin(thetap(options.sbus(SS2))-deltap(options.sbus(SS2))+delta(i));
js16(w1,w2)=V(i)*Vs(options.sbus(SS2))*Ys(options.sbus(SS2))*sin(thetas(options.sbus(SS2))-delta(i)+deltas(options.sbus(SS2)))-V(i)*Vp(options.sbus(SS2))*Yp(options.sbus(SS2))*sin(thetap(options.sbus(SS2))-delta(i)+deltap(options.sbus(SS2)))-V(i)*V(options.ebus(SS2))*Ys(options.sbus(SS2))*sin(thetas(options.sbus(SS2))-delta(i)+delta(options.ebus(SS2)));
js21(w1,w2)=-V(i)*Vp(options.sbus(SS2))*Yp(options.sbus(SS2))*cos(thetap(options.sbus(SS2))-delta(i)+deltap(options.sbus(SS2)))-V(i)*V(options.ebus(SS2))*Ys(options.sbus(SS2))*cos(thetas(options.sbus(SS2))-delta(i)+delta(options.ebus(SS2)))+V(i)*Vs(options.sbus(SS2))*Ys(options.sbus(SS2))*cos(thetas(options.sbus(SS2))-delta(i)+deltas(options.sbus(SS2)));
w2=w2+1;
end
else
if options.ds(j)==1
if i==options.sbus(L2)& j==options.ebus(L2)
j1(l1,l2)=-V(i)*V(j)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j))+V(i)*V(j)*Ys(options.sbus(L2))*sin(thetas(options.sbus(L2))-delta(i)+delta(j));
l2=l2+1;
elseif i==options.ebus(L2)& j==options.sbus(L2)
j1(l1,l2)=-V(i)*V(j)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j))+V(i)*V(j)*Ys(options.sbus(L2))*sin(thetas(options.sbus(L2))-delta(i)+delta(j));
l2=l2+1;
else
j1(l1,l2)=-V(i)*V(j)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j));
l2=l2+1;
end
else
j1(l1,l2)=-V(i)*V(j)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j));
l2=l2+1;
end
if d(j)==2 && options.dp(j)~=1
if options.ds(j)==1
if i==options.sbus(M2) & j==options.ebus(M2)
j2(m1,m2)=V(i)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j))-V(i)*Ys(options.sbus(M2))*cos(thetas(options.sbus(M2))-delta(i)+delta(options.ebus(M2)));
m2=m2+1;
elseif i==options.ebus(M2)& j==options.sbus(M2)
j2(m1,m2)=0;
m2=m2+1;
else
j2(m1,m2)=V(i)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j));
m2=m2+1;
end
else
j2(m1,m2)=V(i)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j));
m2=m2+1;
end
elseif (d(j)==2 && options.dp(j)==1)
j2(m1,m2)=0;
m2=m2+1;
end
if d(i)==2
if options.ds(j)==1
if i==options.sbus(N2) & j==options.ebus(N2)
j3(n1,n2)=-V(i)*V(j)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j))+V(i)*V(options.ebus(N2))*Ys(options.sbus(N2))*cos(thetas(options.sbus(N2))-delta(i)+delta(options.ebus(N2)));
n2=n2+1;
elseif i==options.ebus(N2) & j==options.sbus(N2)
j3(n1,n2)=-V(i)*V(j)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j))+V(i)*V(options.sbus(N2))*Ys(options.sbus(N2))*cos(thetas(options.sbus(N2))-delta(i)+delta(options.sbus(N2)));
n2=n2+1;
else
j3(n1,n2)=-V(i)*V(j)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j));
n2=n2+1;
end
else
j3(n1,n2)=-V(i)*V(j)*Y(i,j)*cos(theta(i,j)-delta(i)+delta(j));
n2=n2+1;
end
end
if d(i)==2
if d(j)==2 && options.dp(j)~=1
if options.ds(j)==1
if i==options.sbus(O2) & j==options.ebus(O2)
j4(o1,o2)=-V(i)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j))+V(i)*Ys(options.sbus(O2))*sin(thetas(options.sbus(O2))-delta(i)+delta(options.ebus(O2)));
o2=o2+1;
elseif i==options.ebus(O2)& j==options.sbus(O2)
j4(o1,o2)=0;
o2=o2+1;
else
j4(o1,o2)=-V(i)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j));;
o2=o2+1;
end
else
j4(o1,o2)=-V(i)*Y(i,j)*sin(theta(i,j)-delta(i)+delta(j));;
o2=o2+1;
end
elseif (d(j)==2 && options.dp(j)==1)
j4(o1,o2)=0;
o2=o2+1;
end
end
if options.dp(j)==1
js5(r1,r2)=0;
r2=r2+1;
if d(i)==2
js8(t1,t2)=0;
t2=t2+1;
end
end
if j==options.sbus(S2)
if options.ds(i)==1 & i==options.ebus(S2)
js6(s1,s2)=V(i)*Vs(options.sbus(S2))*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js7(s1,s2)=-V(i)*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
s2=s2+1;
if d(i)==2
js9(u1,u2)=V(i)*Vs(options.sbus(S2))*Ys(options.sbus(S2))*cos(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
js10(u1,u2)=V(i)*Ys(options.sbus(S2))*sin(thetas(options.sbus(S2))-delta(i)+deltas(options.sbus(S2)));
u2=u2+1;
end
else
js6(s1,s2)=0;
js7(s1,s2)=0;
s2=s2+1;
if d(i)==2
js9(u1,u2)=0;
js10(u1,u2)=0;
u2=u2+1;
end
end
end
if (S2<options.upfcno)
% if options.upfcno==1
% if j==options.sbus(1)
% S2=S2+1;
% end
if options.upfcno==2
if (j==options.sbus(1)|| j==options.sbus(2))
S2=S2+1;
end
end
end
% if (SS2<upfc_num && j==2) && (ds(i)==1 & i==sbus(2))
% SS2=SS2+1;
% end
if (SS2<options.upfcno && j==2)
if options.upfcno==2
if (options.ds(i)==1 & i==options.sbus(2))
SS2=SS2+1;
end
% elseif upfc_num==3
% if (ds(i)==1) & (i==sbus(2)||i==sbus(3))
% SS2=SS2+1;
% end
% elseif upfc_num==4
% if (ds(i)==1) & (i==sbus(2)||i==sbus(3)||i==sbus(4))
% SS2=SS2+1;
% end
end
end
if options.ds(i)==1 & i==options.sbus(SS2)
if options.ds(j)==1 && j==options.ebus(SS2)
js11(w1,w2)=-Vs(options.sbus(SS2))*V(options.ebus(SS2))*Ys(options.sbus(SS2))*sin(thetas(options.sbus(SS2))-deltas(options.sbus(SS2))+delta(options.ebus(SS2)));
js16(w1,w2)=V(i)*V(j)*Ys(options.sbus(SS2))*sin(thetas(options.sbus(SS2))-delta(i)+delta(j));
js21(w1,w2)=V(i)*V(j)*Ys(options.sbus(SS2))*cos(thetas(options.sbus(SS2))-delta(i)+delta(j));
w2=w2+1;
else
js11(w1,w2)=0;
js16(w1,w2)=0;
js21(w1,w2)=0;
w2=w2+1;
end
if options.ds(j)==1 && j==options.ebus(X2)
js12(x1,x2)=Vs(options.sbus(X2))*Ys(options.sbus(X2))*cos(thetas(options.sbus(X2))-deltas(options.sbus(X2))+delta(j));
js17(x1,x2)=-V(i)*Ys(options.sbus(X2))*cos(thetas(options.sbus(X2))-delta(i)+delta(j));
js22(x1,x2)=V(i)*Ys(options.sbus(X2))*sin(thetas(options.sbus(X2))-delta(i)+delta(j));
x2=x2+1;
elseif d(j)==2
js12(x1,x2)=0;
js17(x1,x2)=0;
js22(x1,x2)=0;
x2=x2+1;
end
end
if options.dp(i)==1 && options.dp(j)==1
js13(z1,z2)=0;
js18(z1,z2)=0;
js23(z1,z2)=0;
z2=z2+1;
end
if (options.ds(i)==1 && options.dp(j)==1)
if options.upfcno==2
if i~=options.ebus(1)&& i~=options.ebus(2)
js14(y1,y2)=0;
js15(y1,y2)=0;
js19(y1,y2)=0;
js24(y1,y2)=0;
js25(y1,y2)=0;
y2=y2+1;
end
% elseif upfc_num==3
% if i~=ebus(1)&& i~=ebus(2)&& i~=ebus(3)
% js14(y1,y2)=0;
% js15(y1,y2)=0;
% js19(y1,y2)=0;
% js24(y1,y2)=0;
% js25(y1,y2)=0;
% y2=y2+1;
% end
% elseif upfc_num==4
% if i~=ebus(1)&& i~=ebus(2)&& i~=ebus(3) && i~=ebus(4)
% js14(y1,y2)=0;
% js15(y1,y2)=0;
% js19(y1,y2)=0;
% js24(y1,y2)=0;
% js25(y1,y2)=0;
% y2=y2+1;
% end
end
end
if options.ds(i)==1 && j==busno
S2=S2+1;
end
if (X2<options.upfcno && j==2)
if options.upfcno==2
if (options.ds(i)==1 & i==options.sbus(2))
X2=X2+1;
end
% elseif upfc_num==3
% if (ds(i)==1) & (i==sbus(2)||i==sbus(3))
% X2=X2+1;
% end
% elseif upfc_num==4
% if (ds(i)==1) & (i==sbus(2)||i==sbus(3)||i==sbus(4))
% X2=X2+1;
% end
end
end
end
end
s1=s1+1;r1=r1+1;l1=l1+1;m1=m1+1;
if d(i)==2
t1=t1+1;
u1=u1+1;
n1=n1+1;
o1=o1+1;
% j4
% pause
end
if (options.ds(i)==1)
if options.upfcno==2
if (i==options.sbus(1)||i==options.sbus(2))
x1=x1+1;z1=z1+1;y1=y1+1; w1=w1+1;
end
% elseif upfc_num==3
% if (i==sbus(1)||i==sbus(2)||i==sbus(3))
% x1=x1+1;z1=z1+1;y1=y1+1; w1=w1+1;
% end
% elseif upfc_num==4
% if (i==sbus(1)||i==sbus(2)||i==sbus(3)||i==sbus(4))
% x1=x1+1;z1=z1+1;y1=y1+1; w1=w1+1;
% end
end
end
end
Unrecognized function or variable 'busno'.
  1 Kommentar
John D'Errico
John D'Errico am 21 Aug. 2025
Bearbeitet: John D'Errico am 21 Aug. 2025
It is almost impossible to know what you are doing wrong, or maybe right. The code is unreadable spaghetti. I'm sorry, but it is. Too many variables to count, variables with almost random looking names, that mean nothing to anyone trying to read your code. What does it model? We don't know, because you have given no information as to what it even does. There is not even a single useful comment line in the entire mess of code. (I'll guess it has something to bo with a bus, of some sort. But even at that, it might be the crosstown bus, or something else. We cannot tell form the code and the very little you have said.)
Is the code correct, and maybe you have numerical problems? Numerically singular jacobians are not uncommon. Often it means there may be a problem with your starting values, or possibly it might just be an issue with double preciaion arithmetic not being sufficient. We don't even know what you are doing with that jacobian. Is something diverging? We don't know. Maybe there is just a bug in your code.
I would suggest you need to first find someone willing to work with you in depth. Someone who has knowledge about the relevant area of mathematics. Maybe this is a problem in numerical analysis. I don't even know. You need to explain what the code should do, what it is modeling. You need to spend time documenting what the code does. Or possibly you need to spend some time learning about the numerical methods you are trying to use.
I'm sorry, but undocumented code that does not seem to work can be almost impossible to debug by individuals who don't know what is the purpose of that code. You should understand that I am not trying to be hurtful about this. When you write code as you have done however, it becomes virtually impossible to diagnose problems in it, especially by someone else.

Melden Sie sich an, um zu kommentieren.

Melden Sie sich an, um diese Frage zu beantworten.

Antworten (0)

Kategorien

Mehr zu Newton-Raphson Method finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by