UNIVERSIDAD NACIONAL DEL CENTRO DEL PERÚ FACULTAD DE INGENIERÍA ELÉCTRICA Y ELECTRÓNICA
SOLUCIÓN NUMÉRICA DE LA ECUACIÓN DE OSCILACIÓN MEDIANTE EL MÉTODO DE EULER Y EULER MODIFICADO CURSO:
ESTABILIDAD DE SISTEMAS DE POTENCIA
CATEDRÁTICO:
DR. ING. PEDRO NICOLÁS TORRES MAYTA
INTEGRANTES:
ARAOZ PALOMINO KLINSMAN HIDALGO HILARIO NILO KELVIN HINOSTROZA MEZA ÁNGEL CLIF MALLQUI APOLINARIO LENNIN SANDOVAL OBISPO ERIK OSORIO CHILQUILLO ELIZABETH POMA SUELDO DONATO RAMOS TORRES RYDER RAMOS YANGALI JHON SINCHE BARRA JOSELYN TORRES CLEMENTE POOL VELAZCO FLORES ANDRÉ
SEMESTRE:
X Huancayo – Huancayo – 2018 2018
INTRODUCCIÓN ............................................................................................................................. ............................................................................................................................. 2 I.
INTEGRACIÓN NUMÉRICA ............................................................... ................................................................................................. .................................. 3
II.
MÉTODO DE EULER Y EULER MODIFICADO ........................................................... ...................................................................... ........... 7
III.
EJEMPLO DE APLICACIÓN ............................................................. ............................................................................................ ............................... 13
IV.
SOLUCIÓN MEDIANTE EL USO DE HERRAMIENTAS COMPUTACIONALES .................. 21
CONCLUSIONES .................................................................................................................. ........................................................................................................................... ......... 26
1
En el presente trabajo se muestra 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 poseen 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ó una hoja de cálculo en MS Office Excel y un programa en Matlab para analizar la solución en diferentes puntos mostrando su aproximación al valor verdadero.
2
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:
Donde:
̇ ,, …1
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,, …2
Que deben ser satisfechas a cada instante de tiempo. En general, ecuaciones diferenciales no lineales sólo admiten solución numérica, calculada paso a paso. Diversos algoritmos de integración numérica son conocidos y la opción del “mejor” depende grandemente del problema.
Siendo así, 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:
∆ ∆ ∆ = ∆ = 2! = 3! ⋯ 3 a) PRESICIÓN : Esta característica dada, principalmente, por dos factores:
Error de redondeo
Error de truncamiento 3
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:
∆+ …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. b) 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. 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
∝ + ∝ ∝ − ⋯∝ −+ + ⋯ −+ 0…6
Donde αi y βi son constantes y se desea determinar y n+1, conociendo los
valores anteriores de yi (i < n+1).
4
Esta es una ecuación discreta, lineal, y será estable si todas las raíces de su polinomio característico tuvieran 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 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 estables y los métodos de Euler Reverso y Trapezoidal están en esta categoría. c) 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 mucho 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 estable, puede haber problemas de inestabilidad numérica. 1. 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.
5
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:
0 ̇̇ 0 …7 ℷℷℷ ℷ ℷ ℷℷℷℷℷ …8 { ℷ ℷ ℷ ℷ ℷ ℷ ℷ ℷ ±√ ℷ, 2 4 …9
cuya solución analítica es de la forma
Donde
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.
6
1. 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:
∆ = …10
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:
∝1 ⋯∝1 0
0 ⋯ ∆ 0 ………11
La aplicación del método de Euler para la solución del sistema dado por la
ecuación (7) resulta en:
7
+ ()∆t …12 {+ ( )∆t
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:
0 + ̇̇ ⇒[ ] 0 + [1∆ 1∆∆] []…13
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:
[][] ′ ′ ′ , sin
Sustituyendo los respectivos términos:
8
[++] [] ∆ sin + ∆∗ + ∆∗ sin + ∆∗ + ∆∗ sin
En forma de ecuaciones se tiene:
Se puede suponer:
=
De tal forma que se puede escribir en forma más compacta las ecuaciones:
2. 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 = t 0 y el valor de la derivada calculada para una estimativa del estado en el instante t = t 0+ ∆t. Este es los más simples de los 9
métodos de integración del tipo predictor-corrector, en que surge un paso de 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: Predictor …(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:
…(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):
10
…(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) Cálculo de la potencia acelerante y las derivadas del ángulo y de la velocidad al principio del intervalo, es decir, en el instante t:
…(25)
b) Cálculo del primer valor del ángulo y la velocidad al final del intervalo, es decir, en el instante
∆
:
…(26)
c) Cálculo de la potencia acelerante y de la derivada de la velocidad al final del intervalo, es decir, en el instante
∆
: 11
…(27)
d) Cálculo del valor corregido de la velocidad al final del intervalo, es decir, en el instante
∆
:
…(28)
e) Cálculo de la derivada del ángulo al final del intervalo, es decir, en el instante
∆
:
…(29)
f) Cálculo del valor corregido del ángulo al final del intervalo, es decir, en el instante
∆
:
…(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.
12
Un generador sincrónico a 60 Hz, posee una constante de inercia H igual a 5MJ/MVA y una reactancia transitoria de eje directo X’d = 0.3 p.u, está c onectado
a una barra de potencia infinita a través de un circuito puramente reactivo como se muestra en la siguiente figura:
El generador entrega una potencia real de P elec = 0.8 p.u, y Q elec = 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. Efect úe la simulación numérica hasta t = 0.5.
a) En Pre falla: 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:
0.300.2 0.230 0.65 . .
13
La corriente se calcula:
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:
´ || 1∗1.0.6517 1.80∗ 0.81.80∗ 26.387800.46055 Siendo el ángulo de operación estable inicial:
b) En Falla: 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:
14
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:
c) En Post falla: 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:
ANTES de la perturbación resulta:
15
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. Efectúe la simulación numérica hasta t = 0.5. Las ecuaciones de oscilaciones a considerar resultan ser:
1. MÉTODO DE EULER:
•
[++] [] ∆ sin 0 ;0 [11] [00] 0.05∗60∗0. 8 05 sin26.38 [11] 26.0380.05∗∗.−.0. [11] 0.26.963838 [22] [11] 0.05∗60∗0. 8 51 sin26.38 [22] 20.6.963380.05∗∗.0.963∗180/ −.. [22] 21.9.912647
Para un tiempo de t=0:
=26.38
•
=0
Para un tiempo t=0.05:
;
•
Para un tiempo t=0.1:
;
16
•
Para un tiempo t=0.15:
[33] [22] 0.05∗60∗0. 8 52 sin26.38 [33] 21.9.9126470.05∗∗.1.926∗180/ −.. [33] 32.4.863765 [44] [33] 0.05∗60∗0. 8 53 sin34.665 [44] 32.4.8637650.05∗∗.2.837∗180/ −.. [44] 43.2.674892 [55] [44] 0.05∗60∗0. 8 54 sin42.792 [55] 43.2.6748920.05∗∗.3.648∗180/ −.. [55] 54.3.322443 [66] [55] 0.05∗60∗0. 8 55 sin5 [66] 54.3.3225520.05∗∗.4.3−.25∗180/ 6 6 5. 6 41 [ ] . 6 3. 6 24 ;
•
Para un tiempo t=0.2:
;
•
Para un tiempo t=0.25:
;
•
Para un tiempo t=0.3 (despeje de la falla):
;
Este método se repite hasta poder observar la estabilidad o inestabilidad del sistema mediante la curva angula versus tiempo. Cuadro de resultados:
t (s)
Pm
Pe
0.00
0.8
1.8000
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
(sexag.) 26.38
(rad/s) 0
17
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
……
……
…….
…….
………
2. MÉTODO DE EULER MODIFICADO: Considerando los cálculos previos y las condiciones iniciales de estabilidad:
= 0.46 rad
Primera iteración:
∆
0=
0 rad/s
= 0.1; tdespeje= 0.3s
∆
t1=t0+ = 0+0.1=0.1s Hallando parámetros previos:
л x +. л x −+−. .+. =0- 377= -377
(0.8-0.65 sen 0.46) = 19.27 rad/s²
(0+0.1) = 0.46 + (-377) (0.1) = -37.24 rad (0+0.1) = 0 + (19.27) (0.1) = 1.927 rad/s = 1.927 – 377 = -377.093
=
(0.8-0.65 sen -37.24) = 19.3 rad/s²
(0+0.1) = 0.46 +
(0.1) = -37.14rad (0+0.1) = 0 +
x0.1
(0.1) = 31.83º
x 0.1
(0.1) = 1.9285 rad/s
18
Segunda iteración:
∆
t2=t1+ = 0.1+0.1=0.2s Hallando parámetros previos:
. . л x .+. .+. л x −.+−. .+. =1.9285 - 377 = -375.07
(0.8-0.65 sen -37.14) = 17.16 rad/s²
(0.1+0.1) = -37.14 + (-375.07) (0.1) = -74.65 rad (0.1+0.1) = 1.9285 + (17.16) (0.1) = 3.64 rad/s = 3.64 – 377 = -373.36 =
(0.8-0.65 sen -74.65) = 13.4 rad/s²
(0.1+0.1) = -37.14 +
x0.1
(0.2) = -74.5615rad
(0.2) = 47.94º
(0.1+0.1) = 1.9285 +
x 0.1
(0.2) = 3.46 rad/s
Tercera iteración:
∆
t3=t2+ = 0.2+0.1=0.3s Hallando parámetros previos:
. . л x .+. .+. л x
=3.46 - 377 = -373.54 (0.8-0.65 sen -74.5615) = 11.97 rad/s²
(0.2+0.1) = -74.5615 + (-373.54) (0.1) = 1.18rad (0.2+0.1) = 3.46 + (11.97) (0.1) = 4.657 rad/s = 4.657 – 377 = -372.343 =
(0.8-1.4625 sen 1.18) = -20.8 rad/s²
19
−.+−. .+−.
(0.2+0.1) = 0.837 + (0.3) = -36.46rad (0.2+0.1) = 3.46 +
x0.1
(0.3) = 71º x 0.1
(0.3) = 1.24rad
(0.3) = 3.0175 rad/s
Finalmente se obtiene el cuadro de resultados: Iteración 1
(sexag.) ꞷ (rad/s) 31.83
1.9285
2
47.94
3.46
3
71
3.0175
20
1. CON EL MS OFFICE EXCEL: En la siguiente tabla se muestra los resultados transpuestos, cuando la falla fuera sostenida (sin aclaramiento), y cuando fuera aclarada en los 0.3 segundos como se suponía al principio. Iteración (i) ti (s)
ᵟi(rad)
ᵟi(sexagesimal)
wi (rad/s)
0
0
0.46055
26.3875713
0
1
0.05
0.46055
26.3875713
0.96342613
2
0.1 0.50872131
29.1475838
1.92685226
3
0.15 0.60506392
34.6676089
2.83805945
4
0.2 0.74696689
42.7980503
3.6491002
5
0.25 0.9294219
53.2519523
4.32462941
6
0.3
1.1456534
65.641103
4.8508558
7
0.35 1.3881962
79.537781
5.2426686
8
0.4
1.6503296
94.55692
5.54578141
9
0.45 1.9276187
110.44441
5.8323978
10
0.5
127.153
6.19231597
2.2192386
Ángulo vs. Tiempo 4000 3500 3000 2500 2000 1500 1000 500 0 -500 0
0,5
1
1,5 Falla despejada
2
2,5
3
3,5
Fala sostenida
21
En la figura se visualiza que el ángulo crece indefinidamente cuando la falla persiste en el sistema. La siguiente tabla muestra la evolución del ángulo a través del tiempo hasta los 1.75 segundos una vez se ha despejado la falla a los 0.3 segundos. Paso tiempo (s) Iteración (i) ti (s)
ᵟi(rad)
ᵟi(sexages.)
wi (rad/s) 0
0.05
0
0
0.46055
26.3875713
0.05
1
0.05
0.46055
26.3875713 0.96342613
0.05
2
0.1 0.50872131 29.1475838 1.92685226
0.05
3
0.15 0.60506392 34.6676089 2.83805945
0.05
4
0.2 0.74696689 42.7980503 3.6491002
0.05
5
0.25 0.9294219 53.2519523 4.32462941
0.05
6
0.3 1.14565337 65.641103
4.8508558
0.05
7
0.35 1.38819616 79.537781
3.847479
0.05
8
0.4 1.58057011 90.559997 2.64452723
0.05
9
0.45 1.71279647 98.136009 1.39587582
0.05
10
0.5 1.78259026
0.05
11
0.55 1.79133225 102.63578 -1.01234494
0.05
12
0.6
0.05
13
0.65 1.63099697 93.449243 -3.40344242
0.05
14
0.7 1.46082485 83.699099 -4.64723162
0.05
15
0.75 1.22846327 70.385761 -5.87936182
0.05
16
0.8 0.93449518
53.54263
-6.96818202
0.05
17
0.85 0.58608608 33.580259
-7.6774674
0.05
18
0.9 0.20221271 11.585935
-7.6942725
0.05
19
0.95 -0.18250092 -10.45653 -6.73996616
0.05
20
0.05
21
1.05 -0.75608327 -43.32038 -1.85514151
0.05
22
1.1 -0.84884034 -48.63497
1.54416455
0.05
23
1.15 -0.77163212 -44.21126
5.12110817
1
1.740715
102.1349
0.17483965
99.735623 -2.19436058
-0.51949923 -29.76511 -4.73168088
22
0.05
24
1.2 -0.51557671 -29.54037
8.55136932
0.05
25
1.25 -0.08800824 -5.042501
11.4185114
0.05
26
1.3 0.48291733 27.669125 13.1687793
0.05
27
1.35 1.1413563
0.05
28
1.4 1.81118665 103.77335 12.3981394
0.05
29
1.45 2.43109362
0.05
30
1.5 2.99252493 171.45905 10.9386065
0.05
31
1.55 3.53945526 202.79585 12.0371491
0.05
32
1.6 4.14131271 237.27974 14.6132121
0.05
33
1.65 4.87197331 279.14351 18.4404826
0.05
34
1.7 5.79399744
0.05
35
1.75 6.92750573 396.91684 25.4735512
65.394899 13.3966071 139.2914
331.9716
11.2286262
22.6701658
Ángulo vs. Tiempo 4000 3500 3000 2500 2000 1500 1000 500 0 -500 0
0,5
1
1,5
2
2,5
3
3,5
Falla despejada
En la figura se alcanza notar que el ángulo oscila entre valores positivos y negativos. Pero según transcurre el tiempo, el ángulo crecerá indefinidamente. Lo cual muestra un régimen de inestabilidad, a pesar del despeje de falla a los 0.3 segundos.
23
2. CON MATLAB:
El ingreso de datos arroja los siguientes resultados:
24
El sistema para un tiempo de despeje de 0.3 segundos es inestable:
25
1. Los métodos de Euler y Euler modificado son herramientas eficaces para la solución numérica de la ecuación de oscilación de un sistema, pudiéndose calcular el ángulo y la velocidad en función del tiempo; y de esta manera predecir la estabilidad del sistema en análisis. Para obtener resultados con mayor precisión se debe elegir un Δt
pequeño, aunque incrementa el número y el tiempo de cálculo; éste inconveniente se solucionó con la implementación de una hoja de cálculo en MS Office Excel y un programa en Matlab, ambos permiten graficar la variación de la velocidad y el ángulo en función del tiempo. 2. En el ejercicio desarrollado se hace clara la inestabilidad del sistema a pesar de ser liberada la falla y también como era de esperarse la inestabilidad se hace más grave en el caso de falla sostenida.
26
