Métodos de predicción-corrección de Adams-Bashforth-Moulton •
-
pasos y p
y
h 55 f ( y (t ), u (t )) 59 f ( y (t h), u (t h))
24 37 f ( y (t 2h), u (t 2h)) 9 f ( y (t 3h), u (t 3h))
• Utiliza como corrector el método de Adams-Moulton de 3 pasos y ( t h) y ( t )
h 9 f ( y p (t h), u(t h )) 19 f ( y (t ), u( t ))
,
,
• En ambos métodos el error es proporcional a h 5 ,
Simulación. Master en IPS 11/12
59
Errores en la solución numérica de un modelo dinámico •
Ex sten 3 t pos e errores a cons erar en a ntegrac n e as ecuac ones, a suma e os tres constituye el error global:
•
Errores locales
–
rror e runcam en o epen e e or en e a gor mo e n egrac n y e ama o e paso
de integración (h)
d y (t ) n
n
• Para uno de orden n el error es proporcional a: h dt n • Así un método de orden mayor permite trabajar con un h mayor manteniendo la precisión. – Error de redondeo, está ligado a la precisión de la
representación de los números en el ordenador y depende de: •
e on eo por pa a ras e ong tu
n ta
• Acumulación de operaciones • Tipo de máquina • Tipo de compilador • Tipos de datos: precisión sencilla, doble, etc. – Error de acumulación, es el error que se comete en un
los anteriores pasos de integración. • Un número de pasos de integración grande contribuye a aumentar el efecto de este error
•
Error global es del orden de h n-1 y no de hn Simulación. Master en IPS 11/12
60
Mejoras en la estimación: extrapolación • Existen procedimientos para mejorar la estimación de y(t+h) alternativos a usar un algoritmo de mayor orden. • Estos métodos están basados en el conocimiento del orden del error proporcionado por el método que utilicemos. – Sabemos que si tomamos la serie de Taylor hasta el término de orden n el error es proporcional a hn. – , e h e y2 la solución estimada con un paso h/2. Entonces:
chn h y2 ye c( )n 2
– Eliminando c· hn entre ambas ecuaciones tendremos:
y 2 2 n y1 ye 2 n 1
Simulación. Master en IPS 11/12
61
Métodos de paso variable • En los métodos de paso fijo el error cometido es proporcional a hn, siendo n el orden del algoritmo de integración. –
.
• El número de veces que resuelve el problema de integración de las ecuaciones es tmax /h. – Menor h más número de operaciones, más tiempo de cálculo.
• Solución de compromiso: –
.
• Problema: – Escoger el valor de h. – No se obtienen y soluciones en instantes de tiempo regularmente espaciados, si se necesitan debemos interpolar.
¿Concepto de intervalo
t Simulación. Master en IPS 11/12
La distancia entre dos untos consecutivos en los que se almacena la salida del modelo. 62
• En eneral como arámetros de control de la simulación tenemos dos parámetros: – Paso de integración ( h) – Intervalo de comunicación ( cint).
• Se suele definir un tercer parámetro que es el número de pasos e ntegrac n en un nterva o e comun cac n nstp . • Estos tres parámetros no son independientes, se solicitan sólo os e e os. • En general se solicita cint y nstp, de modo que – e m o o es e paso o: =c n ns p – Si el método es de paso variable se calcula el tamaño del paso inicial como hi=cint/nstp
Simulación. Master en IPS 11/12
63
• Existen diversos métodos de paso variable: – Ex lícitos: Run e-Kutta-Merson Run e-Kutta-Fehlber ... – Implícitos: Gear …
• ¿Cómo estimar el valor del paso de integración? – El error cometido en un paso es: e=c·hn. – Especifico una cota de error máximo en cada paso de integración E. – Selecciono el valor del paso de integración ( h1), de modo que e=E. • Sabiendo que E=c· h1n, entonces: 1
n
E e
• Problema es como estimo el error cometido e. • 2 d y (t ) 2 1 2 1 n h f h y t f h y t f (h y (t )) ( ( )) ( ( )) 2 dt n 2 y t
y t y t 2 f ( y (t h)) ( f ( y (t h))) f ( y (t h)) f ( y (t )) ( f ( y (t )) f ( y (t h))) t h 2· t t h ... Simulación. Master en IPS 11/12 64
• Otro método para estimar el error: – Supongamos, por simplicidad: – La estimación del valor en t+h:
y (t ) f ( y (t )) y (t h) y (t ) hf ( y (t ))
h f ( y (t ))
2!
y
f ( y (t ))
– Si aplicamos la misma fórmula con un paso de integración negativo (-h) podemos re-estimar el valor en t usando el valor calculado en t+h: e
2
h f ( y (t h))
2!
y
– Comparando el valor real en t con el valor estimado en la segunda ecuación - e . • Si e< seguimos la integración con el paso h. • Si no, tomamos el nuevo paso de integración h=h/2 y repetimos la estimación de y(t+h). • Si la condición se cumple durante un determinado número de pasos de integración (por ejemplo 5) el paso de integración pasa a valer el doble h=h*2.
– Excesivos cálculos, tengo que repetir el algoritmo pero hacia atrás.
Simulación. Master en IPS 11/12
65
– Método de Runge-Kutta-Fehlberg: • Estima el valor de t+h usando un método de Run e-Kutta de orden 4 • Estima el error cometido en un paso a partir del valor de y(t+h) calculado mediante un método de Runge-Kutta de orden 5 (y RK5(t+h)) – e= y – yRK5 – La ventaja es que el RK5 utiliza los mismos términos que los calculados en el RK4 más un término adicional
• Se define una tolerancia para el error global: E – E=MAX(Tolabs,Tolrel·|y(t+h)|)
• Se calcula q=[E·h/e]1/4 – Si q≤0.1, se repiten los cálculos con un paso q·h – Si q≥4., el siguiente paso de integración será MAX(4·h,hmax) – Si 0.1 < <4. el si uiente aso de inte ración será MAX ·h h
• Una estrategia más conservadora es considerar q=0.84·[E·h/e] 1/4 • Existen estrategias adicionales para manejar el tamaño del paso de integración pero se escapan del alcance de este curso, en particular para los métodos Simulación. Master en IPS 11/12 66 multipaso.
Simulación. Master en IPS 11/12
67
Estabilidad • El error entre la solución numérica y la analítica debe tender a cero a medida que la integración progresa. paso de integración h • Sin estabilidad las soluciones divergen de la verdadera o son oscilantes
Solución x
oscilante Solución exacta t Solución inestable
Simulación. Master en IPS 11/12
68
d y
·y
y(0) 1; 0.
solución analítica y(t) e - t
• Método de Euler explícito: y (t h) y (t ) h· y (t ) (1 ·h)· y (t )
y( t h ) 1 h 1 y( t )
• Para que y(t): – Sea decreciente: 0< ·h<2
– Sea siempre positiva (negativa, si c.i es menor que cero): 0< ·h<1 Simulación. Master en IPS 11/12
69
• Método de Euler implícito: y( t h ) y( t ) h μ y( t h ) y( t h ) • Para que y(t):
y( t h )
–
1
y( t )
1
– Sea siempre positiva: cualquier h (que siempre es mayor que cero)
• Método de Crank-Nicholson:
h y( t h ) y( t ) (μ y( t h ) μ y( t )) 2
• Para que y(t): – Sea decreciente: cualquier valor de h
y( t h )
1 μ
1 μ
h 2 1 2
– Sea siempre positiva (negativa, si c.i es menor que cero): 0< ·h/2<1
Simulación. Master en IPS 11/12
70
Estabilidad: caso general • •
Se ha estudiado la estabilidad de un modelo matemático formado por una sola ecuación lineal de la forma: d y a· y b·u d t Si en lugar de una ecuación se tiene un modelo matemático consistente en un sistema de ecuaciones diferenciales lineales se puede estudiar la estabilidad global n m f 1 ( y1 , y2 ,..., yn , u1 , u2 ,...um ) a1i · yi b1 j ·u j d t i 1 j 1 n m d y2 f 2 ( y1 , y2 ,..., yn , u1 , u2 ,...um ) a2i · yi b2 j ·u j d t i 1 j 1 ... n m d u n y1 , y2 ,..., yn , u1 , u2 ,...um ani · yi · nj j d t i 1 j 1 y1 a11 a12 ... a1n y1 b11 b12 ... b1m u1 ... ... n m A· y B·u dt ... ... d t ... ... ... ... ... ... ... ... ... yn an1 an 2 ... ann yn bn1 bn 2 ... bnm um
d y1
• •
Si los autovalores de A son reales , la estabilidad de la solución numérica se expresa como una cota de |h| siendo el mayor autovalor de A. Si los autovalores de A son complejos el estudio de la estabilidad resulta más complejo, y se escapa del alcance de este curso. Simulación. Master en IPS 11/12
71
Estabilidad: caso general •
os au ova ores e son rea es, a es a a e a so uc n num r ca se expresa como una cota de | h| siendo el mayor autovalor de A.
• En general los métodos implícitos, aunque conllevan mas cálculo, proporcionan rangos mayores de h para obtener soluciones estables y no oscilantes Simulación. Master en IPS 11/12
72
• •
Si las ecuaciones diferenciales son lineales, hemos visto que existe un criterio para estudiar la . En general los modelos matemáticos son sistemas de ecuaciones diferenciales no lineales de la forma: –
:
f 1 ( y1 , y2 ,..., yn , u1 , u2 ,...um ) d t d y d y1 2
,
,...,
, 1,
2
1
2
n
d y En éste caso no uede onerse el modelo de la forma A· B·u t y no puede estudiarse la estabilidad por el método propuesto. • Así, resulta difícil un estudio de la estabilidad en el caso general de sistemas de ecuaciones diferenciales no lineales. Se suele estudiar la estabilidad local de una aproximación linealizada. Simulación. Master en IPS 11/12
•
73
Linealización de modelos
Simulación. Master en IPS 11/12
,...
... d yn f n ( y1 , y2 ,..., yn , u1 , u2 ,...um ) d t
74
2
m
Ecuación a linealizar
f ( x ' , x , u ) 0 x0 ' , x0 , u 0
Selección del punto de linealización
',
Se suele seleccionar un punto de linealización estac onar o, aunque no es o gator o f x ' 0 f x 0
Se linealiza el modelo:
f x
' '0
x '
f x 0
x
f u
0
0
f u 0
x0
,
0
0
u 0
f d x f f x u 0 x ' 0 dt x 0 u 0 d x dt
f x ' 0
1
f f x x x 0 ' 0
1
f u u 0
Simulación. Master en IPS 11/12
75
Ejemplo linealización dvs (t )
5·vs (t ) 2.5
dt
dv(t ) dt dy(t )
vs (t )
· ·
v(t )
dt
u (t )
dy(t ) dt
·
s
dvs (t ) dt
·
5·vs (t ) 2.5 dv(t )
2· y (t )
dt
u (t ) vs (t )
0.
·
·
·s ·
2· y (t )
v(t ) 0.
Punto de linealización estacionario
u0 , vs 0 , vs 0
0; v0 ; v0 0; y0 ; y0 0
Punto de linealización estacionario, que cumple el modelo matemático vs v0
5·vs 2.5 g
dy (t ) dt
c m
v0
u0 vs 0 c m
0 vs
vs 0 ·e
2· y0
u0
2 c
c
m
m
0 g ·0
vs 0 ·e
2· y0
m·g y0 0.5·log c·vs 0
v (t ) v0 0
Conocido el volta e de alimentación u ueda determinado el flu o de aire v posición de la pelota (y0). No son valores independientes. Simulación. Master en IPS 11/12
la 76
.
Linealización del modelo vs t dt
v
t
5 2.5·
dv(t ) dt
vs (t )
0
u (t )
vs (t )
vs (t ) 5 2.5 m
u t
5·vs (t ) 2.5
u0 vs 0
2
v
t
1 u t 0 vs (t ) 0
2.5·
0
vs (t ) 2.5 u (t ) 0 vs 0
m·g c·v(t ) c·vs (t )·e 2· y (t ) 0
mv(t ) cv(t ) (c·e 2· y (t ) ) vs (t ) (2·c·vs (t )· y (t )·e 2· y (t ) ) y (t ) 0 mv(t ) cv(t ) c·e 2· y0 vs (t ) 2·c·vs 0 · y0 ·e 2·y0 y (t ) 0 dy (t ) v(t ) 0 y(t ) v(t ) 0
Ahora resulta un modelo lineal y podemos estudiar su estabilidad numérica con respecto al paso de integración usado y al m o o e n egrac n
vs (t ) 5 2.5 v(t )
c m
u0 vs 0
v(t )
( ) 2.5 u(t ) 2 vs t
c m
vs 0
e
2· y0
c
vs (t ) 2 ·vs 0 · y0 ·e 2· y0 y (t ) m
Simulación. Master en IPS 11/12
77
Pongamos el modelo linealizado en forma matricial 0 y (t ) 2· y v(t ) 2·c·vs 0 · y0 ·e 0 m vs (t ) 0
1
c m
0
y (t ) 0 2· y0 c·e v(t ) 0 u(t ) . m ( ) v t u s v s0 5 2.5 02 v s 0 0
Si los parámetros del modelo son m=0.1, c=1 y g=9.8 Si se selecciona un punto de trabajo estacionario con u 0=5, entonces vs0=1.5811 e y 0=0. 1 0 y (t ) 0 y (t ) 0 v(t ) - 4.6883 - 10. 6.1977 v(t ) 0 u (t ) vs (t ) 0 0 - 12.90590 vs (t ) 1.5812
al punto de linealización), se necesitan conocer los autovalores de la matriz
- 4.6883 - 10. 6.1977 0 0 - 12.90590
Simulación. Master en IPS 11/12
0.4932, -9.5068 y -12.9059 78
• Ya hemos visto que si los autovalores de son reales, la estabilidad de la solución numérica se ex resa como una cota de h siendo el ma or autovalor.
En este caso los autovalores son: -0.4932, -9.5068 y -12.9059, la cota de la estabilidad de la solución numérica viene dada por 12.9059· h Si el método usado es el de Euler para que la solución numérica sea estable debe elegirse un valor del paso de integración h tal que 12.9059· h < 2, así h<0.1550 Simulemos el ejemplo (ejemplo_EULER_levit_flujo_aire.m) con distintos pasos de integración y veamos lo que pasa ... Simulación. Master en IPS 11/12
79
Solución exacta, la velocidad del flujo del aire y la altura de la pelota aumentan al aumentar el voltaje de alimentación. Observe que la velocidad del flujo del aire no presenta sobrepico.
Simulación. Master en IPS 11/12
80
H=0.1<0.155 La solución estable ero fí ese ue la velocidad del flu o del aire resenta sobrepico, cuando la exacta no lo tiene.
Simulación. Master en IPS 11/12
81
H=0.155 La solución del sistema no lineal aún es estable, pero está en el límite de dejar de serlo. Fíjese en la oscilación que presentan las variables de salida después del cambio en la entrada (la solución exacta no es así).
Simulación. Master en IPS 11/12
82
H=0.2 La solución es oscilante. Fíjese que al linealizar el modelo los autovalores nos dicen que para h>0.155 la solución es inestable, y aún para h=0.2 la solución no se ha inestabilizado. Esto es así porque el modelo lineal es una aproximación del no lineal que es el que se está simulando.
Simulación. Master en IPS 11/12
83
H=0.201 La solución, simulada con el método de Euler explícito, es inestable.
Simulación. Master en IPS 11/12
84
Aquí tiene la simulación del sistema usando el método de Euler explícito, pero con una estrategia de paso variable. Fíjese que la solución es como la exacta y a em s cuan o e s stema se esta za se pue e usar un paso e integración mucho mayor que cuando el sistema está cambiando.
Simulación. Master en IPS 11/12
“
85
”
•
A veces se presentan sistemas con dinámicas mezcladas rápidas y lentas. • Si ambas son muy diferentes, pueden presentarse problemas de integración.
• •
Se dice que un sistema es stiff si su resolución numérica por algunos métodos exige una disminución significativa del paso de integración para evitar inestabilidades. Se consideran “stiff” aquellos problemas en que el cociente entre el mayor y menor
autovalor de la linealización es superior a 3
•
dt
f ( y, u )
dt
Ay Bu
¿El sistema de levitación por flujo de aire es stiff? max e
i
min Re i
•
. 26.35 6 0.493
Moderadamente sti
Fuertement e stiff si
si
max Re i min Re i
max Re i min Re i
6
Fuertemente STIFF
¿El sistema electromecánico de las prácticas es stiff? Simulación. Master en IPS 11/12
86
3
Sistemas stiff: Ejemplo 1
dt dy 2 dt
1000 y1 10
y1 1000 0 y1 10 y 1 y2 0 2 1
Forma matricial
Condiciones iniciales :
0 1
• Autovalores -1000, -1 • san o e m o o e u er exp c o: – Para evitar la inestabilidad debe elegirse un paso de integración h tal que h<2/1000=0.002, ya que 1000 es el mayor autovalor. – Así, si quisiéramos integrar el sistema propuesto entre 0 y 5 segundos deberíamos realizar un mínimo de 5·1000/2=2500 iteraciones. – Esto representa un tiempo de cálculo muy grande la posibilidad de que los errores de redondeo acumulados sean muy grandes. – Ventajas de usar un método implícito, que no tiene esos condicionantes sobre el paso de integración son evidentes. Simulación. Master en IPS 11/12
87
• Si solucionamos el sistema de ecuaciones analíticamente: y 1 t 0.01 1 e 1000 y 2 t 0.01
0.01 999
e 1000t 0.9899 e t
• Un modo de solucionar el problema es: – Transformar la ecuación diferencial rápida en algebraica y usar la ecuación resu an e en a ecuac n en a, con o cua es a pue e n egrarse con un mayor. – Así el sistema de ecuaciones propuesto se transformará en: dy 1 1000 y 1 10 0 1000 y 1 10 y 1 0.01 dt 2
y1 y 2
2
dt dt – Ahora la solución analítica es:
0.01 y 2
y 1 ( t ) 0.01 y 2 ( t ) 0.01 0.99e t – Numéricamente se puede solucionar con el método de Euler explícito pero a ora e m te e esta a a o por e paso e ntegrac n est en h=2>>0.002. Simulación. Master en IPS 11/12 88
• La técnica expuesta es aplicable cuando puedan distinguirse claramente ecuaciones rápidas lentas, no se precise una gran precisión. • En el caso de que las ecuaciones están fuertemente acopladas no es posible utilizar directamente simplificaciones como la anteriormente expuesta. • us r mos o con e s gu en e s s ema: dy 1
500.5 y 1 499.5y 2
dy 2 499.5y 1 500.5 y 2 dt Condiciones iniciales : y 1 (0) 2; y 2 (0) 1
y1 500 .5 499 .5 y1 y 499 .5 500 .5 y 2 2
• Autovalores como en el caso anterior 1 y 1000. , ecuaciones analíticamente: – Dinámica rápidas y lentas no separables.
• Si aplicamos el mismo método que en el caso anterior: – No es válido, se necesitan métodos es ecíficos como el de Gear
2 y(0) 1
y 1 t 1.5e t 0.5e 1000t 2
dy 2
. e . e
0.
2
499.5 .
y 1 ( t ) 1.5e 1.998t 0.5 y 2 . e 1.998 t .
Simulación. Master en IPS 11/12
89
Método de Gear • • • •
Método específico para problemas stiff Formulas implícitas de diversos órdenes variables Ajuste del intervalo de integración variable La elección del orden y el paso h se hace en función de un índice teniendo en cuenta el número de asos el tiem o de cálculo extra asociado • Una vez seleccionado el orden y h, se aplica el método, lo que implica reso ver una ecuac n no- nea en genera un s s ema e ecuac ones no lineales) por medio del método de Newton-Raphson para lo cual se debe evaluar el Jacobiano. Si en tres iteraciones h no converge se disminuye h arbitrariamente y se comienza de nuevo. • El método de Gear permite usar pasos de integración que son varias veces el valor de la menor constante de tiem o del sistema inversa del mayor autovalor), lo cual hace que para sistemas stiff sea un procedimiento más rápido que cualquier otro. Simulación. Master en IPS 11/12
90
1
Lazos algebraicos • Si el modelo presenta una ecuación del tipo: – z=5·(sin(t)-ez) – Esta ecuación debe usarse para calcular z (porque t es el tiempo). – ¿Cómo despejo z? – implícitas para calcular z. • Resolvedores como los usados en los métodos implícitos – zi 1 5·sin( t ) e
z i
– Newton-Raphson zi 1 zi
F z z
1
F ( zi ) zi
i
1 1 5·e
zi
zi 5· sin( t ) e
i
– ... Simulación. Master en IPS 11/12
91
Lazos algebraicos • Si el modelo presenta un lazo algebraico lineal (ej: caudal que circula or dos tubos conectados : pe p ps – pe=2 w – po=1 – w=4·(pe-p) – w=3·(p-po) – El orden de calculo y la manipulación de las ecuaciones (que se puede automatizar) podría ser
• pe=2 • po=1 • · e · • w=3·(p-po)
o
Simulación. Master en IPS 11/12
92
•
Si el modelo presenta un lazo algebraico no lineal: – – – – –
pe=2 o=
w=4·sqrt(p e-p) w=3·sqrt(p-po) El orden de calculo y la manipulación de las ecuaciones podría ser • pe=2 • po=1 • p=(42·p e+ 32·p o)/(42+ 32) • w=4·sqrt(p-p o)
– Esta manipulación no se puede automatizar si el lazo algebraico es no-lineal (como éste), en general no queda más remedio que resolver todo el lazo numéricamente • pe=2 ; po=1 w 4 pe p o. 2 w 0. p p 0 3
Lazo algebraico no lineal a resolver numericamente (ejemplo usando Newton-Raphson) i 1 w1 estimacion inicial del caudal ; p1 estimacion inicial de la presion ; while wi 1 wi y pi 1 pi 1
1
wi 1 wi p p 2·w i i 1 i 3
wi 4 pe pi 2 pe pi wi p p i 0 3 1
end while
Simulación. Master en IPS 11/12
93
“Tearing” o ruptura de lazos: ejemplo ,
1
2
,
.
W (t ) C V 1 P1 (t ) P(t ) W (t ) C V 2 P (t ) P2 (t ) 1
2
.
Si observamos las dos ecuaciones del lazo, vemos que las dos ecuaciones son lineales en W, por lo que no queda más remedio que elegir como variable de Tearing P. Es indiferente elegir cualquiera de las dos ecuaciones como ecuación no lineal a resolver de forma numérica (seleccionemos la segunda), denominada ecuación de residuo. Para resolver el lazo se parte de una estima inicial de P=P i , con dicha estima calculamos W de la primera ecuación: , con el valor calculado de W iremos a la segunda ecuación:
W (t ) C V 2 P (t ) P2 (t )
W (t ) C V 1 P1 (t ) P (t )
0
y calcularemos P por un procedimiento iterativo (por ejemplo Newton-Raphson). Comparamos el valor de P con la estima inicial, si son suficientemente próximos hemos encontrado la solución y si no repetiremos el procedimiento con el nuevo valor de P. De este modo, aplicando el Tearing, en este sencillo ejemplo, hemos reducido un lazo no lineal de tamaño 2, a un lazo lineal de tamaño 1 (una sola ecuación implícita). ap car e earing a azos a ge ra cos e gran mens n no es tr v a e eg r e orma a ecua a as var a es e earing e entre as posibles candidatas (es mejor elegir aquellas que se conozca con mayor certeza su valor inicial, ya que vamos a iterar a partir de él). Lo ideal es que el número de variables de Tearing sea el menor posible, pues ese será el tamaño del lazo no lineal que debamos resolver. Tampoco es trivial elegir bien las ecuaciones de residuo, si es posible elegir deben elegirse aquellas que sean “bien comportadas” para garant zar una uena convergenc a num r ca. Indicar que no existen algoritmos automáticos eficientes para elegir de forma óptima las variables de Tearing y las ecuaciones de residuo, debemos usar la heurística para hacer la mejor elección.
Simulación. Master en IPS 11/12
94
Simulación y Optimización 4º Ingeniería Informática. Tema 2 Simulación
Fuente Original (“Tearing”): . " objetos y simulación de sistemas híbridos en ámbito del control de procesos químicos", Simulación yelOptimización PhD Dissertation, Facultad de Ciencias, 4º Ingeniería Informática. Tema 2 Simulación UNED
95
96
2.4 Problemas DAE •
Muchos problemas se formulan como ecuaciones implícitas
(DAEs implícitos) donde no es posible despejar el modelo matemático •
F (Y , Y , t ) 0
También como conjuntos de ecuaciones diferenciales y algebraicas acopladas (DAE semi-implícitas)
F Y , Z , t 0 G (Y , Z , t )
Y – • Forma secuencial:
– Dado y(t) se resuelve g(y(t), z(t),t) = 0 ==> se obtiene z(t) – Empleando un método explícitode resolución deODE a y' = , , • Forma simultánea
– Resolver y' = f(y,z,t), g(y, z,t)=0 de forma simultánea mediante métodos implícitos (BDF)
• Veamos algún ejemplo, el orden de un DAE y como se resuelven ... Simulación. Master en IPS 11/12
97
, •
Se define el índice/orden de un DAE como el número de derivaciones requeridas .
• • •
Las ODEs son DAEs de orden cero. Los DAEs de orden 2 o más se llaman High-Index Problem E emp o:
x2 (t ) x1 (t )
t t · x3 x3 (t ) x3 (t ) 1 x2 t
• •
t
Veamos cual es el orden del DAE y como reducirlo … (explicado en el aula). Existen algoritmos generales para resolver DAEs de orden 1, pero no existen algoritmos generales para resolver DAEs de orden 2 o superior.
Simulación. Master en IPS 11/12
98
Resolución de DAEs por linealización F (Y , Z ) Y 0 G (Y , Z )
Y F (Y , Z ) ( , ) G Y Z R
Forma matricial
Nomin almente R=0.
G Y G Z
• • • •
Linealizando
Entonces
F Y
Y F y R G y
F z Y
G z Z
F G 1G
Y
Estimamos los jacobianos (Fy, Fz, Gy y Gz). Resolvemos la ecuación diferencial obtenida (Y). Resolvemos la ecuación algebraica G(Y,Z)=0, y obtenemos Z. Ajustamos los valores de Gy y Gz para que los residuos sean cero. F Y Y t
Y
0 F y ·Y F y ·Y F ·t 0 F y y
·Y
F t
·t
y
Simulación. Master en IPS 11/12
99
Resolución de DAEs de orden 1 por método implícitos Caso particul ar de
F (Y , Z ) 0 Y 0 G (Y , Z )
, Y , t ) 0 H (Y
• Reemplazará la derivada de Y por una dy y aproximación en diferencias de primer t t orden (BDF, backward difference formula). Y Y • Resultando: H n 1 n Y t 0 Siendo : h t t hn 1 • Debiendo resolver la ecuación implícita, por ejemplo, por el • DASSL ( Differential Algebraic System Solver ) utiliza otras fórmulas BDF de coeficientes fijos para aproximar la derivada con ordenes mayores y con paso variable para obtener una precisión mayor. • El método de Gear también permite resolver DAEs semiSimulación. Master en IPS 11/12 100 implícitos.
Resolución de DAEs: método de RADAU5 • Si el modelo matemático en forma de DAEs (de orden 1) es: , Y , t ) 0 H (Y • Pue e ponerse e a s gu ente orma:
, •
·
,
ntonces pue e usarse e m to o Kutta implícito) para resolverlo.
m to o unge-
Resolución de DAEs de orden superior • Si el orden del DAE es 2 o su erior no existen métodos generales para su resolución, debe reducirse el orden del DAE a orden 0 o 1, con el procedimiento antes indicado para poder resolverlo como una ODE o con un DAE de orden 1. 101
Simulación. Master en IPS 11/12
2.5 Eventos y Discontinuidades
• Los eventos implican una discontinuidad en y(t) o en su derivada y´(t)=f(y(t),u(t)) • Discontinuidad en la variable y(t) implica detener la integración y reiniciarla en el nuevo valor de y(t) • Discontinuidad en la función derivada: – La aparición del evento hace que f(y,u) pase en ese instante de tiempo de ser f u a ser f u – En estos casos la aplicación directa de los métodos de integración mencionados es incorrecta f(y,u)
t h
y (t h) y (t )
f ( y ( ), u ) d f 1(y,u)
t
d y d t d y d t
Estim ación de f(y,u)
f 1 ( y, u )
ti empo t d
f 2 ( y, u )
t iempo t d -
2
,
-
Simulación. Master en IPS 11/12
102
• La integración correcta exige: – oca zar e ns an e e scon nu a – Reinicializar la integración con las nuevas ecuaciones y condiciones iniciales a partir de ese punto f(y,u) ,
f 2(y,u)
f 1(y,u)
t d
y (t d ) y (t ) f 1 ( y ( ), u ( ))d t h
y (t h) y (t d ) f 2 ( y ( ), u ( ))d t-2h
t-h
t
t+d
t+h
Simulación. Master en IPS 11/12
103
Tipos de eventos • Eventos en el tiempo, ocurren en un instante determinado de tiem o. – Periódicos (ej: periodo de muestreo de un controlador). – Aperiódicos, aquellos que suceden un tiempo después de que haya suce o otro evento.
• Eventos en el estado, ocurren cuando alguna condición en . – Pueden distinguirse tres tipos de eventos en el estado, aquellos en los que la condición se satisface en la dirección positiva, negativa o en ambas direcciones. • x=0 • = • x=0 y dx/dt>0
Simulación. Master en IPS 11/12
104
• Los eventos en el tiempo son fáciles de tratar porque se conoce a priori cuando van a suceder, se ajusta el tamaño del paso de integración al instante en el que va a . • Un posible algoritmo para la detección de los eventos en el estado es el siguiente: – Al finalizar cada paso de integración se comprueba si ha sucedido el evento. – Si el evento ha sucedido, para determinar el instante en el que sucedió, se toma el valor de la variable asociada al evento en el intervalo de integración previo y el actual, y se interpola para conocer el instante del cambio. – e ap ca e nuevo e a gor mo oman o como paso e n egrac n a erenc a en re el instante anterior y el instante en el que se produce el evento.
(x(t)) dt f 1 ( x, u ) si ( x) 0 dx x, u si x 0
t+h
t
(x(t+h)) t+d
Simulación. Master en IPS 11/12
105
• La simulación de modelos híbridos (con eventos) plantea problemas teóricos y computacionales que no se presentan en los modelos continuos. •
“
”
– el número de eventos en el estado es grande comparado con el número de pasos de integración.
•
, – La detección de los eventos. – La determinación de su instante de disparo. – a e ecuc n e os even os. – La solución del problema de re-inicio tras la ejecución de cada evento.
• La única solución al “chattering” es modificar el modelo. • Veamos un ejemplo
Simulación. Master en IPS 11/12
106
dv (t )
• Modelo simplificado del bote de la pelota
dt dy (t )
g
v(t ) dt v (t ) k ·v (t )
•
= Simulación. Master en IPS 11/12
•
¿Por qué la pelota cae por debajo de x=0? –
ons erac ones: • El evento se activa cuando no está activo (false) y se cumple la condición lógica • El evento se desactiva cuando está activo y se deja de cumplir la condición lógica
–
cuando y(t) 0
, pero no podemos representar el cero exactamente. Se dice que se cumple el evento cuando en valor absoluto y(t) sea menor que epsilon.
•
107
A medida que pasan lo sucesivos botes llegamos a valores de y(t) próximos a epsilon y puede suceder que:
– – – –
Detectemos el evento cuando y(t) cruce por el valor inferior de la banda (-epsilon). En ese momento el evento se pone a true (y no se desactivará hasta que y(t) sea mayor que epsilon), la velocidad cambia de signo y la pelota empieza a subir. Pero la velocidad inicial es pequeña y la pelota se va frenando por la fuerza de la gravedad sin llegar a epsilon, y sin desactivar el evento por tanto. En un momento dado la pelota empieza a bajar y vue ve a cruzar por –eps on, pero e even o no se activa (ya que está activo), la velocidad no cambia de signo y la pelota cae de forma ilimitada
Simulación. Master en IPS 11/12
108
g dt 0
dv (t )
• Solución, modificar el modelo – Añadir que v´(t)=0., si y(t)<0.
si y (t ) 0 si y (t ) 0
dy (t )
v(t ) dt v (t ) k ·v (t )
cuando y(t) 0
La simulación se ralentiza
La pelota ya no se cae
Simulación. Master en IPS 11/12
109
Ejemplo Chattering
Simulación. Master en IPS 11/12
110
