Metodo de Euler y Euler ModificadoDescripción completa
Views 235 Downloads 1 File size 5MB
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
INTRODUCCIÓN En el presente trabajo se muestra la 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. 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, y el método de euler modificado cuya precisión es mayor, es decir está más próxima al valor verdadero, sin embargo también cuenta con error. Se implementó un programa (MATLAB), para analizar la solución en diferentes puntos mostrando su aproximación al valor verdadero.
ESTABILIDAD TRANSITORIA
Página 1
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
INTEGRACIÓN NUMÉRICA Para el análisis de la estabilidad transitoria de sistemas de potencia, se debe resolver un conjunto de ecuaciones diferenciales no-lineales de 1 er orden, con condiciones iniciales conocidas, de la forma. ´x =f ( x , r , t ) … … … … … … … … .(1)
Donde: x= vector de variable de estado r= vector de variables algebraicas t= tiempo Las variables algebraicas están relacionadas las restricciones algébricas de la forma 0=f ( x , r , t ) … … … … … … … (2) que deben ser satisfechas a cada instante de tiempo. En general, ecuaciones diferenciales no-lineales sólo admiten solución numérica, i. e., calculada paso a paso. Diversos algoritmos de integración numérica son conocidos y la opción del “mejor” depende grandemente del problema. Así siendo, algunas características inherentes a los métodos de integración deben ser analizadas para llegar a una conclusión. Los principales métodos de integración numérica pueden ser interpretados como aproximaciones basadas en el truncamiento de la expresión de la solución del problema en términos de su serie de Taylor: ESTABILIDAD TRANSITORIA
Página 2
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
x ( t 0 +∆ t )=x ( t 0 ) +
|
2
|
2
3
|
3
dx (t ) d x(t ) (∆ t ) d x (t) ( ∆t ) ∆ t+ + +… (3) 2 3 dt t =t d t t =t 2! d t t=t 3 ! 0
0
0
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 (3), y será proporcional a (∆ t) p +1 … … … … … … … … … … … … … … … ..(4) 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 t n, será dada, por lo tanto, por y ( t n ) =Y n +O ( ∆ t p+1 ) +ε n … … … … … …..(5) ESTABILIDAD TRANSITORIA
Página 3
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
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: Dos tipos de inestabilidad pueden ocurrir en la solución de ecuaciones diferenciales ordinarias:
Inestabilidad inherente; y Inestabilidad inducida.
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
integración
inducida
utilizado
o,
está
relacionada
con
el
método
de
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 ∝0 y n +1+∝1 y n+ ∝2 y n−1+ …+∝k y n−k+1 + β 0
dy n+ 1 dy dy + β 1 n + …+ β k n−k +1 =0 … … … … … … … … … … … … … dt dt dt
Donde αi y βi son constantes y se desea determinar y n+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 ESTABILIDAD TRANSITORIA
Página 4
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
que el método o paso de integración es crítico para la solución, en especial los métodos de Runge-Kutta. 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 diferencias
original
(continua)
(discreta)
y
es aproximada
apenas
algunos
por
una
ecuación
de
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.
ESTABILIDAD TRANSITORIA
Página 5
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
Fig. 1: Discretización de una Función Continua
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:
{
´x 1=x 2 x 1 (0)=x 10 … … … … … … … … .(7) x´ 2=−a 1 x 1−a2 x 2 x 2 ( 0 )=x 20
cuya solución analítica es de la forma
{
ℷ2 x 10−x 20 ℷ t x 20−ℷ1 x 10 ℷ t e + e ℷ2−ℷ1 ℷ2−ℷ1 … … … … … …(8) ℷ2 x 10−x 20 ℷ t x 20−ℷ1 x10 ℷ t x 1 ( t )=ℷ1 e +ℷ2 e ℷ2−ℷ1 ℷ2−ℷ1 x 1 ( t )=
1
1
2
2
Donde −a2 ± √ a22−4 a 1 ℷ1,2= … … … … … … … … … ..(9) 2
ESTABILIDAD TRANSITORIA
Página 6
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
son las raíces del polinomio característico de la ecuación diferencial de 2 do orden asociada a la ecuación (7) o, equivalentemente, son los polos de la función de transferencia de este sistema. 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 la Figura 2. Utilizando el valor de la derivada de la función en el instante de tiempo t = t 0, se puede escribir que: x ( t 0 +∆ t )=x ( t 0 ) +
|
dx (t ) … … … … … … … … … … … … … .(10) dt t =t 0
Fig.2: Interpretación Gráfica del Método de Euler
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
ESTABILIDAD TRANSITORIA
Página 7
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
{
∝0=1∝1=−1 β 0=0 β 1=−∆t … … …(11) α 2=α 3=…=α k =0 β 2= β3 =…=β k =0
La aplicación del método de Euler para la solución del sistema dado por la ecuación (7) resulta en
{
( |) ( |) d x1
( x1 )n +1=( x 1 )n + dt ( x2 )n +1=( x 2 )n +
d x2 dt
∆t
n
… … … … … … … … … … … … … … … … … .. ( 12 ) ∆t
n
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:
[ ][
][ ]
(x ) ( x 1)n x´ 1=x 2 x 1 (0)=x 10 ∆t ⇒ 1 n+1 = 1 …(13) x´ 2=−a 1 x 1−a2 x 2 x 2 ( 0 )=x 20 ( x 2)n+1 −a1 ∆ t 1−a2 ∆ t ( x 2)n
{
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: ESTABILIDAD TRANSITORIA
Página 8
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
[ ][ ]
Xk=
k δ (t k ) x1 = k ω (t k ) x2
X ' K = x 1'k x2
[ ]
∂(δ ( t k ) ) ∂t = ∂(ω ( t k ) ) ∂t
[ ] 'k
[ ][
ω ( tk ) ω (t k ) F ( X K , t K ) = ∂ (δ ( t k ) ) = πf (Pmec + Pmax sin δ ( t k ) ) 2 H ∂ t 2
]
Sustituyendo los respectivos términos:
[
ω (t k ) δ ( t k +1) δ (t k ) = + ∆ t πf ( Pmec + P max sin δ ( t k ) ) ω ( t k+1 ) ω (t k ) H
[ ][ ]
]
En forma de ecuaciones se tiene: δ ( t k+ 1) =δ ( t k ) + ∆ t∗ω ( t k ) ω ( t k+1 ) =ω ( t k ) + ∆t∗(
πf ( P mec + P max sin δ ( t k ) ) ) H
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
ESTABILIDAD TRANSITORIA
Página 9
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
Pmec + Pmax sin δ k πf (¿ ¿ H ) ¿ ω k+1=ω k +∆ t∗¿
MÉTODO DE EULER REVERSO El método de Euler reverso (“backward Euler”) modifica el método de Euler original, utilizando el valor de la derivada en el instante de tiempo t = t 0+ ∆t en la ecuación (10): x ( t 0 +∆ t )=x ( t 0 ) +
|
dx ( t ) … … … … … … … … … … … … ( 14 ) dt t =t 0
La discretización resultante de la aplicación de este método también puede ser expresada por la ecuación (6), donde, en este caso,
ESTABILIDAD TRANSITORIA
Página 10
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
{
∝0=1∝1=−1 β 0=−∆t β1 =0 … … … … … …(15) α 2=α 3=…=α k =0 β 2= β3 =…=β k =0
Este método es más estable que el método de Euler convencional (“forward Euler”), pero su implementación requiere el cálculo del valor de la derivada en un instante de tiempo para el cual aún no se conoce el estado. La implementación explícita, en este caso, requiere una extrapolación inicial para la determinación del valor del estado en el instante t = t 0+ ∆t para permitir el cálculo de la derivada. La aplicación de este método para el sistema ejemplo de la ecuación (7) resulta en
{
( |) ( |) d x1
( x1 )n +1=( x 1 )n + dt ( x2 )n +1=( x 2 )n +
d x2 dt
∆t
n+1
… … … … … … … … … … … … … … … … …(16) ∆t
n +1
Por otro lado, la implementación implícita continua siendo muy simples, debido a la sustitución de las derivadas por sus expresiones analíticas:
[
{
][ ] [ ][ ]
( x 1 )n+1=( x1 )n + ( x 2 )n +1 ∆t 1 −∆ t ( x 1)n+1 1 0 ( x 1)n ⇒ = ………………………… ( x2 )n +1=( x 2 )n ±a 1 ( x 1 )n+1 −a2 ( x 2) n+1 ∆ t a1 ∆t 1+ a2 ∆ t ( x 2)n+1 0 1 ( x 2)n
MÉTODO DE EULER MODIFICADO El método de Euler modificado intenta mejorar el desempeño del método convencional a través de la mejoría de la estimativa de la derivada. Para eso, la ecuación (10) será modificada haciendo la derivada igual a la media entre la derivada en el instante t = t0 y el valor de la derivada calculada para una estimativa del estado en el instante t = t0+ ∆t. Este es los más simples de los métodos de integración del tipo preditor-corrector, en que surge un paso de ESTABILIDAD TRANSITORIA
Página 11
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
predicción del valor del estado en el instante t = t 0+ ∆t y este valor es utilizado para la corrección y obtención de un nuevo valor del estado en ese instante. El método, por lo tanto, puede ser descrito por los siguientes pasos: Paso 1: Preditor
…………………………………(18) Paso 2: Corrector
……………………(19) Se debe notar que la ecuación (18) corresponde a un paso del método de Euler convencional, mostrado en la ecuación (10). Ya la ecuación (19) es muy parecida con el método trapezoidal, siendo la única diferencia el cálculo aproximado (a partir de la estimativa x p) de la derivada en el instante t = t0+ ∆t. La implementación de este método para la solución del sistema ejemplo resulta en
……………………………………(20)
Y, para el paso corrector,
ESTABILIDAD TRANSITORIA
Página 12
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
……………………… (21) Este
método
también
puede
ser
implementado
de
forma
implícita,
sustituyendo las expresiones de las derivadas en las ecuaciones (20) y (21):
……………………(22)
…(23) Se puede, aún, sustituir los valores de x 1p y x2p, obtenidos en las ecuaciones (22), en la ecuación (23), resultando en
……(24) Aplicando este método a la solución a la ecuación de oscilación, se obtiene el siguiente procedimiento de solución: a) Calculo de la potencia acelerante y las derivadas del ángulo y de la velocidad al principio del intervalo, es decir, en el instante t:
ESTABILIDAD TRANSITORIA
Página 13
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
………..(25) b) Cálculo del primer valor del ángulo y la velocidad al final del intervalo, es decir, en el instante
t+ ∆ t :
……….(26) c) Cáculo de la potencia acelerante y de la derivada de la velocidad al final del intervalo, es decir, en el instante
t+ ∆ t :
…(27) d) Cáculo del valor corregido de la velocidad al final del intervalo, es decir, en el instante
t+ ∆ t :
…….(28) e) Cálculo de la derivada del ángulo al final del intervalo, es decir, en el instante t+ ∆ t : ESTABILIDAD TRANSITORIA
Página 14
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
…..(29) f) Cálculo del valor corregido del ángulo al final del intervalo, es decir, en el instante t+ ∆ t :
……….(30) El proceso iterativo se continúa hasta que se presenten nuevos cambios en el sistema o hasta que se concluya el tiempo de estudio.
EJEMPLO: 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 como se muestra en la siguiente figura.
ESTABILIDAD TRANSITORIA
Página 15
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
El generador entrega una potencia real de Pelec = 0.8 p.u, y Qelec = 0.074 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.
RESOLUCIÓN. La impedancia de transferencia entre el voltaje generador y la barra de potencia infinita Antes de que ocurra la falla, la reactancia entre los puntos A y B puede ser encontrada por la combinación serie paralelo:
ESTABILIDAD TRANSITORIA
Página 16
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
La potencia aparente es dada por:
El voltaje interno de la maquina E’ es dado por:
La ecuación de potencia eléctrica antes de la perturbación queda dada por:
Siendo el ángulo de operación estable inicial:
Durante la falla se tiene, que la línea fallada se divide en dos, resultando el diagrama de impedancias como en la Figura siguiente.
Se efectúa una transformación de estrella a delta y se logra:
ESTABILIDAD TRANSITORIA
Página 17
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
La reactancia equivalente entre el generador y la barra de potencia infinita resulta ser:
La ecuación de potencia eléctrica antes de la perturbación queda dada por:
Cuando la falla es despejada, la línea fallada es puesta fuera de servicio, resultando el diagrama de reactancias como sigue.
Donde la reactancia equivalente entre el generador y la barra de potencia infinita resulta ser:
De tal modo que la potencia eléctrica después de la falla queda dada por:
Se procede a plantear las ecuaciones de oscilación, de la forma:
ESTABILIDAD TRANSITORIA
Página 18
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
ANTES de la perturbación resulta:
DURANTE la perturbación resulta:
DESPUÉS la perturbación resulta:
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:
ESTABILIDAD TRANSITORIA
Página 19
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
Solución analítica método de euler:
[
ω (t k ) δ ( t k +1) δ (t k ) = + ∆ t πf ( Pmec −Pmax sin δ ( t k ) ) ω ( t k+1 ) ω (t k ) H
[ ][ ] •
Para un tiempo de t=0
δ ( 0 ) =26.38 •
]
; ω ( 0 ) =0
Para un tiempo t=0.05
[
0 δ ( 1 ) = δ ( 0 ) +0.05 π∗60∗( 0.8−PmaxII sin δ ( 26.38 )) ω (1) ω (0) 5
[ ][ ]
[
0 δ ( 1 ) = 26.38 +0.05 π∗60∗(0.8−0.65 sin δ ( 26.38 ) ) 0 ω (1) 5
[ ][ ] •
;
[ ][
]
;
[ ][
]
δ ( 1 ) = 26.38 0.9638 ω (1)
Para un tiempo t=0.1
[
ω (1) δ ( 2 ) = δ ( 1 ) +0.05 π∗60∗(0.8−PmaxII sin δ ( 26.38 ) ) ω (2) ω ( 1) 5
[ ][ ]
[
0.963∗180/ pi δ ( 2 ) = 26.38 +0.05 π∗60∗(0.8−0.65 sin δ ( 26.38 )) 0.963 ω (2) 5
[ ][ ] •
]
]
Para un tiempo t=0.15
ESTABILIDAD TRANSITORIA
Página 20
] ]
δ ( 2 ) = 29.147 1.926 ω (2)
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
[
ω (2) δ ( 3 ) = δ ( 2 ) +0.05 π∗60∗(0.8−PmaxII sin δ (26.38 )) ω (3) ω ( 2) 5
[ ][ ]
]
[
1.926∗180 / pi δ ( 3 ) = 29.147 +0.05 π∗60∗( 0.8−0.65sin δ ( 29.147 )) 1.926 ω (3) 5
[ ][
]
[ ][
]
δ ( 3 ) = 34.665 2.837 ω (3)
•
[
[ ][ ]
]
;
]
;
]
[
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
[
ω(4 ) δ ( 5 ) = δ ( 4 ) +0.05 π∗60∗(0.8−PmaxII sin δ ( 42.792 )) ω (5) ω (4 ) 5
[ ][ ]
[
]
1.926∗180 / pi δ ( 5 ) = 42.792 +0.05 π∗60∗( 0.8−0.65sin δ ( 42.792 )) 3.648 ω (5) 5
[ ][
]
[ ][
]
δ ( 5 ) = 53.243 4.324 ω (5)
•
;
Para un tiempo t=0.2
ω (3) δ ( 4 ) = δ ( 3 ) +0.05 π∗60∗(0.8−PmaxII sin δ (34.665 )) ω (4 ) ω (3) 5
•
]
Para un tiempo t=0.3 despeje de la falla
ESTABILIDAD TRANSITORIA
Página 21
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
[
ω ( 5) δ ( 6 ) = δ ( 5 ) + 0.05 π∗60∗(0.8−P maxIII sin δ ( 5 ) ) ω (6) ω (5) 5
[ ][ ]
]
[
4.325∗180/ pi δ ( 6 ) = 53.252 +0.05 π∗60∗(0.8−1.4625 sin δ ( 53.252 )) 4.325 ω (6) 5
[ ][
]
[ ][
]
δ ( 6 ) = 65.641 3.624 ω ( 6)
]
;
Este método se repite hasta poder observar la estabilidad o inestabilidad del sistema mediante la curva angula versus tiempo.
Cuadro de valores: t
Pm
Pe
δ
ω
0.00
0.8
1.8000
26.38
0
0.05
0.8
0.6500
26.38
0.9638
0.10
0.8
0.6500
29.147
1.926
0.15
0.8
0.6500
34.665
2.837
0.20
0.8
0.6500
42.792
3.648
0.25
0.8
0.6500
53.243
4.324
0.30
0.8
1.4625
65.641
3.624
0.35
0.8
1.4625
76.022
2.62
0.40
0.8
1.4625
83.529
1.453
0.45
0.8
1.4625
87.692
0.222
0.50
0.8
1.4625
88.327
-1.024
……
……
…….
…….
………
ESTABILIDAD TRANSITORIA
Página 22
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
EULER MODIFICADO Considerando de los cálculos previos: Falla: PeII= 0.65 sen δ Post-Falla: PeIII= Condiciones iníciales en estabilidad: δ ω
= 0.46 rad 0=
0 rad/s
∆ t = 0.15; t despeje= 0.3s
Primera Iteración: t1=t0+ ∆ t = 0+0.1=0.1s Hallando parámetros previos: ∂(δ ( t 0 ) ) ∂t
=0-0.377= -0.377
∂(ω ( t 0 ) ) л x 60 = (0.8-0.65 sen 0.46) = 19.27 rad/s² ∂t 5 δ 1(0+0.1) = 0.46 + (-377) (0.1) = -37.24 rad ESTABILIDAD TRANSITORIA
Página 23
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
ω 1(0+0.1) = 0 + 19.27)(0.1) = -1.927 rad/s
∂(δ ( 0+ 0.1 )) ∂t
= 1.927 – 377 = -377.093
∂(ω ( t 0 ) ) ∂t
л x 60 (0.8-0.65 sen -37.24) = 19.3 rad/s² 5
=
δ 1(0+0.1) = 0.46 +
(−377 ) +(−377.073) 2
δ 1(0.1) = -37.14rad
ω 1(0+0.1) = 0 +
19.27 19.3 2
x0.1
δ 1(0.1) = 31.83º
x 0.1
ω 1(0.1) = 1.9285 rad/s
Segunda Iteración: t2=t1+ ∆ t = 0.1+0.1=0.2s Hallando parámetros previos: ∂(δ 2 ( 0.1 ) ) ∂t
=1.9285 - 0.377 = -375.07
∂(ω 2 ( 0.1 )) л x 60 = (0.8-0.65 sen -37.14) = 17.16 rad/s² ∂t 5
δ 2 (0+0.1) = -37.14 + (-375.07) (0.1) = -74.65 rad ESTABILIDAD TRANSITORIA
Página 24
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
ω 2(0.1+0.1) = 1.9285 + (17.16) (0.1) = 3.64 rad/s
∂(δ 2 ( 0.1+0.1 ) ) ∂t
= 3.64 – 377 = -373.36
∂(ω 2 ( 0.1+0.1 )) ∂t
=
л x 60 (0.8-0.65 sen -74.65) = 13.4 rad/s² 5
δ 2 (0.1+0.1) = -37.14 +
(−375.07 ) +(−373.36) 2
δ 2(0.2) = -74.5615rad
ω 2(0.1+0.1) = 1.9285 +
x0.1
δ 2 (0.2) = 47.94º
17.16+ 13.488 2
x 0.1
ω 2(0.2) = 3.46 rad/s
Tercera Iteración: t3=t2+ ∆ t = 0.2+0.1=0.3s Hallando parámetros previos: ∂(δ 3 ( 0.2 ) ) ∂t
=3.46 - 377 = -373.54
∂(ω 3 ( 0.2 ) ) л x 60 = (0.8-0.65 sen -74.5615) = 11.97 rad/s² ∂t 5 δ 3 (0.2+0.1) = -74.5615 + (-373.54.07) (0.1) = 1.18rad
ω 3(0.2+0.1) = 3.46 + (11.97)(0.1) = 4.657 rad/s
ESTABILIDAD TRANSITORIA
Página 25
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
∂(δ 2 ( 0.2+0.1 ) ) ∂t
= 4.657 – 377 = -372.343
∂(ω 3 ( 0.2+0.1 ) ) ∂t
=
л x 60 (0.8-1.4625 sen 1.18) = -20.8 rad/s² 5
(−373.54 )+(−372.343) 2
δ 3 (0.2+0.1) = 0.837 + δ 3(0.3) = -36.46rad
x0.1
δ 2 (0.3) = 71º 11.97 +(−20.82) 2
ω 3(0.2+0.1) = 3.46 +
δ 2 (0.3) = 1.24rad
x 0.1
ω 2(0.3) = 3.0175 rad/s
Tabla de Resultados
Iteraciones Parámetros
Grados sexagesimales
δ
Rad/s ¿ ω ¿ velocidad)
(Ángulo)
1
31.83 º
1.9285
2
47.94 º
3.46
3
71 º
3.0175
Se procedió a implementar un programa en Matlab. %FACULTAD DE INGENIERIA ELECTRICA Y ELECTRONICA-UNCP %CURSO DE ESTABILIDAD DE SISTEMAS DE POTENCIA %programa con el metodo de euler clc Pmax1=input('Ingrese la potencia antes de la falla: ');%antes de la falla ESTABILIDAD TRANSITORIA
Página 26
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
Pmax2=input('Ingrese la potencia durante la falla: '); %durante la falla Pmax3=input('Ingrese la potencia despues de la falla: ');%despues de la falla Wr0 = input('Ingrese la velocidad incial: '); %velocidad inicial Dt=input('Ingrese incremento de tiempo: '); t0=input('Ingrese tiempo inicial de simulación: '); tf=input('Ingrese tiempo final de simulación: '); H=input('Ingrese la constante de inercia: '); tdespeje=input('Ingrese tiempo de despeje de la falla: '); Pm=input('Ingrese la potencia mecanica: '); f=input('Ingrese la frecuencia de la red en Hz: '); Delta0=asin(Pm/Pmax1); %angulo incial en radianes N=(tf-t0)/Dt;% calculo de numero de puntos a simular v=[]; for I=1:N ti(I)=t0+I*Dt; % calcular tiempo Delta1(I)=Delta0+Dt*( Wr0); %angulo i+1 Wr1(I)= Wr0+Dt*(pi*f/H*(Pm-Pmax2*sin(Delta0)));%velocidad i+1 if ti(I)>=tdespeje Wr1(I)= Wr0+Dt*(pi*f/H*(Pm-Pmax3*sin(Delta0))); end Delta0= Delta1(I); Wr0= Wr1(I); v=[v,ti(I), Wr1(I), Delta1(I)*180/pi]; end v plot(ti, Delta1*180/pi); xlabel('tiempo'); ylabel('angulo'); title('angulo VS tiempo'); ESTABILIDAD TRANSITORIA
Página 27
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
grid; Resultados de la simulación:
ESTABILIDAD TRANSITORIA
Página 28
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
El sistema para un tiempo de despeje de 0.3 segundos el sistema es inestable
ESTABILIDAD TRANSITORIA
Página 29
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
Para un tiempo de despeje de la falla a 0.1 segundos sistema es inestable.
CONCLUSIÓN: Podemos afirmar, que los programas aquí expuestos (Matlab) resuelven, de primer orden; y probablemente podemos destacar los errores que existen por ESTABILIDAD TRANSITORIA
Página 30
UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ
INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
cada método. De los métodos. Afirmamos que los errores del método de Euler, radica en un intervalo proporcional a ∆t=2, mientras que su error global es proporcional a ∆t; este método podría ser inestable si la Ecuación de oscilación, tiene una constante de tiempo con signo negativo, a menos que se utilice una ∆t pequeña, en cambio en el método modificado si la Ecuación de oscilación, no es lineal, se requiere de un método iterativo para cada intervalo. Su error en un intervalo es proporcional a ∆t=3, mientras que su error global lo es a ∆t=2. En fin podemos afirmar que ambos métodos poseen una desventaja, que consiste en que los órdenes de precisión son bajos. Esta desventaja tiene dos facetas, para mantener una alta precisión se necesita una ∆t pequeña, lo que aumenta el tiempo de cálculo y provoca errores de redondeo.
ESTABILIDAD TRANSITORIA
Página 31