Método DE Euler EN Matlab PDF

Preview

Summary

metod Euler...

Description

ESTABILIDAD EN UN SEP MEDIANTE MÉTODO DE EULER EN MATLAB Chicaiza Quisaguano Gilson Eduardo e-mail: [email protected]

Ramos Borja Jessenia Elizabeth e-mail: [email protected]

In addition, it contains the calculations to analyze the stability behavior of a circuit that suffers a failure, in addition a program (MATLAB) was implemented, to analyze the solution in different points showing its approximation to the true value.

RESUMEN: El presente documento tiene por objetivo el estudio de la estabilidad de los sistemas eléctricos de potencia, analizando la estabilidad angular, partiendo de las formas de solución de una ecuación de oscilación que van más allá de una solución explícita, con la solución de una formula dada y aunque es posible demostrar de forma abstracta que la mayoría de estas ecuaciones posen solución, por lo menos de manera local, en general resulta muy difícil expresar explícitamente de que solución se trata; además de conocimientos básicos de circuitos para conocer la estabilidad del sistema y como poder realizar la compensación del mismo con el método de Euler. Principalmente se debe saber que las ecuaciones que se deben analizar no poseen solución de forma cerrada. Estos métodos son el método de Euler que, aunque presente poca precisión, permite hallar una solución aproximada, es decir está más próxima al valor verdadero, sin embargo, también cuenta con error. Además, contiene los cálculos para analizar el comportamiento de estabilidad de un circuito que sufre una falla, además Se implementó un programa (MATLAB), para analizar la solución en diferentes puntos mostrando su aproximación al valor verdadero.

KEY WORDS: Electric Power System, angular stability, oscillation equation, euler method. I. INTRODUCCIÓN

L

a estabilidad de un sistema de potencia, es la capacidad del mismo de volver a unas condiciones iniciales de operación, retornando a un estado de operación aceptable luego de experimentar una perturbación (cortocircuitos, salida de líneas/ transformadores, perdida de generación, perdida o aumento de demanda), debido a las grandes variaciones que se pueden presentar en los sistemas potencia es necesario el estudio de estabilidad. Las características SEP (sistemas eléctricos de potencia) son variantes en el tiempo y oscilatorios, donde las ecuaciones diferenciales que los rigen son no lineales. El método de análisis es estático en régimen transitorio con análisis en diferentes instantes de tiempo para comprender la estabilidad del sistema. El análisis se complementa con un método gráfico bastante simple para poder observar, en el tiempo, el efecto de una falla en el sistema. La inestabilidad del sistema afecta principalmente los generadores, a pesar de la actuación de los elementos de protección, los cambios que pueden causar la variación del rotor. Existen dispositivos o controladores encargados de realizar la compensación. Los efectos que se pueden tener con la perdida de estabilidad o inestabilidad, en el sistema de potencia o en la máquina síncrona de generación, se puede clasificar de distintas formas. Si se analiza el generador como una máquina síncrona que oscila en un punto de equilibrio, se está hablando de que esta máquina es estable. Pero cuando la máquina sufre un cambio en su carga, puede causar un desequilibrio, dependiendo si la modificación es excesiva o brusca, ya que puede causar la pérdida del sincronismo y sobrepasar los límites de estabilidad. En caso de que el sistema pueda recuperarse ante la presencia de una perturbación mientras el sistema no supere los límites de estabilidad del mismo, el sistema puede oscilar en otro punto de equilibrio, mientras esté dentro de los rangos aceptables de tensión, ángulo y frecuencia.

PALABRAS CLAVE: Sistema Eléctrico de Potencia, estabilidad angular, ecuación de oscilación, método de euler.

ABSTRACT: The present document aims to study the stability of electric power systems, analyzing the angular stability, starting from the solution forms of an oscillation equation that go beyond an explicit solution, with the solution of a given formula and although it is possible to demonstrate in an abstract way that most of these equations have a solution, at least locally, in general it is very difficult to express explicitly what solution is involved; In addition to basic knowledge of circuits to know the stability of the system and how to make the compensation of the same with the Euler method. Mainly, it must be known that the equations that must be analyzed do not have a closed form solution. These methods are the Euler method that, although present little precision, allows finding an approximate solution, that is, it is closer to the true value, however, it also has error.

1

II.

OBJETIVOS

Dos tipos de inestabilidad pueden ocurrir en la solución de ecuaciones diferenciales ordinarias:  Inestabilidad inherente  Inestabilidad inducida.

A. OBJETIVO GENERAL Analizar la Estabilidad de los Sistemas Eléctricos de Potencia, mediante la utilización del método de Euler por medio de iteraciones se puede encontrar la estabilidad del ángulo de la Máquina Síncrona y la para lo cual se realizará una programación en el Software Matlab el cual servirá de ayuda para simplificar los cálculos de dicho sistema.

La instabilidad inherente surge cuando errores numéricos son amplificados, a cada paso de la solución, hasta dominar completamente el cálculo y hacer el método divergencia de la solución real. Ya la inestabilidad inducida está relacionada con el método de integración utilizado o, de manera más precisa, con la discretización resultante de la aplicación del método, que depende, también, del paso de integración escogido. Se puede mostrar que los métodos de integración corresponden a resolver una ecuación algebraica de la forma

B. OBJETIVOS ESPECIFICOS   

Definir el método de Euler y su aplicación en los sistemas eléctricos de potencia Realizar la programación del cálculo de la estabilidad de ángulo de una máquina síncrona mediante el Software Matlab. Interpretar los resultados obtenidos y analizar la gráfica de estabilidad del sistema eléctrico de potencia en el régimen transitorio de la máquina. III.

∝0 y n+1 + ∝1 y n + ∝2 y n−1 +…+∝k y n−k +1+ β 0

Donde αi y βi son constantes y se desea determinar yn+1, conociendo los valores anteriores de yi (i < n+1). Esta es una ecuación discreta, lineal, y será estable si todas las raíces de su polinomio característico tuviera módulo menor que 1. El mayor problema está en el hecho que cuanto mayor es la precisión del método, menor será su estabilidad. Esto es más grave para los métodos en que el método o paso de integración es crítico para la solución, en especial los métodos de RungeKutta. Estos métodos son estables sólo para algunos valores de ∆t y son dichos condicionalmente estables. Métodos con margen de estabilidad infinita (independiente de ∆t) son llamados métodos los- estables y los métodos de Euler Reverso y Trapezoidal están en esta categoría. RIGIDEZ Un sistema de ecuaciones diferenciales es dicho rígido (“stiff”) cuando la razón entre la mayor y la menor constantes de tiempo del problema sea muyo mayor que 1. Métodos de integración tradicionales tienen su intervalo de integración definido por la menor constante de tiempo y el tiempo final de la simulación es definido por la mayor constante de tiempo. De esta forma, un gran número de puntos deberán ser determinados, aumentando sobremanera el tiempo de computación. Además de eso, si el método de integración no fuera La-estable, puede haber problemas de inestabilidad numérica. MÉTODOS DE INTEGRACIÓN NUMÉRICA Los métodos de integración numérica, como fue visto en la sección anterior, representan una discretización de la solución, de forma que la ecuación diferencial original (continua) es aproximada por una ecuación de diferencias (discreta) y apenas algunos valores (correspondientes a determinados instantes de tiempo) son calculados, que corresponderán a la solución aproximada obtenida. La Figura 1 ilustra este efecto, en que una función es aproximada a partir de valores discretos.

FUNDAMENTACIÓN TEÓRICA

A. MARCO TEÓRICO PRESICIÓN Esta característica dada, principalmente, por dos factores:  Error de redondeo.  Error de truncamiento. Los errores de redondeos están relacionados a la representación de los números en un computador (aritmética finita) y pueden ser minimizados utilizando doble precisión y/u otros recursos inherentes al computador siendo utilizado para resolver el problema. Y los errores de truncamiento son debidos a la aproximación de la solución real del problema utilizada por el método de integración escogido. El error de truncamiento puede ser analizado a partir de la serie de Taylor de la función, mostrada en la ecuación , y será proporcional a

(∆ t)

p +1

dy dy n+1 +β 1 n + dt dt

(1)

Donde ∆t es el paso de integración escogido y p es el orden de la serie de Taylor utilizada por el método de integración como aproximación de la solución. La solución verdadera del problema, en un instante dado de tempo tn, será dada, por lo tanto, por

y ( tn ) =Y n +O ( ∆ t p+1 ) + ε n (2) Donde Yn es la aproximación calculada por el método de integración, O ( ∆ t p +1 ) es la orden grandeza de la precisión del método y ɛn representa los demás errores que pueden surgir en el proceso. Estabilidad:

2

Gráfico 1: Discretización de una Función Continua

Gráfico 2: Interpretación Gráfica del Método de Euler

Fuente: (Grainger & Stevenson, 1998)

Fuente: (Grainger & Stevenson, 1998)

Para ejemplificar la aplicación de los diversos métodos de integración numérica, será utilizado el siguiente sistema lineal de las ecuaciones diferenciales:

Este método corresponde a la aplicación de la serie de Taylor, mostrada en la ecuación (1), aproximada sólo por los dos primeros términos. De esta forma, la precisión de este método es del orden O(∆t)2 . La discretización del sistema de ecuaciones, dada por el método de Euler, es equivalente a la ecuación (6) haciendo

{

x´ 1=x 2 x1 (0)= x10 (7) ´x 2=−a1 x 1−a2 x 2 x2 ( 0) =x20

{

cuya solución analítica es de la forma

{

∝ 0=1∝1=−1 β 0=0 β1 =−∆ t (8) α2=α 3=…= α k =0 β 2=β 3=…= β k =0

ℷ 2 x 10 −x 20 ℷt x 20 −ℷ 1 x 10 ℷt e + e ℷ ℷ 2 −ℷ 2 −ℷ 1 1 (4) x 20−ℷ ℷ 2 x 10 −x 20 ℷt 1 x 10 ℷt x 1 ( t ) =ℷ e e +ℷ 1 2 ℷ ℷ 1 1 2 −ℷ 2 −ℷ x1 ( t ) =

2

1

2

1

La aplicación del método de Euler para la solución del sistema dado por la ecuación (7) resulta en

{

Donde

−a2 ± √ a2 −4 a1 (5) ℷ 1,2 = 2

El método de Euler puede ser implementado de forma explícita o implícita. En su forma explícita, el valor de las derivadas mostradas en la ecuación (12) son explícitamente calculadas a cada paso del algoritmo. Ya en el caso de haber expresiones analíticas para las derivadas, estas pueden ser sustituidas, resultando un sistema de ecuaciones que no dependen explícitamente de las derivadas. Para el ejemplo dado, la ecuación (7) suministra expresiones analíticas para las derivadas, resultando en el siguiente sistema de ecuaciones:

MÉTODO DE EULER El método de Euler es el más simple de los métodos de integración numérica y puede ser visualizado en el gráfico 3. Utilizando el valor de la derivada de la función en el instante de tiempo t = t0, se puede escribir que:

|

dx

( x2 )n+1=( x 2 )n + dt 2 ∆ t n

son las raíces del polinomio característico de la ecuación diferencial de 2do orden asociada a la ecuación (4) o, equivalentemente, son los polos de la función de transferencia de este sistema.

x (t 0 +∆ t ) =x ( t 0 ) +

( |) ( |) dx

( x1 )n+1=( x 1 )n + dt 1 ∆ t n (9)

2

dx(t) (6) dt t =t

{

[ ][

x´ 1=x 2 x1 (0)= x10 (x ) 1 ∆t ⇒ 1 n+1 = ) (x −a ∆ t 1−a x´2=−a1 x 1−a2 x 2 x2 ( 0) =x20 2 n+1 1 2∆

0

La formulación implícita permite reducir el error numérico que puede haber en la determinación de las derivadas de las funciones, que son sustituidas por expresiones analíticas equivalentes. Para el caso particular de la ecuación de oscilación se conoce que:

3

k X k = x1k = δ ( t k ) ω ( t k) x2

[ ][ ] [ ]

X' K=

x x

'k 1 'k 2

∂(δ ( tk ) ) ∂t = ∂ (ω (t k ) ) ∂t

[ ][ ω (t k )

[ ]

ω (t k ) F ( X K , t K ) = ∂ (δ ( tk ) ) = πf (Pmec + Pmax sin δ ( t k ) ) H ∂2t 2

]

Gráfico 3: Curvas Real vs Aproximada producida por las Aproximaciones del Método de Euler Fuente: (Grainger & Stevenson, 1998)

Sustituyendo los respectivos términos:

[

IV.

]

ω (t k ) δ ( t k +1 ) δt = ( k ) +∆ t πf (Pmec +Pmax sin δ ( t k )) ω ( t k+1) ω (t k ) H

[ ][ ]

En forma de ecuaciones se tiene:

MARCO EXPERIMENTAL

En el marco experimental se va analizar como se realizo los calculos matemáticamente y la simulación del Método de euler en el Sofware MATLAB. EJEMPLO:

δ (t k +1 ) =δ ( tk ) + ∆ t∗ω (t k ) πf (Pmec+ P max sin δ ( t k) ) ω (t k+1 ) =ω ( t k )+ ∆ t∗( ) H

Un generador sincrónico a 60 Hz, posee una constante de inercia de H = 5MJ/MVA y una reactancia transitoria de eje directo X’d = 0.3 p.u, está conectado a una barra de potencia infinita a través de un circuito puramente resistivo. El generador entrega una potencia real de Pelec = 0.9 p.u, y Qelec = 0.084 p.u a la barra de potencia infinita cuyo voltaje es 1.0 p.u. Una falla trifásica ocurre en la mitad de una línea, y la cual es despejada por la puesta fuera de servicio de la línea, mediante la apertura simultanea de ambos extremos de la línea. La falla es despejada en 0.3 segundos. Obtener la solución numérica de la ecuación de oscilación, empleando el método de Euler, con un paso de tiempo Δt = 0.05 segundos. Efectué la simulación numérica hasta t = 0.5.

Se puede suponer:

δ (t k ) =δ k ω (t k ) = ω k De tal forma que se puede escribir en forma más compacta las ecuaciones:

δ k+1=δk +∆ t∗ωk P mec+ Pmax sin δk πf (¿ ¿ H ) ¿ ω k+1=ωk +∆ t∗¿

RESOLUCIÓN. La impedancia de transferencia entre el voltaje generador y la barra de potencia infinita Antes de que ocurra la falla, la

4

La reactancia equivalente entre el generador y la barra de potencia infinita resulta ser:

reactancia entre los puntos A y B puede ser encontrada por la combinación serie paralelo:

X 2=0,20

0,25∗0,25 =0,47 pu X 1=0.20+0.15 0,25+ 0,25

La ecuación de potencia eléctrica antes de la perturbación queda dada por:

E a∗V a sen(δ ) Xs 1.00∗1,12 2 sen(δ) Pele = 0,20 0 Pele =0,65 senδ Durante 2

Pele =

Cuando la falla es despejada, la línea fallada es puesta fuera de servicio, resultando el diagrama de reactancias como sigue. La potencia aparente es dada por: ¿

I=

S 0.9+0,084 j =0,9+0,084 j p .u ¿= 1 ∠0 V

El voltaje interno de la maquina E’ es dado por:

E' =V + j X 1 I =1.0+(0,47 j)(0.9−0.084 j)=1,12 ∠ 2

La ecuación de potencia eléctrica antes de la perturbación queda dada por: 0

Donde la reactancia equivalente entre el generador y la barra de potencia infinita resulta ser:

E a∗V a

sen(δ ) Xs 1.00∗1,12 0 Pele = sen(δ) 0,47 0 Pele =2,38 sen ( δ ) Antes Pele =

X 2=0,20+ 0.15 + 0,25=0,6 De tal modo que la potencia eléctrica después de la falla queda dada por:

E a∗V a sen(δ ) Xs 1,00∗1,12 3 Pele = sen(δ) 0,6 0 Pele =1,866 Despues 3 Pele =

Siendo el ángulo de operación estable inicial: 0

Pele =sen(δ) 0.9=2,38 sen(δ) δ=22,21 °=0,38779 rad

Se procede a plantear las ecuaciones de oscilación, de la forma:

Durante la falla se tiene, que la línea fallada se divide en dos, resultando el diagrama de impedancias como en la Figura siguiente.

d 2 δ πf = (Pmec−Pelec ) 2 H dt H =5 , f =60 Hz ANTES de la perturbación resulta:

d 2 δ 60 π = (0,9−2,38 senδ) 5 dt2 d2 δ =30,159−89,723 senδ dt2 DURANTE la perturbación resulta: Se efectúa una transformación de estrella a delta y se logra:

5

2

d δ

ω ( 2) δ ( 3 ) = δ ( 2 ) +0.05 π∗60∗(0.8−P maxII sin δ ( 26.38 ) ) ω ( 3) ω ( 2) 5

60 π = 5 (0,9−0 senδ)

[ ][ ] [ ][ ] [ ][ ]

dt2 2 d δ =(0,9−0 senδ) dt2

DESPUÉS la perturbación resulta: 2

d δ 60 π (0,9−1,866 senδ) = 5 dt2 d2 δ =30,159−70,346 senδ dt2

Para un tiempo t=0.2

[

ω (t k ) δ ( t k +1 ) δ( t k ) = +∆ t πf (Pmec−Pmax sin δ ( t k )) ω ( t k+1) ω (t k ) H

[ ][ ]

[ ][ ]

[

[

[ [

]

2.837∗180 / pi δ ( 4 ) = 34.665 +0.05 π∗60∗(0.8−0.65 sin δ ( 34.665) ) 2.837 ω ( 4) 5 δ ( 4 ) = 42.792 ; 3.648 ω ( 4)

[ ][ ] [ ][ ]

]

Para un tiempo t=0.25

Para un tiempo de t=0 δ (0 ) =22.14 ;ω (0 ) =0 Para un tiempo t=0.05

]

ω (3 ) δ ( 4 ) = δ ( 3) + 0.05 π∗60∗(0.8−P maxII sin δ ( 34.665 ) ) ω ( 4) ω( 3 ) 5

Considere que la falla es despejada en 0.3 segundos. Obtener la solución numérica de la ecuación de oscilación, empleando el método de Euler, con un paso de tiempo Δt = 0.05 segundos. Efectué la simulación numérica hasta t = 0.5. Las ecuaciones de oscilaciones a considerar resultan ser: Solución analítica método de euler:

[

1.926∗180 / pi δ ( 3 ) = 29.147 + 0.05 π∗60∗(0.8−0.65 sin δ ( 29.147) ) 1.926 ω ( 3) 5 δ ( 3 ) = 34.665 ; 2.837 ω ( 3)

[

ω (4 ) δ ( 5 ) = δ ( 4 ) + 0.05 π∗60∗(0.8−P maxII sin δ ( 42.792 ) ) ω ( 5) ω (4) 5

[ ][ ]

[

]

] ] ]

0 1.926∗180 / pi δ ( 1) = δ ( 0 ) +0.05 δ ( 5 ) = 42.792 + 0.05 π∗60∗(0.9−P maxII sin δ ( 22.14 π∗60∗(0.8−0.65 sin δ ( 42.792) ) ω ( 1) ω ( 0) 3.648 ω ( 5) 5 5 0 δ ( 5 ) 53.243 δ ( 1) = 22,14 + 0.05 π∗60∗(0.9−0.65sin δ ( 22.14 ; ω ( 5) = 4.324 0 ω ( 1) 5 Para un tiempo t=0.3 despeje de la falla δ ( 1) = 22.14 ; 1,2347 ω ( 1) ω( 5 ) δ ( 6 ) = δ ( 5 ) +0.05 π∗60∗(0.8−PmaxIII sin δ (5 )) Para un tiempo t=0.1 ω ( 6) ω (5) ω (1 ) 5 δ ( 2) = δ ( 1 ) +0.05 π∗60∗(0.9−PmaxII sin δ ( 22.14 4.325∗180/ pi ω ( 2) ω ( 1) δ ( 6 ) = 53.252 +0.05 5 π∗60∗( 0.8−1.4625 sin δ ( 53.252 ) 4.325 ω ( 6) 5 1,2347∗180 / pi δ ( 2) = 26.38 +0.05 π∗60∗( 0.8−0.65 sin δ ( 26.3 δ ( 6 ) = 65.641 1,2347 ; ω ( 2) 3.624 5 ω ( 6)

[ ][ ]

[ ][ ] [ ][ ]

[ ][ ]

;

[ ][ ] [ ][ ]

[

[ ][ ] [ ][ ]

[

[ ][ ]

[ ][ ] [ ][ ]

[

δ ( 2) = 29.147 1.926 ω ( 2)

[

[

]

Este método se repite hasta poder observar la estabilidad o inestabilidad del sistema mediante la curva angula versus tiempo.

Para un tiempo t=0.15

Cuadro de valores:

6

δ

Pm

Pe

0.00

0.9

1.8000

22.14

0

0.05

0.9

0.6500

22,14

1,2347

0.10

0.9

0.6500

29.147

1.926

0.15

0.9

0.6500

34.665

2.837

0.20

0.9

0.6500

42.792

3.648

0.25<...