Post on 19-Sep-2019
i
i
Equation Chapter 1 Section 1
Proyecto Fin de Grado
Grado en Ingeniería de Tecnologías Industriales
Flujo de Carga Factorizado en Coordenadas
Cartesianas
Autor: Jose Manuel Ruiz Oltra
Tutores: Catalina Gómez Quiles
Antonio Gómez Expósito
Dep. de Ingeniería Eléctrica
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Sevilla, 2015
iii
iii
Proyecto Fin de Grado
Grado en Ingeniería de Tecnologías Industriales
Flujo de Carga Factorizado en Coordenadas
Cartesianas
Autor:
Jose Manuel Ruiz Oltra
Tutores:
Catalina Gómez Quiles
Profesor Contratado Doctor
Antonio Gómez Expósito
Catedrático de Universidad
Dep. de Ingeniería Eléctrica
Escuela Técnica Superior de Ingeniería
Universidad de Sevilla
Sevilla, 2015
v
v
Proyecto Fin de Grado: Flujo de Carga Factorizado en Coordenadas Cartesianas
Autor: Jose Manuel Ruiz Oltra
Tutores: Catalina Gómez Quiles
Antonio Gómez Expósito
El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros:
Presidente:
Vocales:
Secretario:
Acuerdan otorgarle la calificación de:
Sevilla, 2015
El Secretario del Tribunal
AGRADECIMIENTOS
Llegado a este punto, me gustaría agradecer primero la gran disponibilidad, ayuda y consejo que he recibido
durante todo el curso por parte de mis tutores, Cati y Antonio, sin olvidar todo lo que he aprendido gracias a
ellos. También agradecerles que me hubieran escogido para llevar a cabo todo este trabajo, del cual estoy muy
satisfecho y contento.
Agradecer a todo el personal docente de la Escuela Superior Técnica de Ingeniería de Sevilla toda la atención y
ayuda prestada durante estos cuatro años de aprendizaje, sin ella nada de esto habría sido posible. Me gustaría
no olvidar dar las gracias a los servicios de la escuela que también han aportado su granito de arena, y por último
y no por ello menos importante, gracias a mi familia y amigos por su gran apoyo y compañía.
Jose Manuel Ruiz Oltra
Sevilla, 2015
vii
vii
RESUMEN
En este proyecto se presenta el ‘Flujo de cargas factorizado en coordenadas cartesianas’ (FLFRC) que es un
método distinto al tradicional flujo de cargas resuelto por el método de Newton-Raphson, siendo su predecesor
el ‘Flujo de cargas factorizado’ (FLF), cuya formulación está escrita en coordenadas polares. El FLFRC se trata
de un método completamente nuevo que nunca se había implementado, hasta ahora, extrayéndose de él varias
conclusiones muy interesantes.
En este trabajo se hará un breve recorrido por los métodos convencionales de resolución del problema del flujo
de cargas, se describirán y se mostrarán las principales ventajas del citado FLF y, como núcleo del proyecto y
novedad, se presentará y se analizará detalladamente el FLFRC.
Para ilustrar todo lo anterior se expondrán breves y detallados ejemplos de aplicación directa de los métodos
factorizados en sistemas muy simples, se resumirán de forma esquemática las fases de programación que se han
llevado acabo y se expondrán y analizarán los resultados obtenidos de aplicar el FLFRC a cinco redes de distintos
tamaños, desde una red pequeña de transporte de 14 nudos, hasta la red de transporte del verano de 2008 de
Polonia constituida por 3012 nudos, obteniendo resultados muy interesantes. También se verá un ejemplo de
una posible aplicación muy prometedora del FLFRC, la cual consiste en determinar los nudos o zonas más
críticas de una red y se concluirá analizando parte del trabajo futuro que ofrece la implantación de este método.
ÍNDICE
Agradecimientos vi
Resumen vii
Índice viii
Índice de Tablas ix
Índice de Figuras x
1 Introducción 1
2 Solución convencional del flujo de cargas 2 Método de Newton-Raphson 5 Formulación 6
2.2.1 Formulación en polares 6 2.2.2 Formulación en cartesianas 8
Límites de reactiva en nudos PV 10
3 Solución factorizada en polares 11 Formulación 11 Algoritmo de resolución 13 La importancia del perfil de tensiones inicial 16
4 Solución factorizada en cartesianas 21 El riesgo de los logaritmos de complejos de parte real negativa 21
4.1.1 Identificación del problema 21 4.1.2 Offset como solución 22
Formulación 23 4.2.1 Formulación con offset 23
Algoritmo de resolución con offset 30 Límites de reactiva en nudos PV generalizada 33
5 Entorno de programación 34
6 Resultados 36 Tablas comparativas para distintas redes 36 Cuenca de atracción del FLFRC para el caso ejemplo de dos nudos 39 Comportamiento del FLFRC ante la variación de distintos parámetros 40 Análisis del número de iteraciones usando mapas de colores 43 Ejemplo de aplicación del FLFRC 46
7 Conclusiones 52
8 Trabajo futuro 53
Anexo 54 FLFRC: Formulación sin offset 54 Demostraciones 61 Otras figuras 64
Referencias 70
ix
ix
ÍNDICE DE TABLAS
Tabla 1. Expresiones para los términos del jacobiano en la formulación polar 7
Tabla 2. Expresiones para los términos del jacobiano en la formulación rectangular 9
Tabla 3. Pimer análisis del FLFRC ante casos bases de distintas redes 36
Tabla 4. Iteraciones y tiempo de cálculo para el caso base 37
Tabla 5. Iteraciones y tiempo de cálculo para el límite de máxima cargabilidad 37
Tabla 6. Iteraciones y tiempo de cálculo en el límite de máxima cargabilidad con límites de reactiva 38
Tabla 7. Comparación de cargas máximas de convergencia entre el FLF y el FLFRC 38
ÍNDICE DE FIGURAS
Figura 1. Estados de un sistema eléctrico según su seguridad 1
Figura 2. Algoritmo de resolución del FLF 15
Figura 3. Sistema de 2 nudos (versión 1) 16
Figura 4. Región de atracción con el método de Newton-Raphson 17
Figura 5. Región de atracción con FLF 17
Figura 6. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso factible,
empezando desde 𝑈0 = 𝐾0 = 1 + 𝑗0.0001 18
Figura 7. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso factible,
empezando desde 𝑈0 = 𝐾0 = 0.5 + 𝑗0.5 19
Figura 8. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso infactible,
empezando desde 𝑈0 = 𝐾0 = 1 + 𝑗0.0001 para 𝜆 = 4.2 19
Figura 9. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso infactible,
empezando desde 𝑈0 = 𝐾0 = 0.5 + 𝑗0.5 para 𝜆 = 4.2 20
Figura 10. Sistema de 2 nudos (versión 2) 27
Figura 11. Algoritmo de resolución del FLFRC 32
Figura 12. Nudo fuera de límites de reactiva 33
Figura 13. Fases de la programación 35
Figura 14. Región de atracción con FLFRC 39
Figura 15. Comportamiento de la componente e para un sistema de 14 nudos 40
Figura 16. Comportamiento de la componente f para un sistema de 14 nudos 40
Figura 17. Evolución del residuo para el caso base de una red de 14 nudos 41
Figura 18. Evolución del residuo para el caso base de una red de 3120 nudos 41
Figura 19. Aplicación de distintos perfiles ante la evolución de la carga para la red de 14 nudos 42
Figura 20. Aplicación de distintos perfiles ante la evolución de la carga para la red de 3120 nudos 42
Figura 21. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 14 nudos 43
Figura 22. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 3120 nudos 44
Figura 23. Iteraciones frente a la variación del offset y la carga en la red de 14 nudos 45
Figura 24. Iteraciones frente a la variación del offset y la carga en la red de 3120 nudos 45
Figura 25. Aplicación del FLFRC para determinar nudos críticos (versión 1) 48
Figura 26. Aplicación del FLFRC para determinar nudos críticos (versión 2) 51
Figura 27. Evolución del residuo para el caso base de una red de 118 nudos 64
Figura 28. Evolución del residuo para el caso base de una red de 300 nudos 64
Figura 29. Evolución del residuo para el caso base de una red de 1354 nudos 65
Figura 30. Aplicación de distintos perfiles ante la evolución de la carga para la red de 118 nudos 65
Figura 31. Aplicación de distintos perfiles ante la evolución de la carga para la red de 300 nudos 66
xi
xi
Figura 32. Aplicación de distintos perfiles ante la evolución de la carga para la red de 1354 nudos 66
Figura 33. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 118 nudos 67
Figura 34. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 300 nudos 67
Figura 35. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 1354 nudos 68
Figura 36. Iteraciones frente a la variación del offset y la carga en la red de 118 nudos 68
Figura 37. Iteraciones frente a la variación del offset y la carga en la red de 300 nudos 69
Figura 38. Iteraciones frente a la variación del offset y la carga en la red de 1354 nudos 69
1
1
1 INTRODUCCIÓN
omo se explica en [1] el problema conocido como flujo de cargas (load flow o power flow en lengua
inglesa) consiste en obtener las condiciones de operación en régimen permanente de un sistema de energía
eléctrica. Más concretamente, dados los consumos en cada nudo, y la potencia generada por los
alternadores, se trata de encontrar las tensiones en los nudos y los flujos de potencia por las líneas y
transformadores.
Sin duda alguna, la rutina del flujo de cargas es la más empleada por los ingenieros involucrados en la
explotación y planificación de los sistemas de potencia, bien como aplicación independiente o como subrutina
de aplicaciones más complejas (estabilidad transitoria, colapso de tensiones, problemas de optimización,
simuladores de entrenamiento, etc.).
En la operación diaria, constituye la base del análisis de seguridad del sistema o análisis de contingencias. Esta
herramienta se ejecuta periódicamente para identificar posibles problemas de sobrecargas o tensiones
inaceptables, como consecuencia de la evolución de la carga, o cuando ocurre algún cambio brusco (inesperado
o programado) en la topología de la red. En la planificación, permite simular el estado en que se encontrarán los
distintos escenarios que se están analizando ante una demanda estimada. En la siguiente figura (Figura 1) se
especifica los posibles estados de un sistema eléctrico en función de la seguridad, en el que la aplicación del
flujo de cargas es vital para analizar cada estado, ya sea de manera preventiva o correctiva.
El flujo de cargas consta básicamente de dos etapas: la primera y más decisiva consiste en obtener las tensiones
complejas en todos los nudos eléctricos. Para este propósito no es posible utilizar las herramientas
convencionales de análisis de circuitos lineales, porque las restricciones de contorno no se especifican en
términos de impedancias (cargas) y fuentes de tensión (generadores) sino de potencias, lo cual conduce a un
sistema no lineal de ecuaciones. La segunda etapa consiste simplemente en el cálculo de todas las magnitudes
de interés, como flujos de potencia activa y reactiva, pérdidas, etc., lo cual es inmediato.
Una de las claves del éxito que tuvieron los programas de flujos de carga a comienzos de los setenta fue la
introducción de técnicas numéricas eficientes para la solución de grandes sistemas de ecuaciones lineales. Estas
técnicas, desarrolladas por ingenieros eléctricos específicamente para resolver este problema, marcaron un hito
fundamental en el análisis de los sistemas de potencia.
En este trabajo, además de explicar con detalle las técnicas más utilizadas para el cálculo del flujo de cargas,
incluyendo aspectos complementarios tales como el tratamiento de los límites de reactiva de los generadores, se
explican y se muestran nuevos métodos de aplicación del flujo de cargas, caracterizados por poseer numerosas
e importantes ventajas respecto a las técnicas tradicionalmente utilizadas (ver [2]). Estos métodos son los
denominados ‘Flujo de carga factorizado en coordenadas polares’ (expuestos en la Sección 3) y ‘Flujo de carga
factorizado en coordenadas cartesianas’ (ver Sección 4), siendo este el que se explica con más detalle, debido
a ciertas ventajas y consecuencias importantes que se explicarán posteriormente.
C
Figura 1. Estados de un sistema eléctrico según su seguridad
2
2
Además, como consecuencia del FLFRC, se propondrá y explicará una nueva forma de tratar los límites de
reactiva en generadores, y se introducirá una aplicación de este nuevo método para averiguar los nudos más
críticos en una red dada, en función de su cargabilidad.
Para ilustrar este último método, se hará uso de cinco redes de tamaños muy diferenciados entre sí. Las tres
primeras son redes de prueba del IEEE, disponibles en [3] y son de 14, 118 y 300 nudos. Las dos últimas redes
usadas son redes de transporte reales y recientes, y constan de 1354 y 3120 nudos [4].
2 SOLUCIÓN CONVENCIONAL DEL FLUJO DE
CARGAS
omo se desarrolla en [5], el estado de una red eléctrica de n nudos queda determinado completamente
mediante las tensiones complejas en todos sus nudos. Las leyes de Kirchhoff y los modelos para cada
componente de la red se condensan en las ecuaciones nodales, que en forma compleja se escriben
como:
ℐ = 𝒴𝒰 (1)
ℐ𝑖 = ∑𝒴𝑖𝑗
𝑛
𝑗=1
𝒰𝑗 𝑗 = 1,2,… , 𝑛 (2)
donde 𝒰 es el vector de tensiones nodales, ℐ el vector de intensidades netas inyectadas en los nudos e 𝒴 la matriz
𝑛 × 𝑛 de admitancias de nudos.
Además, en cada nudo debe cumplirse que:
𝒮𝑖 = 𝒮𝐺𝑖 − 𝒮𝐶𝑖 = 𝒰𝑖ℐ𝑖∗ (3)
siendo 𝒮𝑖 la potencia compleja neta inyectada en el nudo 𝑖, obtenida en el caso general como diferencia entre la
potencia generada y la consumida por la carga en dicho punto. La ecuación anterior, aplicada a todos los nudos,
puede escribirse en forma matricial como:
𝒮 = 𝑑𝑖𝑎𝑔(𝒰)ℐ∗ (4)
donde 𝒮 es el vector de potencias complejas nodales y 𝑑𝑖𝑎𝑔(𝒰) denota la matriz diagonal cuyos elementos son
los del vector 𝒰.
Conocida la matriz de admitancias, las expresiones (1) y (4) constituyen un sistema de 2𝑛 ecuaciones
complejas en términos de las 3𝑛 incógnitas complejas contenidas en 𝒮, 𝒰 e ℐ. En teoría, conociendo 𝑛 de dichas
incógnitas podría resolverse el sistema no lineal resultante para obtener las 2𝑛 restantes. En la práctica, las
intensidades complejas nodales nunca son conocidas o especificadas a priori en un sistema de potencia, por lo
que se prefiere eliminarlas sustituyendo (1) en (4). Esto conduce al sistema no lineal de 𝑛 ecuaciones complejas
siguiente:
𝒮 = 𝑑𝑖𝑎𝑔(𝒰)[𝒴𝒰]∗ (5)
C
3
3
Descomponiendo la potencia compleja en su parte real e imaginaria, 𝒮 = 𝑃 + 𝑗𝑄, y utilizando coordenadas
cartesianas para los elementos de la matriz de admitancias, 𝒴 = 𝐺 + 𝑗𝐵, la ecuación anterior se convierte en:
𝑃 + 𝑗𝑄 = 𝑑𝑖𝑎𝑔(𝒰)[𝐺 − 𝑗𝐵]𝒰∗ (6)
𝑃𝑖 + 𝑗𝑄𝑖 = 𝒰𝑖 ∑[𝐺𝑖𝑗 − 𝑗𝐵𝑖𝑗]𝒰𝑗∗
𝑛
𝑗=1
𝑖 = 1,2,… , 𝑛 (7)
Los métodos iterativos más importantes que se describirán posteriormente no pueden trabajar con las
ecuaciones complejas anteriores, porque la presencia de variables conjugadas impide llevar a cabo derivadas en
forma compleja. Es preciso, por tanto, separar dichas ecuaciones en 2𝑛 ecuaciones reales. Habitualmente, las
tensiones se expresan en coordenadas polares, 𝒰 = 𝑉 ∟𝜃, lo que conduce a:
𝑃𝑖 = 𝑉𝑖 ∑𝑉𝑗(𝐺𝑖𝑗 cos 𝜃𝑖𝑗 + 𝐵𝑖𝑗 sin 𝜃𝑖𝑗)
𝑛
𝑗=1
(8)
𝑄𝑖 = 𝑉𝑖 ∑𝑉𝑗(𝐺𝑖𝑗 sin𝜃𝑖𝑗 − 𝐵𝑖𝑗 cos 𝜃𝑖𝑗)
𝑛
𝑗=1
(9)
𝑖 = 1,2,… , 𝑛
mientras que si se utilizan coordenadas cartesianas o rectangulares, 𝒰 = 𝑉𝑟+𝑗𝑉𝑥 , se obtiene:
𝑃𝑖 = 𝑉𝑟𝑖 ∑(𝐺𝑖𝑗𝑉𝑟𝑗 − 𝐵𝑖𝑗𝑉𝑥𝑗)
𝑛
𝑗=1
+ 𝑉𝑥𝑖 ∑(𝐺𝑖𝑗𝑉𝑥𝑗 + 𝐵𝑖𝑗𝑉𝑟𝑗)
𝑛
𝑗=1
(10)
𝑄𝑖 = 𝑉𝑥𝑖 ∑(𝐺𝑖𝑗𝑉𝑟𝑗 − 𝐵𝑖𝑗𝑉𝑥𝑗)
𝑛
𝑗=1
− 𝑉𝑟𝑖 ∑(𝐺𝑖𝑗𝑉𝑥𝑗 + 𝐵𝑖𝑗𝑉𝑟𝑗)
𝑛
𝑗=1
(11)
𝑖 = 1,2,… , 𝑛
En lo sucesivo, salvo indicación contraria, nos referimos a las ecuaciones en coordenadas polares.
Obsérvese que cada nudo aporta dos ecuaciones y cuatro incógnitas, por lo que deben especificarse dos
magnitudes por nudo para que las ecuaciones anteriores puedan resolverse. En función de las condiciones de
contorno impuestas, pueden distinguirse dos tipos principales de nudos:
Nudos de consumo o nudos PQ: nudos donde se conoce el consumo de potencia activa (𝑃𝐶𝑖𝑒𝑠𝑝
) y reactiva
(𝑄𝐶𝑖𝑒𝑠𝑝
), siendo nula la potencia generada (𝑃𝐺𝑖 = 𝑄𝐺𝑖 = 0).
Las restricciones impuestas son, por tanto,
𝑃𝑖𝑒𝑠𝑝
= −𝑃𝐶𝑖𝑒𝑠𝑝
; 𝑄𝑖𝑒𝑠𝑝
= −𝑄𝐶𝑖𝑒𝑠𝑝 (12)
4
4
quedando como incógnitas las dos componentes de la tensión nodal respectiva. La gran mayoría de
nudos de una red, sobre todo en niveles de menor tensión, son de este tipo.
Nudos de generación o nudos PV: nudos donde un generador regula la tensión a un valor especificado
(𝑉𝑖𝑒𝑠𝑝
) e inyecta una potencia activa (𝑃𝐺𝑖𝑒𝑠𝑝
) determinada previamente por consideraciones económicas.
Las restricciones resultantes, que tienen en cuenta el posible consumo local, son:
𝑃𝑖𝑒𝑠𝑝
= 𝑃𝐺𝑖𝑒𝑠𝑝
− 𝑃𝐶𝑖𝑒𝑠𝑝
; 𝑉𝑖 = 𝑉𝑖𝑒𝑠𝑝 (13)
quedando 𝑄𝑖 y 𝜃𝑖 como incógnitas.
Ahora bien, si sólo se considerasen ambos tipos de nudos, todas las potencias activas inyectadas deberían
especificarse de antemano, lo cual es imposible porque las pérdidas en la red, que también deben ser aportadas
por los generadores, no se conocen hasta que se obtienen los flujos de potencia por cada elemento. Es decir, la
potencia activa de al menos un generador no puede ser especificada y debe calcularse al final del proceso.
Afortunadamente, esta incógnita adicional se compensa con el hecho de que, cuando se trabaja con fasores, uno
de los ángulos de fase puede tomarse libremente como origen de fases. Por simplicidad de cálculo, se toma como
origen de fases precisamente el nudo de generación cuya potencia se deja libre. Este nudo, que suele ser un
generador importante con capacidad para regular frecuencia (ver [6]), o un nudo de interconexión con el exterior,
se denomina nudo de referencia, nudo oscilante o, más comúnmente, nudo slack.
Sea 𝑛𝐷 el número de nudos de consumo. Entonces, el número de nudos de generación, sin contar el nudo
slack, será 𝑛𝐺 = 𝑛 − 𝑛𝐷 − 1. Sin pérdida de generalidad, supondremos que los 𝑛𝐷 primeros nudos son de
consumo y que el nudo de referencia es el último. En base a la clasificación de nudos realizada anteriormente,
las ecuaciones que intervienen en el problema del flujo de cargas son las siguientes:
𝑃𝑖𝑒𝑠𝑝
= 𝑉𝑖 ∑𝑉𝑗(𝐺𝑖𝑗 cos 𝜃𝑖𝑗 + 𝐵𝑖𝑗 sin𝜃𝑖𝑗)
𝑛
𝑗=1
(14)
𝑖 = 1,2,… , 𝑛𝐷 + 𝑛𝐺
𝑄𝑖𝑒𝑠𝑝
= 𝑉𝑖 ∑𝑉𝑗(𝐺𝑖𝑗 sin𝜃𝑖𝑗 − 𝐵𝑖𝑗 cos 𝜃𝑖𝑗)
𝑛
𝑗=1
(15)
𝑖 = 1,2,… , 𝑛𝐷
La solución de este problema consiste en encontrar los desfases 𝜃𝑖 , 𝑖 = 1,2,… , 𝑛𝐷 + 𝑛𝐺 , y los módulos de
tensiones 𝑉𝑖 , 𝑖 = 1,2, … , 𝑛𝐷 , que satisfacen las 2𝑛𝐷 + 𝑛𝐺 ecuaciones (14) y (15).
Nótese que fijar la tensión compleja del nudo oscilante, y dejar libre su potencia compleja, implica
simplemente que las dos ecuaciones respectivas no intervienen en el proceso. Dichas ecuaciones servirán
después, una vez resuelto el problema, para hallar precisamente la potencia compleja de dicho nudo.
Del mismo modo, las 𝑛𝐺 ecuaciones (9) excluidas de (15) permitirán calcular la potencia reactiva que necesita
inyectar o absorber cada generador para mantener su tensión al valor especificado. Como la capacidad de un
generador para absorber o generar reactiva está limitada, es necesario comprobar que no se viola ninguno de los
límites, lo cual complica y alarga normalmente el proceso de solución. Más adelante se explicará cómo se aborda
este problema según la metodología adoptada (Sección 2.3).
5
5
Dado que las ecuaciones resultantes son no lineales, su solución debe ser forzosamente iterativa, por lo que es
necesario adoptar unos valores iniciales para las variables del problema. La búsqueda de valores iniciales
adecuados, que hagan converger el proceso iterativo hacia un punto físicamente viable, de entre las muchas
soluciones matemáticamente posibles, no es un problema trivial en el caso general. Afortunadamente, las
características especiales del problema del flujo de cargas, donde sabemos de antemano que las tensiones se
mueven en una banda relativamente pequeña alrededor de su valor nominal, y que los desfases entre nudos
adyacentes se mueven en márgenes estrechos por motivos de estabilidad, hacen que el denominado perfil plano
sea casi siempre la mejor opción para iniciar el proceso iterativo. Dicho perfil consiste en hacer 𝜃𝑖0 = 0 para
todos los nudos y 𝑉𝑖0 = 1 𝑝𝑢 para los nudos de consumo. Si se ha ejecutado previamente un flujo de cargas, y
los cambios en el estado del sistema han sido menores, puede iniciarse el proceso con la solución del caso
anterior. Esto es especialmente útil cuando se analizan distintas perturbaciones partiendo del mismo caso base
(ver [7]). La experiencia demuestra, sin embargo, que utilizar unos valores aparentemente más próximos a la
solución, pero arbitrarios, suele dar peores resultados que el perfil plano.
Una vez resueltas las ecuaciones (14) y (15) es posible calcular cualquier magnitud deseada (los programas
comerciales de flujos de cargas generan una serie de ficheros o salidas gráficas con toda la información que el
usuario especifique). Los flujos de potencia para un único elemento conectado entre los nudos i y j se pueden
obtener de:
𝑃𝑖𝑗 = 𝑉𝑖𝑉𝑗(𝐺𝑖𝑗 cos 𝜃𝑖𝑗 + 𝐵𝑖𝑗 sin𝜃𝑖𝑗) − 𝐺𝑖𝑗𝑉𝑖2 (16)
𝑄𝑖𝑗 = 𝑉𝑖𝑉𝑗(𝐺𝑖𝑗 sin𝜃𝑖𝑗 − 𝐵𝑖𝑗 cos 𝜃𝑖𝑗) + 𝑉𝑖2(𝐵𝑖𝑗 − 𝑏𝑝𝑖𝑗) (17)
donde 𝑏𝑝 denota la suceptancia paralelo del modelo en π visto en [8]. Cuando existen varios elementos en
paralelo de distintas características, debe recurrirse directamente a dicho modelo para calcular los flujos
individuales.
Análogamente, las pérdidas totales del sistema pueden calcularse, una vez hallada la potencia del nudo slack,
mediante suma de las inyecciones de todos los nudos, o bien como suma de las pérdidas de cada elemento.
Método de Newton-Raphson
Antes de detallar este método, cabe destacar a su predecesor, denominado ‘Método de Gauss-Seidel’ (ver [9]),
cuya convergencia es bastante pobre. Su única utilidad hoy día es su utilización como forma de generar valores
iniciales para el método de Newton-Rapshon, en aquellos raros casos en los que la convergencia de éste sea
problemática partiendo del perfil plano.
En general, estos métodos consisten en encontrar el vector x que satisface el sistema no lineal
𝑓(𝑥) = 0 (18)
Como se explica en [2], el método de Newton-Rapshon obtiene sucesivamente nuevos valores mediante
aproximaciones de primer orden de las funciones no lineales involucradas. La ecuación (18) puede aproximarse
por su desarrollo en serie alrededor del punto 𝑥𝑘:
𝑓(𝑥) ≅ 𝑓(𝑥𝑘) + 𝐹(𝑥𝑘)(𝑥𝑘+1 − 𝑥𝑘) = 0 (19)
6
6
donde 𝐹 =𝜕𝑓
𝜕𝑥 es el jacobiano de 𝑓(𝑥). Partiendo del valor inicial 𝑥0 se obtienen correciones ∆𝑥𝑘 resolviendo
el sistema lineal:
−𝐹(𝑥𝑘)∆𝑥𝑘 = 𝑓(𝑥𝑘) (20)
y nuevos valores 𝑥𝑘+1 de:
𝑥𝑘+1 = 𝑥𝑘 + ∆𝑥𝑘 (21)
El proceso iterativo se detiene cuando se cumple que:
𝑚𝑎𝑥𝑖 |𝑓𝑖(𝑥𝑘)| ≤ 𝜀
para un ε suficientemente pequeño (como por ejemplo, 10−5). Para valores 𝑥0 próximos a la solución, el método
de Newton-Raphson converge cuadráticamente (cuando diverge también lo hace cuadráticamente). En el caso
concreto del flujo de cargas, con independencia del tamaño de la red, el número de iteraciones oscila
habitualmente entre 3 y 5 partiendo de perfil plano [10], aunque la actualización de variables de control descrita
en [11] puede aumentar significativamente este número.
A diferencia de otros métodos (ver [9]), que pueden implementarse directamente en forma compleja, la
presencia del conjugado en las expresiones de la potencia obliga a trabajar en forma real cuando se trata de
calcular derivadas.
Formulación
Según la forma en que se expresen las tensiones, se obtiene la formulación polar o rectangular, siendo con
diferencia la primera de ellas la más extendida.
2.2.1 Formulación en polares
Siguiendo la formulación desarrollada en [12] se obtiene que en este caso el vector x, de dimensión 2𝑛𝐷 + 𝑛𝐺,
viene dado por:
𝑥 = [𝜃|𝑉]𝑇 = [𝜃1, 𝜃2, … , 𝜃𝑛−1|𝑉1, 𝑉2, … , 𝑉𝑛𝐷]𝑇 (22)
y las funciones que queremos anular pueden expresarse, para cada nudo, como el residuo o diferencia entre la
potencia especificada y la potencia calculada en el estado actual, es decir:
𝑓(𝑥) = [∆𝑃|∆𝑄]𝑇 = [∆𝑃1, ∆𝑃2, … , ∆𝑃𝑛−1|∆𝑄1, ∆𝑄2, … , ∆𝑄𝑛𝐷]𝑇 (23)
7
7
siendo:
∆𝑃𝑖 = 𝑃𝑖𝑒𝑠𝑝
− 𝑉𝑖 ∑𝑉𝑗(𝐺𝑖𝑗 cos 𝜃𝑖𝑗 + 𝐵𝑖𝑗 sin 𝜃𝑖𝑗)
𝑛
𝑗=1
(24)
𝑖 = 1,2,… , 𝑛 − 1
∆𝑄𝑖 = 𝑄𝑖𝑒𝑠𝑝
− 𝑉𝑖 ∑𝑉𝑗(𝐺𝑖𝑗 sin𝜃𝑖𝑗 − 𝐵𝑖𝑗 cos𝜃𝑖𝑗)
𝑛
𝑗=1
(25)
𝑖 = 1,2,… , 𝑛𝐷
Con esta notación, y dividiendo el jacobiano en submatrices como se ha hecho con los vectores anteriores, la
ecuación (20), aplicada al problema del flujo de cargas, se convierte en ( [13], [14]):
[𝐻 𝑁𝑀 𝐿
]𝑘
[∆𝜃
∆𝑉/𝑉]𝑘
= [∆𝑃∆𝑄
]𝑘
(26)
y la (21)
[𝜃𝑉]𝑘+1
= [𝜃𝑉]𝑘
+ [∆𝜃∆𝑉
]𝑘
(27)
La utilización de ∆𝑉/𝑉 en lugar de ∆𝑉 no afecta numéricamente al algoritmo, pero simplifica los términos
del jacobiano haciéndolo numéricamente más simétrico (estructuralmente ya lo es). Teniendo en cuenta que
−𝜕(𝑓𝑖𝑒𝑠𝑝
− 𝑓𝑖)
𝜕𝑥𝑗=
𝜕𝑓𝑖𝜕𝑥𝑗
(28)
donde 𝑓 es indistintamente 𝑃 o 𝑄 y x se refiere a 𝑉 o 𝜃, los términos del jacobiano se obtienen de las siguientes
definiciones:
𝐻𝑖𝑗 =𝜕𝑃𝑖
𝜕𝜃𝑗 𝑁𝑖𝑗 = 𝑉𝑗
𝜕𝑃𝑖
𝜕𝑉𝑗 𝑀𝑖𝑗 =
𝜕𝑄𝑖
𝜕𝜃𝑗 𝐿𝑖𝑗 = 𝑉𝑗
𝜕𝑄𝑖
𝜕𝑉𝑗 (29)
Para 𝑖 ≠ 𝑗
𝐻𝑖𝑗 = 𝐿𝑖𝑗 = 𝑉𝑖𝑉𝑗(𝐺𝑖𝑗 sin 𝜃𝑖𝑗 − 𝐵𝑖𝑗 cos 𝜃𝑖𝑗)
𝑁𝑖𝑗 = −𝑀𝑖𝑗 = 𝑉𝑖𝑉𝑗(𝐺𝑖𝑗 cos 𝜃𝑖𝑗 + 𝐵𝑖𝑗 sin𝜃𝑖𝑗)
Para 𝑖 = 𝑗
𝐻𝑖𝑖 = −𝑄𝑖 − 𝐵𝑖𝑖𝑉𝑖2 𝐿𝑖𝑖 = 𝑄𝑖 − 𝐵𝑖𝑖𝑉𝑖
2
𝑁𝑖𝑖 = 𝑃𝑖 + 𝐺𝑖𝑖𝑉𝑖2 𝑀𝑖𝑖 = 𝑃𝑖 − 𝐺𝑖𝑖𝑉𝑖
2
Tabla 1. Expresiones para los términos del jacobiano en la formulación polar
8
8
Los valores respectivos se muestran en la tabla anterior (Tabla 1). Obsérvese que entre las expresiones del
jacobiano y las de ∆𝑃 y ∆𝑄 hay muchos términos comunes, lo cual debe aprovecharse para ahorrar esfuerzo de
cálculo. De este modo, el cálculo del jacobiano es casi un subproducto del cálculo de los residuos de potencia.
La solución del flujo de cargas mediante el método de Newton-Rapshon consta entonces de los siguientes
pasos:
1. Inicializar tensiones con el perfil plano o usar la solución de un caso anterior.
2. Calcular [∆𝑃|∆𝑄], así como los términos del jacobiano. Si todas las componentes de este vector son
menores que ε, detener el proceso. En caso contrario, continuar.
3. Obtener [∆𝜃|∆𝑉/𝑉] resolviendo el sistema (26) con las técnicas explicadas en [15].
4. Actualizar [𝜃|𝑉] y volver al paso 2.
Recuérdese que, por cada nudo PV, nos ahorramos una ecuación en el sistema anterior, lo cual constituye la
principal ventaja de la formulación polar.
En muchos casos se consiguen mejoras de convergencia si se trabaja con ∆𝑄/𝑉 en lugar de ∆𝑄, lo cual se
consigue dividiendo simplemente la fila respectiva por 𝑉𝑖. La razón es que, de ese modo, sólo el término
𝑄𝑖𝑒𝑠𝑝
/𝑉𝑖, numéricamente menor que los demás, permanece no lineal en 𝑉𝑖, y a mayor linealidad mayor
convergencia.
Aunque, por claridad de presentación, el jacobiano se ha dividido en submatrices compatibles con las
particiones realizadas en los vectores de tensiones y potencias, en la práctica, las dos ecuaciones de un nudo PQ
y sus dos variables asociadas se ordenan consecutivamente.
2.2.2 Formulación en cartesianas
A diferencia del apartado anterior y como se ve en [16], la dimensión del vector x es ahora 2𝑛 − 2, estando
constituido por los siguientes elementos:
𝑥 = [𝑉𝑟1, 𝑉𝑟2, … , 𝑉𝑟(𝑛−1)|𝑉𝑥1, 𝑉𝑥2, … , 𝑉𝑥(𝑛−1)]𝑇 (30)
El cálculo de los residuos de potencia se realiza en este caso mediante las ecuaciones
∆𝑃𝑖 = 𝑃𝑖𝑒𝑠𝑝
− [𝑉𝑟𝑖 ∑(𝐺𝑖𝑗𝑉𝑟𝑗 − 𝐵𝑖𝑗𝑉𝑥𝑗)
𝑛
𝑗=1
+ 𝑉𝑥𝑖 ∑(𝐺𝑖𝑗𝑉𝑥𝑗 + 𝐵𝑖𝑗𝑉𝑟𝑗)
𝑛
𝑗=1
] (31)
𝑖 = 1,2,… , 𝑛 − 1
∆𝑄𝑖 = 𝑄𝑖𝑒𝑠𝑝
− [𝑉𝑥𝑖 ∑(𝐺𝑖𝑗𝑉𝑟𝑗 − 𝐵𝑖𝑗𝑉𝑥𝑗)
𝑛
𝑗=1
− 𝑉𝑟𝑖 ∑(𝐺𝑖𝑗𝑉𝑥𝑗 + 𝐵𝑖𝑗𝑉𝑟𝑗)]
𝑛
𝑗=1
(32)
𝑖 = 1,2,… , 𝑛𝐷
9
9
Como se observa claramente, el número de ecuaciones es menor que el de incógnitas. La razón es que aún
falta por imponer en los nudos de generación la restricción:
∆𝑉𝑖2 = (𝑉𝑖
𝑒𝑠𝑝)2 − 𝑉𝑟𝑖
2 − 𝑉𝑥𝑖2 = 0 𝑖 = 1,2, … , 𝑛𝐺 (33)
En forma matricial, el flujo de cargas en coordenadas cartesianas queda formulado entonces como:
[𝑆 𝑇𝑈 𝑊𝐶 𝐷
]
𝑘
[∆𝑉𝑟∆𝑉𝑥
]𝑘
= [∆𝑃∆𝑄
∆𝑉2
]
𝑘
(34)
[𝑉𝑟𝑉𝑥
]𝑘+1
= [𝑉𝑟𝑉𝑥
]𝑘
+ [∆𝑉𝑟∆𝑉𝑥
]𝑘
(35)
Para 𝑖 ≠ 𝑗
𝑆𝑖𝑗 = −𝑊𝑖𝑗 = 𝐺𝑖𝑗𝑉𝑟𝑖 + 𝐵𝑖𝑗𝑉𝑥𝑖
𝑇𝑖𝑗 = 𝑈𝑖𝑗 = 𝐺𝑖𝑗𝑉𝑥𝑖 − 𝐵𝑖𝑗𝑉𝑟𝑖
𝐶𝑖𝑗 = 𝐷𝑖𝑗 = 0
Para 𝑖 = 𝑗
𝑆𝑖𝑖 = 𝐼𝑟𝑖 + 𝐺𝑖𝑖𝑉𝑟𝑖 + 𝐵𝑖𝑖𝑉𝑥𝑖 𝑈𝑖𝑖 = −𝐼𝑥𝑖 − 𝐵𝑖𝑖𝑉𝑟𝑖 + 𝐺𝑖𝑖𝑉𝑥𝑖
𝑊𝑖𝑖 = 𝐼𝑟𝑖 − 𝐺𝑖𝑖𝑉𝑟𝑖 − 𝐵𝑖𝑖𝑉𝑥𝑖 𝑇𝑖𝑖 = 𝐼𝑥𝑖 − 𝐵𝑖𝑖𝑉𝑟𝑖 + 𝐺𝑖𝑖𝑉𝑥𝑖
𝐶𝑖𝑖 = 2𝑉𝑟𝑖 𝐷𝑖𝑖 = 2𝑉𝑥𝑖
Tabla 2. Expresiones para los términos del jacobiano en la formulación rectangular
siendo los elementos del jacobiano los recogidos en la tabla anterior (Tabla 2). En dicha tabla, 𝐼𝑟𝑖 e 𝐼𝑥𝑖 se refieren
a la parte real e imaginaria, respectivamente, de la corriente neta inyectada en el nudo i, que se calcula de la
expresión:
𝐼𝑟𝑖 + 𝑗𝐼𝑥𝑖 = ∑(𝐺𝑖𝑗 + 𝑗𝐵𝑖𝑗)(𝑉𝑟𝑗 + 𝑗𝑉𝑥𝑗)
𝑛
𝑗=1
(36)
El proceso iterativo consta de los mismos pasos descritos para la formulación polar (Sección 2.2.1), cambiando
simplemente las expresiones utilizadas.
10
10
Límites de reactiva en nudos PV
Como se desarrolla en [11], la metodología descrita hasta ahora (Secciones 2.2.1 y 2.2.2) sólo tiene en cuenta
las restricciones de igualdad impuestas por los nudos de consumo y generación. En la práctica, sin embargo, las
soluciones suministradas por un programa industrial de cálculo de flujos de carga deben tener en cuenta
restricciones de contorno adicionales. Estas restricciones pueden ser de desiguladad (límites impuestos a ciertas
variables reguladas y de control) o de igualdad (valores deseados para ciertas variables independientes). Los más
importantes son:
Límites de reactiva en nudos PV. Se aplican a la reactiva generada o absorbida por los alternadores, o
compensadores, que regulan la tensión de un nudo. Cuando se alcanza uno de estos límites el nudo PV
pasa a ser PQ.
Límites de tensión en nudos PQ. Menos comunes que los anteriores. Si se alcanza alguno de esto límites
el nudo se convierte en PV por la activación de algún elemento regulador (normalmente la toma de un
transformador).
Otras restricciones que se suelen tener en cuenta son: transformadores de regulación, transformadores
desfasadores e intercambio de potencia entre áreas, aunque estas no son muy comunes. En todos los casos, lo
que se pretende es controlar una magnitud basándose en una o varias variables de control.
La forma de llevar a cabo estos ajustes depende del método numérico utilizado. Básicamente, los
procedimientos existentes pueden agruparse en dos categorías [17]:
1. El vector de estado se adapta dinámicamente, ampliando o disminuyendo el número de variables, para
tener en cuenta en cada caso los límites alcanzados.
2. Se utiliza un proceso realimentado de ajuste de la variable de control, para conocerlo con más detalle,
consultar [11].
La implantación de cualquiera de estos mecanismos incrementa notablemente el número de iteraciones, así
como la complejidad del código.
Centrándonos en límites de reactiva en nudos PV, sea cual sea el método utilizado, la reactiva que tiene que
inyectar un generador o compensador para mantener su tensión constante debe calcularse en cada iteración y
compararse con sus límites. Si se viola algún límite, 𝑄𝑙𝑖𝑚, la tensión del nudo regulado no puede mantenerse al
valor 𝑉𝑒𝑠𝑝, pasando a ser éste un nudo de consumo (PQ) con 𝑄𝑒𝑠𝑝 = 𝑄𝑙𝑖𝑚. A partir de ese momento hay que
monitorizar la tensión de este nudo PQ virtual. Si 𝑉𝑘 > 𝑉𝑒𝑠𝑝 y 𝑄𝑒𝑠𝑝 = 𝑄𝑚𝑎𝑥, o si 𝑉𝑘 < 𝑉𝑒𝑠𝑝 y 𝑄𝑒𝑠𝑝 = 𝑄𝑚𝑖𝑛,
entonces dicho nudo vuelve a ser PV (esto puede ocurrir por interacción entre varios nudos PV).
En el método de Gauss-Seidel este mecanismo se lleva a cabo de forma trivial.
La implantación de este proceso tampoco es difícil en el método de Newton-Raphson en coordenadas polares,
puesto que en cada iteración hay que calcular el jacobiano. Basta con incluir ∆𝑄𝑖 en el vector de residuos y ∆𝑉𝑖
en el vector de estado cuando el nudo i pasa de PV a PQ, y excluirlos en caso contrario, actualizando
acordemente las filas y columnas del jacobiano.
11
11
3 SOLUCIÓN FACTORIZADA EN POLARES
aciendo referencia al artículo [18], las herramientas que existen actualmente del flujo de cargas hacen
uso del método iterativo de Newton-Raphson, ya sea en forma polar (Sección 2.2.1) o recurriendo a
coordenadas rectangulares (Sección 2.2.2). Para casos donde la velocidad de la solución es más
importante que la fiabilidad de la misma, se desarrollaron nuevos métodos computacionalmente menos costosos,
como son el ‘Flujo de cargas en continua’ [19] y el ‘Método desacoplado rápido’ [20] (Fast decoupled load
flow method en inglés), siendo éste último el más popular. La principal diferencia entre estos dos métodos y el
método convencional de Newton-Raphson reside en el tratamiento del jacobiano, el cuál se mantiene constante
en las dos metodologías anteriormente citadas y es calculado únicamente al principio.
La principal ventaja del método de Newton-Raphson sobre los métodos anteriormente descritos, es que éste
posee convergencia cuadrática cuando la solución está lo suficientemente cerca. Esto explica por qué el número
de iteraciones suele resultar indepediente del tamaño del sistema cuando se adopta un perfil inicial plano. Sin
embargo, si la suposición inicial no está lo suficientemente cerca de la solución, como ocurre cuando la red está
muy cargada, el método de Newton-Raphson diverge a la misma velocidad (divergencia cuadrática), aunque el
caso que se esté resolviendo se encuentre en zona factible. Es por ello por lo que las herramientas comerciales
suelen permitir al usuario adoptar distintas alternativas como perfil de inicio, basadas en métodos iterativos más
simples ( [19] y [20]).
En los últimos años ha emergido una aproximación muy prometedora para la solución de sistemas
sobredeterminados gracias a la idea de desglosar el sistema no lineal original en dos subsistemas, uno de ellos
lineal, el cual es resuelto en dos pasos [21].
Recientemente ha sido demostrado que la misma idea anterior también puede ser aplicada en el problema del
flujo de cargas (ver [22]), denomidado ‘Flujo de cargas factorizado’ (FLF), cuya convergencia es más fiable
que el método convencional de Newton-Raphson, particularmente en casos cercanos al límite de máxima
cargabilidad del sistema. El FLF implica la solución de dos sistemas de ecuaciones por cada iteración, siendo el
primero un problema lineal de mínima distancia, y el segundo, un paso similar al de Newton-Raphson.
Como novedad de este método (FLF) se demuestra que, con una inicializacíon acertada puede converger
fiablemente a dos tipos de soluciones, una real y otra compleja. La solución real se alcanza cuando el sistema se
encuentra en zona factible, y la solución compleja cuando se trata de un caso infactible. En [18] se analizan y se
prueban distintas estrategias de inicialización, con el objetivo de reducir el número de iteraciones en cuanto sea
posible, ante distintos estados de carga en la red.
Formulación
Una forma compacta de escribir el sistema de ecuaciones del flujo de cargas por el método de Newton-Raphson
es:
ℎ(𝑥) = 𝑝 (37)
donde,
𝑝 es el vector de potencias activa y reactiva especificadas. En la expresión polar contiene la potencia
activa inyectada en todos los nudos, menos en el nudo slack, y la potencia reactiva inyectada en todos
los nudos de consumo (nudos PQ). Luego, el vector 𝑝 viene dado por:
𝑝 = [𝑃1, 𝑃2, … , 𝑃𝑁−1|𝑄1, 𝑄2, … , 𝑄𝑁𝑃𝑄]𝑇 (38)
H
12
12
𝑥 es el vector de estado. En la formulación polar convencional contiene: ángulos de desfase de todos
los nudos, menos el slack, y módulos de las tensiones de todos los nudos de consumo (nudos PQ).
Luego este vector es de la forma:
𝑥 = [𝑉1, 𝑉2, … , 𝑉𝑁𝑃𝑄|𝜃1, 𝜃2, … , 𝜃𝑁−1]
𝑇 (39)
Partiendo de un valor inicial, 𝑥0, el método de Newton-Raphson devuelve una solución resolviendo de forma
iterativa el siguiente sistema:
𝐻𝑘∆𝑥𝑘 = ∆𝑝𝑘 = 𝑝 − ℎ(𝑥𝑘) (40)
donde ∆𝑥𝑘 = 𝑥𝑘+1 − 𝑥𝑘, 𝐻𝑘 es el jacobiano de ℎ(∙), calculado en el punto actual 𝑥𝑘, y ∆𝑝𝑘 es la diferencia de
potencias o vector de residuos.
La solución del FLF recurre a un vector de estado ligeramente modificado, 𝑥, y requiere la introducción de
dos vectores auxiliares, 𝑢 e 𝑦, como se ve en [22]:
El vector de estado 𝑥, compuesto por 𝑛 = 2𝑁 − 1 elementos, en el que los módulos de las tensiones
han sido reemplazados por una expresión equivalente:
𝑥 = [𝛼1, 𝛼2, … , 𝛼𝑁|𝜃1, 𝜃2, … , 𝜃𝑁−1]𝑇 (41)
donde 𝛼𝑖 = ln𝑉𝑖2.
Vector 𝑦:
𝑦 = {𝑈𝑖 , 𝐾𝑖𝑗, 𝐿𝑖𝑗} (42)
donde 𝑈𝑖 = 𝑉𝑖2 y, para cada rama que conecta el nudo 𝑖 con el nudo 𝑗, la pareja de variables 𝐾𝑖𝑗, 𝐿𝑖𝑗
viene dada por:
𝐾𝑖𝑗 = 𝑉𝑖𝑉𝑗𝑐𝑜𝑠𝜃𝑖𝑗 (43)
𝐿𝑖𝑗 = 𝑉𝑖𝑉𝑗𝑠𝑖𝑛𝜃𝑖𝑗 (44)
Este vector 𝑦 está compuesto por 𝑚 variables, donde 𝑚 = 2𝑏 + 𝑁, siendo 𝑏 el número de ramas y
𝑁 el número de nudos.
Vector 𝑢, compuesto también por 𝑚 variables:
𝑢 = {𝛼𝑖, 𝛼𝑖𝑗 , 𝜃𝑖𝑗} (45)
donde
𝛼𝑖𝑗 = 𝛼𝑖 + 𝛼𝑗 (46)
𝜃𝑖𝑗 = 𝜃𝑖 − 𝜃𝑗 (47)
13
13
El FLF desglosa el sistema original (37) en los siguientes problemas más simples [23]:
𝐸𝑦 = 𝑝 (48)
𝑢 = 𝑓(𝑦) (49)
𝐶𝑥 = 𝑢 (50)
donde 𝐸 y 𝐶 son matrices rectangulares de tamaño 𝑛 × 𝑚 y 𝑚 × 𝑛, respectivamente, y el vector 𝑓(∙) está
definido en términos de funciones elementales no lineales, que pueden ser fácilmente invertidas, obteniéndose
así:
𝑦 = 𝑓−1(𝑢) (51)
El sistema aumentado de arriba, puede ser expresado de manera más compacta, eliminando el vector 𝑢:
𝐸𝑦 = 𝑝 (52)
𝐶𝑥 = 𝑓(𝑦) (53)
Resaltar que (52) es un sistema lineal subdeterminado, mientras que (53) es sobredeterminado. Entre las
infinitas soluciones de (52), sólo la que satisface (53) constituye una solución para el problema original no lineal
(37). Cuando el vector auxiliar restante 𝑦 es eliminado, se obtiene el sistema original:
𝐸𝑓−1(𝐶𝑥) = 𝑝 (54)
obteniendo una solución parecida al sistema original de Newton-Raphson:
[𝐸𝐹𝑘−1𝐶]∆𝑥𝑘 = 𝑝 − 𝐸𝑓−1(𝐶𝑥) = ∆𝑝𝑘 (55)
donde 𝐹 es el jacobiano de 𝑓(∙).
Algoritmo de resolución
Como vimos anteriormente, el FLF puede dividirse en dos problemas, uno lineal y otro no lineal. El primero
de ellos es un problema de mínima distancia, y el último es un procedimiento muy similar al de Newton-
Raphson. Las dos fases del procedimiento surgen de la representación factorizada de sistemas no lineales, luego
puede ser formalmente ejecutado como sigue, asumiendo la no singularidad del jacobiano ( [22], [23]):
Paso 0: inicialización
Para un valor inicial, 𝑥0, construir 𝑦0 = 𝑓−1(𝐶𝑥0). A diferencia del método de Newton-Raphson, en el FLF
es posible inicializar directamente con 𝑦𝑘 = 𝑦0, pudiendo resultar ventajoso a la hora de enfrentarse a casos mas
complejos.
Paso 1:
1.1. Calcular 𝜆 resolviendo el sistema lineal,
(𝐸𝐸𝑇)𝜆 = 𝑝 − 𝐸𝑦𝑘 (56)
14
14
1.2. Obtener �̃� de:
�̃�𝑘 = 𝑦𝑘 + 𝐸𝑇𝜆 (57)
1.3. Calcular �̃�𝑘 = 𝑓(�̃�𝑘)
Paso 2:
2.1. Obtener 𝑥𝑘+1 resolviendo:
�̃�𝑥𝑘+1 = 𝐸�̃�−1�̃�𝑘 (58)
donde �̃� = 𝐸�̃�−1𝐶, calculado en �̃�𝑘 = 𝑓(�̃�𝑘).
2.2. Actualizar 𝑦𝑘+1 = 𝑓−1(𝐶𝑥𝑘+1).
2.3. Comprobar la convergencia: si ||𝑝 − 𝐸𝑦𝑘+1||∞ es lo suficientemente pequeño, entonces parar. En
caso contrario, volver al paso 1.
En la figura de la siguiente página (Figura 2) se muestra un diagrama de flujo que recoge de manera mas visual
toda la información anterior.
Cabe destacar que el coste computacional del paso 1 es muy alto. Aún así, este esfuerzo extra es compensado
por el paso 2, siendo este menos costoso que el paso convencional de Newton-Raphson. En general, cada
iteración (que incluye la resolución de los dos pasos anteriores) ha mostrado ser más rápida que el
procedimiento de Newton-Raphson para problemas reales complejos y grandes [22].
Además, se conoce que la matriz del jacobiano de Newton-Raphson se convierte singular cuando se alcanza
el punto crítico o “curva de la nariz” (bifurcación). Después de ese punto, no existen soluciones físicas y el
proceso iterativo de Newton-Raphson deja de funcionar [22].
Resultados mostrados en [22] demuestran que para sistemas muy grandes, el FLF es más robusto en términos
de convergencia que el método tradicional de Newton-Raphson, particularmente en escenarios muy cargados,
en el borde del límite de cargabilidad.
Otra gran ventaja, además novedosa, del FLF es que las soluciones físicamente posibles (caso factible) pueden
ser siempre alcanzadas en el dominio real, mientras que cuando se va superando el límite de máxima cargabilidad
del sistema, empiezan a surgir componentes imaginarias en la solución (ver [18]).
15
15
Figura 2. Algoritmo de resolución del FLF
16
16
La importancia del perfil de tensiones inicial
Como se explicó en el Capítulo 2, al tratarse de un problema iterativo, es necesario empezar suponiendo una
solución (perfil inicial), generalmente en las variables del problema (tensiones y fases). Existen muchas
posibilidades de perfiles iniciales, dando lugar a varias posibilidades como solución. La primera opción posible
es que el proceso converja a una solución real, tras más o menos iteraciones. Otra posibilidad es que el proceso
iterativo converja a soluciones físicamente imposibles o inviables, o incluso que no converja a ninguna solución,
debido a la lejanía entre la suposición inicial y una posible solución del problema. Es por ello por lo que la
elección de un perfil inicial correcto es muy importante, marcando la diferencia entre resolver el problema o no.
La elección del perfil no es trivial, debido a la cantidad de soluciones matemáticas posibles. Afortunadamente,
en la explotación normal del sistema eléctrico, las tensiones se mantienen muy próximas al valor en por unidad
(oscilan en una banda muy pequeña) y los desfases entre nudos adyacentes se mueven en márgenes estrechos.
Por este motivo la elección de un perfil inicial plano es casi siempre la mejor opción para iniciar el proceso
iterativo, y consiste en hacer 𝜃𝑖0 = 0 para todos los nudos y 𝑉𝑖
0 = 1 (pu) para los nudos de consumo (nudos
PQ).
Aunque el perfil plano funciona en la mayoría de los casos, la convergencia nunca está garantizada. Como se
ve en [24], la única forma de saber si un valor inicial es adecuado para obtener una solución, es determinar la
zona de atracción del perfil inicial. Estas zonas pueden ser de tres tipos:
1. El perfil inicial está en el interior de la región de atracción de la solución y el método numérico converge.
2. El perfil inicial está fuera del área de atracción de la solución. El método numérico diverge si se parte
de dicha suposición inicial.
3. Aunque el perfil inicial se encuentre dentro de la región de atracción, el método numérico diverge.
Desgraciadamente no es posible definir analíticamente la región de atracción, por lo que para definirla se usan
métodos numéricos. Una forma de hacerlo es crear una malla en el plano XY, de manera que cada punto sea un
perfil inicial distinto. Por ejemplo, se usará el eje de ordenadas como módulo de tensiones (𝑉𝑖0) y el eje de
abscisas como ángulo de fase (𝜃𝑖0), en grados. Una vez que se tiene la malla con todos los perfiles iniciales, basta
con aplicar en cada punto el método numérico escogido y representar en la misma malla el resultado en términos
de iteraciones para analizar la convergencia.
A continuación, se presenta un ejemplo que pone en práctica lo dicho anteriormente (Figura 3). El sistema de
estudio consta de dos nudos, uno de ellos es el nudo slack, por lo que el problema queda resuelto al obtener la
tensión del otro nudo, siendo ésta de 𝑉2 = 0.9209 − 𝑗0.0913 𝑝𝑢, conociendo que su potencia compleja
demanda es de 𝒮2 = 0.9 + 𝑗0.6 𝑝𝑢. La malla de valores iniciales que se toma es bastante amplia, el eje de
abscisas (ángulo de fases) está comprendido entre -150 y 150 grados, y el eje de ordenadas (módulo de tensiones)
comprende los valores desde 0 a 10 en pu.
Figura 3. Sistema de 2 nudos (versión 1)
Esta idea se ha aplicado a dos métodos numéricos distintos, el primero (Figura 4) es el método convencional
de Newton-Raphson, y el segundo (Figura 5), el FLF. Cabe destacar que la región de atracción depende del
algoritmo de resolución o del método numérico que se aplique.
17
17
Figura 4. Región de atracción con el método de Newton-Raphson
Figura 5. Región de atracción con FLF
Nota: para ambos métodos se ha aplicado la misma tolerancia o criterio de convergencia, 10−5, y las dos
figuras han sido extraídas de [25] .
Los diferentes colores significan el número de iteraciones necesarias para obtener una solución, mientras que
la zona más oscura o marrón significa que el método no converge.
En la figura donde se aplica el método de Newton-Raphson (Figura 4), vemos que la región de atracción es
bastante amplia, siendo el perfil plano (𝑉𝑖0 = 1 + 0𝑗) el área que requiere menos iteraciones. Si la comparamos
con la figura en la que se aplica el FLF (Figura 5) notamos que la región de atracción es enormemente más
amplia que en el caso anterior y que la zona de no convergencia se ha reducido notablemente. También merece
destacar que en este último caso, en general, son pocas las iteraciones que se usan para alcanzar la solución (tres
iteraciones en una zona bastante amplia) por lo que, como conclusión, podríamos decir que el FLF tiene una
región de atracción mucho más amplia que la del método de Newton-Raphson, y que además necesita menos
iteraciones para alcanzar la solución. Esto último lo convierte, de nuevo, en un método más robusto y fiable que
el método de Newton-Raphson.
18
18
Más adelante (Sección 6), se aplicará esta misma idea al mismo sistema, pero haciendo uso del ‘Flujo de
Cargas Factorizado en Coordenadas Cartesianas’ (FLFRC), obteniendo distintas conclusiones muy
interesantes. El método anteriormente nombrado, FLFRC, es un algoritmo de resolución algo distinto que el
anterior FLF, el cual se definirá y se explicará detalladamente en la próxima sección (Sección 4).
Hasta ahora se ha hablado sobre la estrategia de inicialización para el flujo de cargas, concretamente para el
caso base del sistema. Sin embargo, a medida que la cargabilidad del sistema va aumentando y usamos el método
convencional de Newton-Raphson, el perfil inicial plano va perdiendo las ventajas que indicamos anteriormente,
debido a que las tensiones van disminuyendo desde su valor nominal. Incluso cuando nos acercamos al límite
de máxima cargabilidad, si empezamos con perfil plano, es muy posible que no converja a una solución.
Como se dijo, el FLF ofrece más rubostez que el método convencional de Newton-Raphson. Gracias a esta
nueva técnica, adoptar un perfil plano como perfil de inicio sigue siendo una buena opción a la hora de ir
aumentando la carga del sistema, incluso funciona bien si nos acercamos al límite de máxima cargabilidad, tal
y como se demuestra en [18].
Por otra parte, cuando superamos el límite de máxima cargabilidad del sistema, el método de Newton-Raphson
no converge a ninguna solución, a diferencia del FLF que sí converge, obteniendo resultados con componentes
imaginarias. En este caso, si se sigue haciendo uso del perfil inicial plano, se alcanza un número alto de
iteraciones; sin embargo, si se usa un perfil inicial con una cierta componente imaginaria, como se demuestra en
[18], mejora bastante el comportamiento del FLF ante casos infactibles (mas allá del límite de máxima
cargabilidad), reduciendo el número de iteraciones. A continuación se exponen algunos ejemplos recogidos en
[18] sobre lo mencionado anteriormente:
Figura 6. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso factible, empezando
desde 𝑈0 = 𝐾0 = 1 + 𝑗0.0001
En la Figura 6 se ha fijado el nivel de carga en 𝜆 = 4, justo antes del límite de máxima cargabilidad de esta
red (éste es de 𝜆 = 4.03). Se observa que adoptando un perfil inicial plano con una componente imaginaria casi
nula, el número de iteraciones es bastante bajo, a pesar de que nos encontramos casi en el límite de máxima
cargabilidad, lo que demuestra la robustez del método.
En la Figura 7 se analiza el mismo caso anterior, pero esta vez se ha tomado un perfil inicial distinto, en el que
la componente imaginaria es mayor, siendo este 𝑈0 = 𝐾0 = 0.5 + 𝑗0.5.
19
19
Figura 7. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso factible, empezando
desde 𝑈0 = 𝐾0 = 0.5 + 𝑗0.5
En la figura anterior observamos que a pesar de empezar con una cierta componente imaginaria más elevada,
el método converge, aunque en una iteración más que en el caso anterior (Figura 6). A medida que se va
avanzando en las iteraciones, notamos como las componentes imaginarias de los módulos de tensiones se van
anulando poco a poco, ya que nos encontramos ante un caso factible. Cabe destacar que aunque el perfil inicial
es distinto al anterior, se alcanza la misma solución en ambos casos.
Ahora, analizemos de nuevo los dos casos anteriores, pero esta vez cambiando la carga a 𝜆 = 4.2, carga que
supera el límite de máxima cargabilidad para esta red de 14 nudos:
Figura 8. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso infactible, empezando
desde 𝑈0 = 𝐾0 = 1 + 𝑗0.0001 para 𝜆 = 4.2
20
20
Figura 9. Evolución de la tensión en nudos PQ en un sistema de 14 nudos ante caso infactible, empezando
desde 𝑈0 = 𝐾0 = 0.5 + 𝑗0.5 para 𝜆 = 4.2
Observamos que en la Figura 8 se alcanza una solución compleja en 16 iteraciones, notando un
comportamiento oscilatorio, mientras que en la Figura 9 se consigue la misma solución compleja, pero esta vez
en bastantes menos iteraciones, 5. También se puede apreciar como en este caso la componente imaginaria de
cada módulo de las tensiones de los nudos es distinta de cero, ya que nos encontramos ante un caso infactible.
Como conclusión de lo anterior, si se recurre a un inicio complejo, aumenta significativamente el área de
convergencia del FLF ante casos infactibles. En cambio, si se usa este mismo tipo de perfil ante casos factibles,
se empeora un poco la convergencia del método, pudiendo tardar algo más en converger el FLF (normalmente
una iteración más). Para más información acerca de la estrategia de inicialización, consultar [18].
21
21
4 SOLUCIÓN FACTORIZADA EN CARTESIANAS
n este apartado se introducirá, se explicará y analizará detalladamente una variante del FLF visto en el
apartado anterior, denominado ‘Flujo de Cargas Factorizado en Coordenadas Cartesianas’ (FLFRC).
En cuanto a la formulación, la diferencia entre el FLF y FLFRC es que, en el primero, se usan coordenadas
polares para definir el vector tensión de cada nudo, mientras que en el segundo, se hace uso de coordenadas
cartesianas:
𝒰𝑖 = 𝑉𝑖𝑒𝑗𝜃𝑖 = 𝑒𝑖 + 𝑗𝑓𝑖 (59)
Esto implica varios cambios importantes en la formulación del problema y en el método de resolución,
obteniéndose distintas conclusiones respecto a las que se vieron del FLF, las cuales se especificarán en apartados
posteriores.
El riesgo de los logaritmos de complejos de parte real negativa
4.1.1 Identificación del problema
Anteriormente, en el FLF, se vió que una vez pasado el límite de máxima cargabilidad del sistema empiezan
a surgir componentes imaginarias en la solución (𝑉𝑖 y 𝜃𝑖). Sin embargo, llegado a un cierto valor muy alto de la
carga del sistema (ya en caso infactible), el método no converge, dando resultados ‘NaN’ (‘Not a Number’).
El problema anterior aparece a la hora de calcular logaritmos de números complejos cuya parte real es negativa.
Esto implica que pueda no cumplirse la propiedad ln(𝑎 ∗ 𝑏) = ln(𝑎) + ln (𝑏), lo que influye en el cálculo de
la matriz 𝐶, matriz que se mantiene numéricamente y dimensionalmente constante mientras no se cambie la
topología de la red, haciendo que el método no converja a ninguna solución.
A continuación se detalla el porqué cuando la parte real de 𝑎 o 𝑏, o ambas, se hace negativa, la propiedad
ln(𝑎 ∗ 𝑏) = ln(𝑎) + ln (𝑏) puede no cumplirse:
Como se ve en [26], el logaritmo natural de un número complejo 𝑎 es otro número complejo 𝑥 = ln (𝑎) que
sea solución de la ecuación:
𝑎 = 𝑒𝑥 (60)
La ecuación anterior no tiene solución única, de hecho, tiene un número infinito de soluciones. Si 𝑎 = 𝑣 + 𝑗𝑤,
entonces:
ln(𝑎) = ln|𝑎| + 𝑗𝐴𝑛𝑔(𝑎) + 𝑗2𝜋𝑘 (61)
donde |𝑎| es el módulo de 𝑎 (|𝑎| = √𝑣2 + 𝑤2), y 𝐴𝑛𝑔(𝑎) es el argumento principal de 𝑎 (𝐴𝑛𝑔(𝑎) =
tan−1 𝑤
𝑣), que está comprendido entre (−𝜋, 𝜋].
Lo anterior hace que la propiedad ln(𝑎 ∗ 𝑏) = ln(𝑎) + ln (𝑏) pueda no ser correcta, ya que ambos lados de
la igualdad pueden diferir por un entero múltiplo de 𝑗2𝜋. A continuación se exponen algunos ejemplos:
E
22
22
Ejemplos
ln(−𝑗) = −𝜋
2𝑗 ≠
𝜋
2𝑗 + 𝜋𝑗 = ln(𝑗) + ln (−1)
Sea 𝑎 un entero negativo:
ln(𝑎 ∗ 𝑎) = ln|𝑎|2 ≠ ln(𝑎) + ln (𝑎)
ln((−2) ∗ (−2)) = ln((−2)2) = 1.3863 ≠ 1.3863 + 𝑗2𝜋 = ln(−2) + ln (−2)
Sea 𝑎 = −0.5 + 𝑗 y 𝑏 = 2 + 𝑗0.5
ln(𝑎 ∗ 𝑏) = 0.835 + 𝑗2.2794 = ln(𝑎) + ln (𝑏), luego se cumple la propiedad.
Sea 𝑎 = −0.5 + 𝑗 y 𝑏 = −2 + 𝑗0.5
ln(𝑎 ∗ 𝑏) = 0.835 − 𝑗1.3521 ≠ 0.835 + 𝑗4.9311 = ln(𝑎) + ln (𝑏)
Fin de ejemplo
Como conclusión, cuando aparecen partes reales negativas, puede no cumplirse la propiedad ln(𝑎 ∗ 𝑏) =ln(𝑎) + ln (𝑏), por lo que la construcción de la matriz 𝐶 no sería válida, lo que hace que el problema diverja.
4.1.2 Offset como solución
Como solución al problema anterior, se introduce la idea de sumar un offset que se basa en añadir un número
natural positivo 𝑚 a las variables del problema. Es decir, consiste en sumar una constante a todas las 𝑒 y 𝑓 de
cada nudo, obteniendo un vector de tensiones ampliado 𝒰𝑖′:
𝒰𝑖 = 𝑒𝑖 + 𝑗𝑓𝑖
𝒰𝑖′ = 𝑒𝑖
′ + 𝑗𝑓𝑖′ = (𝑒𝑖 + 𝑚) + 𝑗(𝑓𝑖 + 𝑚)
Esta suma se realiza en el momento en el que se construye el perfil inicial de tensiones, y cuando el problema
converja a una solución, se obtendrán las componentes 𝑒𝑖′ y 𝑓𝑖
′ a las que se tendrán que restar el offset introducido
para obtener la solución real:
𝑒𝑖 = 𝑒𝑖′ − 𝑚
𝑓𝑖 = 𝑓𝑖′ − 𝑚
Esta idea permite que a la hora de calcular logaritmos, la parte real nunca sea negativa. Lo anterior implica que
la matriz 𝐶 (que siempre es constante durante el problema, a menos que cambie la topología de la red) sea
siempre válida, al cumplirse siempre la relación ln(𝑎 ∗ 𝑏) = ln(𝑎) + ln (𝑏).
Todo esto permite, una vez pasado el límite de máxima cargabilidad, conseguir la convergencia del método
para cargas mayores que con el FLF, ya que éste último se quedaba en una cierta carga como límite, debido a la
aparición de logaritmos de complejos con parte real negativa. La solución del offset no es extensible al FLF, ya
que se complica bastante la formulación matemática.
23
23
La introducción del offset como solución implica dos problemas:
1. Complica y cambia bastante la formulación y el algoritmo de resolución respecto al FLFRC sin offset,
tal y como se explicará y analizará posteriormente (Sección 4 y Sección 4.3).
2. El offset se trata de un parámetro libre del problema, es decir, se podrá variar como convenga,
dependiendo de la situación. Este offset no afecta al valor numérico de la solución, pero, el problema
radica en que cada red necesita distintos offsets, dependiendo del valor de la carga (una vez pasado el
límite de máxima cargabilidad, cuando aparecen logaritmos de números negativos). Como se verá
posteriormente, para casos factibles e infactibles sin riesgo de logaritmos de números negativos, el valor
del offset no importa (siempre que sea distinto de 0).
Formulación
A continuación, se verá la formulación del FLFRC con la introducción del offset, para ayudar a comprenderla
mejor, se introducirá un ejercicio resuelto paso a paso a modo de ejemplo aplicado a un caso muy sencillo de
una red de dos nudos. Antes de empezar con la formulación anterior, cabe destacar que anterior a ella se implantó
y se estudió otra formulación distinta, sin offset. Esta formulación fue la primera en analizarse, y tenía algunas
ventajas tales como el reducido número de variables y su similitud con el FLF en cuanto a la formulación. Sin
embargo, fue aquí donde nos dimos cuenta del problema citado anteriormente, concretamente se localizó en dos
variables del problema, observando que la parte real del logaritmo complejo se hacía negativa. Esta primera
formulación sin offset está recogida en el apartado Anexo, y al igual que en la formulación con offset, se ha
añadido el mismo ejemplo de ejercicio, con el objetivo de ayudar a comprender mejor esta formulación. Todo
lo anterior dió lugar a la introducción del offset del cuál se habló antes, y cuya formulación se detalla a
continuación.
4.2.1 Formulación con offset
Al introducir el offset en el vector complejo de tensiones, la formulación del FLFRC cambia drásticamente
con respecto a la formulación sin offset recogida en Anexo.
Como ya se dijo, el vector de tensiones 𝒰𝑖 se sustituye por 𝒰𝑖′:
𝒰𝑖 = 𝑒𝑖 + 𝑗𝑓𝑖 → 𝒰𝑖′ = 𝑒𝑖
′ + 𝑗𝑓𝑖′ = (𝑒𝑖 + 𝑚) + 𝑗(𝑓𝑖 + 𝑚)
donde 𝑚 es el valor del offset. Luego, las variables nuevas son ahora 𝑒𝑖′ = 𝑒𝑖 + 𝑚 y 𝑓𝑖
′ = 𝑓𝑖 + 𝑚
Recordemos que la introducción del offset se realiza cuando se construye el perfil inicial de tensiones, y que
cuando el problema converge, basta con obtener [𝑒1′ , 𝑒2
′ , … , 𝑒𝑁′ |𝑓1
′, 𝑓2′, … , 𝑓𝑁
′ ] y deshacer el cambio de la
siguiente forma:
𝑒𝑖 = 𝑒𝑖′ − 𝑚 (62)
𝑓𝑖 = 𝑓𝑖′ − 𝑚 (63)
En este caso, el sistema de ecuaciones a resolver sigue siendo el mismo visto en la Sección 3.1:
𝐸𝑦 = 𝑝 (64)
𝑢 = 𝑓(𝑦) (65)
𝐶𝑥 = 𝑢 (66)
24
24
La introducción del offset hace que vuelva a cambiar el vector de estado 𝑥, y con ello todas las variables del
problema, las cuales se definen de la siguiente forma:
El vector de estado 𝑥, compuesto por 2𝑁 elementos:
𝑥 = [𝛼1, 𝛼2, … , 𝛼𝑁|𝛽1, 𝛽2, … , 𝛽𝑁]𝑇 (67)
donde
𝛼𝑖 = ln 𝑒𝑖′ = ln(𝑒𝑖 + 𝑚) (68)
𝛽𝑖 = ln𝑓𝑖′ = ln(𝑓𝑖 + 𝑚) (69)
Vector 𝑦, cuya dimensión es ahora 4(𝑁 + 𝑏):
𝑦 = {𝑒𝑖′2, 𝑓𝑖
′2, 𝑒𝑖′, 𝑓𝑖
′, 𝑒𝑖′𝑒𝑗
′, 𝑓𝑖′𝑓𝑗
′, 𝑒𝑖′𝑓𝑗
′, 𝑓𝑖′𝑒𝑗
′} (70)
El vector 𝑢, también de dimensión 4(𝑁 + 𝑏), se obtiene de una forma más sencilla de la que se vió en
la formulación anterior:
𝑢 = ln𝑦 = {ln(𝑒𝑖′2) , ln(𝑓𝑖
′2) , ln 𝑒𝑖′ , ln 𝑓𝑖
′ , ln(𝑒𝑖′𝑒𝑗
′) , ln(𝑓𝑖′𝑓𝑗
′) , ln(𝑒𝑖′𝑓𝑗
′) , ln(𝑓𝑖′𝑒𝑗
′)} (71)
volviéndose a cumplir la relación
𝑦 = exp(𝑢) = 𝑓−1(𝑢)
Vector 𝑝, de 2𝑛 − 1 componentes:
𝑝 = [𝑈𝑠𝑙𝑎𝑐𝑘 , 𝑈1, 𝑈2, … , 𝑈𝑁𝑃𝑉|𝑃1, 𝑃2, … , 𝑃𝑁−1|𝑄1, 𝑄2, … , 𝑄𝑁𝑃𝑄
]𝑇 (72)
Sin embargo, a la hora de resolver 𝑝 = 𝐸𝑦, no se usa el vector 𝑝. En su lugar, se recurre a un vector 𝑝
modificado, �̂�, como sigue a continuación :
�̂� = 𝑝 − 𝑝𝐼 (73)
donde 𝑝𝐼 es un vector que contiene los términos independientes que resultan al aplicar las siguientes
ecuaciones:
𝑈𝑖𝑆𝑃 = 𝑈𝑖 (𝑝𝑎𝑟𝑎 𝑛𝑢𝑑𝑜𝑠 𝑃𝑉) (74)
𝑃𝑖𝑆𝑃 = 𝐺𝑖𝑖𝑈𝑖 + ∑[𝐺𝑖𝑗𝑈𝑖𝑗 + 𝐵𝑖𝑗𝑊𝑖𝑗]
𝑁−1
𝑗≠𝑖
(75)
𝑄𝑖𝑆𝑃 = −𝐵𝑖𝑖𝑈𝑖 + ∑[𝐺𝑖𝑗𝑊𝑖𝑗 − 𝐵𝑖𝑗𝑈𝑖𝑗]
𝑁𝑃𝑄
𝑗≠𝑖
(76)
25
25
donde
𝑈𝑖 = 𝑉𝑖2 = 𝑒𝑖
2 + 𝑓𝑖2 (77)
𝑈𝑖𝑗 = 𝑒𝑖𝑒𝑗 + 𝑓𝑖𝑓𝑗 (78)
𝑊𝑖𝑗 = −𝑒𝑖𝑓𝑗 + 𝑓𝑖𝑒𝑗 (79)
Los términos independientes citados anteriormente son aquellos que no dependen de ninguna variable
contenida en el vector 𝑦, definido en (70). Todo esto se verá mas claramente en un ejemplo de aplicación
del FLFRC con offset.
Al cambiar el vector auxiliar 𝑦, la construcción de la matriz 𝐸 también cambia drásticamente, pues esta
expresa la relación del vector 𝑦 con el vector 𝑝. Esta relación viene dada por las ecuaciones (74), (75)
y (76) para las cuales es necesario transformar las variables {𝑈𝑖 , 𝑈𝑖𝑗 ,𝑊𝑖𝑗}, para que estén dadas en
función de 𝑒𝑖′ y 𝑓𝑖
′, ya que originalmente están escritas en función de 𝑒𝑖 y 𝑓𝑖:
Sustituyendo (62) y (63) en (77), (78) y (79), se obtienen las siguientes expresiones en función de las
variables 𝑒𝑖′ y 𝑓𝑖
′ de {𝑈𝑖 , 𝑈𝑖𝑗 ,𝑊𝑖𝑗}:
𝑈𝑖 = 𝑒𝑖′2 + 𝑓𝑖
′2 − 2𝑚𝑒𝑖′ − 2𝑚𝑓𝑖
′ + 2𝑚2 (80)
𝑈𝑖𝑗 = 𝑒𝑖′𝑒𝑗
′ + 𝑓𝑖′𝑓𝑗
′ − 𝑚𝑒𝑖′ − 𝑚𝑒𝑗
′ − 𝑚𝑓𝑖′ − 𝑚𝑓𝑗
′ + 2𝑚2 (81)
𝑊𝑖𝑗 = −𝑒𝑖′𝑓𝑗
′ + 𝑓𝑖′𝑒𝑗
′ + 𝑚𝑒𝑖′ − 𝑚𝑒𝑗
′ − 𝑚𝑓𝑖′ + 𝑚𝑓𝑗
′ (82)
Aplicando las ecuaciones anteriores a (74), (75) y (76) ya es posible construir la nueva matriz 𝐸, de
dimensiones 𝑛 × 4(𝑁 + 𝑏).
Matriz 𝐶, obtenida igual que anteriormente, de dimensiones 4(𝑁 + 𝑏) × 2𝑁 y definida por las
siguientes igualdades:
ln 𝑒𝑖′2 = 2 ln 𝑒𝑖
′ = 2𝛼𝑖
ln 𝑓𝑖′2 = 2 ln 𝑓𝑖
′ = 2𝛽𝑖
ln 𝑒𝑖′ = 𝛼𝑖
ln 𝑓𝑖′ = 𝛽𝑖
ln 𝑒𝑖′𝑒𝑗
′ = ln 𝑒𝑖′ + ln 𝑒𝑗
′ = 𝛼𝑖 + 𝛼𝑗
ln 𝑓𝑖′𝑓𝑗
′ = ln𝑓𝑖′ + ln𝑓𝑗
′ = 𝛽𝑖 + 𝛽𝑗
ln 𝑒𝑖′𝑓𝑗
′ = ln 𝑒𝑖′ + ln𝑓𝑗
′ = 𝛼𝑖 + 𝛽𝑗
ln 𝑓𝑖′𝑒𝑗
′ = ln𝑓𝑖′ + ln 𝑒𝑗
′ = 𝛽𝑖 + 𝛼𝑗
Cabe destacar que las igualdades anteriores sólo son ciertas si no existen logaritmos con parte real
negativa, tal y como se vio en el apartado 4.1.
26
26
Una vez obtenidas estas expresiones, la relación entre 𝑢 y 𝑥 es inmediata:
[ ln 𝑒𝑖
′2
ln 𝑓𝑖′2
ln 𝑒𝑖′
ln 𝑓𝑖′
ln 𝑒𝑖′𝑒𝑗
′
ln 𝑓𝑖′𝑓𝑗
′
ln 𝑒𝑖′𝑓𝑗
′
ln 𝑓𝑖′𝑒𝑗
′]
= 𝐶
[ 𝛼1
⋮𝛼𝑁
𝛽1
⋮𝛽𝑁]
Finalmente, se obtiene el mismo sistema de ecuaciones que se obtuvo anteriormente, pero este presenta algunas
diferencias:
[𝐸�̃�−1𝐶′]𝑥𝑘+1′ = 𝐸�̃�−1�̂� (83)
1) 𝐶′ es la matriz reducida de 𝐶 y se usa para que el producto [𝐸�̃�−1𝐶′] sea invertible. Esta matriz se
obtiene suprimiendo la columna correspondiente a 𝛽𝑠𝑙𝑎𝑐𝑘 de 𝐶, pero, al ser no nulo 𝛽𝑠𝑙𝑎𝑐𝑘 debido a que
𝛽𝑠𝑙𝑎𝑐𝑘 = ln 𝑓𝑠𝑙𝑎𝑐𝑘′ = ln(𝑓𝑠𝑙𝑎𝑐𝑘 + 𝑚) = ln𝑚 (suponiendo que 𝑓𝑠𝑙𝑎𝑐𝑘 es nulo), hay que hacer una
corrección en el vector 𝑢.
2) �̂� es una versión corregida del vector 𝑢, producto del problema anterior. Este se obtiene restando a 𝑢 la
columna correspondiente a 𝛽𝑠𝑙𝑎𝑐𝑘 de la matriz 𝐶:
�̂� = 𝑢 − 𝐶(: , 𝛽𝑠𝑙𝑎𝑐𝑘) ln(𝑚) (84)
En este caso, el jacobiano de 𝑓(∙) es de una forma más simple que la vista en el FLFRC sin offset:
�̃�−1 = 𝑑𝑖𝑎𝑔(𝑦) (85)
donde 𝑑𝑖𝑎𝑔(𝑦) es una matriz diagonal, de tamaño 4(𝑁 + 𝑏) × 4(𝑁 + 𝑏), cuyos elementos son los del vector
auxiliar 𝑦. La demostración de la expresión (85) queda recogida en el apartado ‘Anexo’.
Por último, decir que del sistema de ecuaciones final (83) se obtiene 𝑥𝑘+1′ , siendo 𝑥′ un vector reducido del
vector de estado 𝑥. Para obtener éste último vector, habría que añadir 𝛽𝑠𝑙𝑎𝑐𝑘 a 𝑥′, siendo 𝛽𝑠𝑙𝑎𝑐𝑘 = ln (𝑚) (con
la suposición que se vió anteriormente).
27
27
Ejemplo. Aplicación del FLFRC con offset a un sistema de dos nudos
Al igual que en el ejemplo desarrollado en el Anexo, se analiza el siguiente sistema:
Figura 10. Sistema de 2 nudos (versión 2)
El objetivo del problema es obtener la tensión compleja del nudo dos de la figura, 𝑉 ∟𝜃, sabiendo que la
potencia demandada en ese punto es de 200 MW y la reactiva demandada de 50 MW. Sabemos que la solución
es 𝒰2 = 0.9001 − 𝑗0.2 .
Empezamos definiendo el valor del offset y el criterio de convergencia:
Valor del offset: 𝑚 = 2
Criterio de convergencia, ∆𝑝𝑘 ≤ 0.001.
Datos de la red:
𝒴 = [𝐺11 + 𝑗𝐵11 𝐺12 + 𝑗𝐵21
𝐺21 + 𝑗𝐵21 𝐺22 + 𝑗𝐵22] = [
−𝑗10 𝑗10𝑗10 −𝑗10
], donde 𝒴 es la matriz de admitancias.
𝑈1 = 𝑒12 + 𝑓1
2 = 1 𝑝𝑢
𝑝 = [𝑈1
𝑃2
𝑄2
] = [1
−2−0.5
] 𝑃2 =−200 𝑀𝑊
100 𝑀𝑊= −2 𝑝𝑢
𝑄2 =−50 𝑀𝑊
100 𝑀𝑊= −0.5 𝑝𝑢
Vector de estado 𝑥 y vectores auxiliares, 𝑢 e 𝑦:
𝑥 = [
𝛼1
𝛼2
𝛽1
𝛽2
] → 𝑥′ = [
𝛼1
𝛼2
𝛽2
] ; 𝑦 =
[ 𝑒1
′2
𝑒2′ 2
𝑓1′2
𝑓2′2
𝑒1′
𝑒2′
𝑓1′
𝑓2′
𝑒1′𝑒2
′
𝑓1′𝑓2
′
𝑒1′𝑓2
′
𝑓1′𝑒2
′ ]
; 𝑢 = ln𝑦 =
[
2𝛼1
2𝛼2
2𝛽1
2𝛽2
𝛼1
𝛼2
𝛽1
𝛽2
𝛼1 + 𝛼2
𝛽1 + 𝛽2
𝛼1 + 𝛽2
𝛽1 + 𝛼2]
28
28
Matriz 𝐸 y vector �̂�:
𝑈1 = 𝑒1′2 + 𝑓1
′2 − 2𝑚𝑒1′ − 2𝑚𝑓1
′ + 2𝑚2
𝑈21 = 𝑒2′𝑒1
′ + 𝑓2′𝑓1
′ − 𝑚𝑒2′ − 𝑚𝑒1
′ − 𝑚𝑓2′ − 𝑚𝑓1
′ + 2𝑚2
𝑊21 = −𝑒2′𝑓1
′ + 𝑓2′𝑒1
′ + 𝑚𝑒2′ − 𝑚𝑒1
′ + 𝑚𝑓1′ − 𝑚𝑓2
′
𝑈1 = 𝑈1
𝑃2 = 𝐺22𝑈2 + 𝐺21𝑈21 + 𝐵21𝑊21
𝑄2 = −𝐵22𝑈2 + 𝐺21𝑊21 − 𝐵21𝑈21
Sustituyendo las expresiones anteriores en las expresiones de 𝑈1, 𝑃2 y 𝑄2 , y ordenando los términos,
obtenemos:
𝑈1 = 𝑒1′2 + 𝑓1
′2 − 2𝑚𝑒1′ − 2𝑚𝑓1
′ + 2𝑚2
𝑃2 = 𝐺22𝑒2′ 2 + 𝐺22𝑓2
′2 + 𝑒2′(−2𝑚𝐺22 − 𝑚𝐺21 + 𝑚𝐵21) + 𝑓2
′(−2𝑚𝐺22 − 𝑚𝐺21 − 𝑚𝐵21)+ 𝑒1
′(−𝑚𝐺21 − 𝑚𝐵21) + 𝑓1′(−𝑚𝐺21 + 𝑚𝐵21) + 𝑒1
′𝑒2′(𝐺21) + 𝑓1
′𝑓2′(𝐺21)
+ 𝑒2′𝑓1
′(−𝐵21) + 𝑓2′𝑒1
′(𝐵21) + 2𝑚2(𝐺22 + 𝐺21)
𝑄2 = −𝐵22𝑒2′ 2 − 𝐵22𝑓2
′2 + 𝑒2′(2𝑚𝐵22 + 𝑚𝐺21 + 𝑚𝐵21) + 𝑓2
′(2𝑚𝐵22 − 𝑚𝐺21 + 𝑚𝐵21)+ 𝑒1
′(𝑚𝐵21 − 𝑚𝐺21) + 𝑓1′(𝑚𝐵21 + 𝑚𝐺21) + 𝑒1
′𝑒2′(−𝐵21) + 𝑓1
′𝑓2′(−𝐵21)
+ 𝑒2′𝑓1
′(−𝐺21) + 𝑓2′𝑒1
′(𝐺21) + 2𝑚2(−𝐵22 − 𝐵21)
Con esto, la matriz queda definida como:
𝐸 = [1 0 1 0 −4 0 −4 0 0 0 0 00 0 0 0 −20 20 20 −20 0 0 10 −100 10 0 10 20 −20 20 −20 −10 −10 0 0
]
Y por último, el vector �̂�:
𝑝𝐼 = [
2𝑚2
2𝑚2(𝐺22 + 𝐺21)
−2𝑚2(𝐵22 + 𝐵21)
] = [800]
�̂� = 𝑝 − 𝑝𝐼 = [1
−2−0.5
] − [800] = [
−7−2
−0.5]
29
29
Cálculo de la matriz 𝐶 y de �̂�:
𝑢 = 𝐶𝑥 → ln(𝑦) = 𝐶 [
𝛼1
𝛼2
𝛽1
𝛽2
] =
[
2𝛼1
2𝛼2
2𝛽1
2𝛽2
𝛼1
𝛼2
𝛽1
𝛽2
𝛼1 + 𝛼2
𝛽1 + 𝛽2
𝛼1 + 𝛽2
𝛽1 + 𝛼2]
luego 𝐶 =
[ 2 0 0 00 2 0 00 0 2 00 0 0 21 0 0 00 1 0 00 0 1 00 0 0 11 1 0 00 0 1 11 0 0 10 1 1 0]
→ 𝐶′ =
[ 2 0 00 2 00 0 00 0 21 0 00 1 00 0 00 0 11 1 00 0 11 0 10 1 0]
�̂� = 𝑢 − 𝐶(: , 𝛽1) ∗ ln(𝑚) = 𝑢 −
[ 002000100101]
ln 2
�̃�−1 = 𝑑𝑖𝑎𝑔(𝑦)
Una vez definidas y construidas todas las variables del problema, se ejecuta el mismo proceso iterativo que se
ha desarrollado en el Ejemplo del apartado Anexo, hasta que se cumple la condición:
max(|∆𝑝𝑘|) ≤ 0.001
30
30
Una vez cumplida, se obtienen los valores de 𝛼2 y 𝛽2, por lo que hay que deshacer el cambio:
𝑒2′ = exp(𝛼2) → 𝑒2 = 𝑒2
′ − 𝑚 = 0.9000
𝑓2′ = exp(𝛽2) → 𝑓2 = 𝑓2
′ − 𝑚 = −0.2000
obteniéndose como solución:
𝒰2 = 0.9000 − 𝑗0.2000 ≈ 0.9001 − 𝑗0.2
Fin de ejemplo
Algoritmo de resolución con offset
De forma muy similar al algoritmo de resolución del FLF, se definen los siguientes pasos de resolución:
Paso 0: inicialización
Se escoge un perfil inicial de tensiones, normalmente un perfil plano, además es necesario fijar un valor para
el offset distinto de 0, ya que sin él, el proceso no podría ejecutarse debido a que el producto 𝐸𝐸′ se convierte
en singular. Con él se obtiene el vector de estado 𝑥0′ = [ln 𝑒1
′ , ln 𝑒2′ , … , ln 𝑒𝑁
′ | ln 𝑓1′, ln 𝑓2
′, … , ln 𝑓𝑁′ ]𝑇.
Una vez se tiene 𝑥0′ , se construye 𝑦0 = 𝑓−1(𝐶𝑥0
′ ).
Cabe destacar que, si vamos a analizar un caso que se encuentre en zona infactible, hay que inicializar el
método con un perfil con una cierta componente imaginaria en las variables 𝑒′ y 𝑓′ (es decir, 𝑒′ y 𝑓′ ambas
complejas), pues en caso contrario, el método no funciona. El valor óptimo de esta semilla imaginaria se verá
en la Sección 6, aunque podemos adelantar que con un valor muy pequeño (por ejemplo, 𝑗0.001), ya es
suficiente para que converja. En caso de que no sepamos si se trata de un caso factible o infactible, podemos
inicializar siempre con un perfil complejo. Esto último puede implicar que en zona factible se tarde un poco más
en converger (una iteración más, normalmente), tal y como se verá en la Sección 6.
Paso 1:
1.1. Construir �̂� = 𝑝 − 𝑝𝐼
1.2. Calcular ∆�̂�𝑘 = �̂� − �̂�𝑘, donde �̂�𝑘 = 𝐸𝑦𝑘
1.3. Obtener ∆𝑦𝑘, resolviendo:
∆�̂�𝑘 = 𝐸∆𝑦𝑘
1.4. Obtener �̃�𝑘 de:
�̃�𝑘 = 𝑦𝑘 + ∆𝑦𝑘
31
31
Paso 2:
2.1. Calcular 𝑢𝑘, transformarlo en �̂�𝑘 y calcular el jacobiano �̃�−1:
𝑢𝑘 = ln �̃�𝑘
�̃�−1 = 𝑑𝑖𝑎𝑔(�̃�𝑘)
Paso 3:
3.1. Obtener 𝑥𝑘+1′ , resolviendo:
[𝐸�̃�−1𝐶′]𝑥𝑘+1′ = 𝐸�̃�−1�̂�𝑘
3.2. Completar 𝑥𝑘+1′ con ln𝑚 , para obtener 𝑥𝑘+1
3.3. Obtener 𝑦𝑘+1:
𝑦𝑘+1 = ln𝑢𝑘+1
donde, 𝑢𝑘+1 = 𝐶𝑥𝑘+1
Paso 4:
4.1. Comprobar la convergencia, si ||�̂� − 𝐸𝑦𝑘+1||∞ es lo suficientemente pequeño, entonces seguir con
el paso 5. En caso contrario, volver al paso 1.2.
Paso 5:
5.1. Obtener:
[𝑒1′ , 𝑒2
′ , … , 𝑒𝑁′ |𝑓1
′, 𝑓2′, … , 𝑓𝑁
′ ]𝑇 = exp (𝑥′)
5.2. Deshacer el cambio, restando el valor del offset 𝑚 a cada componente del vector anterior,
obteniendo:
[𝑒1, 𝑒2, … , 𝑒𝑁|𝑓1, 𝑓2, … , 𝑓𝑁]𝑇
5.3. Construir el vector de tensiones en forma rectangular:
𝒰 = [
𝑒1 + 𝑗𝑓1𝑒2 + 𝑗𝑓2
⋮𝑒𝑁 + 𝑗𝑓𝑁
]
En la Figura 11 se muestra un diagrama de flujo que resume todo el proceso de manera más gráfica, al igual
que se hizo con el FLF.
32
32
Figura 11. Algoritmo de resolución del FLFRC
33
33
Límites de reactiva en nudos PV generalizada
Al igual que se vió en la Sección 2.3, la formulación del FLF y del FLFRC sólo tienen en cuenta las
restricciones de igualdad impuestas por los nudos de consumo y de generación. Sin embargo, en la realidad se
tienen en cuenta otros tipos de restricciones.
La restricción adicional más habitual es
la de límite de reactiva en nudos PV,
aplicándose a la potencia reactiva
generada o absorbida por los generadores
que regulan la tensión de un nudo.
Cuando se sobrepasa unos de estos
límites, el nudo de generación (nudo PV)
pasa a ser un nudo de consumo ficticio,
PQ (inyectando o absorbiendo potencia).
La implantación de un algoritmo que
tenga en cuenta la restricción anterior
hace que se incremente el número de
iteraciones del método numérico. El
algoritmo que se usaba anteriormente
implicaba calcular en cada iteración la
reactiva que un generador tiene que
inyectar para mantener su tensión
constante, y compararla con los límites. Si
se viola algún límite, la tensión del nudo
regulado no puede mantenerse en el valor
𝒰𝑒𝑠𝑝, pasando a ser éste un nudo de
consumo con 𝑄𝑒𝑠𝑝 = 𝑄𝑙𝑖𝑚. En las
ilustraciones de la derecha se expone un
ejemplo muy sencillo de un generador
que, debido a un aumento en la demanda,
supera los límites de reactiva,
convirtiéndose en un nudo PQ ficticio.
Como novedad, se propone un nuevo
algoritmo de tratamiento de límites de
reactiva en nudos PV. Éste está basado en
la idea de no tener que calcular la reactiva
inyectada por los generadores ni tener que
compararla con los límites en cada
iteración, sino sólamente cuando estamos
los suficientemente cerca de una solución
fiable.
Para implementar este algoritmo, basta con iniciar el método con un perfil complejo. Al traterse de zona
factible (ya que el control de reactiva en zona infactible carece de sentido), se convergerá siempre a una solución
en el dominio real, lo que implica que la parte imaginaria introducida en el perfil inicial de tensiones se irá
haciendo en cada iteración más pequeña, hasta hacerse nula.
Aprovechando la propiedad anterior, basta con incluir la condición de que se realice el control de reactiva sólo
cuando se está lo suficientemente cerca de una solución, es decir, cuando la parte imaginaria que se mencionó
anteriormente sea menor que un valor lo suficientemente pequeño, como por ejemplo de 10−4. Todo esto
ayudaría a rebajar la parte de cálculo del método numérico completo, ya que no se calcularía la reactiva ni se
compararía con sus límites en todas las iteraciones. Por último cabe destacar que al modificarse la topología de
un nudo (pasar de PV a PQ) implica un cambio importante en la matriz 𝐸 y en los vectores 𝑝 y 𝑝𝐼, lo que
complica bastante la programación del método, ya que hay que añadir y quitar dinámicamente filas y columnas
de los elementos anteriormente nombrados.
Figura 12. Nudo fuera de límites de reactiva
34
34
5 ENTORNO DE PROGRAMACIÓN
in entrar en detalle en el código del programa y en sus respectivas partes, el objetivo de este apartado
es definir e ilustrar las fases en las que se ha dividido la programación completa del FLFRC. El siguiente
esquema resumiría todo lo anterior:
S
35
35
Figura 13. Fases de la programación
Como se observa, consta de cuatro fases, pero la primera fase de lectura de párametros de la red solo sería
necesario que se ejecutara una única vez, mientras la topología de la red no cambie. De esta forma se guardan
los datos obtenidos en esa primera fase y empezaríamos directamente desde la segunda. La definición de los
parámetros de entrada dependen del usuario y del tipo de análisis que quiera realizar, y a ello le seguiría la
ejecución del FLFRC, en el cual se ha implementado la medición del tiempo de ejecución para extraer
conclusiones y compararlo con otros programas. Por último vendría la presentación de aquellos resultados en
los que se esté interesado, ya que resuelto el problema de tensiones, todas las demás magnitudes eléctricas son
conocidas.
36
36
6 RESULTADOS
n este apartado se aplicará el FLFRC a cinco redes muy distintas con el objetivo de analizar el método
y extraer conclusiones. Para ello se ha realizado una serie de tablas comparativas, que recogen
información de todas las redes (cinco en total). Posteriormente se desarrollará un estudio más detallado
en dos de éstas redes, de manera más gráfica. La primera red de las dos nombradas anteriormente consta de 14
nudos, y es una red de prueba propuesta por la IEEE [3], la segunda red consta de 3120 nudos, y representa la
red polaca de transporte del verano de 2008 [4]. En el apartado Anexo se realizará el mismo estudio a otras redes
de diversos tamaños, dos de ellas son otras redes de prueba, propuestas por la IEEE, de 118 y 300 nudos [3] y
otra de 1354 nudos que representa la red de trasnporte europea de alto voltaje [4].
Uno de los objetivos que se persigue, entre otros, es el de minimizar el número de iteraciones y el tiempo de
cálculo, aunque, para este último el programa desarrollado no ha sido optimizado. Sin embargo, se obtienen
unos tiempos bastante buenos en comparación con otra plataforma, MatPower V5.1. [4], cuyos resultados
veremos a continuación en forma de tablas.
Antes de analizar los resultados cabe destacar que el FLFRC contiene hasta dos parámetros libres, los cuales
hay que fijar antes de iniciar el método, siendo estos: el valor del offset (𝑚) y la semilla imaginaria (S.I.). El
hecho de tener esta libertad no es un problema si nuestro único objetivo es el de obtener el estado de la red, pero
si queremos reducir el número de iteraciones y el coste computacional, habrá que sintonizar las variables
adecuadas al valor óptimo.
El nivel de carga del sistema se identifica con la letra 𝜆, y es un valor constante que multiplica a todas las
potencias de todos los nudos, ya sean generadas o demandadas, lo que aumenta o disminuye el nivel de carga
original del sistema. Este valor pertenece al modelo de la red, por lo que se establece antes de aplicar el FLFRC.
A partir de un cierto valor de 𝜆 el caso se convierte en infactible, obteniéndose resultados con parte imaginaria
en cada componente, este valor es denominado como límite de máxima cargabilidad.
Nota: los estudios se han realizado en un ordenador portátil con un procesador Intel Celeron P4600, 2.00 GHz
2.00 GHz, con una memoria RAM de 3GB y con Windows 8.1 como sistema operativo.
Tablas comparativas para distintas redes
A continuación se recoge en forma de tabla un primer análisis del FLFRC aplicado al caso base de cada red
citada anteriormente. Consideramos el caso base como aquel en el que el nivel de carga es el original que traen
las redes. En este primer estudio se ha utilizado un perfil plano, usando una semilla imaginaria (𝑆. 𝐼.) nula y el
valor del offset (𝑚) que minimiza el número de iteraciones y tiempo de cálculo. El límite de convergencia se ha
establecido en 10−4, valor a partir del cual (valores menores del límite de convergencia) se considera que se ha
obtenido una solución lo suficientemente fiable.
Tabla 3. Pimer análisis del FLFRC ante casos bases de distintas redes
E
37
37
Ahora se analizan los mismos casos anteriores realizando algunas modificaciones. En las siguientes tablas que
se muestran, las primeras dos columnas, se trata el método variando la semilla imaginaria. La primera de ellas
es con una 𝑆. 𝐼. de 𝑗0.3 y la otra con una 𝑆. 𝐼. nula que son sumadas al perfil plano, usando siempre el mismo
valor de offset (𝑚 = 50), valor genérico que permite la convergencia en todos estos casos. Todo lo anterior se
compara con una herramienta de flujo de cargas, denominado MatPower. La primera tabla que se observa se ha
realizado para el caso base de cada sistema, y la siguiente a ella, para el límite de máxima cargabilidad. En ellas
se analiza también las iteraciones y el tiempo de cálculo que el método ha necesitado.
Tabla 4. Iteraciones y tiempo de cálculo para el caso base
Tabla 5. Iteraciones y tiempo de cálculo para el límite de máxima cargabilidad
En la Tabla 4 se ha ejecutado el FLFRC ante el caso base de cada red. Se observa que se obtienen unos
resultados muy parecidos si lo comparamos con MatPower, a pesar de no estar optimizado el código, y que la
principal diferencia entre inicializar el método con una semilla imaginaria o no, reside en el tiempo de cálculo,
siendo este menor cuando no se incluye una semilla imaginaria.
Si comparamos la Tabla 3 con la segunda columna de la Tabla 4 observamos que el hecho de tomar un offset
genérico que sirva ante todos los casos, en vez de usar un offset óptimo para cada red, implica a veces una
iteración más. En concreto lo anterior sucede para las redes de 118 y 300 nudos. Como ventaja se observa que
para la red más grande, de 3120 nudos, el FLFRC converge en una iteración menos que el método de Newton-
Raphson, sin embargo, para las redes de 300 y 1354 nudos, ocurre lo contrario.
En la Tabla 5, se ha realizado el mismo estudio que en la Tabla 4 con la diferencia de que en ésta se ha variado
el estado de la carga, el cuál se ha fijado en el máximo posible (límite de máxima cargabilidad), justo antes de
que aparezcan componentes imaginarias en la solución. Se extraen las mismas conclusiones que anteriormente
y cabe destacar que es posible reducir el número de iteraciones y el tiempo de cálculo escogiendo los valores
óptimos de los parámetros libres para cada red.
En la siguiente tabla se analizará el mismo caso que vimos en la tabla anterior, es decir, se analizará el límite
de máxima cargabilidad, pero esta vez implementando el control de límites de reactiva en nudos PV. Este control
se ha realizado de dos formas distintas:
38
38
FLFRC 1: se calcula la reactiva inyectada por los generadores y se compara con los límites cuando se
está lo suficientemente cerca de una solución que, en este caso, se considera cuando se alcanzan tres
iteraciones. Es decir, cuando el programa ha realizado tres iteraciones, se implementa el control de
reactiva.
FLFRC 2: el algoritmo que se utiliza es el límite de reactiva para nudos PV generalizada, visto en el
Apartado 4.4. Como en este caso es necesario una semilla imaginaria, se ha escogido una de valor 𝑗0.01
para todas las redes, aunque éste no es el valor óptimo, valor que se podría sintonizar para cada red con
el objetivo de reducir el número de iteraciones.
Tabla 6. Iteraciones y tiempo de cálculo en el límite de máxima cargabilidad con límites de reactiva
En la Tabla 6 se puede observar cómo el primer método (FLFRC 1) realiza un número de iteraciones mayor
que en el segundo método (FLFRC 1) para redes pequeñas, aunque para redes grandes ocurre lo contrario.
Además también se ve cómo entre ambos métodos varía ligeramente el límite máxima cargabilidad, como por
ejemplo para el caso de 118 y 1354 nudos. Los resultados obtenidos son muy similares a los de MatPower,
aunque hay que resaltar que este programa no puede realizar este estudio para la red de 3120 nudos.
Como novedad del FLFRC se muestra la Tabla 7, ofreciéndose una comparativa entre el FLF y el FLFRC,
la cual justifica por qué se formuló esta versión del FLFRC con offset, ya que el resultado de introducir el offset
en la formulación en polares (FLF) resulta un modelo matemático mas complejo, y si no se introduce el offset
se corre el peligro de ejecutar un caso en polares que no converja debido al problema de los logaritmos de
complejos de parte real negativa.
En esta tabla podemos apreciar la carga máxima que soporta la fomulación en polares (FLF), a partir de la cual
el método no converge. Con estos niveles de carga nos referimos al 𝜆𝑚𝑎𝑥 para que el método converja, ya en
zona infactible con resultados complejos. Estos valores se comparan con las cargas máximas que se han podido
alcanzar en el FLFRC de dos formas distintas, la primera usando el valor típico del offset, y la segunda regulando
los parámetros para lograr la convergencia. Sin embargo este límite mostrado es superable si se escogen los
parámetros libres óptimos para cada carga, lo que implicaría un estudio más complejo sobre estrategias de
inicialización que se deja pendiente como trabajo futuro (ver Sección 8).
Tabla 7. Comparación de cargas máximas de convergencia entre el FLF y el FLFRC
39
39
Se observa que se alcanzan cargas mucho mayores de las que se consiguen en el FLF, lo que implica una
novedad ventajosa. Se aprecia como usando un valor del offset genérico puede obtenerse un ligero aumento en
la carga máxima de convergencia para cada red respecto al FLF. Sin embargo, si se regula éste valor del offset
y la semilla imaginaria, observamos cómo la carga máxima de convergencia aumenta notablemente para cada
red, en especial para la red de 14 nudos en la que se ha alcanzado un valor muy extremo, lo que demuestra la
robustez del método. A pesar de todo esto, aunque resulte una ventaja respecto al FLF, estos valores de carga
tan extremadamente altos no tienen sentido en la realidad, ya que nunca se alcanzarían estos valores.
Una vez visto resultados del FLFRC en forma de tablas, se analizarán otros temas de manera mas visual,
gráficamente. Principalmente se verá la relación que existe entre los parámetros libres y el número de iteraciones,
concluyéndose con una aplicación muy interesante del FLFRC.
Cuenca de atracción del FLFRC para el caso ejemplo de dos nudos
Como se comentó en el Apartado 3.3, se ha aplicado la idea de visualizar la región de atracción del FLFRC
para el sistema de la Figura 3. El resultado es el siguiente:
Figura 14. Región de atracción con FLFRC
Podemos apreciar que se obtienen bastantes diferencias respecto al método convencional de Newton-Raphson
(Figura 4) y el FLF (Figura 5). En este caso se observa que la zona de convergencia es del 100%, a diferencia
que en los anteriores métodos, lo que le convierte en un método más robusto. La zona amarilla situada en la
parte inferior, significa que se ha alcanzado una solución en 20 iteraciones, pero esta solución no es la solución
real del problema (recordemos que ésta era de 𝒰2 = 0.9209 − 𝑗0.0913 𝑝𝑢), resultando |𝒰2| = 0.118 𝑝𝑢,
solución que es totalmente inviable. A pesar de la ventaja de mayor robustez, observamos que el FLFRC necesita
algunas iteraciones más para converger si lo comparamaos con el FLF, para ciertas zonas (ver Figura 5).
40
40
Comportamiento del FLFRC ante la variación de distintos parámetros
En este apartado se analizará el comportamiento del FLFRC ante la variación de distintos parámetros, como
por ejemplo el nivel de carga 𝜆, el perfil inicial, la semilla imaginaria, etc.
En la Figura 15 y Figura 16 se muestran la parte real e imaginaria de las componentes e y f, respectivamente,
para cada nudo de la red de 14 nudos. El estudio se ha realizado variando el nivel de la carga, de manera que se
diferencien el caso factible del infactible, con el objetivo de visualizar cuando empiezan a aparecer componentes
imaginarias en la solución. Se ha tomado un valor genérico del offset 𝑚 = 50, y una semilla imaginaria de 𝑗0.5.
Figura 15. Comportamiento de la componente e para un sistema de 14 nudos
Figura 16. Comportamiento de la componente f para un sistema de 14 nudos
41
41
Se observa que conforme aumenta la carga, la tensión de cada nudo va disminuyendo y que, una vez superado
el límite de máxima cargabilidad (𝜆 = 4.0602), empiezan a aparecer componentes imaginarias en la solución.
En las próximas gráficas se muestra como evoluciona el máximo residuo conforme se avanza en las iteraciones
necesarias para alcanzar la solución, con dos métodos distintos. El primero de ellos es el FLFRC y el segundo
el método de Newton-Raphson, obtenido desde MatPower. Para este estudio se ha tomado como offset 𝑚 = 10,
siendo nula la semilla imaginaria, partiéndose desde un perfil plano.
Figura 17. Evolución del residuo para el caso base de una red de 14 nudos
Figura 18. Evolución del residuo para el caso base de una red de 3120 nudos
42
42
Observamos que la curva del FLFRC está por encima de la curva de Newton-Raphson para la red de 14 nudos,
demostrando una desventaja. Sin embargo, en la red de 3120 nudos ocurre lo contrario, lo cual resulta ventajoso
ya que se converge antes.
A continuación se expone como afecta el uso de distintos perfiles iniciales ante un aumento progresivo del
nivel de la carga, desde zona factible a zona infactible. Ante zona factible, se usan dos perfiles iniciales distintos,
el primero de ellos es el perfil plano, y el segundo un perfil con una cierta componente imaginaria (1 + 𝑗0.5).
Ante zona infactible siempre se usa el perfil inicial complejo, 1 + 𝑗0.5.
Figura 19. Aplicación de distintos perfiles ante la evolución de la carga para la red de 14 nudos
Figura 20. Aplicación de distintos perfiles ante la evolución de la carga para la red de 3120 nudos
43
43
Observamos que comenzar con un perfil complejo requiere a veces una iteración más que si se empezara con
el perfil plano. A partir del límite de máxima cargabilidad (marcado con una línea roja vertical) se usa siempre
un perfil inicial complejo, con parte imaginaria de 𝑗0.5. Cabe destacar que las puntas que se observan
(iteraciones) pueden reducirse si se usan los valores óptimos del offset y de la semilla imaginaria.
Análisis del número de iteraciones usando mapas de colores
Para concluir con las gráficas, en este apartado se analizará el comportamiento del método ante la variación de
distintos parámetros para analizar cómo varía el número de iteraciones necesario para converger.
Las gráficas siguientes se tratan de mapas de colores, cada color significa un número distinto de iteraciones,
cuya leyenda se muestra en el margen derecho. Los tonos más cálidos indican un mayor número de iteraciones,
mientras que los azulados significan un número de iteraciones más reducido.
La Figura 21 y la Figura 22 se han realizado manteniendo constante el valor del offset en 𝑚 = 50, variando
la semilla imaginaria (eje vertical) y el nivel de carga 𝜆 (eje horizontal), apreciándose las distintas iteraciones
con distintos colores. La parte real del perfil inicial de tensiones es siempre 1 𝑝𝑢.
Figura 21. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 14 nudos
44
44
Figura 22. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 3120 nudos
En zona factible (parte izquierda de la franja negra vertical) se observa que conforme menos semilla imaginaria
se use, menor número de iteraciones necesita para converger, aunque en el rango de semilla imaginaria que se
muestra siempre se converge en pocas iteraciones. Si nos fijamos en la parte derecha de la franja, es decir, en
zona infactible, se observa que conforme nos acercamos al límite superior e inferior del eje de ordenadas el
número de iteraciones crece bastante, empeorando el comportamiento del método. Es en la zona intermedia de
este eje donde se alcanza el menor número de iteraciones, es decir, con una semilla imaginaria de 𝑗0.5
aproximadamente, aunque 𝑗0.3 también es un buen valor.
Como conclusión podemos extraer que es posible establecer un perfil general complejo para todo tipo de
situaciones, cuya parte imaginaria está entorno a 𝑗0.3 y 𝑗0.5, lo que ayudaría al método a agilizar los cálculos,
reduciendo el número de iteraciones necesario en zona infactible para cualquier red, aunque el precio a pagar es
que en zona factible puede aumentar levemente el número de iteraciones.
Todo esto no significa que nos encontremos en el valor óptimo de la semilla imaginaria, pero si que sería una
posible aproximación. Destacar por último que este estudio se ha realizado fijando el offset en 70, luego el
resultado mostrado puede variar si se varía este valor, mejorando o empeorando el método.
En los mapas de colores que se exponen a continuación se ha aplicado la misma idea que en la Figura 21 y
Figura 22, con la diferencia de que en este caso se ha fijado la semilla imaginaria en 𝑗0.3, variando el offset (eje
vertical) y la carga del sistema (eje horizontal).
45
45
Figura 23. Iteraciones frente a la variación del offset y la carga en la red de 14 nudos
Figura 24. Iteraciones frente a la variación del offset y la carga en la red de 3120 nudos
Ante zona factible, comprobamos que con cualquier valor del offset el método converge, y que este valor no
influye en el número de iteraciones ante una carga dada, ya que se observa que el cambio del número de
iteraciones se produce una vez pasado un determinado valor de la carga.
Ante zona infactible observamos que el offset sí que puede influir en el número de iteraciones ante una carga
constante. Como conclusión se puede decir que conforme aumenta la carga del sistema, habría que aumentar el
valor del offset.
Cabe destacar que si el offset es nulo, el método no converge, de ahí ese color marrón oscuro en la parte
inferior. Por último decir de nuevo que, el valor de la semilla imaginaria se ha fijado en 𝑗0.3 para realizar el
estudio pero, si la variamos, los resultados podrían cambiar a peor o a mejor.
46
46
Ejemplo de aplicación del FLFRC
Para concluir con los resultados, se ilustrará una de las posibles aplicaciones del FLFRC. En concreto veremos
cómo determinar los nudos más críticos de una red, basándonos en la solución compleja que ofrece el FLFRC
ante casos infactibles.
A medida que va aumentando la carga, aquel nudo donde empiecen a aparecer componentes imaginarias en la
solución, será el más crítico, pudiéndose definir zonas críticas atendiendo a este criterio. Para ello es necesario
establecer un criterio que nos permita decidir cuando considerar si un nudo es crítico o no. Se ha escogido el
siguiente:
Descomponiendo las componentes del vector tensión:
𝒰 = 𝑒 + 𝑗𝑓 = (𝑒𝑟 + j𝑒𝑖) + (𝑓𝑟 + j𝑓𝑖)𝑗 = (𝑒𝑟 + 𝑓𝑟𝑗) + (−𝑓𝑖 + 𝑒𝑖𝑗) = 𝒰𝑟 + 𝒰𝑖
donde 𝑒𝑟 y 𝑒𝑖 es la parte real e imaginaria de la componente 𝑒, respectivamente (al igual ocurre con la
componente 𝑓). 𝒰𝑟 y 𝒰𝑖 representarían la contribución al vector 𝒰 por parte de componentes reales e
imaginarias, respectivamente.
Cuando estamos en caso factible, 𝒰 = 𝒰𝑟, ya que las componentes imaginarias de 𝑒 y 𝑓 son nulas.
Cuando pasamos al caso infactible la componente 𝒰𝑖 = (−𝑓𝑖 + 𝑒𝑖𝑗) empieza a ser distinto de cero, por
lo que el criterio que se escoge para determinar si un nudo es crítico o no, consiste en estudiar la siguiente
inecuación:
√𝑓𝑖2 + 𝑒𝑖
2 ≥ 𝑘
donde 𝑘 es un umbral a partir del cual el resultado se considera complejo, es decir, un resultado
infactible. Este umbral podemos fijarlo, por ejemplo, en 0.1, aunque la elección de este parámetro es
algo libre, ya que podríamos haber escogido 0.001 en vez de 0.1.
Los primeros nudos que sobrepasen éste límite 𝑘 aumentando la carga poco a poco, lo consideraremos nudos
críticos, ya que serán los nudos cuya componente imaginaria será la de mayor valor. El primer nudo que
sobrepase ese valor será el más crítico, pero conforme aumentamos la carga podemos ir definiendo zonas
críticas, catalogándolas como zonas más vulnerables.
Para ilustrar esta aplicación se ha escogido una red de 118 nudos, red de prueba propuesta por IEEE [3]. En el
siguiente diagrama se muestra como va evolucionando la red aplicando toda la teoría desarrollada anteriormente,
observando como poco a poco ciertos nudos se van convirtiendo en nudos críticos.
47
47
48
48
Figura 25. Aplicación del FLFRC para determinar nudos críticos (versión 1)
49
49
En este primer mapa observamos que los nudos más críticos son los nudos 34 y 36, entre los que se encuentran
dos unidades de generación y otras dos de consumo. Si seguimos aumentando la carga poco a poco, empiezan
a surgir nuevas zonas críticas, todas ellas caracterizadas por tener bastantes nudos de generación.
Visto esto, nos formulamos la siguiente pregunta, ¿qué pasaría si en vez de cargar toda la red, cargamos todos
los nudos a excepción de los más críticos (nudos 34 y 36)? La respuesta se ilustra a continuación:
50
50
51
51
Figura 26. Aplicación del FLFRC para determinar nudos críticos (versión 2)
En este caso vemos que los nudos más críticos vuelven a ser el 34 y 36, obteniendo casi las mismas zonas
críticas que se obtuvieron anteriormente. Como principal diferencia vemos que la carga necesaria para que
aparezcan los nudos críticos (𝜆 = 3.4919) es mayor de la que se obtenía en el mapa anterior (𝜆 = 3.1949), por
lo que el límite de máxima cargabilidad ha aumentado. Como conclusión podemos decir que, efectivamente, los
nudos 34 y 36 son los más críticos de la red, por lo que si se “refuerzan” estos nudos, como por ejemplo
aumentando la capacidad de generación, reforzando las líneas o añadiendo elementos compensadores, la red
estaría preparada para una mayor carga.
Nota: el esquema de la red de 118 nudos usado en esta aplicación ha sido obtenido de [27].
52
52
7 CONCLUSIONES
na vez analizado el comportamiento del FLFRC se concluye este trabajo mencionando algunas de las
conclusiones más importantes que se han extraído a lo largo del desarrollo del mismo. Antes de
continuar, hay que aclarar que la aparición de componentes imaginarias en el módulo y argumento
de la tensión de un nudo carece de sentido, ya que se tratan de magnitudes que no se podrían dar en la realidad.
A pesar de todo ello, tiene bastante utilidad, ya que gracias a éste fenómeno podríamos averiguar los nudos más
críticos de una red, o determinar ante casos infactibles como de alejados estamos del límite de máxima
cargabilidad.
Como se explicó al principio del apartado de Resultados, el FLFRC tiene una mayor región de atracción que
el método convencional de Newton-Raphson y que el FLF. Además siempre converge ante una gama muy
amplia de perfiles iniciales para el sistema de dos nudos que se analizó, luego esto lo convierte en un método
más robusto. Como desventaja presenta que en determinadas regiones el FLFRC necesita más iteraciones que
el FLF para alcanzar la solución.
Fijándonos en las tablas mostradas en el capítulo de Resultados, se observa que el número de iteraciones es
bastante parecido si lo comparamos con la aplicación MatPower V5.1., aunque el FLFRC conlleva a veces
alguna iteración más. Además, el tiempo de cálculo es bastante similar en ambos casos a pesar de no estar
optimizado el código del FLFRC. No hay que olvidar que es posible reducir el número de iteraciones de las
tablas donde se estudia cada red en su límite de máxima cargabilidad, para ello hay que sintonizar los parámetros
libres en los valores óptimos para cada red. Como ventaja del FLFRC también se ha visto como aumenta bastante
la carga máxima sin que el método no converja (carga máxima de convergencia) respecto al FLF, gracias a la
introducción del offset para solucionar el problema con los logaritmos complejos de parte real negativa.
En cuanto a la evolución del residuo en cada iteración, observamos que el FLFRC saca ventaja al método de
Newton-Raphson cuando se trata de redes grandes, a pesar de que en redes pequeñas a veces ocurre lo contrario,
lo que no resulta ventajoso.
Si nos fijamos en las demás gráficas podemos afirmar que si sintonizamos adecuadamente los parámetros
libres del FLFRC podemos minimizar el número de iteraciones necesario para alcanzar una solución ante un
determinado caso. Además, es posible definir un perfil general de parámetros libres que sirva para abordar
cualquier caso, ya sea factible o infactible. Este se puede definir estableciendo el valor del offset entre 50 o 70,
por ejemplo, y una semilla imaginaria entre 𝑗0.3 y 𝑗0.5. La desventaja de iniciar el proceso con el perfil general
anterior es que, ante casos factibles, el método suele tardar más en converger, normalmente una iteración más,
ya que se debe anular la componente imaginaria que se añadió en el perfil inicial de tensiones.
También se ha desarrollado una aplicación muy interesante del FLFRC para detemrinar los nudos más críticos
de una red, pudiendo aumentar el límite de máxima cargabilidad de una red si los nudos más críticos de la misma
se refuerzan, tal y como se vió en el final de la Sección 6.
Por último hay que destacar que todo esto no es un tema cerrado, ya que quedan muchos aspectos que mejorar
y problemas interesantes que resolver en un futuro próximo. En la siguiente sección se incluyen algunos de éstos
propósitos.
U
53
53
8 TRABAJO FUTURO
n este apartado se indican algunas ideas que han ido surgiendo durante la realización de este trabajo, y que
podrían desarrollarse en un futuro próximo para tratar de mejorar ciertos aspectos del FLFRC. Algunas
de ellas son:
Optimización del código: el objetivo es reducir el tiempo de ejecución del programa.
Carga máxima sin que no converja: consiste en averiguar el valor λmax después del límite y compararla
con la que se obtuvo en el FLF.
Offset personalizado en cada nudo: como actualmente se presenta el problema de que se alcanza un
determinado estado de carga máximo (tras sobrepasar el límite de máxima cargabilidad), se pretende
aumentar este límite utilizando distintos offset para los distintos nudos.
Estudio de la cargabilidad de cada nudo: se trata de profundizar más en la última aplicación mostrada
en el Capítulo 6, es decir, buscar el 𝜆𝑚𝑎𝑥 por nudo.
Otras estrategias de inicialización: consiste en encontrar nuevas estrategias para inicializar el proceso
iterativo ante casos factibles e infactibles, como por ejemplo empezando directamente con el vector
auxiliar ‘𝑦’ como perfil inicial o encontrar otros valores óptimos. Con esto se pretende mejorar el
comportamiento del FLFRC, además de aumentar aún más la carga máxima de convergencia.
E
54
54
ANEXO
FLFRC: Formulación sin offset
De forma similar a la formulación vista en el FLF (Sección 3.1), tenemos el sistema de ecuaciones original del
flujo de carga:
ℎ(𝑥) = 𝑝 (86)
El FLFRC hace desglosar el sistema anterior (86) en el mismo sistema de tres ecuaciones que vimos en la
Sección 3.1:
𝐸𝑦 = 𝑝 (87)
𝑢 = 𝑓(𝑦) (88)
𝐶𝑥 = 𝑢 (89)
La solución del FLFRC implica una modificación en el vector de estado 𝑥, pues ahora se hace uso de
coordenadas cartesianas (𝒰𝑖 = 𝑒𝑖 + 𝑗𝑓𝑖 ), lo que conlleva un cambio en todas las variables del problema:
El vector de estado 𝑥, compuesto por 𝑛 = 2𝑁 elementos (𝑁 es el número de nudos), viene dado por:
𝑥 = [𝑎1, 𝑎2, … , 𝑎𝑁|𝑏1, 𝑏2, … , 𝑏𝑁]𝑇 (90)
donde,
𝛼𝑖 = ln𝒰𝑖 = 𝑎𝑖 + 𝑗𝑏𝑖 (91)
Vector y:
𝑦 = {𝑈𝑖, 𝑈𝑖𝑗 ,𝑊𝑖𝑗} (92)
donde
𝑈𝑖 = 𝑉𝑖2 = 𝑒𝑖
2 + 𝑓𝑖2 (93)
y para cada rama que conecta el nudo 𝑖 con el nudo 𝑗, la pareja de variables 𝑈𝑖𝑗 y 𝑊𝑖𝑗 viene dada
por:
𝑈𝑖𝑗 = 𝑒𝑖𝑒𝑗 + 𝑓𝑖𝑓𝑗 (94)
𝑊𝑖𝑗 = −𝑒𝑖𝑓𝑗 + 𝑓𝑖𝑒𝑗 (95)
Dándose las siguientes propiedades:
𝑈𝑖𝑗 = 𝑈𝑗𝑖
𝑊𝑖𝑗 = −𝑊𝑗𝑖
55
55
Al igual que anteriormente, el vector 𝑦 está compuesto por 𝑚 = 2𝑏 + 𝑁 variables, siendo 𝑏 el
número de ramas y 𝑁 el número de nudos.
El vector 𝑢, compuesto también por 𝑚 variables, obedece a la relación no lineal 𝑢 = 𝑓(𝑦) (88), por
lo que viene dado por:
𝑢 = [
ln𝑈𝑖
𝑅𝑒(ln𝑤𝑖𝑗)
𝐼𝑚(ln𝑤𝑖𝑗)] (96)
donde,
𝑤𝑖𝑗 = 𝒰𝑖𝒰𝑗∗ = 𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗 (97)
𝑈𝑖 = 𝒰𝑖𝒰𝑖∗ = 𝑒𝑖
2 + 𝑓𝑖2 (98)
Verificándose la siguiente igualdad:
𝑦 = exp(𝑢) = [
𝑈𝑖
𝑈𝑖𝑗
𝑊𝑖𝑗
] = 𝑓−1(𝑢) (99)
𝑢 = ln 𝑦 = [
ln𝑈𝑖
𝑅𝑒(ln𝑤𝑖𝑗)
𝐼𝑚(ln𝑤𝑖𝑗)] (100)
Antes de continuar, cabe destacar que fue aquí donde se identificó el problema citado en el primer
párrafo de la Sección 4.2, el cuál dió lugar a la formulación con offset. Concretamente se observó
que la parte real de 𝑤𝑖𝑗, es decir 𝑈𝑖𝑗, se hacía negativa, luego aparecían logaritmos complejos de
parte real negativa, dando lugar al problema citado en la Sección 4.1.1.
El vector 𝑝 también es diferente al del FLF. Este contiene el módulo de la tensión al cuadrado de
cada nudo de generación (nudos PV), la potencia activa inyectada en todos los nudos, menos en el
nudo slack, y la potencia reactiva inyectada en todos los nudos de consumo (nudos PQ). Luego viene
dado por:
𝑝 = [𝑈1, 𝑈2, … , 𝑈𝑁𝑃𝑉|𝑃1, 𝑃2, … , 𝑃𝑁−1|𝑄1, 𝑄2, … , 𝑄𝑁𝑃𝑄
]𝑇 (101)
A la hora de calcular la relación entre el vector 𝑝 con el vector 𝑦, es decir, la matriz 𝐸 (87), es
necesario conocer la expresión analítica de la potencia activa y reactiva inyectada en cada nudo,
además de la tensión de cada nudo de generación. Al tratarse de coordenadas cartesianas, estas
expresiones vienen dadas de la siguiente forma, ya escritas en variables del vector 𝑦, lo que conduce
a una relación directa entre los vectores 𝑝 e 𝑦:
𝑈𝑖𝑆𝑃 = 𝑈𝑖 (𝑝𝑎𝑟𝑎 𝑛𝑢𝑑𝑜𝑠 𝑃𝑉) (102)
𝑃𝑖𝑆𝑃 = 𝐺𝑖𝑖𝑈𝑖 + ∑[𝐺𝑖𝑗𝑈𝑖𝑗 + 𝐵𝑖𝑗𝑊𝑖𝑗]
𝑁−1
𝑗≠𝑖
(103)
56
56
𝑄𝑖𝑆𝑃 = −𝐵𝑖𝑖𝑈𝑖 + ∑[𝐺𝑖𝑗𝑊𝑖𝑗 − 𝐵𝑖𝑗𝑈𝑖𝑗]
𝑁𝑃𝑄
𝑗≠𝑖
(104)
siendo la dimensión de 𝐸 de 𝑛 × 𝑚, donde 𝑛 = 2𝑁 − 1, en cuya matriz no se incluye la relación
con la tensión del nudo slack, 𝑈𝑠𝑙𝑎𝑐𝑘 .
Al variar la forma del vector de estado 𝑥 y el vector auxiliar 𝑢, la matriz 𝐶 que los relaciona también
se ve modificada, cumpliéndose la siguiente ecuación:
𝑢 = 𝐶′𝑥′ (105)
𝑢 = [
ln𝑈𝑖
𝑅𝑒(ln𝑤𝑖𝑗)
𝐼𝑚(ln𝑤𝑖𝑗)] = [
2𝑎𝑖
𝑎𝑖 + 𝑎𝑗
𝑏𝑖 − 𝑏𝑗
] (106)
Gracias a la anterior relación (106), se obtienen los coeficientes de la matriz 𝐶. La dimensión de 𝐶′ es de 𝑚 × (𝑛 − 1), luego la dimensión de 𝐶 será de 𝑚 × (𝑛 + 1).
Notar que 𝑥′ es una expresión reducida del vector de estado 𝑥, en la cual no se tiene en cuenta ni la
componente 𝑎𝑠𝑙𝑎𝑐𝑘 ni 𝑏𝑠𝑙𝑎𝑐𝑘 ya que al tratarse del nudo slack, la tensión y fase ya son conocidas,
permitiendo que las ecuaciones en las que interviene sean dimensionalmente correctas.
Al igual que anteriormente, 𝐶′ es una matriz reducida de la matriz completa 𝐶 obtenida al eliminar
las columnas referentes a 𝑎𝑠𝑙𝑎𝑐𝑘 y 𝑏𝑠𝑙𝑎𝑐𝑘 de la matriz 𝐶.
Eliminando el vector 𝑢, se puede escribir de forma más compacta el sistema de ecuaciones desglosado:
𝐸𝑦 = 𝑝 (107)
𝐶𝑥 = 𝑓(𝑦) (108)
Cuando se elimina el vector 𝑦, se obtiene el sistema:
𝐸𝑓−1(𝐶𝑥) = 𝑝 (109)
a partir del cual, al igual que en el FLF, se obtiene una solución parecida a la del sistema original de Newton-
Raphson:
[𝐸𝐹𝑘−1𝐶]∆𝑥𝑘 = 𝑝 − 𝐸𝑓−1(𝐶𝑥) = ∆𝑝𝑘 (110)
Por lo que, finalmente, el vector de estado reducido 𝑥′ se obtiene de la siguiente forma:
�̃�𝑥𝑘+1′ = 𝐸�̃�−1�̃�′ (111)
[𝐸�̃�−1𝐶′]𝑥𝑘+1′ = 𝐸�̃�−1�̃�′ (112)
siendo �̃� el jacobiano de 𝑓(∙); Por tanto que �̃�−1 viene dado por:
57
57
�̃�−1 = [𝜕𝑦
𝜕𝑢] = [
𝑈𝑖 0 00 𝑈𝑖𝑗 −𝑊𝑖𝑗
0 𝑊𝑖𝑗 𝑈𝑖𝑗
] (113)
Cabe destacar que la relación anterior no es trivial, y en el Anexo se da su demostración.
Por último, no hay que olvidar que las matrices involucradas en este problema para redes grandes, son muy
dispersas, debido a que tienen mucha cantidad de elementos nulos.
Ejemplo. Aplicación del FLFRC sin offset a un sistema de dos nudos
Para ilustrar mejor la formulación vista anteriormente, se resolverá el siguiente problema:
El objetivo del problema es obtener la tensión compleja del nudo dos, 𝑉 ∟𝜃, sabiendo que la potencia
demandada en ese punto es de 200 MW y la reactiva demandada de 50 MW. Sabemos que la solución es 𝒰2 =0.9001 − 𝑗0.2 (problema ya resuelto en [18]).
Primero, establecemos un criterio de convergencia, ya que se trata de un proceso iterativo. Se escoge:
∆𝑝𝑘 ≤ 0.001
A continación, obtenemos los datos de la red:
𝒴 = [𝐺11 + 𝑗𝐵11 𝐺12 + 𝑗𝐵21
𝐺21 + 𝑗𝐵21 𝐺22 + 𝑗𝐵22] = [
−𝑗10 𝑗10𝑗10 −𝑗10
], donde 𝒴 es la matriz de admitancias.
Ahora, definamos el vector de estado 𝑥 y los vectores auxiliares 𝑢 e 𝑦:
𝒰1 = 1 + 𝑗0 → 𝛼1 = ln𝒰1 = 𝑎1 + 𝑗𝑏1 = 0
𝒰2 = 𝑒2 + 𝑗𝑓2 → 𝛼2 = ln𝒰2 = 𝑎2 + 𝑗𝑏2
𝑃2 =−200 𝑀𝑊
100 𝑀𝑊= −2 𝑝𝑢
𝑝 = [𝑃2
𝑄2] = [
−2−0.5
]
𝑄2 =−50 𝑀𝑊
100 𝑀𝑊= −0.5 𝑝𝑢
58
58
𝑥 = [
𝑎1
𝑎2
𝑏1
𝑏2
] → 𝑥′ = [𝑎2
𝑏2] ; 𝑦 = [
𝑈1
𝑈2
𝑈12
𝑊12
] ; 𝑢 =
[
ln 𝑈1
ln 𝑈2
𝑅𝑒(ln(𝑤12))
𝐼𝑚(ln(𝑤12))] = [
2𝑎1
2𝑎2
𝑎1 + 𝑎2
𝑏1 − 𝑏2
]
siendo 𝑤12 = 𝑈12 + 𝑗𝑊12
Cálculo de las matrices 𝐸 y 𝐶, matrices que se mantienen constante durante todo el problema:
𝑃2 = 𝐺22𝑈2 + 𝐺21𝑈21 + 𝐵21𝑊21 = {𝑈21 = 𝑈12} = 𝐺22𝑈2 + 𝐺21𝑈12 + 𝐵21𝑊12
𝑄2 = −𝐵22𝑈2 + 𝐺21𝑊21 − 𝐵21𝑈21 = {𝑊21 = −𝑊12} = −𝐵22𝑈2 − 𝐺21𝑊12 − 𝐵21𝑈12
𝑝 = 𝐸𝑦 → [𝑃2
𝑄2] = 𝐸 [
𝑈1
𝑈2
𝑈12
𝑊12
] = [0 𝐺22 𝐺21 −𝐵21
0 −𝐵22 −𝐵21 −𝐺21] [
𝑈1
𝑈2
𝑈12
𝑊12
]
luego 𝐸 = [0 0 0 −100 10 −10 0
]
Para el cálculo de la matriz 𝐶:
𝑢 = 𝐶𝑥 →
[
ln 𝑈1
ln 𝑈2
𝑅𝑒(ln(𝑤12))
𝐼𝑚(ln(𝑤12))] = 𝐶 [
𝑎1
𝑎2
𝑏1
𝑏2
] = [
2𝑎1
2𝑎2
𝑎1 + 𝑎2
𝑏1 − 𝑏2
]
luego 𝐶 = [
2 0 0 00 2 0 01 1 0 00 0 1 −1
] → 𝐶′ = [
0 02 01 00 −1
]
Una vez obtenido todo lo necesario, empieza el proceso iterativo:
Primera iteración:
Suponiendo un perfil inicial plano:
𝒰0 = [𝒰1
𝒰2] = [
1 + 𝑗01 + 𝑗0
] → 𝑦0= [
𝑈1
𝑈2
𝑈12
𝑊12
] = [
1
1
1
0
]
59
59
∆𝑝0 = 𝑝 − 𝑝0 = 𝑝 − 𝐸𝑦0 = [−2
−0.5] − 𝐸 [
1110
] = [−2
−0.5]
Como max(|∆𝑝0|) > 0.001, sigue el proceso iterativo:
∆𝑝0 = 𝐸∆𝑦0 → ∆𝑦0 = [
0−0.05
00.2
] → �̃�0 = 𝑦0 + ∆𝑦0 = [
00.95010.2
]
�̃�0 =
[
ln 𝑈1
ln 𝑈2
𝑅𝑒(ln(𝑤12))
𝐼𝑚(ln(𝑤12))] = [
0−0.05130.01960.1974
]
�̃�−1 =
[ 𝑈1 0 0 00 𝑈2 0 00 0 𝑈𝑖𝑗 −𝑊𝑖𝑗
0 0 𝑊𝑖𝑗 𝑈𝑖𝑗 ] = [
1 0 0 00 0.95 0 00 0 1 −0.20 0 0.2 1
]
Ahora, basta con resolver la siguiente ecuación, obteniendo 𝑥1′ :
[𝐸�̃�−1𝐶′]𝑥1′ = 𝐸�̃�−1�̃�0 → 𝑥1
′ = [𝑎2
1
𝑏21] = [
−0.0804−0.2174
]
Luego el vector tensión será 𝒰21 = exp(𝑎2
1 + 𝑗𝑏21) = 0.9010 − 𝑗0.1990
Una vez obtenida la solución, calculamos el residuo de potencias:
Continuando con la solución anterior:
𝒰1 = [𝒰1
1
𝒰21] = [
1 + 𝑗00.9010 − 𝑗0.1990
] → 𝑦1= [
𝑈1
𝑈2
𝑈12
𝑊12
] = [
1
0.8514
0.901
0.199
]
∆𝑝1 = 𝑝 − 𝑝1 = 𝑝 − 𝐸𝑦1 = [−2
−0.5] − 𝐸 [
10.85140.9010.199
] = [−0.01−0.004
]
Como max(|∆𝑝1|) > 0.001, continuamos con el proceso iterativo:
60
60
Segunda iteración:
∆𝑝1 = 𝐸∆𝑦1 → ∆𝑦1 = [
0−0.0004
00.001
] → �̃�1 = 𝑦1 + ∆𝑦1 = [
10.8510.9010.2
]
�̃�1 =
[
ln 𝑈1
ln𝑈2
𝑅𝑒(ln(𝑤12))
𝐼𝑚(ln(𝑤12))] = [
0−0.1613−0.08020.2184
]
�̃�−1 =
[ 𝑈1 0 0 00 𝑈2 0 00 0 𝑈𝑖𝑗 −𝑊𝑖𝑗
0 0 𝑊𝑖𝑗 𝑈𝑖𝑗 ] = [
1 0 0 00 0.851 0 00 0 0.901 −0.20 0 0.2 0.901
]
Ahora, se resuelve de nuevo la siguiente ecuación, obteniendo 𝑥2′ :
[𝐸�̃�−1𝐶′]𝑥2′ = 𝐸�̃�−1�̃�1 → 𝑥2
′ = [𝑎2
2
𝑏22] = [
−0.0812−0.2186
]
Luego el vector tensión será: 𝒰22 = exp(𝑎2
2 + 𝑗𝑏22) = 0.9001 − 𝑗0.1999
Una vez obtenida la solución, calculamos de nuevo el residuo de potencias:
Continuando con la solución anterior:
𝒰2 = [𝒰1
2
𝒰22] = [
1 + 𝑗00.9001 − 𝑗0.1999
] → 𝑦1= [
𝑈1
𝑈2
𝑈12
𝑊12
] = [
1
0.8501
0.9001
0.1999
]
∆𝑝2 = 𝑝 − 𝑝2 = 𝑝 − 𝐸𝑦2 = [−2
−0.5] − 𝐸 [
10.85010.90010.1999
] = [−0.001
0]
Como max(|∆𝑝2|) = 0.001, finaliza el proceso iterativo, por lo que podemos decir que se ha encontrado la
solución haciendo dos iteraciones completas, obteniendo la solución:
𝒰2 = 0.9001 − 𝑗0.1999 ≈ 0.9001 − 𝑗0.2
Fin de ejemplo
61
61
Demostraciones
Jacobiano para el FLF
𝐹−1 =𝜕𝑦
𝜕𝑢= [
𝑈𝑖 0 00 𝑈𝑖𝑗 −𝑊𝑖𝑗
0 𝑊𝑖𝑗 𝑈𝑖𝑗
]
Partiendo de 𝑦 = [
𝑦1
𝑦2
𝑦3
] = [
𝑈𝑖
𝑈𝑖𝑗
𝑊𝑖𝑗
] y de 𝑢 = [
𝑢1
𝑢2
𝑢3
] = [
ln𝑈𝑖
𝑅𝑒(ln𝑤𝑖𝑗)
𝐼𝑚(ln𝑤𝑖𝑗)] = [
ln𝑈𝑖
𝑅𝑒(ln(𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗))
𝐼𝑚(ln(𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗))] tenemos que:
𝐹−1 =
[ 𝜕𝑦1
𝜕𝑢1
𝜕𝑦1
𝜕𝑢2
𝜕𝑦1
𝜕𝑢3
𝜕𝑦2
𝜕𝑢1
𝜕𝑦2
𝜕𝑢2
𝜕𝑦2
𝜕𝑢3
𝜕𝑦3
𝜕𝑢1
𝜕𝑦3
𝜕𝑢2
𝜕𝑦3
𝜕𝑢3]
𝜕𝑦1
𝜕𝑢1=
𝜕(𝑈𝑖)
𝜕(ln𝑈𝑖) →
𝜕(ln𝑈𝑖)
𝜕(𝑈𝑖)=
1
𝑈𝑖 →
𝜕𝑦1
𝜕𝑢1= 𝑈𝑖
𝜕𝑦1
𝜕𝑢2= 0
𝜕𝑦1
𝜕𝑢3= 0
𝜕𝑦2
𝜕𝑢1= 0
𝜕𝑦2
𝜕𝑢2=
𝜕(𝑈𝑖𝑗)
𝜕(𝑅𝑒(ln(𝑈𝑖𝑗+𝑗𝑊𝑖𝑗)))
Ayudándonos de la siguiente igualdad:
𝑢2 + 𝑗𝑢3 = ln(𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗) → 𝑒𝑢2+𝑗𝑢3 =𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗
e identificando términos (parte real igual a parte real y parte imaginaria igual a parte imaginaria):
62
62
𝜕𝑒𝑢2+𝑗𝑢3
𝜕𝑢2= 𝑒𝑢2+𝑗𝑢3 = 𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗
𝜕𝑈𝑖𝑗
𝜕𝑢2= 𝑈𝑖𝑗 =
𝜕𝑦2
𝜕𝑢2
𝜕𝑊𝑖𝑗
𝜕𝑢2= 𝑊𝑖𝑗 =
𝜕𝑦3
𝜕𝑢2
𝜕𝑒𝑢2+𝑗𝑢3
𝜕𝑢2=
𝜕(𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗)
𝜕𝑢2=
𝜕𝑈𝑖𝑗
𝜕𝑢2+ 𝑗
𝜕𝑊𝑖𝑗
𝜕𝑢2
𝜕𝑦2
𝜕𝑢3=
𝜕(𝑈𝑖𝑗)
𝜕(𝐼𝑚(ln(𝑈𝑖𝑗+𝑗𝑊𝑖𝑗)))
Ayudándonos de nuevo de la igualdad anterior:
𝑢2 + 𝑗𝑢3 = ln(𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗) → 𝑒𝑢2+𝑗𝑢3 =𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗
e identificando términos (parte real igual a parte real y parte imaginaria igual a parte imaginaria),
𝜕𝑒𝑢2+𝑗𝑢3
𝜕𝑢3= 𝑗𝑒𝑢2+𝑗𝑢3 = −𝑊𝑖𝑗 + 𝑗𝑈𝑖𝑗
𝜕𝑈𝑖𝑗
𝜕𝑢3= −𝑊𝑖𝑗 =
𝜕𝑦2
𝜕𝑢3
𝜕𝑊𝑖𝑗
𝜕𝑢3= 𝑈𝑖𝑗 =
𝜕𝑦3
𝜕𝑢3
𝜕𝑒𝑢2+𝑗𝑢3
𝜕𝑢3=
𝜕(𝑈𝑖𝑗 + 𝑗𝑊𝑖𝑗)
𝜕𝑢3=
𝜕𝑈𝑖𝑗
𝜕𝑢3+ 𝑗
𝜕𝑊𝑖𝑗
𝜕𝑢3
𝜕𝑦3
𝜕𝑢1= 0
Por lo que finalmente queda demostrado que:
𝐹−1 =𝜕𝑦
𝜕𝑢=
[ 𝜕𝑦1
𝜕𝑢1
𝜕𝑦1
𝜕𝑢2
𝜕𝑦1
𝜕𝑢3
𝜕𝑦2
𝜕𝑢1
𝜕𝑦2
𝜕𝑢2
𝜕𝑦2
𝜕𝑢3
𝜕𝑦3
𝜕𝑢1
𝜕𝑦3
𝜕𝑢2
𝜕𝑦3
𝜕𝑢3]
= [
𝑈𝑖 0 00 𝑈𝑖𝑗 −𝑊𝑖𝑗
0 𝑊𝑖𝑗 𝑈𝑖𝑗
]
63
63
Jacobiano para el FLFRC
𝐹−1 =𝜕𝑦
𝜕𝑢=
[ 𝑒𝑖
′2 0 0 0 0 0 0 0
0 𝑓𝑖′2 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 𝑓𝑖′𝑓𝑗
′ 0 0
0 0 0 0 0 0 𝑒𝑖′𝑓𝑗
′ 0
0 0 0 0 0 0 0 𝑓𝑖′𝑒𝑗
′]
Partiendo de 𝑦 =
[ 𝑦1
𝑦2
𝑦3
𝑦4
𝑦5
𝑦6
𝑦7
𝑦8]
=
[ 𝑒𝑖
′2
𝑓𝑖′2
𝑒𝑖′
𝑓𝑖′
𝑒𝑖′𝑒𝑗
′
𝑓𝑖′𝑓𝑗
′
𝑒𝑖′𝑓𝑗
′
𝑓𝑖′𝑒𝑗
′]
y de 𝑢 =
[ 𝑢1
𝑢2
𝑢3
𝑢4
𝑢5
𝑢6
𝑢7
𝑢8]
=
[ ln 𝑒𝑖
′2
ln 𝑓𝑖′2
ln 𝑒𝑖′
ln 𝑓𝑖′
ln 𝑒𝑖′𝑒𝑗
′
ln 𝑓𝑖′𝑓𝑗
′
ln 𝑒𝑖′𝑓𝑗
′
ln 𝑓𝑖′𝑒𝑗
′]
tenemos que:
𝜕𝑦𝑖
𝜕𝑢𝑖= (
𝜕𝑢𝑖
𝜕𝑦𝑖)−1 = (
1
𝑦𝑖)−1 = 𝑦𝑖 𝑖 = 1,2, … ,8
La expresión anterior es válida gracias a que todas las derivadas parciales son de la forma 𝜕𝑥
𝜕 ln𝑥
𝜕𝑦𝑖
𝜕𝑢𝑗= 0
La expresión anterior es válida para 𝑖, 𝑗 = 1,2,… ,8 con 𝑖 ≠ 𝑗 .
Por lo que finalmente queda demostrado la expresión final del jacobiano para el FLFRC.
64
64
Otras figuras
A continuación se muestran las mismas figuras expuestas en la Sección 6 pero aplicadas a otras redes de 118,
300 y 1354 nudos:
Figura 27. Evolución del residuo para el caso base de una red de 118 nudos
Figura 28. Evolución del residuo para el caso base de una red de 300 nudos
65
65
Figura 29. Evolución del residuo para el caso base de una red de 1354 nudos
Figura 30. Aplicación de distintos perfiles ante la evolución de la carga para la red de 118 nudos
66
66
Figura 31. Aplicación de distintos perfiles ante la evolución de la carga para la red de 300 nudos
Figura 32. Aplicación de distintos perfiles ante la evolución de la carga para la red de 1354 nudos
67
67
Figura 33. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 118 nudos
Figura 34. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 300 nudos
68
68
Figura 35. Iteraciones frente a la variación de la semilla imaginaria y la carga en la red de 1354 nudos
Figura 36. Iteraciones frente a la variación del offset y la carga en la red de 118 nudos
69
69
Figura 37. Iteraciones frente a la variación del offset y la carga en la red de 300 nudos
Figura 38. Iteraciones frente a la variación del offset y la carga en la red de 1354 nudos
70
70
REFERENCIAS
[1] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Introducción,» de Análisis y operación de sistemas
de energía eléctrica, McGraw-Hill, 2002, pp. 139-140.
[2] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Método de Newton-Raphson,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 147-152.
[3] [En línea]. Available: http://www.ee.washington.edu/research/pstca/. [Último acceso: 2015].
[4] [En línea]. Available: http://www.pserc.cornell.edu//matpower/.
[5] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Formulación del problema,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 140-143.
[6] A. Gómez Expósito y F. L. Alvarado, «Control de frecuencia y de tensiones,» de Análisis y operación de
sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 217-259.
[7] J. L. Martínez Ramos y V. H. Quintana, «Operación del sistema de transporte,» de Análisis y operación
de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 311-381.
[8] A. Fernández Otero y J. Cidrás Pidre, «Elementos de los sistemas de energía eléctrica,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 69-134.
[9] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Métodos iterativos simples,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 144-147.
[10] B. Stott, «Effective Starting Process for Newton Raphson Load Flows,» de IEEE Proceedings, vol.118,
1971, pp. 983-987.
[11] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Ajuste de límites y reguladores,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 158-161.
[12] A. Gómez Expósito y F. L. ALvarado, «Flujo de cargas, Formulación polar,» de Análisis y operación de
sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 148-150.
[13] W. F. Tinney y C. E. Hart, «Power Flow Solution by Newton's Method,» de IEEE Transactions on Power
Apparatus and Systems, vol. PAS-86, 1967, pp. 1449-1460.
[14] J. W. Walker y W. F. Tinney, «Direct Solutions of Sparse Network Equations by Optimally Ordered
Triangular Factorization,» de Proceedings of the IEEE, vol. 55, 1967, pp. 1801-1809.
[15] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Sistema de gran dimensión,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 163-169.
71
71
[16] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Formulación cartesiana,» de Análisis y operación
de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 151-152.
[17] B. Stott, «Review of Load Flow Calculation Methods,» de Proceedings of the IEEE, vol. 62, 1974, pp.
916-929.
[18] A. Gómez Expósito, C. Gómez Quiles y W. Vargas, «Factored Solution of Infeasible Load Flow Cases,»
In Power Systems Computation Conference (PSCCT), (pp. 1-7). IEEE., 2014.
[19] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Flujo de cargas en continua,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 156-157.
[20] A. Gómez Expósito y F. L. Alvarado, «Flujo de cargas, Método desacoplado rápido,» de Análisis y
operación de sistemas de energía eléctrica, McGraw-Hill, 2002, pp. 152-155.
[21] A. Gómez Expósito, A. Abur, A. de la Villa Jaén y C. Gómez Quiles, «A Multilevel State Estimation
Paradigm for Smart Grids,» de Proceedings of the IEEE, vol 99 (6), 2011, pp. 952-976.
[22] A. Gómez Expósito y C. Gómez Quiles, «Factorized Load Flow,» de IEEE Transactions on Power
Systems, vol. 28 (4), 2013, pp. 4607-4614.
[23] A. Gómez Expósito, «Factored Solution of Nonlinear Equation Systems,» Proceedings of the Royal
Society - A..
[24] F. Milano, «Power Flow Analysis,» de Power System Modelling and Scripting, pp. 79-82.
[25] W. Vargas, C. Gómez Quiles y A. Gómez Expósito, «Fast Computation of Maximum Loading Points via
the Factored Load Flow,» enviado a IEEE Transactions on Power Systems, actualmente en revisión.
[26] «Wikipedia,» [En línea]. Available: https://es.wikipedia.org/wiki/Logaritmo_complejo.
[27] [En línea]. Available: http://www.medwelljournals.com/fulltext/?doi=ijepe.2009.276.288.