Análisis de bifurcaciones en sistemas eléctricos de potencia

Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Análisis de bifurcaciones en sistemas eléctricos de potencia Autor: Manuel Cobacho Rodríguez Tutor: Francisco Gordillo Álvarez Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 iii Trabajo Fin de Máster Máster Universitario en Ingeniería Industrial Análisis de bifurcaciones en sistemas eléctricos de potencia Autor: Manuel Cobacho Rodríguez Tutor: Francisco Gordillo Álvarez Catedrático de Universidad Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 v Trabajo Fin de Máster: Análisis de bifurcaciones en sistemas eléctricos de potencia Autor: Manuel Cobacho Rodríguez Tutor: Francisco Gordillo Álvarez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2022 El Secretario del Tribunal vii A mi familia A mis seres queridos ix Agradecimientos Quiero comenzar agradeciendo a mi familia y seres queridos por todo el apoyo y el cariño mostrados siempre, y en especial a mi madre y mi hermano, que han estado siempre a mi lado apoyándome en los momentos más delicados. Agradecerle a Francisco, mi tutor del TFM, su ayuda y sus consejos a lo largo de la realización de este trabajo, y por todos los conocimientos adquiridos gracias a él. Por último, dar las gracias también a todos los compañeros y amigos que he conocido y me han acompañado a lo largo de mi trayectoria académica. Manuel Cobacho Rodríguez Sevilla, 2022 ÍNDICE DE TABLAS Tabla 3-1. Valores de los parámetros del sistema. 21 xvii ÍNDICE DE FIGURAS Figura 2-1. Microrred CC típica con fuentes y cargas [2]. 3 Figura 2-2. Representación de f(x) asociada a un sistema dinámico y retrato de estados correspondiente [3]. 5 Figura 2-3. Ejemplo de diagrama de bifurcaciones [3]. 6 Figura 2-4. Caso Bifurcación Silla-Nodo para C<0 [3]. 7 Figura 2-5. Caso Bifurcación Silla-Nodo para C=0 [3]. 8 Figura 2-6. Caso Bifurcación Silla-Nodo para C>0 [3]. 8 Figura 2-7. Caso Bifurcación Tridente para C>0 [3]. 9 Figura 2-8. Caso Bifurcación Tridente para C=0 [3]. 9 Figura 2-9. Caso Bifurcación Tridente para C<0 [3]. 9 Figura 2-10. Diagrama de fases de la bifurcación de Hopf en función del parámetro λ [3]. 11 Figura 2-11. Diagrama de bifurcaciones de la bifurcación de Hopf supercrítica [3]. 11 Figura 2-12. Diagrama de bifurcaciones de la bifurcación de Hopf subcrítica [3]. 12 Figura 2-13. Ejemplo de rutina FUNC. 13 Figura 2-14. Ejemplo de rutina STPNT y otras rutinas vacías. 14 Figura 2-15. Ejemplo de rutina BCND. 14 Figura 2-16. Ejemplo de rutina ICND. 14 Figura 2-17. Ejemplo de rutina FOPT. 15 Figura 2-18. Ejemplo de fichero de constantes. 16 Figura 2-19. Ejemplo de fichero principal. 17 Figura 2-20. Diagrama de bifurcaciones asociado al ejemplo de la figura 19. 17 Figura 3-1. Esquema microrred de partida [4]. 19 Figura 3-2. Modelado de cargas activas [4]. 20 Figura 3-3. Modelo de 2 fuentes (izquierda) y modelo equivalente (derecha) [4]. 21 Figura 4-1. Representación de las soluciones estacionarias del modelo 1. 12 Figura 4-2. Detección de la bifurcación silla-nodo (LP). 13 Figura 4-3. Detección de la bifurcación de Hopf (HB). 13 Figura 4-4. Diagrama de bifurcaciones del modelo 1. 14 Figura 4-5. Plano de fases para P=30kW. 15 Figura 4-6. Evolución de las variables de estado con el tiempo para P=30kW. 15 Figura 4-7. Plano de fases para P=45kW. 16 Figura 4-8. Evolución de las variables de estado con el tiempo para P=45kW. 16 Figura 4-9. Trayectoria de un punto situado fuera del ciclo límite. 17 Figura 4-10. Evolución de las variables de estado con el tiempo para un punto situado fuera del ciclo límite. 17 Figura 4-11. Plano de fases para una potencia entre PII y PI. 18 Figura 4-12. Plano de fases para una potencia superior a PI. 19 xix Figura 4-13. Plano de fases para potencias bajas (P=500W). 19 Figura 4-14. Evolución de las variables de estado con el tiempo para P=500W. 20 Figura 5-1. Representación de las soluciones estacionarias del modelo 2. 22 Figura 5-2. Detección de la bifurcación de Hopf (HB). 22 Figura 5-3. Representación del ciclo límite inestable del modelo 2. 23 Figura 5-4. Diagrama de bifurcaciones del modelo 2. 23 Figura 5-5. Plano de fases para P=25kW (región estable). 24 Figura 5-6. Evolución de las variables de estado del sistema para P=25kW. 25 Figura 5-7. Trayectoria del estado del sistema sometido a una pequeña perturbación en la región PU. 25 Figura 5-8. Evolución de las variables de estado para una pequeña perturbación en la región PU. 26 Figura 5-9. Evolución del estado del sistema tras cruzar el ciclo límite inestable. 26 Figura 5-10. Evolución de las variables de estado tras cruzar el ciclo límite inestable. 27 Figura 5-11. Plano de fases para P=75kW (entre PII y PI). 27 Figura 5-12. Evolución del estado del sistema al ser empujado fuera del equilibrio para P=75kW. 28 Figura 5-13. Plano de fases para P=160kW (entre PII y PI). 28 Figura 5-14. Evolución del estado del sistema para P=160kW. 29 Figura 5-15. Plano de fases para una potencia superior a PI. 29 1 1. INTRODUCCIÓN n este capítulo introductorio se expondrá el contexto del trabajo junto a los objetivos marcados en la ejecución de este. 1.1 Contexto En este trabajo se va a llevar a cabo un estudio del comportamiento de modelos simplificados de una microrred eléctrica mediante la conocida teoría de bifurcaciones, un campo de las matemáticas que analiza el cambio de comportamiento que tiene lugar en sistemas dinámicos no lineales cuando se producen pequeñas variaciones en los valores de los parámetros que los definen. De esta manera, se analizarán los distintos modos de comportamiento del sistema a través de modelos matemáticos sencillos, utilizando lo que se conoce como diagramas de bifurcaciones para posteriormente, en base a estos, realizar simulaciones que nos permitan explicar lo que ocurre físicamente en cada región diferenciada de dichos diagramas. 1.2 Motivación y objetivos En este apartado vamos a desarrollar cuales son los objetivos que se quieren alcanzar con la realización de este trabajo fin de máster. Por un lado, se pretende llevar a la práctica la teoría de bifurcaciones, un campo de las matemáticas que nos resulta muy útil a la hora de estudiar el comportamiento de los sistemas dinámicos no lineales. Con ello, se pretende analizar algunas de las bifurcaciones que se presentan de forma práctica más frecuentemente en los sistemas dinámicos no lineales mediante su representación en el denominado ‘diagrama de bifurcaciones’. Por otro lado, una vez representado el diagrama de bifurcaciones, el objetivo será analizar el cambio de comportamiento que se produce en el sistema como consecuencia de dichas bifurcaciones, familiarizándonos con conceptos muy interesantes como equilibrio estable e inestable, ciclo límite, retrato de estados, etc. 1.3 Organización del documento Una vez que se ha puesto en contexto el trabajo y se han establecido los objetivos de este, es importante clarificar el contenido de los siguientes capítulos del documento: • En el capítulo 2 se va a desarrollar una serie de conceptos básicos para entender de forma clara el contexto del trabajo y se explicará de forma superficial el funcionamiento de la herramienta de software utilizada. • A continuación, en el capítulo 3 se describirá el sistema objeto de estudio. • Posteriormente, los capítulos 4 y 5 estarán dedicados a los dos modelos utilizados para representar el sistema del capítulo 3. E Introducción 2 • Finalmente, en el capítulo 6 se expondrán las principales conclusiones de los resultados obtenidos para los modelos de los capítulos 4 y 5. Cabe destacar que se han incluido unos anexos para plasmar el código utilizado en AUTO y los diagramas de bloques construidos en simulink para la realización de simulaciones. 3 2. CONCEPTOS BÁSICOS Y HERRAMIENTAS na vez que se ha dado una visión global del propósito de este trabajo y cómo se va a organizar, en este capítulo se va a exponer una serie de conceptos que es de obligado cumplimiento conocer para proceder a la realización de dicho trabajo. Para comenzar veremos lo que es una microrred CC y porqué están jugando un papel cada vez más importante en lo referente a la generación de energía. A continuación, explicaremos de forma básica pero suficiente una serie de conceptos relativos a la teoría de bifurcaciones, ya que la iremos aplicando continuamente a lo largo de este trabajo. Para terminar, comentaremos de manera superficial el funcionamiento de la herramienta de software que utilizaremos de apoyo para representar los diagramas de bifurcaciones que iremos graficando según avanza el trabajo. 2.1 Las microrredes de corriente continua (CC) Las redes de corriente continua (CC) están integradas en una multitud de sistemas autónomos y han ido ganando terreno en el campo de la distribución de energía en forma de pequeñas microrredes. Este resurgimiento se ha producido debido, por una parte, a la fácil integración de tecnología renovable, y por otra, al creciente uso de cargas electrónicas en el sistema, lo cual hace de la distribución en continua una opción interesante para crear sistemas más eficientes [1]. Figura 2-1. Microrred CC típica con fuentes y cargas [2]. U Conceptos básicos y herramientas 4 Una microrred CC (véase Figura 1) es un sistema generación de energía local cercano al punto de consumo que puede funcionar de manera independiente de la red de suministro de energía principal, generando y utilizando su propia energía CC. Las formas más comunes de generar energía usando microrredes CC son las placas solares, las celdas de combustible y las turbinas de viento. Utilizándose de manera correcta, todas estas unidades de una microrred CC pueden almacenar o generar energía CC, siendo directamente compatible con ordenadores, teléfonos móviles, LEDs y muchos otros dispositivos finales [1]. Un artículo de la plataforma igor-tech.com [1] recoge una serie de ventajas que tiene el uso de las microrredes CC: o Sostenibilidad: desde el punto de vista de la sostenibilidad, las microrredes CC producen un menor impacto sobre el medio ambiente que las redes de energía alterna. Un aspecto que nos permite mejorar la sostenibilidad es la autonomía del sistema de energía, ya que, al ser capaz de controlar la generación, la distribución y el uso de la energía, los propietarios de edificios pueden perseguir sus objetivos de sostenibilidad con una mayor flexibilidad. Llevando un control en generación de energía, los propietarios de edificios pueden decidir usar placas solares o turbinas eólicas para producir energía CC renovable y limpia en lugar de depender de centrales eléctricas situadas a grandes distancias que están probablemente quemando combustibles fósiles para producir energía. Al utilizar una microrred CC, la generación de energía es local, lo que supone una minimización de pérdidas a lo largo de las líneas de distribución de alrededor del 6%, disminuyendo las emisiones globales para producir energía. o Fiabilidad: una microrred CC es una herramienta indispensable para suministrar energía de forma fiable a instalaciones y edificios. Una funcionalidad determinante que hace que las microrredes DC sean extraordinariamente útiles es su modo de funcionamiento en isla, lo que significa que pueden ser desconectadas de la red principal y seguir alimentando cargas locales por sí mismas. Esta funcionalidad es útil típicamente en situaciones de emergencia, tales como catástrofes o bajo condiciones meteorológicas adversas, que dejan sin suministro a la red principal. Además de su capacidad para funcionar en isla, las microrredes CC son también más resilientes que las infraestructuras de corriente alterna. En el caso de las infraestructuras CA, hay miles de kilómetros que no están adecuadamente aclimatados. Sin embargo, una microrred CC, al ser una estructura más pequeña, puede ser adecuadamente acondicionada para soportar catástrofes o situaciones meteorológicas severas. o Integración: una microrred de continua se puede integrar perfectamente para dar energía a un edificio inteligente. Debido a que cada vez más las compañías están poniendo el foco en crear entornos de trabajo que mejoren la productividad, salud y bienestar de sus empleados, los edificios inteligentes se han convertido en una solución emergente. Dentro de estos edificios inteligentes hay muchos sensores, luces, pantallas, y otros dispositivos que forman la denominada tecnología “IoT” (Internet of Things, por sus siglas en inglés), que interactúan mutuamente para mejorar el lugar de trabajo. Con el fin de distribuir energía y datos para conectar estos dispositivos se puede establecer una infraestructura PoE (Power over Ethernet, en inglés). Con el carácter de bajo voltaje CC de PoE, la mejor manera de generar energía es con una microrred CC, ya que es un proceso de uso de corriente continua exclusivamente y puede prescindir de conversiones CA/CC innecesarias dentro del edificio. Otra integración de una microrred CC es su capacidad de convertirse en una parte del edificio inteligente. Conectando generadores de energía a una plataforma del edificio inteligente, los 5 Análisis de bifurcaciones en sistemas eléctricos de potencia administradores de las instalaciones serán capaces de distribuir energía eficientemente a los dispositivos minimizando el uso innecesario de energía. Un administrador de instalaciones será capaz también de monitorizar la generación de energía de una microrred CC, permitiendo la toma de decisiones energéticas más inteligentes en cumplimiento de los objetivos energéticos de la compañía. 2.2 La teoría de bifurcaciones Antes de seguir hay que señalar que a lo largo de este trabajo vamos a tratar con sistemas dinámicos no lineales, que serán definidos a partir de una serie de ecuaciones diferenciales, las cuales contendrán tanto parámetros como variables de estado de dicho sistema. Nuestro objetivo será realizar un análisis cualitativo de nuestro sistema, es decir, llevaremos a cabo un estudio de los modos de comportamiento de este a la vez que efectuaremos un análisis de estabilidad, indicando regiones estables e inestables y puntos de bifurcación. Antes de pasar a explicar lo que son las bifurcaciones vamos a desarrollar una serie de conceptos que es necesario conocer. Sabemos que un sistema dinámico se define por un espacio de estados X y un campo vectorial f definido sobre X. Sea 𝑥󰇗=𝑓(𝑥) un sistema dinámico tal que x  (a, b)  . En la figura 2-2 podemos observar un ejemplo de f(x) y del retrato de estados asociado. Figura 2-2. Representación de f(x) asociada a un sistema dinámico y retrato de estados correspondiente [3]. A cualquier estado inicial x(0) le corresponde una trayectoria x(t). En función del lugar que ocupe el estado inicial en el retrato de estados, este seguirá una trayectoria u otra. Las trayectorias se componen de una parte transitoria y de una permanente, a la que llamamos atractor o conjunto límite. Los atractores representan el comportamiento a largo plazo de las trayectorias. El conjunto de estados iniciales que conducen a un mismo atractor forma su cuenca de atracción. Además, las cuencas de atracción están separadas por separatrices. Debido a las no-linealidades, un sistema dinámico puede tener múltiples atractores. El espacio de estados se descompone en las cuencas de atracción de sus distintos atractores. La representación gráfica del mapa de todas las cuencas de atracción en el espacio de estados se denomina retrato de estados. Este aporta una perspectiva global de los modos de comportamiento de un sistema. Conceptos básicos y herramientas 12 la bifurcación de Hopf supercrítica, el símbolo ○ representa la amplitud del ciclo límite inestable. Figura 2-12. Diagrama de bifurcaciones de la bifurcación de Hopf subcrítica [3]. Aquí conviene analizar la estabilidad del sistema para ambos tipos de bifurcaciones de Hopf. En ambos casos existe una pérdida de estabilidad del equilibrio para λ = 0 si el parámetro crece. Para la bifurcación supercrítica el equilibrio estable se transforma en un ciclo límite estable de pequeña amplitud. El sistema se mantendría en un entorno de equilibrio y, por tanto, tendríamos una pérdida de estabilidad suave o no catastrófica. En el caso de la bifurcación subcrítica la región de atracción del punto de equilibrio se encuentra acotada por un ciclo límite inestable que se va comprimiendo conforme el parámetro λ se acerca al valor crítico, el cual desaparece al alcanzar dicho valor. A partir de ese valor, al crecer el parámetro, el sistema es lanzado fuera del equilibrio, teniendo lugar una pérdida de estabilidad dura o catastrófica. La diferencia entre los dos tipos de pérdidas de estabilidad es que, en la primera, al devolver al parámetro a valores negativos, el sistema es capaz de volver al equilibrio, y, por tanto, el sistema es controlable. No ocurre lo mismo en el segundo caso, donde si intentamos devolver al sistema a su equilibrio anterior, al haber tenido lugar una pérdida de estabilidad catastrófica, muy probablemente el sistema ya no se encuentre dentro de la cuenca de atracción de dicho equilibrio, alejándose por tanto las trayectorias de dicho equilibrio a medida que aumenta t. 2.3 El software Auto Pasamos ahora a presentar el software que utilizaremos para realizar los análisis de estabilidad y representar los diagramas de bifurcaciones de los sistemas con los que trabajaremos. Auto es un software diseñado para resolver problemas de bifurcaciones y continuación de ecuaciones diferenciales ordinarias. Este software es capaz de efectuar un análisis de bifurcaciones limitado para sistemas algebraicos pero los principales algoritmos de Auto, sin embargo, están dirigidos a la continuación de soluciones de sistemas de ecuaciones diferenciales ordinarias sujetas a determinadas condiciones de contorno (incluidas las condiciones iniciales) y a limitaciones integrales. 13 Análisis de bifurcaciones en sistemas eléctricos de potencia 2.3.1 Ficheros de entrada A la hora de comunicarse con este software el usuario debe generar varios ficheros: • Fichero de ecuaciones xxx.f90 El fichero xxx.f90 contiene las rutinas de Fortran FUNC, STPNT, BCND, ICND, FOPT y PVLS. Habrá que sustituir el nombre utilizado aquí para el archivo (xxx) por el nombre que el usuario decida. Si alguna de las rutinas a las que hemos hecho referencia no tuvieran utilidad para el problema concreto bajo estudio, se procedería dejando el cuerpo de dichas rutinas vacío de contenido. En este fichero es donde se definen las ecuaciones del sistema en estudio junto con las variables y parámetros que las gobiernan, además de condiciones de contorno, funciones objetivo y demás que formaran parte del problema a resolver. ➢ Rutina FUNC: esta rutina es utilizada para la definición de las ecuaciones que gobiernan el sistema. Además, en el cuerpo de dicha rutina se declararán las variables y los parámetros asociados a dichas ecuaciones. Figura 2-13. Ejemplo de rutina FUNC. ➢ Rutina STPNT: en ella se asignará un determinado valor a los parámetros fijos de las ecuaciones del sistema y se dará un valor inicial a las variables para que el programa inicie los cálculos a partir de un punto inicial. Conceptos básicos y herramientas 14 Figura 2-14. Ejemplo de rutina STPNT y otras rutinas vacías. ➢ Rutina BCND: en el cuerpo de esta rutina se definirán las condiciones de contorno del problema. Figura 2-15. Ejemplo de rutina BCND. ➢ Rutina ICND: la rutina ICND define las condiciones integrales del problema. Figura 2-16. Ejemplo de rutina ICND. 15 Análisis de bifurcaciones en sistemas eléctricos de potencia ➢ Rutina FOPT: esta rutina define la función objetivo para casos de problemas de optimización. Figura 2-17. Ejemplo de rutina FOPT. ➢ Rutina PVLS: esta rutina sirve para definir “medidas de solución”. • Fichero de constantes c.xxx En caso de que nuestro fichero de ecuaciones fuera xxx.f90, las constantes que definen el modo de ejecutar los cálculos deben ser guardadas en un fichero denominado c.xxx. El formato de este fichero es libre, con valores asignados a constantes (constant=value) separados por comas o espacios. Como curiosidad, los comentarios deben comenzar con el carácter “#” o “!” y terminan cuando se acaba la línea donde estamos escribiendo dicho comentario. Existen diferentes tipos de constantes en función de su finalidad, entre las que se encuentran las siguientes: o Constantes del problema Estas constantes definen las características principales del problema a resolver, tales como la dimensión del sistema de ecuaciones, el número de condiciones de contorno o el máximo número de parámetros que pueden ser utilizados en todas las rutinas del fichero de ecuaciones. o Constantes de discretización Las constantes de discretización configuran el modelo de mallado utilizado para la resolución de ecuaciones diferenciales ordinarias. Se definen parámetros como el número de intervalos de mallado o el número de puntos de colocación de Gauss por intervalo de mallado. Estos parámetros serán fijados por el usuario en función de las características de problema concreto a resolver. Conceptos básicos y herramientas 16 o Tolerancias Estas constantes definen características del procedimiento de resolución numérica de ecuaciones no lineales (ya sea el método de Newton, Chord, Muller, etc) tales como los criterios de convergencia o el número máximo de iteraciones permitidas para la localización de soluciones especiales (bifurcaciones, puntos de silla, etc). o Parámetros libres Para que un problema este bien planteado necesita de la existencia de un número determinado de parámetros libres, es decir, parámetros que tengan la posibilidad de variar conforme tiene lugar la resolución numérica de dicho problema. Por ello, se definen una serie de constantes con el fin de indicar al programa el número de parámetros libres que se van a utilizar. AUTO guarda los parámetros en un vector denominado PAR(). o Constantes de computación Son constantes que sirven para indicar al programa qué tipo de soluciones especiales queremos que detecte, el tipo de problema que el usuario quiere resolver (cálculo de soluciones estacionarias de ecuaciones diferenciales ordinarias, cálculo de soluciones periódicas, problemas de optimización, etc). Figura 2-18. Ejemplo de fichero de constantes. • Fichero principal xxx.auto Este fichero es un ‘script’ en cuyo cuerpo se ejecuta una secuencia de comandos de Python que el usuario introduce. Es el fichero principal, y desde él se ejecutan los ficheros de ecuaciones y de constantes explicados anteriormente. En él, el usuario suministra paso por paso el programa que resuelve el sistema de ecuaciones del problema. Además, el usuario puede usar el comando ‘plot’ en el fichero principal para obtener el respectivo diagrama de bifurcaciones del problema. 17 Análisis de bifurcaciones en sistemas eléctricos de potencia Figura 2-19. Ejemplo de fichero principal. Figura 2-20. Diagrama de bifurcaciones asociado al ejemplo de la figura 19. 2.3.2 Ficheros de salida Una vez ejecutado el fichero principal, AUTO nos devuelve una salida standard por consola con los puntos que forman las soluciones especiales del problema y puntos elegidos por el usuario, y, además, tres ficheros de salida, que son los siguientes: • Fichero b.xxx Este fichero contiene la información necesaria a partir de la cual se genera el diagrama de bifurcaciones del sistema en estudio. Además de soluciones especiales como puntos de bifurcación o de inicio y fin, el fichero b.xxx nos proporciona puntos de equilibrio del sistema generados a partir de pequeñas variaciones del parámetro o los parámetros libres que estemos utilizando. Conceptos básicos y herramientas 18 • Fichero s.xxx El fichero s.xxx contiene gráficos completos e información de reinicio para determinadas soluciones etiquetadas. • Fichero d.xxx El último fichero de salida nos proporciona mensajes de diagnóstico, así como el historial de convergencia, autovalores y multiplicadores de Floquet. Podemos controlar la cantidad de información que queremos que nos proporcione este fichero de diagnóstico a través de la constante de AUTO ‘IID’. La teoría de bifurcaciones es un campo de las matemáticas que se centra en el estudio de los cambios cualitativos o topológicos que se dan en los sistemas cuando al variar cuasiestáticamente uno o varios parámetros de estos, se produce un cambio significativo en el comportamiento de tales sistemas. Para la resolución de dichas ecuaciones diferenciales ordinarias (EDO), el software Auto se apoya en un conjunto de algoritmos internos que se encargan de realizar las iteraciones pertinentes y que para ello se ayuda de una serie de constantes en función de las cuales el método de iteración será uno u otro. 19 3. DESCRIPCIÓN DEL SISTEMA omo se ha comentado en el capítulo introductorio, este trabajo va a centrarse en el análisis cualitativo del comportamiento de una microrred de corriente continua. Para ello, vamos a partir de unos modelos que han sido objeto de estudio en un trabajo realizado por 4 miembros de la Universidad Federal de Santa Catarina (André Pires Nóbrega Tahim, Daniel J. Pagano, Eduardo Lenz y Vinicius Stramosk) [4]. La microrred a partir de la cual vamos a obtener los modelos con los que vamos a trabajar está representada en la figura 3-1. Se trata de una microrred que opera aislada de la red general de distribución, lo que se conoce como operación en ‘modo isla’. Figura 3-1. Esquema microrred de partida [4]. Esta microrred se suele dar en situaciones de emergencia, cuando un fallo en la red principal deja aisladas ciertas cargas del sistema, teniendo la microrred el papel de generar y consumir su propia energía, garantizando la continuidad del servicio de dichas cargas. Podemos observar cómo la microrred se divide en dos subsistemas bien diferenciados. Por un lado, el subsistema formado por todas las fuentes del sistema, y, por otro, el subsistema que reúne a todas las cargas de este. • En el subsistema formado por las fuentes, tenemos por un lado las placas fotovoltaicas junto a los reguladores de carga MPPT. Estos últimos sirven para adaptar la entrada fotovoltaica a la tensión de las baterías, lo que permite hacer funcionar la placa solar en su punto de máxima potencia para así obtener la máxima eficiencia de producción. Por otro lado, están las baterías y los convertidores de potencia. Estos últimos son utilizados para desacoplar las cargas de las fuentes en caso de perturbaciones, además de para ajustar los niveles de tensión a los valores requeridos por cada dispositivo de la red. En este subsistema encontramos también las líneas de distribución que son modeladas por una resistencia y una inductancia equivalente. C Descripción del sistema 20 • En el subsistema formado por las cargas encontramos las cargas existentes en el sistema aislado. Estas cargas necesitan unos niveles de tensión apropiados para poder funcionar, y es por ello por lo que se utilizan unos convertidores llamados ‘de punto de carga’, que lo que hacen es fijar una tensión constante a la salida. Además, en este subsistema encontramos cargas resistivas y condensadores. 3.1 Obtención de los modelos Los modelos que vamos a utilizar para llevar a cabo nuestros análisis de estabilidad y nuestras simulaciones parten del modelo del sistema de la figura 3-1. Haremos una serie de consideraciones para la construcción de nuestros modelos: ➢ Asumiendo la hipótesis de que existen condiciones climáticas constantes, el conjunto formado por las fuentes distribuidas de generación (placas solares) y sus convertidores de potencia pueden ser modelados como fuentes de potencia constante, vistos desde los terminales del bus. ➢ Las baterías y sus convertidores de potencia serán modelados como una fuente ideal de tensión en serie con una resistencia ➢ Las líneas de distribución de la energía serán modeladas como una resistencia en serie con una inductancia ➢ Los condensadores de entrada a los convertidores de punto de carga junto con un banco de condensadores conectado al bus son modelados como un único condensador equivalente. ➢ Las cargas activas que llevan aparejadas un convertidor de ‘punto de carga’ serán modeladas de dos maneras distintas: o Por un lado, el primer modelo corresponderá a una carga de potencia constante ideal. o Por otra parte, el segundo modelo de cargas activas vendrá definido como un convertidor de ‘punto de carga’ a trozos. En este caso, para tensiones por debajo de la tensión de diseño de salida del convertidor, las cargas se comportarán como cargas pasivas. Figura 3-2. Modelado de cargas activas [4]. 21 Análisis de bifurcaciones en sistemas eléctricos de potencia 3.2 Parámetros del sistema El circuito equivalente que obtenemos realizando las consideraciones del apartado anterior se encuentra reflejado en la figura 3-3. Figura 3-3. Modelo de 2 fuentes (izquierda) y modelo equivalente (derecha) [4]. Del modelo equivalente mostrado en la figura 3-3 resalta el hecho de que las cargas de potencia constante y las fuentes de potencia constante han sido reagrupadas en una sola carga de potencia equivalente P, siendo: 𝑃=𝑃𝐿+𝑃𝑆 , 𝑃𝐿≥0 𝑦 𝑃𝑆≤0 Tomaremos como parámetros del sistema aquellos que vienen definidos en el artículo que hemos tomado de apoyo para la realización de este trabajo. Rd (Ω) Ld (mH) C (mF) R (Ω) Vref (V) P (W) 0.11 0.291 1 500 380 Variable Tabla 3-1. Valores de los parámetros del sistema. Como se puede deducir de la tabla, el parámetro libre será la potencia P. 3.3 Ecuaciones y variables de estado Las ecuaciones diferenciales que definen nuestro sistema se obtienen realizando los correspondientes balances de tensión e intensidad al circuito equivalente de la figura 3-3. Realizando el balance de tensión se obtiene lo siguiente: 𝑉𝑟𝑒𝑓=𝑅𝑑𝑖𝑠+𝐿𝑑𝑑𝑖𝑠 𝑑𝑡+𝑣𝑏𝑢𝑠 Modelo Nº1 16 16 Figura 4-7. Plano de fases para P=45kW. Representando de nuevo la evolución de las variables de estado con respecto al tiempo, partiendo otra vez de un punto inicial situado en el interior del ciclo límite, observamos cómo para tiempos negativos la amplitud de las oscilaciones y, por tanto, la amplitud del ciclo límite inestable, ha disminuido con respecto a la situación anterior (para P=30kW). Además, como era de esperar, para tiempos positivos la trayectoria termina estabilizándose en el punto de equilibrio existente para P=45kW. Figura 4-8. Evolución de las variables de estado con el tiempo para P=45kW. Por otro lado, si partimos de un punto inicial situado fuera del ciclo límite podemos ver cómo las variables de estado siguen una trayectoria que se aleja del ciclo límite inestable (figura 4-9). 17 Análisis de bifurcaciones en sistemas eléctricos de potencia Figura 4-9. Trayectoria de un punto situado fuera del ciclo límite. Para visualizar esto último más claramente representaremos la evolución de las variables de estado con el tiempo (figura 4-10) y observaremos cómo para tiempos positivos la amplitud de las oscilaciones de las variables de estado aumenta, indicándonos que la trayectoria que siguen tiende a alejarse del ciclo límite. Figura 4-10. Evolución de las variables de estado con el tiempo para un punto situado fuera del ciclo límite. Las simulaciones realizadas hasta ahora se han circunscrito a la zona del diagrama de bifurcaciones donde coexisten un equilibrio estable y otro inestable. Podemos plantearnos ahora realizar simulaciones para valores de P superiores a PII. Estudiaremos primero lo que ocurre en la región situada entre PII y PI, teniendo en cuenta que en esta zona coexisten dos equilibrios inestables, como se puede ver en la figura 4-4. Modelo Nº1 18 18 Figura 4-11. Plano de fases para una potencia entre PII y PI. Examinando la figura 4-11, se puede apreciar la existencia de dos puntos de equilibrio inestable: un equilibrio de tipo fuente espiral (e2) y otro de tipo punto de silla (e1). Si partiendo del equilibrio de tipo fuente espiral empujamos al sistema fuera de dicho equilibrio, observamos cómo las trayectorias son atraídas por el punto de silla para, una vez cerca de este equilibrio, ser expulsado lejos de ambos equilibrios. Veamos ahora lo que sucede para valores de P superiores a PI, valor a partir del cual sabemos de antemano, examinando la figura 4-4, que no existen puntos de equilibrio. Para comprobar esto último vamos a utilizar una herramienta disponible en pplane8 que nos permite encontrar fácilmente los puntos de equilibrio, en caso de que los hubiera. Gracias a esta herramienta podemos representar tanto las curvas donde las derivadas de la primera variable de estado con respecto al tiempo son nulas, como las curvas donde las derivadas de la segunda variable de estado con respecto al tiempo son nulas. De esto último podemos deducir que, donde dos de estas curvas se cruzan existe un punto de equilibrio, ya que las derivadas de ambas variables son cero. 19 Análisis de bifurcaciones en sistemas eléctricos de potencia Figura 4-12. Plano de fases para una potencia superior a PI. Podemos observar en la figura 4-12 que no hay cruce entre las curvas mencionadas anteriormente y, por tanto, concluimos que no existen puntos de equilibrio en esta región. Por último, comprobaremos si para potencias bajas el ciclo límite ha desaparecido, en cuyo caso podremos deducir que esto se ha debido a la colisión de dicho ciclo límite con el equilibrio de tipo punto de silla e1. En este caso la simulación se ha realizado en Simulink, mediante la construcción del correspondiente diagrama de bloques, debido a que el programa pplane8 de Matlab presentaba inconvenientes a consecuencia del paso de la variable vbus por cero. Aplicando un valor de P=500W, obtenemos el siguiente plano de fases: Figura 4-13. Plano de fases para potencias bajas (P=500W). Modelo Nº1 20 20 Representando una vez más la evolución de las variables de estado con respecto al tiempo, se puede observar cómo, partiendo de un estado muy alejado del punto de equilibrio, la trayectoria seguida por las variables de estado tiende hacia dicho equilibrio, lo cual confirma la existencia de una bifurcación homoclina que ha hecho desaparecer el ciclo límite como consecuencia de la colisión de este último con el equilibrio inestable (punto de silla). Figura 4-14. Evolución de las variables de estado con el tiempo para P=500W. 21 Análisis de bifurcaciones en sistemas eléctricos de potencia 5. MODELO Nº2 5.1 Definición del modelo ara este segundo caso se modelarán las cargas activas como un convertidor de punto de carga, siendo la diferencia con el primer modelo que para tensiones bajas la carga se comporta como una resistencia. En este caso representaremos la carga como una función definida a trozos de la forma siguiente: 𝒊(𝒗𝒃𝒖𝒔)= { 𝑷 𝒗𝒃𝒖𝒔 𝑠𝑖 𝑣𝑏𝑢𝑠>𝑉𝑡ℎ 𝑷 𝑽𝒕𝒉 𝟐·𝒗𝒃𝒖𝒔 𝑠𝑖 𝑣𝑏𝑢𝑠≤𝑉𝑡ℎ } , siendo Vth la tensión de salida de diseño del convertidor (150V). Por tanto, las ecuaciones que definen este modelo son las siguientes: 𝒅𝒊𝒔 𝒅𝒕=𝟏 𝑳𝒅(𝑽𝒓𝒆𝒇−𝒗𝒃𝒖𝒔−𝑹𝒅𝒊𝒔) (𝟑) 𝒅𝒗𝒃𝒖𝒔 𝒅𝒕 =𝟏 𝑪(𝒊𝒔−𝒗𝒃𝒖𝒔 𝑹−𝒊(𝒗𝒃𝒖𝒔)) (𝟒), debiendo sustituir la expresión de i(vbus), definida a trozos anteriormente, en la ecuación (4) según sea el valor de vbus. 5.2 Diagrama de bifurcaciones 5.2.1 Cálculo de soluciones estacionarias Al igual que en el primer modelo, en este otro necesitamos proporcionarle a AUTO un punto de equilibrio inicial para que proceda a calcular el resto de puntos de equilibrio del sistema. En este caso no es posible utilizar el mismo punto de equilibrio que usamos en el primer modelo debido a que no se cumplen las condiciones de la función a trozos. Sin embargo, si utilizamos el equilibrio e2 del modelo anterior para P = 100kW sí se cumple la condición de la función a trozos para tensiones mayores que Vth, y por tanto podemos utilizar las mismas ecuaciones que en el modelo 1, siendo el punto de equilibrio utilizado el siguiente: 𝒆𝟐(𝑷=𝟏𝟎𝟎𝒌𝑾)=(𝐼𝑠2,𝑉 𝑏𝑢𝑠2)=(𝟐𝟖𝟕.𝟕𝟔𝟖𝑨 ,𝟑𝟒𝟖.𝟑𝟒𝟓𝑽) Comprobamos que se cumple la condición de la función a trozos 𝑉 𝑏𝑢𝑠2>𝑉𝑡ℎ. P Modelo Nº2 22 22 Figura 5-1. Representación de las soluciones estacionarias del modelo 2. Sabemos por el artículo [4] que existen dos zonas donde el sistema presenta equilibrios estables: la primera de ellas es la zona situada a la izquierda de la bifurcación de Hopf (véase figura 5-1), y la otra es la curva de color rojo que aparece para potencias altas en la figura 5-1. Si abrimos el fichero de diagnóstico generado por AUTO podemos observar que la bifurcación de Hopf que obtenemos para este segundo modelo es la misma que obtuvimos para el primero. Esto se debe a que las ecuaciones que rigen el comportamiento del sistema para tensiones por encima de Vth son las mismas que las utilizadas para el primer modelo. Figura 5-2. Detección de la bifurcación de Hopf (HB). 5.2.2 Cálculo de soluciones periódicas Como ya hicimos con el modelo anterior, en este otro vamos a realizar un cálculo de las soluciones periódicas que aparecen a raíz de la bifurcación de Hopf. Para ello, vamos a dibujar el ciclo límite que parte de dicha bifurcación. Le indicaremos al programa principal de AUTO que realice un nuevo “run” a partir de la bifurcación de Hopf, suministrándole como parámetros libres la potencia P y el período de las órbitas periódicas. Con ello, obtenemos el ciclo límite que se puede visualizar en la figura 5-3. 23 Análisis de bifurcaciones en sistemas eléctricos de potencia Figura 5-3. Representación del ciclo límite inestable del modelo 2. En este caso, al igual que ocurría con el primer modelo, aparece un ciclo límite inestable que parte de la bifurcación de Hopf. Aparentemente, el ciclo límite desaparece, lo cual resulta extraño al no haber indicios de que haya una bifurcación u otro fenómeno que lo haga desaparecer. Si ojeamos el artículo de la referencia [1], lo que sucede es que tiene lugar una bifurcación de tipo silla-nodo de órbitas periódicas, lo cual resuelve el misterio del final extraño de nuestro ciclo límite. Y que el programa no reconozca la bifurcación se debe a que AUTO no es capaz de detectar todos los tipos de bifurcaciones existentes, siendo esta una de ellas. Esta bifurcación, según podemos ver en [4], es el origen de un ciclo límite estable. Para poder representarlo hemos tenido que pasarle al programa AUTO un fichero de datos con la información de un período completo de una de las órbitas generadas dentro de dicho ciclo límite. Con ello, nos queda el diagrama de bifurcaciones final representado en la figura 5-4. Figura 5-4. Diagrama de bifurcaciones del modelo 2. Hay varias cosas a destacar en el diagrama de bifurcaciones de la figura 5-4. Para empezar, el diagrama posee una región estable para potencias bajas (potencias por debajo de la bifurcación silla-nodo de órbitas periódicas). Esto es, si una perturbación saca al sistema de su estado de equilibrio, este siempre vuelve al equilibrio. Por otra parte, si nos encontramos en la zona nombrada en el diagrama como “PU”, aunque el equilibrio es estable localmente, una perturbación en la tensión del bus puede hacer que el estado del sistema cruce el ciclo límite inestable, provocando que la trayectoria sea atraída por el ciclo límite estable y comience a oscilar en Modelo Nº2 24 24 torno a él. Como un detalle menor, se puede observar que la curva que representa el ciclo límite estable parte de la bifurcación silla-nodo con un color verde y a partir de una potencia de 150kW se vuelve de color rojo. Esto sucede porque el fichero de datos que le proporcionamos a AUTO para que calculara las soluciones periódicas encerradas por el ciclo límite estable utilizaba los datos de una simulación para un valor de la potencia de 150kW. 5.3 Simulaciones En este apartado, al igual que se procedió en el apartado de simulaciones del modelo 1, se va a representar en el plano de fases los distintos comportamientos que presenta el sistema para este segundo modelo. Lo primero que se hará será reflejar el comportamiento que presenta el sistema dentro de la región estable (véase figura 5-4). Como resulta evidente, en esta región, cualquier perturbación en el estado de la tensión devolverá al sistema de nuevo al punto de equilibrio. Figura 5-5. Plano de fases para P=25kW (región estable). Como se puede apreciar en la figura 5-5, para dos puntos iniciales distintos alejados del punto de equilibrio, la trayectoria que siguen sus órbitas termina confluyendo en el punto de equilibrio existente para P=25kW. En la figura 5-6 se puede observar de forma más evidente la trayectoria que siguen las variables de estado con el tiempo para uno de los puntos iniciales anteriores. 25 Análisis de bifurcaciones en sistemas eléctricos de potencia Figura 5-6. Evolución de las variables de estado del sistema para P=25kW. A continuación, se va a estudiar la región del diagrama de bifurcaciones cuyo análisis resulta más interesante. Esta zona, donde coexisten el ciclo límite estable y el ciclo límite inestable, se encuentra etiquetada en la figura 5-4 como PU. Partiendo de una situación en la cual el sistema se encuentra en un estado de equilibrio estable, si este es sometido a una pequeña perturbación de manera que el estado del sistema no cruza el ciclo límite inestable, este convergerá de nuevo al punto equilibrio estable. Figura 5-7. Trayectoria del estado del sistema sometido a una pequeña perturbación en la región PU. En la figura 5-7 se nos muestra la trayectoria que seguiría el estado del sistema una vez que ha sido sometido a una pequeña perturbación, retornando de nuevo al punto de equilibrio estable. Trazando la evolución de las variables de estado con el tiempo para el caso anterior, se puede apreciar de manera más clara el retorno del estado del sistema al punto de equilibrio (figura 5-8). Referencias 32 32 33 Análisis de bifurcaciones en sistemas eléctricos de potencia ANEXO 1: CODIGO UTILIZADO EN AUTO Modelo 1 ➢ Fichero principal #Calculo de las soluciones estacionarias limcyc = run(e='limitcycle',c='limitcycle') limcyc = limcyc + run(e='limitcycle',c='limitcycle',DS='-') limcyc = rl(limcyc) save(limcyc,'limcyc') #Calculo de las soluciones periodicas a partir de la bifurcacion de Hopf for solution in limcyc('HB'): limcyc2=run(solution,c='limitcycle.2') save(limcyc2,'limcyc2') limcyc = limcyc + limcyc2 limcyc = rl(limcyc) save(limcyc,'limcyc') #Representacion del diagrama de bifurcaciones p = plot("limcyc") p.config(top_title='P frente a MAX Vbus (curva verde)') cl() ➢ Fichero de ecuaciones !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- ! Model 1 : Ideal CPL -> i(vbus) = P/vbus !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) ! ---------- ---- ! Evaluates the algebraic equations or ODE right hand side ! Input arguments : ! NDIM : Dimension of the ODE system ! U : State variables ! ICP : Array indicating the free parameter(s) ! PAR : Equation parameters ! Values to be returned : ! F : ODE right hand side values ! Normally unused Jacobian arguments : IJAC, DFDU, DFDP (see manual) IMPLICIT NONE INTEGER NDIM, IJAC, ICP(*) DOUBLE PRECISION U(NDIM), PAR(*), F(NDIM), DFDU(*), DFDP(*) DOUBLE PRECISION Vbus,Is,P,LD,VREF,RD,C,R Vbus=U(1) Is=U(2) P=PAR(1) R=PAR(2) VREF=PAR(3) RD=PAR(4) C=PAR(5) LD=PAR(6) F(1)=(Is - (Vbus/R) - (P/Vbus)) / C F(2)=(VREF - Vbus - RD*Is) / LD END SUBROUTINE FUNC !---------------------------------------------------------------------- !---------------------------------------------------------------------- Anexo 1: Codigo utilizado en auto 34 34 SUBROUTINE STPNT(NDIM,U,PAR,T) ! ---------- ----- ! Input arguments : ! NDIM : Dimension of the ODE system ! Values to be returned : ! U : A starting solution vector ! PAR : The corresponding equation-parameter values ! T : Not used here IMPLICIT NONE INTEGER NDIM DOUBLE PRECISION U(NDIM), PAR(*), T ! Initialize the equation parameters PAR(1)=100000.0d0 PAR(2)=500.0d0 PAR(3)=380.0d0 PAR(4)=0.11d0 PAR(5)=0.001d0 PAR(6)=0.000291d0 ! Initialize the solution U(1)=31.57089686d0 U(2)=3167.537301d0 END SUBROUTINE STPNT !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! The following subroutines are not used here, ! but they must be supplied as dummy routines SUBROUTINE BCND END SUBROUTINE BCND SUBROUTINE ICND END SUBROUTINE ICND SUBROUTINE FOPT END SUBROUTINE FOPT SUBROUTINE PVLS END SUBROUTINE PVLS !---------------------------------------------------------------------- !---------------------------------------------------------------------- ➢ Fichero de constantes utilizado para el cálculo de las soluciones estacionarias parnames = {1:'P', 2:'R', 3:'VREF', 4:'RD', 5:'C', 6:'LD'} unames = {1:'Vbus', 2:'Is'} NDIM= 2, IPS = 1, IRS = 0, ILP = 1 ICP = [1,11] NTST= 10, NCOL= 4, IAD = 3, ISP = 1, ISW = 1, IPLT= 0, NBC= 0, NINT= 0 NMX= 250, NPR= 20, MXBF= 10, IID = 2, ITMX= 8, ITNW= 5, NWTN= 3, JAC= 0 EPSL= 1e-06, EPSU = 1e-06, EPSS =1e-4 DS = 0.01, DSMIN= 0.001, DSMAX= 10000, IADS= 1 NPAR = 6, THL = {11: 0.0}, THU = {} UZSTOP = {1 : [1000.0, 350000.0]} ➢ Fichero de constantes utilizado para el cálculo de las soluciones periódicas parnames = {1:'P', 2:'R', 3:'VREF', 4:'RD', 5:'C', 6:'LD'} unames = {1:'Vbus', 2:'Is'} NDIM= 2, IPS = 2, IRS = 2, ILP = 1 ICP = [1, 11] NTST= 20, NCOL= 4, IAD = 3, ISP = 2, ISW = 1, IPLT= 0, NBC= 0, NINT= 0 NMX= 400, NPR= 20, MXBF= 10, IID = 2, ITMX= 8, ITNW= 5, NWTN= 3, JAC= 0 EPSL= 1e-06, EPSU = 1e-06, EPSS =1e-4 DS = 0.01, DSMIN= 0.001, DSMAX= 10000, IADS= 1 NPAR = 6, THL = {}, THU = {} UZSTOP = {1 : [300.0, 85000.0]} 35 Análisis de bifurcaciones en sistemas eléctricos de potencia Modelo 2 ➢ Fichero principal #Calculo de las soluciones estacionarias limcyc = run(e='polconvbif',c='polconvbif') limcyc = limcyc + run(e='polconvbif',c='polconvbif',DS='-') limcyc = rl(limcyc) save(limcyc,'limcyc') #Calculo de las soluciones estacionarias (ciclo limite inestable) a partir de la bifurcacion de Hopf for solution in limcyc('HB'): limcyc2=run(solution,c='polconvbif.2') limcyc2 = rl(limcyc2) save(limcyc2,'limcyc2') #Calculo del ciclo limite estable a partir de los datos de un periodo obtenidos por simulacion limcyc3 = run('polconvbif_2',c='polconvbif.3') limcyc3 = limcyc3 + run(e='polconvbif_2',c='polconvbif.3', DS='-') save(limcyc3,'limcyc3') limcyc = limcyc + limcyc2 + limcyc3 limcyc = rl(limcyc) save(limcyc,'limcyc') #Representacion del diagrama de bifurcaciones p = plot("limcyc") p.config(top_title='P frente a Vbus') cl() ➢ Fichero de ecuaciones utilizado para el cálculo de las soluciones estacionarias y el ciclo límite inestable !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! POL converter represented as a piecewise (VCCS) !---------------------------------------------------------------------- !---------------------------------------------------------------------- SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) ! ---------- ---- ! Evaluates the algebraic equations or ODE right hand side ! Input arguments : ! NDIM : Dimension of the ODE system ! U : State variables ! ICP : Array indicating the free parameter(s) ! PAR : Equation parameters ! Values to be returned : ! F : ODE right hand side values ! Normally unused Jacobian arguments : IJAC, DFDU, DFDP (see manual) IMPLICIT NONE INTEGER NDIM, IJAC, ICP(*) DOUBLE PRECISION U(NDIM), PAR(*), F(NDIM), DFDU(*), DFDP(*) DOUBLE PRECISION Vbus,Is,P,LD,VREF,RD,C,R,Vth Vbus=U(1) Is=U(2) P=PAR(1) R=PAR(2) VREF=PAR(3) RD=PAR(4) C=PAR(5) LD=PAR(6) Vth=PAR(7) !Caso Piecewise ! F(1)=(Is - (Vbus/R) - (P/Vbus)) / C ! F(2)=(VREF - Vbus - RD*Is) / LD Anexo 1: Codigo utilizado en auto 36 36 IF(Vbus.GT.Vth) THEN F(1)=(Is - (Vbus/R) - (P/Vbus)) / C F(2)=(VREF - Vbus - RD*Is) / LD ELSE F(1)=(Is - (Vbus/R) - P*Vbus/(Vth*Vth)) / C F(2)=(VREF - Vbus - RD*Is) / LD END IF END SUBROUTINE FUNC !---------------------------------------------------------------------- !---------------------------------------------------------------------- SUBROUTINE STPNT(NDIM,U,PAR,T) ! ---------- ----- ! Input arguments : ! NDIM : Dimension of the ODE system ! Values to be returned : ! U : A starting solution vector ! PAR : The corresponding equation-parameter values ! T : Not used here IMPLICIT NONE INTEGER NDIM DOUBLE PRECISION U(NDIM), PAR(*), T ! Initialize the equation parameters PAR(1)=100000.0d0 PAR(2)=500.0d0 PAR(3)=380.0d0 PAR(4)=0.11d0 PAR(5)=0.001d0 PAR(6)=0.000291d0 PAR(7)=150.0d0 ! Initialize the solution U(1)=348.3455215d0 U(2)=287.7679d0 END SUBROUTINE STPNT !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! The following subroutines are not used here, ! but they must be supplied as dummy routines SUBROUTINE BCND END SUBROUTINE BCND SUBROUTINE ICND END SUBROUTINE ICND SUBROUTINE FOPT END SUBROUTINE FOPT SUBROUTINE PVLS END SUBROUTINE PVLS !---------------------------------------------------------------------- !---------------------------------------------------------------------- ➢ Fichero de ecuaciones utilizado para el cálculo del ciclo límite estable !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! Equation file for the computation of the stable limit cycle !---------------------------------------------------------------------- !---------------------------------------------------------------------- SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP) ! ---------- ---- ! Evaluates the algebraic equations or ODE right hand side ! Input arguments : ! NDIM : Dimension of the ODE system ! U : State variables ! ICP : Array indicating the free parameter(s) ! PAR : Equation parameters ! Values to be returned : ! F : ODE right hand side values ! Normally unused Jacobian arguments : IJAC, DFDU, DFDP (see manual) 37 Análisis de bifurcaciones en sistemas eléctricos de potencia IMPLICIT NONE INTEGER NDIM, IJAC, ICP(*) DOUBLE PRECISION U(NDIM), PAR(*), F(NDIM), DFDU(*), DFDP(*) DOUBLE PRECISION Vbus,Is,P,LD,VREF,RD,C,R,Vth Vbus=U(1) Is=U(2) P=PAR(1) R=PAR(2) VREF=PAR(3) RD=PAR(4) C=PAR(5) LD=PAR(6) Vth=PAR(7) !Caso Piecewise ! F(1)=(Is - (Vbus/R) - (P/Vbus)) / C ! F(2)=(VREF - Vbus - RD*Is) / LD IF(Vbus.GT.Vth) THEN F(1)=(Is - (Vbus/R) - (P/Vbus)) / C F(2)=(VREF - Vbus - RD*Is) / LD ELSE F(1)=(Is - (Vbus/R) - P*Vbus/(Vth*Vth)) / C F(2)=(VREF - Vbus - RD*Is) / LD END IF END SUBROUTINE FUNC !---------------------------------------------------------------------- !---------------------------------------------------------------------- SUBROUTINE STPNT(NDIM,U,PAR,T) ! ---------- ----- ! Input arguments : ! NDIM : Dimension of the ODE system ! Values to be returned : ! U : A starting solution vector ! PAR : The corresponding equation-parameter values ! T : Not used here IMPLICIT NONE INTEGER NDIM DOUBLE PRECISION U(NDIM), PAR(*), T ! Initialize the equation parameters PAR(1)=150000.0d0 PAR(2)=500.0d0 PAR(3)=380.0d0 PAR(4)=0.11d0 PAR(5)=0.001d0 PAR(6)=0.000291d0 PAR(7)=150.0d0 END SUBROUTINE STPNT !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! The following subroutines are not used here, ! but they must be supplied as dummy routines SUBROUTINE BCND END SUBROUTINE BCND SUBROUTINE ICND END SUBROUTINE ICND SUBROUTINE FOPT END SUBROUTINE FOPT SUBROUTINE PVLS END SUBROUTINE PVLS !---------------------------------------------------------------------- !---------------------------------------------------------------------- ➢ Fichero de constantes utilizado para el cálculo de las soluciones estacionarias parnames = {1:'P', 2:'R', 3:'VREF', 4:'RD', 5:'C', 6:'LD',7:'Vth'} unames = {1:'Vbus', 2:'Is'} NDIM= 2, IPS = 1, IRS = 0, ILP = 1 ICP = [1] NTST= 10, NCOL= 4, IAD = 3, ISP = 1, ISW = 1, IPLT= 0, NBC= 0, NINT= 0 NMX= 250, NPR= 20, MXBF= 10, IID = 2, ITMX= 8, ITNW= 5, NWTN= 3, JAC= 0 Anexo 1: Codigo utilizado en auto 38 38 EPSL= 1e-06, EPSU = 1e-06, EPSS =1e-4 DS = 0.01, DSMIN= 0.001, DSMAX= 7000, IADS= 1 NPAR = 7, THL = {11: 0.0}, THU = {} UZSTOP = {1 : [1000.0, 350000.0]} ➢ Fichero de constantes utilizado para el cálculo de las soluciones periódicas (ciclo límite inestable) parnames = {1:'P', 2:'R', 3:'VREF', 4:'RD', 5:'C', 6:'LD', 7:'Vth'} unames = {1:'Vbus', 2:'Is'} NDIM= 2, IPS = 2, IRS = 0, ILP = 1 ICP = [1, 11] NTST= 10, NCOL= 4, IAD = 3, ISP = 1, ISW = 1, IPLT= 0, NBC= 0, NINT= 0 NMX= 400, NPR= 20, MXBF= 10, IID = 2, ITMX= 8, ITNW= 5, NWTN= 3, JAC= 0 EPSL= 1e-06, EPSU = 1e-06, EPSS =1e-4 DS = 0.01, DSMIN= 0.0001, DSMAX= 7000, IADS= 1 NPAR = 7, THL = {}, THU = {} UZSTOP = {1 : [1000.0, 200000.0]} ➢ Fichero de constantes utilizado para el cálculo del ciclo límite estable parnames = {1:'P', 2:'R', 3:'VREF', 4:'RD', 5:'C', 6:'LD', 7:'Vth'} unames = {1:'Vbus', 2:'Is'} dat = 'datos2' NDIM= 2, IPS = 2, IRS = 0, ILP = 0 ICP = [1,11] NTST= 10, NCOL= 4, IAD = 3, ISP = 0, ISW = 1, IPLT= 0, NBC= 0, NINT= 0 NMX= 600, NPR= 20, MXBF= 10, IID = 2, ITMX= 8, ITNW= 7, NWTN= 3, JAC= 0 EPSL= 1e-06, EPSU = 1e-06, EPSS =1e-4 DS = 0.001, DSMIN= 0.0000001, DSMAX= 2000 , IADS= 1 NPAR = 7, THL = {11: 0.0}, THU = {} #UZSTOP = {1 : [10000.0, 85000.0]} 39 Análisis de bifurcaciones en sistemas eléctricos de potencia ANEXO 2: DIAGRAMAS DE BLOQUES (SIMULINK) Modelo 1 Modelo 2 Código implementado en el bloque ‘MATLAB function’ (fcn): function y = fcn(u1,u2,u3) if(u2>u3) y=u1/u2; else y=u1*u2/(u3^2); end