function flujopot_sol6n
clc
z=rand(12,1);
while 1
clc
%Solución del sistema de ecuaciones de la función flujopot
options = optimoptions('fsolve','Display','iter'); % Option to display output
[x,fval,exitflag,output] =fsolve(@flujopot6n,z,options)
if exitflag==1
resultado=funcion(x(1), x(2), x(3), x(4), x(5), x(6),x(7), x(8), x(9), x(10), x(11), x(12));
for i=1:length(x)
fprintf('x(%2d)=%11.7f\n',i, x(i)) end
for i=1:length(x)
fprintf('%2d f%2d(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x(8),x(9),x(10),x(11),x(12))=%20.18f\n',i,i,resultado(i))
end
return;
end
%Ordenación de datos y soluciones
u1=1;%Debe aparecer aquí porque es distinta u1 que en la función flujopot
d1=0;%Idem
z1=x(1);
z2=x(2);
z3=x(3);
z4=x(4);
p2=-0.19;
q2=-0.075;
z5=x(5);
z6=x(6);
p3=-1.7;
q3=-0.59;
u4=0.95;
z8=x(7);
z7=x(8);
p4=-0.15;
z9=x(9);
z10=x(10);
p5=-0.98;
q5=-0.32;
z11=x(11);
z12=x(12);
p6=0;
q6=0;
%Expresión en complejos
U1=u1*cos(d1)+j*u1*sin(d1);
S1=z1+j*z2;
U2=z3*cos(z4)+j*z3*sin(z4);
S2=p2+j*q2;
U3=z5*cos(z6)+j*z5*sin(z6);
S3=p3+j*q3;
U4=u4*cos(z7)+j*u4*sin(z7);
S4=p4+j*z8;
U5=z9*cos(z10)+j*z9*sin(z10);
S5=p5+j*q5;
U6=z11*cos(z12)+j*z11*sin(z12);
S6=p6+j*q6;
%Cálculo del flujo de potencia en la línea 12
Z12=2.5216e-3+j*0.02171+j*0.01917;%en pu (dato del problema)
y12=1/Z12;%no es igual que el elemento 12 de la matriz Ybus
I12=y12*(U1-U2);%Corriente de 1 a 2
I21=y12*(U2-U1);
S12=U1*conj(I12);%flujo de potencia en la línea 12
S21=U2*conj(I21);
SL12=S12+S21;%pérdidas en la línea 12
%Cálculo del flujo de potencia en la línea 23
Z23=3.3393e-3+j*0.02875;%dato del problema
y23=1/Z23;
I23=y23*(U2-U3);%corriente de 2 a 3
I32=y23*(U3-U2);
S23=U2*conj(I23);%flujo de potencia en la línea 23
S32=U3*conj(I32);
SL23=S23+S32;%pérdidas en la línea 23
%Cálculo del flujo de potencia en la línea 34
Z34=0.1235+j*0.396+j*0.03233;
y34=1/Z34;
I34=y34*(U3-U4);%Corriente de 3 a 4
I43=y34*(U4-U3);
S34=U3*conj(I34);%flujo de potencia en la línea 34
S43=U4*conj(I43);
SL34=S34+S43;%pérdidas en la línea 34
%Cálculo del flujo de potencia en la línea 16 (transformador T1)
Z16=j*0.01917;
y16=1/Z16;
I16=y16*(U1-U6);
I61=y16*(U6-U1);
S16=U1*conj(I16);
S61=U6*conj(I61);
%Cálculo del flujo de potencia en la línea 56
Z56=8.9955e-3+j*0.07746;
y60=j*0.1267;
y56=1/Z56;
I56=y56*(U5-U6)+y60*U5;
I65=y56*(U6-U5)+y60*U6;
S56=U5*conj(I56);
S65=U6*conj(I65);
%Interpretación de los resultados de flujo
%Los flujos presentados son S12, S23 y S13
%Se define positivo en la dirección 12 (23 y 13)
%Es decir, si sale positivo en la parte real e imaginaria, el flujo de
%activa y reactiva va de 1 a 2 (o de 2 a 3; o de 1 a 3)
%Si sale positivo en la parte real y negativo en la parte imaginaria, el
%flujo de activa va de 1 a 2, pero el de reactiva de 2 a 1.
z=x;
end
end
function f=flujopot6n(z)
% if nargin<1
% z=rand(1,12);
% end
%Datos
u1=1;%Módulo tensión nudo 1 (slack) en pu
d1=0;%Argumento tensión nudo 1 (slack)
PD2=0.19;%en pu
QD2=0.075;%en pu
PD3=1.7;%en pu
QD3=0.59;%en pu
PD4=0.15;
u4=0.95;
PD5=0.98;
QD5=0.32;
PD6=0;
QD6=0;
PG2=0;
QG2=0;
PG3=0;
QG3=0;
PG4=0;
QG4=0;
PG5=0;
QG5=0;
PG6=0;
QG6=0;
%Cálculo de las potencias inyectadas en cada nudo
p2=PG2-PD2;
q2=QG2-QD2;
p3=PG3-PD3;
q3=QG3-QD3;
p4=PG4-PD4;
p5=PG5-PD5;
q5=QG5-QD5;
p6=PG6-PD6;
q6=QG6-QD6;
%La construcción de la matriz de admitancias Ybus es manual en este ejemplo
Ybus=[1.5031-76*5079*j -1.5031+24.3695*j 0 0 0 52.1739*j;
-1.5031+24.3695*j 5.480-58.6011*j -3.9849+34.3142*j 0 0 0;
0 -3.9849+34.3142*j 4.6111-36.4277*j -0.6262+2.1635*j 0 0;
0 0 -0.6262+2.1635*j 0.6262-2.1605*j 0 0;
0 0 0 0 1.4792-12.6111*j -1.4792+12.7378*j;
52.1739*j 0 0 0 -1.4792+j*12.7378 1.4792-64.7850*j];
%Matriz de módulos de Ybus
y=abs(Ybus);
%Matriz de argumentos de Ybus
g=angle(Ybus);
%Definición del vector solución
%Para el nudo 1 (slack) las incógnitas son z1 y z2
z1=z(1);
z2=z(2);
%Para el nudo 2 (PQ) las incógnitas son z3 y z4
z3=z(3);
z4=z(4);
%Para el nudo 3 (PQ) las incógnitas son z5 y z6
z5=z(5);
z6=z(6);
z7=z(7);
z8=z(8);
z9=z(9);
z10=z(10);
z11=z(11);
z12=z(12);
%Sistema de ecuaciones
%Como son 3 nudos, deben aparecer 6 ecuaciones
f=[u1*(u1*y(1,1)*cos(d1-d1-g(1,1))+z3*y(1,2)*cos(d1-z4-g(1,2))+z5*y(1,3)*cos(d1-z6-g(1,3))+u4*y(1,4)*cos(d1-z7-g(1,4))+z9*y(1,5)*cos(d1-z10-g(1,5))+z11*y(1,6)*cos(d1-z12-g(1,6)))-z1;
z3*(u1*y(2,1)*cos(z4-d1-g(2,1))+z3*y(2,2)*cos(z4-z4-g(2,2))+z5*y(2,3)*cos(z4-z6-g(2,3))+u4*y(2,4)*cos(z4-z7-g(2,4))+z9*y(2,5)*cos(z4-z10-g(2,5))+z11*y(2,6)*cos(z4-z12-g(2,6)))-p2;
z5*(u1*y(3,1)*cos(z6-d1-g(3,1))+z3*y(3,2)*cos(z6-z4-g(3,2))+z5*y(3,3)*cos(z6-z6-g(3,3))+u4*y(3,4)*cos(z6-z7-g(3,4))+z9*y(3,5)*cos(z6-z10-g(3,5))+z11*y(3,6)*cos(z6-z12-g(3,6)))-p3;
u4*(u1*y(4,1)*cos(z7-d1-g(4,1))+z3*y(4,2)*cos(z7-z4-g(4,2))+z5*y(4,3)*cos(z7-z6-g(4,3))+u4*y(4,4)*cos(z7-z7-g(4,4))+z9*y(4,5)*cos(z7-z10-g(4,5))+z11*y(4,6)*cos(z7-z12-g(4,6)))-p4;
z9*(u1*y(5,1)*cos(z10-d1-g(5,1))+z3*y(5,2)*cos(z10-z4-g(5,2))+z5*y(5,3)*cos(z10-z6-g(5,3))+u4*y(5,4)*cos(z10-z7-g(5,4))+z9*y(5,5)*cos(z10-z10-g(5,5))+z11*y(5,6)*cos(z10-z12-g(5,6)))-p5;
z11*(u1*y(6,1)*cos(z12-d1-g(6,1))+z3*y(6,2)*cos(z12-z4-g(6,2))+z5*y(6,3)*cos(z12-z6-g(6,3))+u4*y(6,4)*cos(z12-z7-g(6,4))+z9*y(6,5)*cos(z12-z10-g(6,5))+z11*y(6,6)*cos(z12-z12-g(6,6)))-p6;
u1*(u1*y(1,1)*sin(d1-d1-g(1,1))+z3*y(1,2)*sin(d1-z4-g(1,2))+z5*y(1,3)*sin(d1-z6-g(1,3))+u4*y(1,4)*sin(d1-z7-g(1,4))+z9*y(1,5)*sin(d1-z10-g(1,5))+z11*y(1,6)*sin(d1-z12-g(1,6)))-z2;
z3*(u1*y(2,1)*sin(z4-d1-g(2,1))+z3*y(2,2)*sin(z4-z4-g(2,2))+z5*y(2,3)*sin(z4-z6-g(2,3))+u4*y(2,4)*sin(z4-z7-g(2,4))+z9*y(2,5)*sin(z4-z10-g(2,5))+z11*y(2,6)*sin(z4-z12-g(2,6)))-q2;
z5*(u1*y(3,1)*sin(z6-d1-g(3,1))+z3*y(3,2)*sin(z6-z4-g(3,2))+z5*y(3,3)*sin(z6-z6-g(3,3))+u4*y(3,4)*sin(z6-z7-g(3,4))+z9*y(3,5)*sin(z6-z10-g(3,5))+z11*y(3,6)*sin(z6-z12-g(3,6)))-q3;
u4*(u1*y(4,1)*sin(z7-d1-g(4,1))+z3*y(4,2)*sin(z7-z4-g(4,2))+z5*y(4,3)*sin(z7-z6-g(4,3))+u4*y(4,4)*sin(z7-z7-g(4,4))+z9*y(4,5)*sin(z7-z10-g(4,5))+z11*y(4,6)*sin(z7-z12-g(4,6)))-z8;
z9*(u1*y(5,1)*sin(z10-d1-g(5,1))+z3*y(5,2)*sin(z10-z4-g(5,2))+z5*y(5,3)*sin(z10-z6-g(5,3))+u4*y(5,4)*sin(z10-z7-g(5,4))+z9*y(5,5)*sin(z10-z10-g(5,5))+z11*y(5,6)*sin(z10-z12-g(5,6)))-q5;
z11*(u1*y(6,1)*sin(z12-d1-g(6,1))+z3*y(6,2)*sin(z12-z4-g(6,2))+z5*y(6,3)*sin(z12-z6-g(6,3))+u4*y(6,4)*sin(z12-z7-g(6,4))+z9*y(6,5)*sin(z12-z10-g(6,5))+z11*y(6,6)*sin(z12-z12-g(6,6)))-q6];
end
function f=funcion(z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12)
% if nargin<1
% z=rand(1,12);
% end
%Datos
u1=1;%Módulo tensión nudo 1 (slack) en pu
d1=0;%Argumento tensión nudo 1 (slack)
PD2=0.19;%en pu
QD2=0.075;%en pu
PD3=1.7;%en pu
QD3=0.59;%en pu
PD4=0.15;
u4=0.95;
PD5=0.98;
QD5=0.32;
PD6=0;
QD6=0;
PG2=0;
QG2=0;
PG3=0;
QG3=0;
PG4=0;
QG4=0;
PG5=0;
QG5=0;
PG6=0;
QG6=0;
%Cálculo de las potencias inyectadas en cada nudo
p2=PG2-PD2;
q2=QG2-QD2;
p3=PG3-PD3;
q3=QG3-QD3;
p4=PG4-PD4;
p5=PG5-PD5;
q5=QG5-QD5;
p6=PG6-PD6;
q6=QG6-QD6;
%La construcción de la matriz de admitancias Ybus es manual en este ejemplo
Ybus=[1.5031-76*5079*j -1.5031+24.3695*j 0 0 0 52.1739*j;
-1.5031+24.3695*j 5.480-58.6011*j -3.9849+34.3142*j 0 0 0;
0 -3.9849+34.3142*j 4.6111-36.4277*j -0.6262+2.1635*j 0 0;
0 0 -0.6262+2.1635*j 0.6262-2.1605*j 0 0;
0 0 0 0 1.4792-12.6111*j -1.4792+12.7378*j;
52.1739*j 0 0 0 -1.4792+j*12.7378 1.4792-64.7850*j];
%Matriz de módulos de Ybus
y=abs(Ybus);
%Matriz de argumentos de Ybus
g=angle(Ybus);
%Definición del vector solución
%Para el nudo 1 (slack) las incógnitas son z1 y z2
% z1=z(1);
% z2=z(2);
% %Para el nudo 2 (PQ) las incógnitas son z3 y z4
% z3=z(3);
% z4=z(4);
% %Para el nudo 3 (PQ) las incógnitas son z5 y z6
% z5=z(5);
% z6=z(6);
% z7=z(7);
% z8=z(8);
% z9=z(9);
% z10=z(10);
% z11=z(11);
%z12=z(12);
%Sistema de ecuaciones
%Como son 3 nudos, deben aparecer 6 ecuaciones
f=[u1*(u1*y(1,1)*cos(d1-d1-g(1,1))+z3*y(1,2)*cos(d1-z4-g(1,2))+z5*y(1,3)*cos(d1-z6-g(1,3))+u4*y(1,4)*cos(d1-z7-g(1,4))+z9*y(1,5)*cos(d1-z10-g(1,5))+z11*y(1,6)*cos(d1-z12-g(1,6)))-z1;
z3*(u1*y(2,1)*cos(z4-d1-g(2,1))+z3*y(2,2)*cos(z4-z4-g(2,2))+z5*y(2,3)*cos(z4-z6-g(2,3))+u4*y(2,4)*cos(z4-z7-g(2,4))+z9*y(2,5)*cos(z4-z10-g(2,5))+z11*y(2,6)*cos(z4-z12-g(2,6)))-p2;
z5*(u1*y(3,1)*cos(z6-d1-g(3,1))+z3*y(3,2)*cos(z6-z4-g(3,2))+z5*y(3,3)*cos(z6-z6-g(3,3))+u4*y(3,4)*cos(z6-z7-g(3,4))+z9*y(3,5)*cos(z6-z10-g(3,5))+z11*y(3,6)*cos(z6-z12-g(3,6)))-p3;
u4*(u1*y(4,1)*cos(z7-d1-g(4,1))+z3*y(4,2)*cos(z7-z4-g(4,2))+z5*y(4,3)*cos(z7-z6-g(4,3))+u4*y(4,4)*cos(z7-z7-g(4,4))+z9*y(4,5)*cos(z7-z10-g(4,5))+z11*y(4,6)*cos(z7-z12-g(4,6)))-p4;
z9*(u1*y(5,1)*cos(z10-d1-g(5,1))+z3*y(5,2)*cos(z10-z4-g(5,2))+z5*y(5,3)*cos(z10-z6-g(5,3))+u4*y(5,4)*cos(z10-z7-g(5,4))+z9*y(5,5)*cos(z10-z10-g(5,5))+z11*y(5,6)*cos(z10-z12-g(5,6)))-p5;
z11*(u1*y(6,1)*cos(z12-d1-g(6,1))+z3*y(6,2)*cos(z12-z4-g(6,2))+z5*y(6,3)*cos(z12-z6-g(6,3))+u4*y(6,4)*cos(z12-z7-g(6,4))+z9*y(6,5)*cos(z12-z10-g(6,5))+z11*y(6,6)*cos(z12-z12-g(6,6)))-p6;
u1*(u1*y(1,1)*sin(d1-d1-g(1,1))+z3*y(1,2)*sin(d1-z4-g(1,2))+z5*y(1,3)*sin(d1-z6-g(1,3))+u4*y(1,4)*sin(d1-z7-g(1,4))+z9*y(1,5)*sin(d1-z10-g(1,5))+z11*y(1,6)*sin(d1-z12-g(1,6)))-z2;
z3*(u1*y(2,1)*sin(z4-d1-g(2,1))+z3*y(2,2)*sin(z4-z4-g(2,2))+z5*y(2,3)*sin(z4-z6-g(2,3))+u4*y(2,4)*sin(z4-z7-g(2,4))+z9*y(2,5)*sin(z4-z10-g(2,5))+z11*y(2,6)*sin(z4-z12-g(2,6)))-q2;
z5*(u1*y(3,1)*sin(z6-d1-g(3,1))+z3*y(3,2)*sin(z6-z4-g(3,2))+z5*y(3,3)*sin(z6-z6-g(3,3))+u4*y(3,4)*sin(z6-z7-g(3,4))+z9*y(3,5)*sin(z6-z10-g(3,5))+z11*y(3,6)*sin(z6-z12-g(3,6)))-q3;
u4*(u1*y(4,1)*sin(z7-d1-g(4,1))+z3*y(4,2)*sin(z7-z4-g(4,2))+z5*y(4,3)*sin(z7-z6-g(4,3))+u4*y(4,4)*sin(z7-z7-g(4,4))+z9*y(4,5)*sin(z7-z10-g(4,5))+z11*y(4,6)*sin(z7-z12-g(4,6)))-z8;
z9*(u1*y(5,1)*sin(z10-d1-g(5,1))+z3*y(5,2)*sin(z10-z4-g(5,2))+z5*y(5,3)*sin(z10-z6-g(5,3))+u4*y(5,4)*sin(z10-z7-g(5,4))+z9*y(5,5)*sin(z10-z10-g(5,5))+z11*y(5,6)*sin(z10-z12-g(5,6)))-q5;
z11*(u1*y(6,1)*sin(z12-d1-g(6,1))+z3*y(6,2)*sin(z12-z4-g(6,2))+z5*y(6,3)*sin(z12-z6-g(6,3))+u4*y(6,4)*sin(z12-z7-g(6,4))+z9*y(6,5)*sin(z12-z10-g(6,5))+z11*y(6,6)*sin(z12-z12-g(6,6)))-q6];
end
