PROYECTO SISTEMAS ELECTRICOS DE POTENCIA
Yulieth Beltran, Samuel Rojas
Objetivo general:
• Desarrollar una metodología para el análisis de flujo de carga, considerando las incertidumbres
presentes tanto en la generación como en la demanda.
Objetivo Especifico:
• Realizar una primera aproximación al problema de flujo de carga mediante la selección, formulación y
resolución eficiente de un sistema IEEE de prueba utilizando ecuaciones no lineales en Matlab®.
El propósito de este proyecto es abordar una primera aproximación al problema de Flujo de Carga en sistemas
de potencia, fundamental para la operación y planificación de redes eléctricas. Mediante la selección de un
sistema de prueba IEEE, se buscará establecer y resolver el sistema de ecuaciones no lineales que gobierna
las variables de estado del sistema eléctrico. La resolución de este sistema mediante las herramientas
de Matlab® permitirá analizar el comportamiento del sistema bajo diferentes condiciones de operación,
proporcionando una base sólida para futuras aplicaciones y estudios más avanzados en el campo de los
sistemas eléctricos.
METODOLOGIA
Para este proyecto se escogio un sistema de nueve nodos haciendo uso de las datos proporcionados por el
sistema, realizando los respectivos calculos que se requieren para esta primera entrega.
1
clc
clear all
close all
format long
Datos del sistema
Tensiones nodales
En la tabla 1 se muestran las tensiones nodales del sistema en magnitud con su respectivo angulo:
Tabla 1. Tensiones nodales
Sb=100e6;
V1=1.04*exp(i * 0 * pi / 180);
V2=1.03*exp(i * 9.28 * pi / 180);
2
V3=1.03*exp(i * 4.66 * pi / 180);
V4=1.03*exp(i * -2.22 * pi / 180);
V5=1*exp(i * -3.99 * pi / 180);
V6=1.01*exp(i * -3.69 * pi / 180);
V7=1.03*exp(i * 3.72 * pi / 180);
V8=1.02*exp(i * 0.73 * pi / 180);
V9=1.03*exp(i * 1.97 * pi / 180);
1.1 Linea de trasmision
En la tabla 1.1 se muestran los respectivos datos de la linea en este caso los datos 41, 27 y 39 representan los
transformadores del sistema pero estos al interfierir con el resultado final no se tendran encuenta.
Tabla 1.1. Parametros de la linea
%Para la linea 1/4
X14=0.0576;
%Para la linea 2/7
X27=0.0625;
%Para la linea 3/9
X39=0.0586;
%Para la linea 4/5
R45=0.0100;
X45=0.0850;
B45=(0.0880)*i;
%Para la linea 4/6
R46=0.0170;
X46=0.0920;
B46=(0.0790)*i;
%Para la linea 5/7
R57=0.0320;
X57=0.1610;
B57=(0.1530)*i;
3
%Para la linea 6/9
R69= 0.0390;
X69=0.1700;
B69=(0.1790)*i;
%Para la linea 7/8
R78=0.0085;
X78=0.0720;
B78=(0.0745)*i;
%Para la linea 8/9
R89=0.0119;
X89=0.1008;
B89=(0.1045)*i;
Calculo de [Y]b
Con los datos de la tabla 1.1 se realizara el calculo de la matriz Ybus, para este caso las Y1, Y2, Y3 no se
tendran encuenta puesto que afectan al resultado final del sistema.
Y11=inv(i*X14);
Y12=0; Y13=0; Y15=0; Y16=0; Y17=0; Y18=0; Y19=0;
Y14=-(inv(i*X14));
Y22=inv(i*X27);
Y21=0; Y23=0; Y25=0; Y26=0; Y24=0; Y28=0; Y29=0;
Y27=-(inv(i*X27));
Y33=inv(i*X39);
Y31=0; Y32=0; Y34=0; Y35=0; Y36=0; Y37=0; Y38=0;
Y39=-(inv(i*X39));
Y44=inv(i*X14)+inv(R45+i*X45)+inv(R46+i*X46)+B45+B46;
Y42=0; Y43=0; Y47=0; Y48=0; Y49=0;
Y41=-(inv(i*X14));
Y45=-(inv(R45+i*X45));
Y46=-(inv(R46+i*X46));
Y55=inv(R45+i*X45)+inv(R57+i*X57)+B45+B57;
Y51=0; Y52=0; Y53=0; Y56=0; Y58=0; Y59=0;
Y54=-(inv(R45+i*X45));
Y57=-(inv(R57+i*X57));
Y66=inv(R46+i*X46)+inv(R69+i*X69)+B46+B69;
Y61=0; Y62=0; Y63=0; Y65=0; Y67=0; Y68=0;
Y64=-(inv(R46+i*X46));
Y69=-(inv(R69+i*X69));
Y77=inv(i*X27)+inv(R57+i*X57)+inv(R78+i*X78)+B57+B78;
Y71=0; Y73=0; Y74=0; Y76=0; Y79=0;
4
Y72=-(inv(i*X27));
Y75=-(inv(R57+i*X57));
Y78=-(inv(R78+i*X78));
Y88=inv(R78+i*X78)+inv(R89+i*X89)+B78+B89;
Y81=0; Y82=0; Y83=0; Y84=0; Y85=0; Y86=0;
Y87=-(inv(R78+i*X78));
Y89=-(inv(R89+i*X89));
Y99=inv(i*X39)+inv(R69+i*X69)+inv(R89+i*X89)+B69+B89;
Y91=0; Y92=0; Y94=0; Y95=0; Y97=0;
Y93=-(inv(i*X39));
Y96=-(inv(R69+i*X69));
Y98=-(inv(R89+i*X89));
Matriz Ybus
Se arma la matriz 6x6 puesto que se eliminaron los valores de los transformadores que se encuntran el los bus
1, 2 y 3, esta matriz se empleara para todos los calculos posteriores.
Ybus=[Y11 Y12 Y13 Y14 Y15 Y16 Y17 Y18 Y19; Y21 Y22 Y23 Y24 Y25 Y26 Y27 Y28 Y29; Y31
Y32 Y33 Y34 Y35 Y36 Y37 Y38 Y39; Y41 Y42 Y43 Y44 Y45 Y46 Y47 Y48 Y49; Y51 Y52 Y53
Y54 Y55 Y56 Y57 Y58 Y59; Y61 Y62 Y63 Y64 Y65 Y66 Y67 Y68 Y69; Y71 Y72 Y73 Y74 Y75
Y76 Y77 Y78 Y79; Y81 Y82 Y83 Y84 Y85 Y86 Y87 Y88 Y89; Y91 Y92 Y93 Y94 Y95 Y96 Y97
Y98 Y99];
Vector de corrientes
I=Ybus*[V1;V2;V3;V4;V5;V6;V7;V8;V9];
Vector de perdidas de potencia aparente
S=[V1;V2;V3;V4;V5;V6;V7;V8;V9].*conj(I);
V_=[V1;V2;V3;V4;V5;V6;V7;V8;V9];
Stotal=sum(S);
Una vez obtenido la Ybus, debemos definir del sistema que son incognitas y que son valores reales para este
caso tenemos un nodo slack el cual es fijo, los generadores con su respectivo valor de tension y las cargas de
los nodos 5, 6 y 8 con el valor de su potencia activa como reactiva.
%parametros de entrada
V2S=1.03;
O2=9.28*pi/180;
V1PV=1.04;
V3PV=1.03;
P1=0.72;
P3=0.85;
P4=0;
Q4=0;
P7=0;
5
Q7=0;
P9=0;
Q9=0;
P5=-1.25;
Q5=-0.5;
P6=-0.9;
Q6=-0.3;
P8=-1;
Q8=-0.35;
Luego declaramos las variables que no sabemos y realizamos una serie de ecuaciones en total un total de 18
ya que debemos tener el valor real como imaginario.
%Parametros de ecuaciones no lineales en variable real:
syms P2 Q2 O1 O3 Q1 Q3 V5PQ O5 V6PQ O6 V8PQ O8 V4PQ O4 V7PQ O7 V9PQ O9
Eq1 =real(conj((P1 +
i*Q1)/(V1PV*exp(i*O1))))==real(Ybus(1,1)*(V1PV*exp(i*O1))+Ybus(1,2)*(V2S*exp(i*O2))
+Ybus(1,3)*(V3PV*exp(i*O3))+Ybus(1,4)*(V4PQ*exp(i*O4))+Ybus(1,5)*(V5PQ*exp(i*O5))
+Ybus(1,6)*(V6PQ*exp(i*O6))+Ybus(1,7)*(V7PQ*exp(i*O7))+Ybus(1,8)*(V8PQ*exp(i*O8))
+Ybus(1,9)*(V9PQ*exp(i*O9)));
Eq2=imag(conj((P1+i*Q1)/(V1PV*exp(i*O1))))==imag(Ybus(1,1)*(V1PV*exp(i*O1))
+Ybus(1,2)*(V2S*exp(i*O2))+Ybus(1,3)*(V3PV*exp(i*O3))+Ybus(1,4)*(V4PQ*exp(i*O4))
+Ybus(1,5)*(V5PQ*exp(i*O5))+Ybus(1,6)*(V6PQ*exp(i*O6))+Ybus(1,7)*(V7PQ*exp(i*O7))
+Ybus(1,8)*(V8PQ*exp(i*O8))+Ybus(1,9)*(V9PQ*exp(i*O9)));
Eq3=real(conj((P2+i*Q2)/(V2S*exp(i*O2))))==real(Ybus(2,1)*(V1PV*exp(i*O1))
+Ybus(2,2)*(V2S*exp(i*O2))+Ybus(2,3)*(V3PV*exp(i*O3))+Ybus(2,4)*(V4PQ*exp(i*O4))
+Ybus(2,5)*(V5PQ*exp(i*O5))+Ybus(2,6)*(V6PQ*exp(i*O6))+Ybus(2,7)*(V7PQ*exp(i*O7))
+Ybus(2,8)*(V8PQ*exp(i*O8))+Ybus(2,9)*(V9PQ*exp(i*O9)));
Eq4=imag(conj((P2+i*Q2)/(V2S*exp(i*O2))))==imag(Ybus(2,1)*(V1PV*exp(i*O1))
+Ybus(2,2)*(V2S*exp(i*O2))+Ybus(2,3)*(V3PV*exp(i*O3))+Ybus(2,4)*(V4PQ*exp(i*O4))
+Ybus(2,5)*(V5PQ*exp(i*O5))+Ybus(2,6)*(V6PQ*exp(i*O6))+Ybus(2,7)*(V7PQ*exp(i*O7))
+Ybus(2,8)*(V8PQ*exp(i*O8))+Ybus(2,9)*(V9PQ*exp(i*O9)));
Eq5=real(conj((P3+i*Q3)/(V3PV*exp(i*O3))))==real(Ybus(3,1)*(V1PV*exp(i*O1))
+Ybus(3,2)*(V2S*exp(i*O2))+Ybus(3,3)*(V3PV*exp(i*O3))+Ybus(3,4)*(V4PQ*exp(i*O4))
+Ybus(3,5)*(V5PQ*exp(i*O5))+Ybus(3,6)*(V6PQ*exp(i*O6))+Ybus(3,7)*(V7PQ*exp(i*O7))
+Ybus(3,8)*(V8PQ*exp(i*O8))+Ybus(3,9)*(V9PQ*exp(i*O9)));
Eq6=imag(conj((P3+i*Q3)/(V3PV*exp(i*O3))))==imag(Ybus(3,1)*(V1PV*exp(i*O1))
+Ybus(3,2)*(V2S*exp(i*O2))+Ybus(3,3)*(V3PV*exp(i*O3))+Ybus(3,4)*(V4PQ*exp(i*O4))
+Ybus(3,5)*(V5PQ*exp(i*O5))+Ybus(3,6)*(V6PQ*exp(i*O6))+Ybus(3,7)*(V7PQ*exp(i*O7))
+Ybus(3,8)*(V8PQ*exp(i*O8))+Ybus(3,9)*(V9PQ*exp(i*O9)));
Eq7=real(conj((P4+i*Q4)/(V4PQ*exp(i*O4))))==real(Ybus(4,1)*(V1PV*exp(i*O1))
+Ybus(4,2)*(V2S*exp(i*O2))+Ybus(4,3)*(V3PV*exp(i*O3))+Ybus(4,4)*(V4PQ*exp(i*O4))
+Ybus(4,5)*(V5PQ*exp(i*O5))+Ybus(4,6)*(V6PQ*exp(i*O6))+Ybus(4,7)*(V7PQ*exp(i*O7))
+Ybus(4,8)*(V8PQ*exp(i*O8))+Ybus(4,9)*(V9PQ*exp(i*O9)));
Eq8=imag(conj((P4+i*Q4)/(V4PQ*exp(i*O4))))==imag(Ybus(4,1)*(V1PV*exp(i*O1))
+Ybus(4,2)*(V2S*exp(i*O2))+Ybus(4,3)*(V3PV*exp(i*O3))+Ybus(4,4)*(V4PQ*exp(i*O4))
+Ybus(4,5)*(V5PQ*exp(i*O5))+Ybus(4,6)*(V6PQ*exp(i*O6))+Ybus(4,7)*(V7PQ*exp(i*O7))
+Ybus(4,8)*(V8PQ*exp(i*O8))+Ybus(4,9)*(V9PQ*exp(i*O9)));
Eq9=real(conj((P5+i*Q5)/(V5PQ*exp(i*O5))))==real(Ybus(5,1)*(V1PV*exp(i*O1))
+Ybus(5,2)*(V2S*exp(i*O2))+Ybus(5,3)*(V3PV*exp(i*O3))+Ybus(5,4)*(V4PQ*exp(i*O4))
6
+Ybus(5,5)*(V5PQ*exp(i*O5))+Ybus(5,6)*(V6PQ*exp(i*O6))+Ybus(5,7)*(V7PQ*exp(i*O7))
+Ybus(5,8)*(V8PQ*exp(i*O8))+Ybus(5,9)*(V9PQ*exp(i*O9)));
Eq10=imag(conj((P5+i*Q5)/(V5PQ*exp(i*O5))))==imag(Ybus(5,1)*(V1PV*exp(i*O1))
+Ybus(5,2)*(V2S*exp(i*O2))+Ybus(5,3)*(V3PV*exp(i*O3))+Ybus(5,4)*(V4PQ*exp(i*O4))
+Ybus(5,5)*(V5PQ*exp(i*O5))+Ybus(5,6)*(V6PQ*exp(i*O6))+Ybus(5,7)*(V7PQ*exp(i*O7))
+Ybus(5,8)*(V8PQ*exp(i*O8))+Ybus(5,9)*(V9PQ*exp(i*O9)));
Eq11=real(conj((P6+i*Q6)/(V6PQ*exp(i*O6))))==real(Ybus(6,1)*(V1PV*exp(i*O1))
+Ybus(6,2)*(V2S*exp(i*O2))+Ybus(6,3)*(V3PV*exp(i*O3))+Ybus(6,4)*(V4PQ*exp(i*O4))
+Ybus(6,5)*(V5PQ*exp(i*O5))+Ybus(6,6)*(V6PQ*exp(i*O6))+Ybus(6,7)*(V7PQ*exp(i*O7))
+Ybus(6,8)*(V8PQ*exp(i*O8))+Ybus(6,9)*(V9PQ*exp(i*O9)));
Eq12=imag(conj((P6+i*Q6)/(V6PQ*exp(i*O6))))==imag(Ybus(6,1)*(V1PV*exp(i*O1))
+Ybus(6,2)*(V2S*exp(i*O2))+Ybus(6,3)*(V3PV*exp(i*O3))+Ybus(6,4)*(V4PQ*exp(i*O4))
+Ybus(6,5)*(V5PQ*exp(i*O5))+Ybus(6,6)*(V6PQ*exp(i*O6))+Ybus(6,7)*(V7PQ*exp(i*O7))
+Ybus(6,8)*(V8PQ*exp(i*O8))+Ybus(6,9)*(V9PQ*exp(i*O9)));
Eq13=real(conj((P7+i*Q7)/(V7PQ*exp(i*O7))))==real(Ybus(7,1)*(V1PV*exp(i*O1))
+Ybus(7,2)*(V2S*exp(i*O2))+Ybus(7,3)*(V3PV*exp(i*O3))+Ybus(7,4)*(V4PQ*exp(i*O4))
+Ybus(7,5)*(V5PQ*exp(i*O5))+Ybus(7,6)*(V6PQ*exp(i*O6))+Ybus(7,7)*(V7PQ*exp(i*O7))
+Ybus(7,8)*(V8PQ*exp(i*O8))+Ybus(7,9)*(V9PQ*exp(i*O9)));
Eq14=imag(conj((P7+i*Q7)/(V7PQ*exp(i*O7))))==imag(Ybus(7,1)*(V1PV*exp(i*O1))
+Ybus(7,2)*(V2S*exp(i*O2))+Ybus(7,3)*(V3PV*exp(i*O3))+Ybus(7,4)*(V4PQ*exp(i*O4))
+Ybus(7,5)*(V5PQ*exp(i*O5))+Ybus(7,6)*(V6PQ*exp(i*O6))+Ybus(7,7)*(V7PQ*exp(i*O7))
+Ybus(7,8)*(V8PQ*exp(i*O8))+Ybus(7,9)*(V9PQ*exp(i*O9)));
Eq15=real(conj((P8+i*Q8)/(V8PQ*exp(i*O8))))==real(Ybus(8,1)*(V1PV*exp(i*O1))
+Ybus(8,2)*(V2S*exp(i*O2))+Ybus(8,3)*(V3PV*exp(i*O3))+Ybus(8,4)*(V4PQ*exp(i*O4))
+Ybus(8,5)*(V5PQ*exp(i*O5))+Ybus(8,6)*(V6PQ*exp(i*O6))+Ybus(8,7)*(V7PQ*exp(i*O7))
+Ybus(8,8)*(V8PQ*exp(i*O8))+Ybus(8,9)*(V9PQ*exp(i*O9)));
Eq16=imag(conj((P8+i*Q8)/(V8PQ*exp(i*O8))))==imag(Ybus(8,1)*(V1PV*exp(i*O1))
+Ybus(8,2)*(V2S*exp(i*O2))+Ybus(8,3)*(V3PV*exp(i*O3))+Ybus(8,4)*(V4PQ*exp(i*O4))
+Ybus(8,5)*(V5PQ*exp(i*O5))+Ybus(8,6)*(V6PQ*exp(i*O6))+Ybus(8,7)*(V7PQ*exp(i*O7))
+Ybus(8,8)*(V8PQ*exp(i*O8))+Ybus(8,9)*(V9PQ*exp(i*O9)));
Eq17=real(conj((P9+i*Q9)/(V9PQ*exp(i*O9))))==real(Ybus(9,1)*(V1PV*exp(i*O1))
+Ybus(9,2)*(V2S*exp(i*O2))+Ybus(9,3)*(V3PV*exp(i*O3))+Ybus(9,4)*(V4PQ*exp(i*O4))
+Ybus(9,5)*(V5PQ*exp(i*O5))+Ybus(9,6)*(V6PQ*exp(i*O6))+Ybus(9,7)*(V7PQ*exp(i*O7))
+Ybus(9,8)*(V8PQ*exp(i*O8))+Ybus(9,9)*(V9PQ*exp(i*O9)));
Eq18=imag(conj((P9+i*Q9)/(V9PQ*exp(i*O9))))==imag(Ybus(9,1)*(V1PV*exp(i*O1))
+Ybus(9,2)*(V2S*exp(i*O2))+Ybus(9,3)*(V3PV*exp(i*O3))+Ybus(9,4)*(V4PQ*exp(i*O4))
+Ybus(9,5)*(V5PQ*exp(i*O5))+Ybus(9,6)*(V6PQ*exp(i*O6))+Ybus(9,7)*(V7PQ*exp(i*O7))
+Ybus(9,8)*(V8PQ*exp(i*O8))+Ybus(9,9)*(V9PQ*exp(i*O9)));
Ya declaradas las 18 ecuaciones hacemos uso de la funcion vpasolve para obtner como resultado final las
variables de estado del sistema:
%% Solucion
sol=vpasolve(Eq1,Eq2,Eq3,Eq4,Eq5,Eq6,Eq7,Eq8,Eq9,Eq10,Eq11,Eq12,Eq13,Eq14,Eq15,Eq16,
Eq17,Eq18,P2,Q2,O1,O3,Q1,Q3,V5PQ,O5,V6PQ,O6,V8PQ,O8,V4PQ,O4,V7PQ,O7,V9PQ,O9,
[1;1;0;0;1;1;1;0;1;0;1;0;1;0;1;0;1;0]);
%sol=vpasolve(Eq1,Eq2,Eq3,Eq4,Eq5,Eq6,Eq7,Eq8,Eq9,Eq10,Eq11,Eq12,Eq13,Eq14,Eq15,Eq16
,Eq17,Eq18,P2,Q2,O1,O3,Q1,Q3,V5PQ,O5,V6PQ,O6,V8PQ,O8,V4PQ,O4,V7PQ,O7,V9PQ,O9);