Programa de flujos de potencia

%PROGRAMA DE FLUJOS DE POTENCIA PARA LA RED DE 115 KV
clear all;
clc;
sbase=100;
nodos=32;
nlineas=36;
Z=[0.0186+0.0703i 0.0133+0.0502i 0.0020+0.0072i 0.0002+0.0007i 0.0066+0.0249i 0.0441+0.1608i 0.0577+0.2103i 0.0577+0.2103i 0.0285+0.2251i 0.0284+0.1038i 0.0166+0.0598i 0.0190+0.0682i 0.0391+0.1426i 0.0009+0.0031i 0.0010+0.0037i 0.0040+0.0304i 0.0232+0.0489i 0.0284+0.0598i 0.0568+0.1196i0.0665+0.1421i 0.0203+0.0427i 0.0291+0.1061i 0.0282+0.0770i 0.0043+0.0163i 0.0381+0.1041i 0.0481+0.1755i 0.0128+0.0467i 0.0048+0.0104i 0.0107+0.0383i 0.0178+0.0471i 0.0068+0.0249i 0.0125+0.0326i 0.0014+0.0037i 0.0471+0.1718i 0.0144+0.0380i 0.0485+0.1282i];
ZR=[real(Z)]*1.02;
ZXl=[imag(Z)*j];
zlinea=[ZR+ZXl];
btot=[0.0064i 0.0046i 0.0027i 0.0001i 0.0023i 0.0156i 0.0204i 0.0231i 0.0101i 0.0063i0.0072i 0.0138i 0.0003i 0.0004i 0.053i 0.0048i 0.0058i 0.0116i 0.0127i 0.0041i 0.0110i 0.0067i 0.0015i 0.0091i 0.0170i 0.0045i 0.0009i 0.0040i 0.0043i 0.0152i 0.0024i 0.0033i 0.0004i 0.0166i 0.0035i 0.0118i];
v=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
o=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

Pg=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0]/sbase;
Qg=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]/sbase;

k1=0.60*[0 0 19 38 19 19 27.90 8.9063 8.9063 38 9.5 19 14.25 19 9.5 9.5 5.9375 8.9063 19 19 8.9063 0 27.90 8.9063 9.5 19 38 38 14.25 38 19 0];
%K1 es la potencia activa demandada.
k2=0.60*[0 0 6.2449 12.4898 6.2449 6.2449 9.1723 2.9273 2.9273 12.4898 3.1225 6.2449 4.6837 6.2449 3.12253.1225 1.9515 2.9273 6.2449 6.2449 2.9273 0 9.1723 2.9273 3.1225 6.2449 12.4899 12.4899 4.6837 12.4899 6.2449 0];
%K2 es la potencia reactiva demandada.

Pd=[0.0*k1+k1]/sbase;
Qd=[0.0*k2+k2]/sbase;
Ybus=zeros(nodos,nodos);
send=[3 2 7 2 3 3 5 7 5 8 10 7 11 13 32 15 15 16 18 19 21 27 28 22 12 19 19 27 23 24 29 24 24 2 30 2];
rece=[2 4 2 14 4 5 6 5 8 9 9 12 7 12 1 16 17 18 22 18 22 22 22 2322 21 20 28 24 26 24 30 25 29 31 31];
tipo=[1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3];
tol=1e-12;
error=1;
Qmin=-100/sbase;
Qmax=140/sbase;
ntransf=2;
sendt=[1 32];
recet=[2 22];
ztransf=[0+0.02083i 0+0.02083i];
t=[1 1];
it=1;
ncapacit=0;

for(i=1:1:nlineas)
Ybus(send(i),send(i))=(1/zlinea(i))+Ybus(send(i),send(i))+(btot(i)/2);Ybus(rece(i),rece(i))=(1/zlinea(i))+Ybus(rece(i),rece(i))+(btot(i)/2);
Ybus(send(i),rece(i))=-1/zlinea(i);
Ybus(rece(i),send(i))=-1/zlinea(i);
end
%Modificacion de Ybus con transformadores
for i=1:1:ntransf
A(i)=1/(ztransf(i)*t(i));
B(i)=A(i)*(1/t(i)-1);
C(i)=A(i)*(t(i)-1);
Ybus(sendt(i),sendt(i))= Ybus(sendt(i),sendt(i))+C(i)+A(i);
Ybus(recet(i),recet(i))=Ybus(recet(i),recet(i))+A(i)+B(i);
Ybus(sendt(i),recet(i))=-A(i);
Ybus(recet(i),sendt(i))= -A(i);
end

% Nueva Modificacion de la Ybus incluyendo capacitores en derivacion
for i=1:1:ncapacit
Ybus(nodoshunt(i),nodoshunt(i))=Ybus(nodoshunt(i),nodoshunt(i))+shuntB(i);
end

YbusR=Ybus;
v2=v;
o2=o;
for ik=1:1:nlineas+1
Ybus=YbusR;
v=v2;
o=o2;
error=1;
it=0;if(ik==nlineas+1)
Ybus=YbusR;
else
Ybus(send(ik),send(ik))=Ybus(send(ik),send(ik))-(1/zlinea(ik));
Ybus(rece(ik),rece(ik))= Ybus(rece(ik),rece(ik))-(1/zlinea(ik));
Ybus(send(ik),rece(ik))=0.0+0.0j;
Ybus(rece(ik),send(ik))=0.0+0.0j;
end
while (error>tol)
for(i=1:1:nodos)
Paux=0;
Qaux=0;
for(j=1:1:nodos)Paux=Paux+v(i)*v(j)*abs(Ybus(i,j))*cos(o(i)-o(j)-angle(Ybus(i,j)));
Qaux=Qaux+v(i)*v(j)*abs(Ybus(i,j))*sin(o(i)-o(j)-angle(Ybus(i,j)));
end
P(i)=Paux;
Q(i)=Qaux;
dP(i)=Pg(i)-Pd(i)-P(i);
dQ(i)=Qg(i)-Qd(i)-Q(i);
if(tipo(i)==1)
dP(i)=0;
dQ(i)=0;
end...