Estudio del efecto de incertidumbres en la trayectoria de una aeronave

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Estudio del efecto de incertidumbres en la trayectoria de una aeronave Autor: Rubén García Díaz Tutor: Rafael Vázquez Valenzuela Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Estudio del efecto de incertidumbres en la trayectoria de una aeronave Autor: Rubén García Díaz Tutor: Rafael Vázquez Valenzuela Profesor Titular Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Trabajo Fin de Grado: Estudio del efecto de incertidumbres en la trayectoria de una aeronave Autor: Rubén García Díaz Tutor: Rafael Vázquez Valenzuela El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Índice Índice I Índice de Figuras III Índice de Tablas V 1 Introducción 1 1.1 Motivación 1 1.2 Objetivo y metodología 1 1.3 Estructura del documento 2 2 Método de GPC 5 2.1 Descripción 5 2.1.1 El caos polinómico de Wiener 5 2.1.2 El caos polinómico de Wiener-Askey 6 2.1.3 Correspondencia entre las variables aleatorias y el tipo de polinomio de Wiener-Askey 6 2.1.4 Propiedad de los polinomios ortogonales 7 2.2 Ejemplo: Aplicación del caos polinómico a la ecuación de crecimiento exponencial 8 2.2.1 Resultados numéricos 9 3 Modelo de variación de masa en crucero 13 4 Incertidumbre en m015 4.1 Distribución uniforme 15 4.1.1 Función de densidad 15 4.1.2 Cálculo de los coeficientes del GPC 15 4.1.3 Esperanza 20 4.1.4 Desviación típica 21 4.1.5 Comparación con el método de Montecarlo 22 4.2 Distribución gamma 23 4.2.1 Función de densidad 23 4.2.2 Cálculo de los coeficientes del GPC 24 4.2.3 Esperanza 25 4.2.4 Desviación típica 26 5 Incertidumbre en CD029 5.1 Cálculo de los coeficientes del método GPC 30 5.2 Esperanza 31 5.3 Desviación típica 32 5.4 Comparación con el método de Montecarlo 33 6 Incertidumbre en c 35 6.1 Cálculo de los coeficientes del GPC 36 6.2 Esperanza 37 I II Índice 6.3 Desviación típica 39 7 Incertidumbre en m0yCD043 7.1 Cálculo de los coeficientes del GPC 43 7.2 Esperanza 45 7.3 Desviación típica 47 7.4 Comparación con el Método de Montecarlo 48 8 Incertidumbre en m0,CD0yc51 8.1 Cálculo de los coeficientes del GPC 51 8.2 Esperanza 52 8.3 Desviación típica 54 9 Incertidumbre en m0,CD0,cyCD259 9.1 Cálculo de los coeficientes del GPC 60 9.2 Esperanza 62 9.3 Desviación típica 63 9.4 Comparación con el método de Montecarlo 64 10Conclusiones 67 Bibliografía 69 Índice de Figuras 2.1 Esquema de Askey de polinomios ortogonales 7 2.2 Valor de los coeficientes h para el caso ejemplo 10 2.3 Evolución de la media y la solución determinista a lo largo del tiempo para el caso ejemplo 10 2.4 Evolución de la desviación típica a lo largo del tiempo para el caso ejemplo 11 4.1 Función de densidad para distribución uniforme m016 4.2 Coeficientes para distribución uniforme con incertidumbre en m0para P=317 4.3 Coeficientes para distribución uniforme con incertidumbre en m0para P=518 4.4 Tiempo necesario para resolver el sistema de ecuaciones diferenciales en función de P, para el caso de incertidumbre en m0, distribución uniforme 19 4.5 Espacio que ocupa en memoria la variable que guarda los coeficientes h en función de P, para el caso de incertidumbre en m0, distribución uniforme 19 4.6 Valor esperado de la masa lo largo del tiempo. Distribución uniforme m020 4.7 Desviación típica de la masa lo largo del tiempo. Distribución uniforme m022 4.8 Cociente entre la desviación típica y el combustible consumido. Distribución uniforme m023 4.9 Función de densidad para distribución gamma m024 4.10 Coeficientes para distribución gamma con incertidumbre en m0para P=325 4.11 Valor esperado de la masa lo largo del tiempo. Distribución gamma m026 4.12 Desviación típica de la masa lo largo del tiempo. Distribución gamma m027 4.13 Cociente entre la desviación típica y el combustible consumido. Distribución gamma m028 5.1 Función de densidad para distribución uniforme CD029 5.2 Coeficientes para distribución uniforme con incertidumbre en CD031 5.3 Valor esperado de la masa lo largo del tiempo. Distribución uniforme CD032 5.4 Desviación típica de la masa lo largo del tiempo. Distribución uniforme CD033 5.5 Cociente entre la desviación típica y el combustible consumido. Distribución uniforme CD034 6.1 Función de densidad para distribución uniforme c35 6.2 Coeficientes para distribución uniforme con incertidumbre en c37 6.3 Valor esperado de la masa lo largo del tiempo. Distribución uniforme c38 6.4 Desviación típica de la masa lo largo del tiempo. Distribución uniforme c39 6.5 Comparación de la desviación típica para los casos de incertidumbre en m0,CD0yc40 6.6 Cociente entre la desviación típica y el combustible consumido. Distribución uniforme c 40 6.7 Comparación entre el cociente de la desviación típica y el combustible consumido. Para los casos de incertidumbre en m0,CD0yc41 7.1 Coeficientes para distribución uniforme con incertidumbre en m0yCD045 7.2 Tiempo necesario para resolver el sistema de ecuaciones diferenciales en función de P, distribución uniforme m0yCD045 7.3 Valor esperado de la masa lo largo del tiempo. Distribución uniforme m0yCD046 7.4 Desviación típica de la masa lo largo del tiempo. Distribución uniforme m0yCD048 7.5 Comparación de la desviación típica de la masa lo largo del tiempo. Distribución uniforme m0yCD049 III 2 Método de GPC Un proceso estocástico es una variable aleatoria que varía en el tiempo. Tal y como puede ser la solución de una ecuación diferencial cuya condición inicial es aleatoria o que contiene alguna variable aleatoria, son ejemplos de procesos estocásticos.Este será el tipo de procesos considerado en este trabajo. El uso de procesos estocásticos permite modelar la existencia de incertidumbres en diferentes aplicaciones, pero al mismo tiempo complica la obtención de soluciones, ya que en principio habría una solución diferente para cada posible valor de condiciones iniciales y/o variables. Por tanto, hay que realizar un tratamiento estadístico de los procesos y hablar de valores tales como distribución, media o covarianza del proceso. El método “Generalized Polynomial Chaos”, también llamado GPC en la literatura, es un procedimiento que permite analizar la evolución en el tiempo de este tipo de funciones. En concreto, es posible mediante GPC obtener estimaciones de media y varianza del proceso a lo largo de su evolución. El método se basa en un desarrollo en serie (multidimensional, si hubiera más de un parámetro aleatorio afectando la ecuación) del proceso, descomponiéndolo en una suma de funciones aleatorias. Estas funciones son polinomios ortogonales, elegidos en función de las distribuciones implicadas, dando parte del nombre al método. Se denomina generalizado porque pueden considerarse diferentes distribuciones (originalmente, sólo podía tratar variables gaussianas). Finalmente, el término caos hace referencia a la aleatoriedad de las variables y no tiene relación alguna con la posible dinámica caótica de un sistema. 2.1 Descripción 2.1.1 El caos polinómico de Wiener Con el título de "The Wiener polynomial chaos" fue conocido originalmente el método de GPC [ 4 ], como se ha dicho anteriormente, solo era aplicable a variables gaussiana. Según este método, un proceso estocástico X(θ) en función de una variable aleatoria θ puede ser escrito como el desarrollo en serie de un coeficiente aipor el polinomio ortogonal de Hermite Hn(ξi1,...,ξin): X(θ) = a0H0 + ∞ ∑ i1=1 ai1H1(ξi1)(θ)) + ∞ ∑ i1=1 i1 ∑ i2=1 ai1i2H2(ξi1(θ),ξi2(θ)) + ∞ ∑ i1=1 i1 ∑ i2=1 i2 ∑ i3=1 ai1i2i3H3(ξi1(θ),ξi2(θ),ξi3(θ))+... (2.1) Donde Hn(ξi1,...,ξin) denota el polinomio ortogonal de Hermite de orden n en las variables (ξi1,...,ξin) y donde el término Hn hace referencia al polinomio de Hermite en términos de la variable Gaussiana ξ con media cero y varianza unidad. La fórmula anterior es para el caso discreto; en caso de tener variables continuas, sólo hay que sustituir los sumatorios por integrales. 5 6Capítulo 2. Método de GPC Por comodidad, la ecuación (2.1) puede ser reescrita como: X(θ) = ∞ ∑ j=0 ajΨj(ξ)(2.2) Donde las funciones Hn(ξi1,...,ξin)yΨj(ξ)están relacionadas uno a uno. Este desarrollo en serie es posible, debido a la semejanza que existe entre la función de distribución Gaussiana y la función de peso de los polinomios de Hermite. En caso de que la variable no fuese Gaussiana, no se tendría una convergencia óptima. 2.1.2 El caos polinómico de Wiener-Askey Posteriormente, este método que solo era válido para variables gaussianas, se extendió a otros tipos de distribuciones estadísticas, dando lugar a la generalización del método, también conocido como " the Wiener-Askey polynomial chaos". Esta extensión a otro tipo de variables estadísticas y por tanto, a otros tipos de polinomios ortogonales, dan lugar a lo que se conoce como el esquema de Askey. Como se observa en la figura 2.1 este esquema no es mas que una agrupación de los tipos de polinomios ortogonales, en los laterales de la figura, se muestra la nomenglatura de las funciones hipergeométricas. Como en la sección anterior, un proceso X(θ) en función de una variable aleatoria θ puede escribirse como un desarrollo en serie de un coeficiente cipor una función que es un polinomio ortogonal In(ζi1,...,ζin): X(θ) = a0I0+ ∞ ∑ i1=1 ci1I1(ζi1(θ))+ ∞ ∑ i1=1 i1 ∑ i2=1 ci1i2I2(ζi1(θ),ζi2(θ))+ ∞ ∑ i1=1 i1 ∑ i2=1 i2 ∑ i3=1 ci1i2i3I3(θ),ζi2(θ),ζi3(θ))+... (2.3) Donde In(ζi1,...,ζin) denota el caos polinomial de orden n en términos del vector aleatorio ζ= (ζi1,...,ζin) .La diferencia con el caso anterior, es que In no está limitado a los polinomios de Hermite, sino que puede ser cualquier tipo de polinomio del esquema de Askey (Figura 2.1). Por comodidad, la expresión (2.3) puede ser reescrita como: X(θ) = ∞ ∑ j=0 cjΦj(ζ)(2.4) Sin embargo, en la práctica no se emplea este desarrollo, pues a la hora de introducirlo en un programa de cálculo numérico como puede ser Matlab, no se puede hacer un cálculo con infinitos términos, de manera que se realiza una aproximación, quedando la ecuación (2.4) de la siguiente forma: X(θ)≈ P ∑ j=0 cjΦj(ζ)(2.5) Donde P, debe ser un número de tal forma que en la aproximación, el error que se cometa al aproximar la variable, se encuentre por debajo de una cierta tolerancia. 2.1.3 Correspondencia entre las variables aleatorias y el tipo de polinomio de Wiener-Askey En el apartado "El caos polinomial de Wiener-Askey", se vio la relación existente entre el tipo de variables aleatorias y el tipo de polinomio ortogonal. Ésto se debe a que algunos tipos de polinomios ortogonales procedentes del esquema de Askey tienen funciones pesos con la misma forma que la función de densidad de ciertos tipos de distribuciones aleatorias. La relación entre estas variables aleatorias y el tipo de polinomio ortogonal que mejor aproxima su función 2.1 Descripción 7 Figura 2.1 Esquema de Askey de polinomios ortogonales. de densidad, se muestra en la tabla 2.1. Donde la primera columna indica si la variable aleatoria es discreta o continua. La segunda columna indica el tipo de distribución aleatoria. La tercera columna muestra el tipo de polinomio ortogonal que mejor aproxima a esa variable aleatoria. La cuarta columna indica el intervalo en los que puede tomar valores esos polinomios. Es importante destacar que, en el caso de los polinomios de Legendre, que es un caso particular de los polinomios de Jacobi (con α=β=0 ) se corresponden con una variable aleatoria muy importante, la uniforme y por este motivo, se considera aparte en la tabla 2.1. Además, será la variable aleatoria mas usada en este trabajo Por lo tanto, en la práctica, se selecciona el tipo del polinomio, en función de la variable aleatoria que se desea estudiar, de acuerdo con lo mostrado en la tabla 2.1. 2.1.4 Propiedad de los polinomios ortogonales Cualquier sistema de polinomios ortogonales {Qn(x),n∈N} , donde Qn(x) es un polinomio de orden n y con N=0,1,2,...,Ncumple la siguiente propiedad: Z D Qn(x)Qm(x)w(x)dx =h2 nδnm,n,m∈N(2.6) Siendo el dominio D, el soporte que figura en la tabla 2.1. La expresión (2.6) es válida para los casos en que 8Capítulo 2. Método de GPC Tabla 2.1 Relación entre el caos polinomial de Wiener-Askey y el tipo de variable aleatoria. Variable aleatoria Caos Wiener-Askey Soporte Continua Gaussiana Caos Hermite (−∞,∞) Gamma Caos Laguerre [0,∞) Beta Caos Jacobi [a,b] Uniforme Caos Legendre [a,b] Discreta Poisson Caos Charlier 0,1,2,... Binomial Caos Krawtchouck 0,1,2,...,N Binomial Negativa Caos Meixner 0,1,2,... Hipergeométrica Caos Hahn 0,1,2,...,N las funciones sean continuas, para el caso en que las funciones sean discretas, se aplica la expresión (2.7), donde es posible que M=∞. M ∑ i=0 Qn(xi)Qm(xi)w(xi) = h2 nδnm,n,m∈N(2.7) La importancia de esta propiedad, radica en que la integral de la multiplicación de dos polinomios ortogonales entre sí, vale 0, excepto en el caso de que esos dos polinomios sean del mismo grado. Esta propiedad, se podrá utilizar en las ecuaciones que tienen que cumplir los coeficientes para simplificarlas y poder calcularlos de forma menos costosa. 2.2 Ejemplo: Aplicación del caos polinómico a la ecuación de crecimiento exponencial Con este sencillo ejemplo, extraido de [ 5 ] se pretende ilustrar el procedimiento a seguir para obtener la media y la desviación típica de una ecuación diferencial cuando sus coeficientes son aleatorios. En este caso, resolveremos la ecuación diferencial de crecimiento exponencial, en la cual y(t) representa la población de una especie determinada a lo largo del tiempo y res la tasa de crecimiento: dy dt =ry(t)(2.8) La solución determinista a esta ecuación, para una condición inicial y(0) = y0es: y(t) = y0ert (2.9) En este ejemplo se quiere modelar la tasa de crecimiento r como una variable aleatoria uniforme de media r0 y semiancho δr . Haciendo uso del caos polinomial, los pasos a seguir para resolver (2.8) son los siguientes: Se expresa la variable y(t) en función de unos coeficientes a calcular y del caos polinomial correspondiente. En este caso, en función de los polinomios de Legendre: y(t;r)≈ P ∑ i=0 hiLi(∆)(2.10) Donde hi tiene que ser calculado mediante la ecuación (2.8) y Li representa el polinomio i de Legendre.La variable uniforme rpuede ser reescrita en función de los polinomios ortogonales de Legendre: r=r0L0(∆)+δrL1(∆)(2.11) 2.2 Ejemplo: Aplicación del caos polinómico a la ecuación de crecimiento exponencial 9 Sustituyendo las ecuaciones (2.10) y (2.11) en (2.8) se obtiene: P ∑ i=0 ˙ hi(t)Li(∆)≈ P ∑ i=0 (r0L0(∆)+δrL1(∆))hi(t)Li(∆)(2.12) El siguiente paso consiste en multiplicar la ecuación (2.12) por Ll(∆) para l=0,...,P , tomar esperanza respecto a ∆y aplicar la propiedad de los polinomios ortogonales anteriormente descrita, obteniendo P+1 ecuaciones: ˙ hi(t)E[L2 l(∆)] = P ∑ i=0 (r0E[L0(∆)Li(∆)Ll(∆)]+δrE[L1(∆)Li(∆)Ll(∆)])hi(t)(2.13) Reordenando y llamando C0 i,l=E[LiLlL0] E[L2 l] y C1 i,l=E[LiLlL1] E[L2 l] queda el siguiente sistema de ecuaciones con P+1 ecuaciones diferenciales, que permite calcular los coeficientes hi . Como se observa en (2.14), los coeficientes verifican ecuaciones lineales, pues el problema inicial también es lineal. ˙ hl= P ∑ i=0 hi(r0C0 i,l+δrC1 i,l)(2.14) siendo las condiciones iniciales: h0(0) = y0,hl(0) = 0,l=1,...,P(2.15) Una vez obtenidos los coeficientes, la media y la desviación típica se calculan como sigue: E[y(t;r)] = P ∑ i=0 hi(t)E[Li(∆)] = P ∑ i=0 hi(t)E[Li(∆)L0(∆)] = h0(t)E[L2 0(∆)] = h0(t)(2.16) Donde se ha usado la propiedad de los polinomios ortogonales. Para el caso de la varianza, quedaría de la siguiente forma: Var[y(t;r)] = E[y2(t;r)] −E[y(t;r)]2= P ∑ i=0 P ∑ j=0 hi(t)hj(t)E[Li(∆)Lj(∆)]−h2 0= P ∑ i=1 h2 i(t)E[L2 i(∆)] (2.17) 2.2.1 Resultados numéricos Una vez descrito el procedimiento para obtener la media y la varianza, se completará el ejemplo con resultados numéricos. Para ello, se van emplear los siguientes valores: Como condición inicial, y0=1 . La tasa de crecimiento, como se dijo en el apartado anterior, va a seguir una distribución uniforme de media r0=2 y semiancho δr=0.2 . Además, P=3 , para garantizar una precisión adecuada. Esta precisión es adecuada, debido a que los coeficientes hi , disminuyen su orden de magnitud en dos unidades respecto al anterior, por lo que con P=3, se tiene una diferencia de 6 ordenes de magnitud entre h0yh3 Teniendo ya los datos iniciales y empleando la ecuación diferencial (2.14) junto a las condiciones iniciales de (2.15), se obtienen fácilmente los coeficientes que permiten modelar la variable y , como se muestra en la figura 2.2. Para obtener el valor de la media en todo instante de tiempo, se aplica (2.16). Además, se representa sobre la misma gráfica el valor de la solución determinista (sin considerar r aleatoria) para cada instante de tiempo, con objeto de comparar ambas soluciones. Los resultados se pueden ver en la figura 2.3. Por ultimo, con (2.17) se obtiene la evolución de la varianza a lo largo del tiempo: basta con hacerle la raíz cuadrada a la varianza para tener la evolución de la desviación típica, como se muestra en la figura 2.4. 10 Capítulo 2. Método de GPC Figura 2.2 Valor de los coeficientes h para el caso ejemplo. Figura 2.3 Evolución de la media y la solución determinista a lo largo del tiempo para el caso ejemplo. 2.2 Ejemplo: Aplicación del caos polinómico a la ecuación de crecimiento exponencial 11 Figura 2.4 Evolución de la desviación típica a lo largo del tiempo para el caso ejemplo. 3 Modelo de variación de masa en crucero En este capítulo, se muestra el modelo que permite calcular la variación de la masa de un vuelo en crucero a lo largo del tiempo. Dicho modelo será el que se emplee en los siguientes capítulos. Se usará el modelo empleado en [ 2 ] y [ 3 ] con el objeto de poder reproducir los mismos resultados que se obtienen en dichos documentos. Las ecuaciones empleadas son las del vuelo simétrico en el plano vertical (con rumbo constante), con las simplificaciones de tierra plana para altitud y velocidad constante: dx dt =V(3.1) dm dt =−cT (3.2) T=D(3.3) L=mg (3.4) Siendo x la distancia recorrida, t el tiempo, V es la velocidad de la aeronave, T es el empuje, D es la resistencia aerodinámica, L es la sustentación, m es la masa de la aeronave, g es la aceleración de la gravedad y c es el consumo específico, el cual puede ser una función de la altitud y la velocidad; al ser éstas dos constantes en nuestro modelo, ctambién lo es. La resistencia aerodinámica puede escribirse como: D=1 2ρV2SCD(3.5) Siendo ρ la densidad, S es la superficie alar y CD es el coeficiente de resistencia aerodinámica. Para modelarlo se usará el modelo de polar parabólica de coeficientes constantes CD=CD0+CD2C2 L , siendo CD0 y CD2 constante y CLes el coeficiente de sustentación, cuya expresión es: CL=2L ρV2S(3.6) Combinando las ecuaciones anteriores, se obtiene la siguiente expresión que permite calcular la evolución de la masa en función del tiempo: dm dt =−c(1 2ρV2SCD0+m22CD2g2 ρV2S)(3.7) Que puede ser reescrita como: dm dt =−(A+Bm2)(3.8) 13 20 Capítulo 4. Incertidumbre en m0 4.1.3 Esperanza Una vez conocido los coeficientes del desarrollo en serie, se expresa la evolución de la esperanza a lo largo del tiempo en función de dichos coeficientes. Para ello, se hace uso de la definición de esperanza, aplicada a la ecuación (4.2). Como resultado, se obtiene la siguiente expresión: E[m(t;m0)] = P ∑ i=0 hi(t)E[Li(∆)] = P ∑ i=0 hi(t)E[Li(∆)L0(∆)] = h0(t)E[L2 0(∆)] = h0(t)(4.7) Para llegar a este resultado, se ha multiplicado por L0(∆) y se ha usado la propiedad de los polinomios ortogonales, dando lugar a la simplificación de la expresión. Por último, se ha usado que L0(∆) = 1 , llegando así a la expresión final. En la tabla 4.2. se muestran la evolución de la esperanza para diferentes valores de tiempo, siendo la primera columna los instantes de tiempo en los que se muestra la esperanza. La segunda columna representa los valores de la esperanza en dicho instantes de tiempo mediante el método de GPC, en la tercera columna se tiene los valores obtenidos en el artículo de Rafael Vázquez y Damian Rivas [ 2 ] en el que se emplea este mismo método. Por último, en la cuarta columna se muestran los resultados obtenidos en el Trabajo Fin de Grado de Manuel Ángel Zapata Habas [ 3 ] usando el método de Montecarlo. En la figura 4.6 se representa la evolución de la esperanza de la masa a lo largo del tiempo. Como es de esperar, disminuye, debido al consumo de combustible a lo largo del vuelo en crucero. Figura 4.6 Valor esperado de la masa lo largo del tiempo. Distribución uniforme m0. Como se puede observar, los valores obtenidos para la media, son muy similares a los obtenidos en el artículo y en el Trabajo Fin de Grado, por lo que se puede concluir que P=3 proporciona un excelente grado 4.1 Distribución uniforme 21 Tabla 4.2 Valores de la esperanza para distintos instantes de tiempos, Distribución uniforme m0 . Siendo GPC*, los resultados obtenidos por Rafael Vázquez Valenzuela y Damian Rivas Rivas usando el método de GPC. tiempo (s) E[m(t;m0)] (Kg).GPC E[m(t;m0)] (Kg).GPC* E[m(t;m0)], (Kg) MC 2000 77485.6 77485 77485.1 4000 73477.1 73477 73477.2 6000 69595.9 69596 69595.1 8000 65831.7 65831 65831.9 10000 62174.8 62175 62174.8 12000 58616.5 58616 58616.6 de precisión. No obstante, para justificar su elección, en la tabla 4.3. se compara los resultados obtenidos para diferente P en los mismos instantes de tiempo. Siendo la primera columna los instantes de tiempo, la segunda la esperanza para P=3 y la tercera la esperanza para P=5 . La variación de P=3 a P=5 se produce en el último decimal, esto significa que con P=3se tiene una precisión de 15 cifras significativas. Tabla 4.3 Valores de la esperanza para P=3 y P=5 , Distribución uniforme m0. tiempo (s) E[m(t;m0)](Kg) P=3E[m(t;m0)] (Kg)P=5 2000 77485.5991161437577485.59911614374 4000 73477.07185930982 73477.07185930982 6000 69595.9313042597469595.93130425975 8000 65831.7096090836765831.70960908369 10000 62174.8363654595862174.83636545959 12000 58616.53319849474 58616.53319849474 4.1.4 Desviación típica En primer lugar, se obtiene la expresión de la varianza en función de los coeficientes hi tal y como se muestra en la ecuación (4.8). En esta ecuación, se vuelve a hacer uso de la propiedad de los polinomios ortogonales y de la relación entre la masa y los coeficientes hide la ecuación (4.2). Var[m(t;m0)] = E[m2(t;m0)] −E[m(t;m0)]2= P ∑ i=0 P ∑ j=0 hi(t)hj(t)E[Li(∆)Lj(∆)]−h2 o= P ∑ i=0 h2 i(t)E[L2 i(∆)] (4.8) Para obtener la evolución de la desviación típica a lo largo del tiempo, basta con hacer la raiz cuadrada a la ecuación (4.8). En la figura 4.7. se puede observar dicha evolución, que tal y como se puede observar, tiende a disminuirse con el tiempo. En la tabla 4.4. se muestran los valores de la desviación típica para diferentes valores de tiempo, siendo la primera columna dichos instantes de tiempo. La segunda columna representa los valores de la desviación típica obtenidos mediante GPC. La tercera columna son los valores obtenidos en el artículo de Rafael Vázquez Valenzuela y Damian Rivas Rivas. La cuarta columna son los valores obtenidos por Manuel Ángel Zapata Habas empleando el método de Montecarlo. Como se pueden ver, los valores obtenidos por GPC son bastante precisos, existiendo un error menor que 0.1 Kg con respecto al método de Montecarlo. Por último, se representa el coeficiente de variación que no es más que el cociente entre la desviación típica y el consumo de combustible. Como se puede observar en la figura 4.8. este coeficiente tiende a disminuir con 22 Capítulo 4. Incertidumbre en m0 Figura 4.7 Desviación típica de la masa lo largo del tiempo. Distribución uniforme m0. Tabla 4.4 Valores de la desviación típica para distintos instantes de tiempos, Distribución uniforme m0 . Siendo GPC*, los resultados obtenidos por Rafael Vázquez Valenzuela y Damian Rivas Rivas usando el método de GPC. tiempo (s) σ(t;m0)] (Kg).GPC σ[m(t;m0)] (Kg).GPC* σ[m(t;m0), (Kg) MC 2000 2787.7 2787 2787.6 4000 2696.8 2696 2696.6 6000 2613.5 2613 2613.3 8000 2536.9 2536 2536.8 10000 2466.6 2467 2466.5 12000 2402.1 2402 2402.2 el tiempo, lo cual significa que la incertidumbre en la masa inicial tiende a reducirse con el tiempo. 4.1.5 Comparación con el método de Montecarlo En este apartado, se compara el coste de implementar computacionalmente el método de Montecarlo y el método de GPC. Para realizar esta comparación, se tendrá en cuenta el tiempo que emplea cada método en obtener la desviación típica y el coste que tiene en memoria. Para el método de Montecarlo [ 3 ], tiene un tiempo de ejecución de 3.5 segundos. En cuanto al coste en memoria, este método necesita generar 33554432 muestras para calcular la esperanza lo que tiene un coste de almacenamiento de 134217728 bytes. Para el cálculo de la desviación típica se usan el mismo número de muestras, sin embargo, en cada iteración tiene que guardar las muestras de las iteraciones anteriores, siendo el coste total de almacenamiento de 268435440 bytes. El coste total de almacenar estas dos variables es de 402653168 bytes. 4.2 Distribución gamma 23 Figura 4.8 Cociente entre la desviación típica y el combustible consumido. Distribución uniforme m0. Para el método de GPC, se ha calculado la esperanza y desviación típica cada 10 segundos en un intervalo temporal de 15000 segundos. El tiempo de ejecución para estas condiciones es de 0.35 segundos. La variable que almacena los coeficientes hi tiene 1501xP valores, en este caso P=3 por lo que se tiene un total de 4503 valores, ocupando un espacio en memoria de 48032 bytes. Como se ha visto anteriormente la esperanza es igual al coeficiente h0 por lo que el coste de esta variable ya está incluido en la variable que almacena los coeficientes. Para el caso de la desviación típica, esta variable tiene 1501 valores y ocupa un espacio en memoria de 12008 bytes. Por tanto, el coste total es de 60040 bytes. Para este caso, se ha demostrado que el método de GPC es más eficiente que el de Montecarlo ya que el tiempo de ejecución es 10 veces mas rápido mientras que el espacio en memoria es unas 2500 veces menor. Además otro aspecto a tener en cuenta es que usando el método de Montecarlo, se tiene una expresión analítica de la evolución de la masa a lo largo del tiempo. Si no se tuviese esta expresión, tendría que resolver un número de ecuaciones diferenciales igual al número de muestras. Mientras que en este trabajo no se ha hecho uso de esa expresión analítica para encontrar la solución. 4.2 Distribución gamma 4.2.1 Función de densidad Para este apartado, se ha tomado una distribución gamma para la masa inicial. Esta función G(k,θ) viene definida por dos parámetros. El parámetro k es el parámetro de forma y θ es el parámetro de escala. La función de densidad viene dada por la siguiente expresión: f(x;k,θ) = xk−1e−x/θ θkΓ(k)(4.9) donde Γ es la función gamma de Euler. Para este trabajo, se va a emplear los valores de k=8.5 y θ=1 , con el objetivo de reproducir los mismos resultados que en el artículo [ 2 ]. La función de densidad de la masa inicial se representa en la figura 4.9. 24 Capítulo 4. Incertidumbre en m0 Figura 4.9 Función de densidad para distribución gamma m0. 4.2.2 Cálculo de los coeficientes del GPC En este apartado, se calcularan los coeficientes del método GPC para el caso de la distribución gamma. La principal diferencia es que en este caso, se usarán los polinomios de Laguerre en vez de los polinomios de Legendre. Siguiendo el mismo procedimiento que para el caso uniforme, se llega al mismo sistema de P+1 ecuaciones diferenciales que se obtuvo para el caso uniforme tal y como se muestra en la expresión (4.10.). ˙ hl(t) = −Aδ0l−B P ∑ i=0 P ∑ j=0 hi(t)hj(t)Ci jl,l=0,...,P(4.10) Sin embargo, el coeficiente Ci jl es diferente al del caso uniforme. La expresión para este caso, se muestra en la ecuación (4.11.). Ci jl =E[φk−1 iφk−1 jφk−1 l] E[(φk−1 l)2](4.11) Donde φk−1 i denota el polinomio asociado de Laguerre de grado i . Para la distribución gamma, la masa inicial puede ser reescrita en función de los polinomios de Laguerre como: m0=¯m0φk−1 0(G)−δm √3kφk−1 1(G) . Por lo que se puede extraer que las condiciones iniciales para resolver el sistema de ecuaciones diferenciales son las mostradas en la expresión (4.12.). h0(0) = ¯m0,h1(0) = −δm √3khl(0) = 0,para l =2,...,P(4.12) Resolviendo para P=3 y usando los datos definidos en el capítulo 3, se obtienen los valores de los coeficientes mostrados en la figura 4.10. Como se puede observar, se sigue cumpliendo que cada coeficiente se encuentra dos ordenes de magnitud por debajo del anterior. 4.2 Distribución gamma 25 Figura 4.10 Coeficientes para distribución gamma con incertidumbre en m0para P=3. 4.2.3 Esperanza Una vez obtenido los coeficientes hi , el cálculo de la esperanza es inmediato, pues es la misma expresión que (4.7) y por tanto, la evolución de la esperanza con el tiempo es el coeficiente h0 . En la figura 4.11 se muestra la evolución de la esperanza a lo largo del tiempo, como se observa, mantiene la tendencia descendiente debido al consumo de combustible durante el vuelo. En la tabla 4.5. se muestra los valores de la esperanza para diferentes instantes de tiempos, además se compara estos resultados con los obtenidos en el artículo de Rafael Vázquez Valenzela y Damian Rivas Rivas [ 2 ] usando este mismo método y los resultados obtenidos en el Trabajo Fin de Grado de Manuel Ángel Zapata Habas [ 3 ]. Por tanto, la primera columna representa los instante de tiempo en los que se muestra la esperanza. La segunda columna son los resultados obtenidos en este trabajo usando GPC. La tercera, los valores obtenidos en el artículo usando también el método de GPC y la cuarta columna representa los valores obtenidos empleando el método de Montecarlo. Tabla 4.5 Valores de la esperanza para distintos instantes de tiempos, Distribución gamma m0. tiempo (s) E[m(t;m0)] (Kg).GPC E[m(t;m0)] (Kg).GPC* E[m(t;m0), (Kg) MC 2000 77485.5 77485 77485.8 4000 73477.1 73477 73477.8 6000 69595.9 69596 69596.1 8000 65831.7 65831 65831.6 10000 62174.8 62175 62174.3 12000 58616.5 58616 58616.2 Como se puede observar se obtienen los mismos resultados que en el trabajo Fin de Grado y en el artículo, con errores inferiores a un kilogramo. A continuación, se va a comprobar para este tipo de distribución aleatoria, el efecto que tiene pasar de P=3 a P=5. Estos resultados, se muestran en la tabla 4.6. Siendo la primera columna los instantes de tiempo en los que se toman valores, la segunda el valor de la esperanza para P=3 y la tercera columna, los valores de la esperanza para P=5. Se han marcado en rojo los valores que 26 Capítulo 4. Incertidumbre en m0 Figura 4.11 Valor esperado de la masa lo largo del tiempo. Distribución gamma m0. son diferentes para que sea mas fácil localizar la diferencia. Tabla 4.6 Valores de la esperanza para P=3 y P=5 , Distribución gamma m0. tiempo (s) E[m(t;m0)](Kg) P=3E[m(t;m0)] (Kg)P=5 2000 77485.59985630068 77485.59985630068 4000 73477.0746238157673477.07462381573 6000 69595.9313047811 69595.9313047811 8000 65831.7193400454365831.71934004541 10000 62174.85069152920 62174.85069152919 12000 58616.526906453758616.5269064536 Como se puede observar, los cambios vuelven a tener una precisión de 15 cifras significativas, por lo que se puede intuir que con P=3 es una buena aproximación para los casos en los que se tiene una sola variable aleatoria. 4.2.4 Desviación típica La expresión que determina la variación de la varianza a lo largo del tiempo en este caso, es análoga a la expresión (4.8), basta con cambiar los polinomios de Legendre por los polinomios generalizados de Laguerre, tal y como se muestra en la expresión (4.13). 4.2 Distribución gamma 27 Var[m(t;m0)] = E[m2(t;m0)] −E[m(t;m0)]2 = P ∑ i=0 P ∑ j=0 hi(t)hj(t)E[φk−1 i(G)φk−1 j(G)]−h2 o = P ∑ i=0 h2 i(t)E[φk−1 i(G)2] (4.13) Para obtener la desviación típica, basta hacer la raíz cuadrada a la expresión anterior. En la figura 4.12. se puede ver la evolución de la desviación típica a lo largo del tiempo. Figura 4.12 Desviación típica de la masa lo largo del tiempo. Distribución gamma m0. En la tabla 4.7. se muestran los valores de la desviación típica para diferentes valores de tiempo, siendo la primera columna dichos instantes de tiempo. La segunda columna representa los valores de la desviación típica obtenidos mediante GPC. La tercera columna son los valores obtenidos en el artículo de Rafael Vázquez Valenzuela y Damian Rivas Rivas. La cuarta columna son los valores obtenidos por Manuel Ángel Zapata Habas empleando el método de Montecarlo. Como se pueden ver, los valores obtenidos por GPC son bastante precisos, existiendo un error menor que un Kilogramo con respecto al método de Montecarlo. Finalmente, en la figura 4.13. se representa el cociente entre la desviación típica y la esperanza del combustible consumido a lo largo del tiempo. La tendencia es decreciente y muy similar al caso uniforme. 28 Capítulo 4. Incertidumbre en m0 Tabla 4.7 Valores de la desviación típica para distintos instantes de tiempos, Distribución gamma m0. tiempo (s) σ(t;m0)] (Kg).GPC σ[m(t;m0)] (Kg).GPC* σ[m(t;m0)], (Kg) MC 2000 2786.5 2786 2786.2 4000 2694.6 2695 2694.8 6000 2610.2 2610 2610.2 8000 2532.8 2533 2532.5 10000 2461.8 2462 2461.9 12000 2396.5 2397 2396.5 Figura 4.13 Cociente entre la desviación típica y el combustible consumido. Distribución gamma m0. 5 Incertidumbre en CD0 En este capítulo, se estudia el caso de tener incertidumbre en el parámetro CD0 . Esta incertidumbre puede venir por un fallo en la medición de este coeficiente o por el deterioro de las superficies de la aeronave. Esta incertidumbre se modelará como una distribución uniforme centrada en el valor nominal y un semiancho del 10% sobre el valor inicial. Esta función de distribución es la mostrada en la figura 5.1. Figura 5.1 Función de densidad para distribución uniforme CD0. A diferencia del caso anterior en el que m0 solo aparecía en la condición inicial y por tanto se iba propagando la incertidumbre a lo largo del tiempo, en este caso, la incertidumbre afecta en todo instante de tiempo, por lo que se espera un efecto totalmente diferentes en los resultados comparándolos con el capítulo anterior. Por otro lado, la masa inicial m0 en este caso es fija e igual al valor nominal, es decir m0=81633Kg . En los siguientes apartados se verá como afecta el caso de tener incertidumbre en CD0 al sistema de ecuaciones que hay que resolver para obtener los coeficientes, se calcularán dichos coeficientes y por último, se 29 36 Capítulo 6. Incertidumbre en c 6.1 Cálculo de los coeficientes del GPC Como ya se ha mencionado, en este caso, c se encuentra tanto en el coeficiente A y en el B, por lo que las expresiones de estas constantes que se va a emplear en este capítulo son las mostradas en la ecuación (6.1) y (6.2). A=1 2ρV2SCD0(6.1) B=2CD2g2 ρV2S(6.2) Por otro lado, se escribe cen función de los polinomios de Legendre: c=¯cL0(∆)+δcL1(∆)(6.3) La masa, en función del desarrollo en serie del GPC queda: m(t;CD0) = P ∑ i=0 hi(t)Li(∆)(6.4) Sustituyendo (6.1-6.4) en la ecuación diferencial de la masa (3.8), se obtiene la siguiente expresión: P ∑ i=0 ˙ hi(t)Li(∆) = −A(¯cL0(∆)+δcL1(∆))−B P ∑ i=0 P ∑ j=0 hi(t)hj(t)Li(∆)Lj(∆)(¯cL0(∆)+δcL1(∆)) (6.5) Ahora, se opera de forma similar a como se operó en los capítulos anteriores, es decir, se multiplica por Ll(∆) , se toma esperanza con respecto a Ll(∆) , se aplica lo propiedad de los polinomios ortogonales y por último, se reordena la expresión. Quedando un sistema de P+1 ecuaciones diferenciales. ˙ hl(t) = −A(cδ0l+δcδ1l)−B P ∑ i=0 P ∑ j=0 hi(t)hj(t)(¯cX0 i jl +δcX1 i jl),l=0,...,P(6.6) Siendo los coeficientes Xk i jl constante y su expresión se muestra en la ecuación (6.7) Xk i jl =E[LiLjLlLk] E[L2 i]k=0,1(6.7) Las condiciones iniciales para este caso, sigue siendo las mismas que en el apartado anterior, se muestran en la ecuación (6.8) h0(0) = m0,hl(0) = 0,l=1,...,P(6.8) Por lo que ya se puede resolver el sistema de ecuaciones diferenciales. Se ha vuelto a seleccionar el valor de P=3 para este caso. La evolución de estos coeficientes se muestran en la figura 6.2. Como se observa, los coeficientes siguen la tendencia de disminuir 2 ordenes de magnitud con respecto al orden de magnitud anterior. 6.2 Esperanza 37 Figura 6.2 Coeficientes para distribución uniforme con incertidumbre en c. 6.2 Esperanza Para el cálculo de la esperanza, se hace uso de la ecuación 6.9, que como se puede ver, es la misma expresión que la empleada para los capítulos anteriores. E[m(t;c)] = P ∑ i=0 hi(t)E[Li(∆)] = P ∑ i=0 hi(t)E[Li(∆)L0(∆)] = h0(t)E[L2 0(∆)] = h0(t)(6.9) En la figura 6.3 se representa gráficamente la evolución de la esperanza con respecto al tiempo, como se puede observar sigue una tendencia descendente debido al consumo de combustible durante el vuelo y muy similar a los casos anteriores. En la tabla 6.1 se recogen el valor de la esperanza para diferentes valores de tiempo. La primera columna, representa los instantes de tiempo en los que se calcula la esperanza y la segunda columna, representa el valor que toma la esperanza en dicho instantes de tiempo. Tabla 6.1 Valores de la esperanza para distintos instantes de tiempos, Distribución uniforme c. tiempo (s) E[m(t;c)] (Kg) 2000 77487.5 4000 73481.2 6000 69602.4 8000 65840.5 10000 62186.1 12000 58630.3 38 Capítulo 6. Incertidumbre en c Figura 6.3 Valor esperado de la masa lo largo del tiempo. Distribución uniforme c. En la tabla 6.2 se compara los valores obtenido de la esperanza en diferentes instantes de tiempo para los tres casos estudiados en este trabajo, considerando incertidumbres en m0 , CD0 y c . Por tanto, la primera columna representa los instantes de tiempo. La segunda columna muestra el valor esperado de la masa considerando incertidumbres en m0. La tercera considerando incertidumbre en CD0y la cuarta, en c. Tabla 6.2 Comparación de los diferentes valores esperados de la masa para diferentes instantes de tiempo, considerando los casos de incertidumbres en m0,CD0yc. tiempo (s) E[m(t;m0)] (Kg) E[m(t;CD0)](Kg) E[m(t;c) (Kg) 2000 77485.6 77487.3 77487.5 4000 73477.1 73480.3 73481.2 6000 69595.9 69600.6 69602.4 8000 65831.7 65837.6 65840.5 10000 62174.8 62181.8 62186.1 12000 58616.5 58624.4 58630.3 Como se puede apreciar, considerando incertidumbre únicamente en m0 da lugar a valores esperados mas bajos que considerando las otras dos incertidumbres (por separado). Por otro lado, la incertidumbre en c es la que da lugar a valores esperados de la masa mas altos. 6.3 Desviación típica 39 6.3 Desviación típica La evolución de la desviación típica viene dada por la raiz cuadrada de la siguiente expresión: Var[m(t;c] = E[m2(t;c)] −E[m(t;c)]2= P ∑ i=0 P ∑ j=0 hi(t)hj(t)E[Li(∆)Lj(∆)]−h2 o= P ∑ i=0 h2 i(t)E[L2 i(∆)] (6.10) Como se puede observar, es la misma expresión que en los casos anteriores. En la figura 6.4, se muestra la evolución de la desviación típica a lo largo del tiempo. Como se puede observar, tiende a aumentar con el tiempo. Figura 6.4 Desviación típica de la masa lo largo del tiempo. Distribución uniforme c. En la tabla 6.3 se muestra los valores de la desviación típica para ciertos instantes de tiempo, siendo la primera columna dichos instantes de tiempo y la segunda columna, los valores obtenidos mediante el método de GPC. En la figura 6.5 se compara las desviaciones típicas para los diferentes caso estudiados hasta ahora. Como se observa, el caso con una mayor desviación típica es el caso en el que se tiene incertidumbre en m0 , mientras que en el que se tiene menor desviación típica es en el caso de CD0. En la figura 6.6 se representa el cociente entre la desviación típica y la masa consumida, como se puede observar, sigue una tendencia decreciente, por lo que la incertidumbre tiende a disminuir a medida que se consume combustible. Finalmente, en la figura 6.7 se compara con el cociente obtenido en los otros apartados. Como se puede observar, m0 es el caso en el que la incertidumbre es mayor, tomando un valor alto al principio y posteriormente siguiendo una tendencia decreciente a lo largo del tiempo. Esto es debido a que la incertidumbre se tiene en 40 Capítulo 6. Incertidumbre en c Figura 6.5 Comparación de la desviación típica para los casos de incertidumbre en m0,CD0yc. Figura 6.6 Cociente entre la desviación típica y el combustible consumido. Distribución uniforme c. 6.3 Desviación típica 41 Tabla 6.3 Valores de la desviación típica para distintos instantes de tiempos, Distribución uniforme CD0. tiempo (s) σ[m(t;c)] (Kg).GPC 2000 235.2 4000 455.1 6000 661.6 8000 856.3 10000 1040.8 12000 1216.3 la condición inicial. En los otros dos casos, la incertidumbre en c es mayor que en CD0 . Por lo que se puede concluir que el caso en el que menos incertidumbre se tiene es en CD0 ya que tanto la desviación típica como el cociente entre la desviación típica y la masa consumida es menor al resto. Figura 6.7 Comparación entre el cociente de la desviación típica y el combustible consumido. Para los casos de incertidumbre en m0,CD0yc. 7 Incertidumbre en m0yCD0 En este capítulo, se pretende estudiar el caso en el que existan incertidumbres tanto en la masa inicial como en CD0 . Ambas variables se van a modelar como una distribución uniforme centradas en el valor nominal y semiancho del 10% del valor nominal. La función de densidad de estas funciones se pueden observar en las figuras 4.1 y 5.1. En los siguientes capítulos, se calcularán los coeficientes necesarios para aplicar el método. Posteriormente, se obtendrá la evolución de la esperanza y la desviación típica de la masa a lo largo del tiempo. Por último, se comparará la eficiencia de este método con el de Montecarlo. 7.1 Cálculo de los coeficientes del GPC En este capítulo se aborda el paso de tener una sola variable aleatoria, a tener dos. En primer lugar, el desarrollo en serie de la masa, ya no depende de una sola variable hi , sino que depende también de una variable hj. Por tanto, el desarrollo en serie queda de la siguiente forma: m(t;m0,CD0) = P ∑ i=0 P ∑ j=0 hi j(t)Li(∆1)Lj(∆2)(7.1) Donde ∆1 hace referencia a la incertidumbre en m0 y ∆2 hace referencia a la incertidumbre en CD0 . Además, la masa inicial y el coeficiente de resistencia sin sustentación CD0pueden ser reescritas como: m0=¯m0L0(∆1)+δmL1(∆1)(7.2) CD0=¯ CD0L0(∆2)+δCD0L1(∆2)(7.3) Sustituyendo (7.1), (7.2) y (7.3) en (3.8) y usando la expresión de A descrita en (5.2), se obtiene el sistema de (P+1)2ecuaciones diferenciales, tal y como se muestra en (7.4). P ∑ i=0 P ∑ j=0 ˙ hi j(t)Li(∆1)Lj(∆2) = −A(¯ CD0L0(∆2)+δCD0L1(∆2)) −B P ∑ i1=0 P ∑ i2=0 P ∑ j1=0 P ∑ j2=0 hi1j1(t)hi2j2(t)Li1(∆1)Li2(∆1)Lj1(∆2)Lj2(∆2) (7.4) Donde los subindices i , i1 e i2 hacen referencia a la variable m0 y los subindices j , j1 y j2 hacen referencia a la variable CD0 . La forma de proceder para simplificar el sistema de ecuaciones es tratar cada variable por separado. En primer lugar, se multiplica por Ll1(∆1) , se toma esperanzas con respecto a ∆1 y se aplica la propiedad de los polinomios ortogonales. Los polinomios ortogonales de Legendre referidos a la variable CD0 , son independientes a los polinomios asociados a la variable m0 , por lo que pueden salir fuera de la esperanza. 43 44 Capítulo 7. Incertidumbre en m0yCD0 P ∑ j=0 ˙ hl1j(t)E[Ll1(∆1)2]Lj(∆2) = −Aδ0l1(¯ CD0L0(∆2)+δCD0L1(∆2)) −B P ∑ i1=0 P ∑ i2=0 P ∑ j1=0 P ∑ j2=0 hi1j1(t)hi2j2(t) xE[Li1(∆1)Li2(∆1)Ll1(∆1)]Lj1(∆2)Lj2(∆2) (7.5) Reordenando la ecuación anterior queda: P ∑ j=0 ˙ hl1j(t)Lj(∆2) = −Aδ0l(¯ CD0L0(∆2)+δCD0L1(∆2)) −B P ∑ i1=0 P ∑ i2=0 P ∑ j1=0 P ∑ j2=0 hi1j1(t)hi2j2(t)Ci1i2l1Lj1(∆2)Lj2(∆2) (7.6) Finalmente, se realiza los mismos pasos para la variable CD0 , multiplicando en este caso por Lm(∆2) , dando lugar a un sistema de (P+1)2ecuaciones diferenciales: ˙ hlm(t)Lj(∆2) = −Aδ0l(¯ CD0δ0m+δCD0δ1m)−B P ∑ i1=0 P ∑ i2=0 P ∑ j1=0 P ∑ j2=0 hi1j1(t)hi2j2(t)Ci1i2l1Cj1j2m(7.7) Para poder obtener los coeficientes hlm queda establecer las condiciones iniciales, estas son las siguientes: h0m(0) = ¯m0,h1m(0) = δmhlm(0) = 0,para l =2,...,P;m=0,...,P(7.8) Para este caso, se vuelve a usar el valor de P=3, dando lugar a 16 coeficientes aleatorios. Estos coeficientes, se muestran en la siguiente expresión:     h00 h10 h20 h30 h01 h11 h21 h31 h02 h12 h22 h32 h03 h13 h23 h33     (7.9) En esta matriz, los coeficientes diagonales tienen el mismo orden de magnitud y van disminuyendo diagonalmente el orden de magnitud en dos unidades. Es decir, los coeficientes h01 y h01 tienen dos ordenes de magnitud menos que el coeficiente h00 . Los coeficientes h02 , h11 y h20 sus ordenes de magnitud son dos unidades menor que h01 y h01 . En la figura 7.1 se representa estos coeficientes agrupados por sus ordenes de magnitud. En la figura 7.2 se muestra el tiempo necesario para resolver el sistema de ecuaciones en función de P, como se observa, sigue una evolución cuadrática. Por lo cual con P=3 se tiene una muy buena aproximación y el tiempo que tarda en resolver el sistema de ecuaciones diferenciales, unos 0.6091 segundos, no es excesivo. En cuanto al espacio que ocupa en memoria, con P=3, la variable que almacena los coeficientes aleatorios ocupa un espacio de 1920128 bytes. El espacio que ocupa esta variable en memoria también sigue una ley cuadrática, por lo que al aumentar P, aumenta cuadráticamente el espacio que ocupa dicha variable en memoria. 7.2 Esperanza 45 Figura 7.1 Coeficientes para distribución uniforme con incertidumbre en m0yCD0. Figura 7.2 Tiempo necesario para resolver el sistema de ecuaciones diferenciales en función de P, distribución uniforme m0yCD0. 7.2 Esperanza Como en capítulos anteriores, la esperanza a lo largo del tiempo se expresa como: E[m(t;m0,CD0)] = P ∑ i=0 P ∑ j=0 hi j(t)E[Li(∆1)Lj(∆2)] = P ∑ i=0 P ∑ j=0 hli(t)E[Li(∆1)L0(∆1)]E[Lj(∆2)L0(∆2)] = =h00(t)E[L2 0(∆1)]E[L2 0(∆2)] =h00(t) (7.10) 52 Capítulo 8. Incertidumbre en m0,CD0yc se obtiene el sistema de (P+1)3ecuaciones diferenciales, tal y como se muestra en (8.7). P ∑ i=0 P ∑ j=0 P ∑ k=0 ˙ hi jk(t)Li(∆1)Lj(∆2)Lk(∆3) = −A(¯ CD0L0(∆2)+δCD0L1(∆2))(¯cL0(∆3)+δcL1(∆3)) −B P ∑ i1,i2=0 P ∑ j1,j2=0 P ∑ k1,k2=0 (¯cL0(∆3)+δcL1(∆3))hi1j1k1(t)hi2j2k2(t) xLi1(∆1)Li2(∆1)Lj1(∆2)Lj2(∆2)Lk1(∆3)Lk2(∆3) (8.7) Donde los subindices i , i1 e i2 hacen referencia a la variable m0 , los subindices j , j1 y j2 hacen referencia a la variable CD0 y los subindices k , k1 y k2 hacen referencia a la variable c . La forma de proceder para simplificar el sistema de ecuaciones es tratar cada variable por separado. A continuación se realizan los mismos pasos que en el capítulo anterior, se multiplica por Ll(∆1) y se toma esperanza respecto a ∆1 , se reordena términos y se repite la misma operación con Lm(∆2) y Ln(∆3) . Dando lugar a un sistema de (P+1)3 ecuaciones diferenciales. ˙ hlmn(t) = −Aδ0l(¯ CD0δ0m+δCD0δ1m)(¯cδ0n+δcδ1n) −B P ∑ i1=0 P ∑ i2=0 P ∑ j1=0 P ∑ j2=0 P ∑ k1=0 P ∑ k2=0 hi1j1k1(t)hi2j2k2(t)Ci1i2l1Cj1j2m(¯cX0 k1k2n+δcX1 k1k2n)(8.8) Para poder obtener los coeficientes hlmn queda establecer las condiciones iniciales, estas son las siguientes: h0mn(0) = ¯m0,h1mn(0) = δmhlmn(0) = 0,para l =2,...,P;m=0,...,P;n=0,...,P(8.9) Para este caso, se vuelve a usar el valor de P=3, dando lugar a 64 coeficientes aleatorios. Estos coeficientes siguen la misma tendencia de ir disminuyendo de magnitud. En concreto, a medida que aumenta la suma de los subindices de los coeficientes en una unidad con respecto al anterior, el orden de magnitud cae en dos unidades. En la figura 8.1 se representa los cuatro grupos de coeficientes con mayor orden de magnitud, que son los mas representativos. Es decir, la primera gráfica representan aquellos coeficientes que sus subindices suman 0, en la segunda aquellos cuyos subindices suman 1, en la tercera los que suman 2 y en la cuarta, los que suman 3. En la figura 8.2 se muestra el tiempo necesario para resolver el sistema de ecuaciones en función de P, como se observa, sigue una evolución cúbica. Por lo cual con P=3 se tiene una muy buena aproximación y el tiempo que tarda en resolver el sistema de ecuaciones diferenciales, unos 53.38 segundos, que ya empieza a ser un tiempo considerable. En cuanto al espacio que ocupa en memoria, con P=3, la variable que almacena los coeficientes aleatorios ocupa un espacio de 768512 bytes. Como era de esperar, a mayor número de coeficientes, mas espacio en memoria se necesita. 8.2 Esperanza Como en capítulos anteriores, la esperanza a lo largo del tiempo vuelve a ser el coeficiente h000 , para llegar a este resultado, basta con extender la ecuación (7.11) a tres variables. En la figura 8.3 se muestra la evolución de la esperanza a lo largo del tiempo. Como se observa, esta gráfica es muy similar a los casos en los que solo se tenía una variable aleatoria. Por lo que el efecto de las incertidumbres en m0,CD0yces muy poco significativo sobre la esperanza. En la tabla 8.1 se muestra el valor de la esperanza en diferentes instantes de tiempo. Por tanto, la primera columna representa los instantes de tiempo. La segunda columna el valor de la esperanza para dichos instantes de tiempo usando el método de GPC. Comparando estos resultados con los obtenidos en los capítulos anteriores, se puede observar como el efecto 8.2 Esperanza 53 Figura 8.1 Coeficientes para distribución uniforme con incertidumbre en m0,CD0yc. Figura 8.2 Tiempo necesario para resolver el sistema de ecuaciones diferenciales en función de P, distribución uniforme m0,CD0yc. 54 Capítulo 8. Incertidumbre en m0,CD0yc Figura 8.3 Valor esperado de la masa lo largo del tiempo. Distribución uniforme m0,CD0yc. Tabla 8.1 Valores de la esperanza para distintos instantes de tiempos, Distribución uniforme m0,CD0yc. tiempo (s) E[m(t;m0)] (Kg).GPC 2000 77485.8 4000 73477.9 6000 69597.6 8000 65834.4 10000 62178.7 12000 58621.6 de considerar estas incertidumbres no tiene demasiado impacto sobre la evolución del valor esperado de la masa, pues apenas hay 10kg de diferencia entre los diferentes casos. 8.3 Desviación típica Usando la definición de varianza y expresando la masa en función de los coeficientes hi jk se tiene: Var[m(t;m0,CD0,c)] = E[m2(t;m0,CD0,c)]−E[m(t;m0,CD0,c)] = P ∑ i=0 P ∑ i0=0 P ∑ j=0 P ∑ j0=0 P ∑ k=0 P ∑ k0=0 hi jk(t)hi0j0k0(t) xE[Li(∆1)Li0(∆1)Lj(∆2)Lj0(∆2)Lk(∆3)Lk0(∆3)]−h2 00(t) (8.10) Aplicando la propiedad de los polinomios ortogonales se pueden simplificar los sumatorios, quedando la 8.3 Desviación típica 55 expresión final de la varianza como se muestra en la ecuación (7.10). Var[m(t;m0,CD0,c)] = P ∑ i=0 P ∑ j=0 P ∑ k=0 h2 i jk(t)E[Li(∆1)2]E[Lj(∆2)2]E[Lk(∆3)2]−h2 00(t)(8.11) Para obtener la desviación típica basta con hacer la raíz cuadrada de le expresión (8.11). En la figura 8.4 se muestra la evolución de la desviación típica a lo largo del tiempo.Como se observa, sigue una tendencia decreciente, hasta un cierto instante de tiempo en el cual empieza a crecer debido a que las incertidumbres en CD0yctiende a crecer con el tiempo. Figura 8.4 Desviación típica de la masa lo largo del tiempo. Distribución uniforme m0,CD0yc. En la tabla 8.2 se muestra los valores de la desviación típica para ciertos instantes de tiempo, siendo la primera columna dichos instantes de tiempo. La segunda columna, los valores obtenidos mediante el método de GPC. Dado que se tiene la evolución de la masa a lo largo del tiempo, es posible calcular de forma teórica la desviación típica, esta solución teórica se muestra en la tercera columna. Tabla 8.2 Valores de la desviación típica para distintos instantes de tiempos, Distribución uniforme m0yCD0. tiempo (s) σ[m(t;m0,CD0,c)] (Kg).GPC σ[m(t;m0,CD0,c), (Kg) 2000 2802.0 2802.0 4000 2752.4 2752.4 6000 2734.3 2734.3 8000 2744.1 2744.1 10000 2778.0 2778.0 12000 2832.5 2832.5 56 Capítulo 8. Incertidumbre en m0,CD0yc Un aspecto a tener en cuenta es que con P=3, se tiene un error absoluto en la esperanza de 1x10−8 , este mismo resultado se repite para P=4. Por lo que P=3 garantiza que el método converge rápido y garantiza la misma precisión que P=4 para este caso. En la figura 8.5 se compara la desviación típica en este caso, con la desviación típica obtenida en el caso de tener incertidumbre en m0 y con la desviación típica de tener incertidumbre en m0 y CD0 . Como se observa, el efecto de la incertidumbre en c da lugar a que en este caso, la desviación típica sea siempre superior al resto de casos. Además, el efecto en la desviación típica de c es mayor que el de CD0 , ya que modifica mas la gráfica. Figura 8.5 Comparación de la desviación típica de la masa lo largo del tiempo. Distribución uniforme m0 , CD0yc. Por otro lado, en la figura 8.6 se representa el cociente entre la desviación típica y el consumo de combustible, como se puede observar decrece con el tiempo, por lo que la incertidumbre tiende a reducirse con el tiempo. 8.3 Desviación típica 57 Figura 8.6 Cociente entre la desviación típica de la masa y el consumo de combustible a lo largo del tiempo. Distribución uniforme m0,CD0yc. 9 Incertidumbre en m0,CD0,cyCD2 En este último capítulo, se pretende estudiar el caso en el que existan incertidumbres tanto en la masa inicial m0 , en el consumo específico c , en el coeficiente de resistencia parásita CD0 y en el coeficiente de resistencia inducida CD2 . Ambas variables se van a modelar como una distribución uniforme centradas en el valor nominal y semiancho del 10% del valor nominal. La función de densidad de estas funciones se pueden observar en las figuras 4.1, 5.1 y 6.1 para el caso de m0 , CD0 y c . Mientras que para el caso de CD2 , al ser la primera vez que se considera esta variable, su función de densidad se muestra en la figura 9.1. Figura 9.1 Función de densidad para CD2, distribución uniforme. En los siguientes apartados, se calcularán los coeficientes necesarios para aplicar el método. Posteriormente, se obtendrá la evolución de la esperanza y la desviación típica de la masa a lo largo del tiempo. Por último, se comparará la eficiencia del método GPC frente al de Montecarlo. 59 60 Capítulo 9. Incertidumbre en m0,CD0,cyCD2 9.1 Cálculo de los coeficientes del GPC En este capítulo se deducirá la expresión que permite calcular los coeficientes del método. En primer lugar, el desarrollo en serie de la masa queda como se expresa en la expresión (9.1). m(t;m0,CD0,c,CD2) = P ∑ i=0 P ∑ j=0 ∑ k=0 P ∑ q=0 hi jkq(t)Li(∆1)Lj(∆2)Lk(∆3)Lq(∆4)(9.1) Donde ∆1 hace referencia a la incertidumbre en m0 , ∆2 hace referencia a la incertidumbre en CD0 , ∆3 hace referencia a la incertidumbre en c y ∆4 hace referencia a la incertidumbre en CD2 . Además, la masa inicial, el coeficiente de resistencia aerodinámico CD0 parásito, el consumo específico c y el coeficiente de resistencia inducido CD2se expresan como: m0=¯m0L0(∆1)+δmL1(∆1)(9.2) CD0=¯ CD0L0(∆2)+δCD0L1(∆2)(9.3) c=¯cL0(∆3)+δcL1(∆3)(9.4) CD0=¯ CD2L0(∆4)+δCD2L1(∆4)(9.5) En ese caso, al ser c y CD0 y CD2 variables aleatorias, es necesario volver a definir las constantes A y B para que vuelvan a ser constantes, esta nueva expresión de A y B se muestran en (9.6) y (9.7). A=1 2ρV2S(9.6) B=2g2 ρV2S(9.7) Sustituyendo (9.1), (9.2), (9.3) y (9.4) en (3.8) y usando las expresiones de A y B descritas en (9.7) y (9.8), se obtiene el sistema de (P+1)4ecuaciones diferenciales, tal y como se muestra en (9.8). P ∑ i=0 P ∑ j=0 P ∑ k=0 P ∑ q=0 ˙ hi jkq(t)Li(∆1)Lj(∆2)Lk(∆3)Lq(∆4) = −A(¯ CD0L0(∆2)+δCD0L1(∆2))(¯cL0(∆3)+δcL1(∆3)) −B P ∑ i1,i2=0 P ∑ j1,j2=0 P ∑ k1,k2=0 P ∑ q1,q2=0 (¯cL0(∆3)+δcL1(∆3)) ×(¯ CD2L0(∆4)+δCD2L1(∆4))hi1j1k1q1(t)hi2j2k2q2(t) ×Li1(∆1)Li2(∆1)Lj1(∆2)Lj2(∆2)Lk1(∆3)Lk2(∆3)Lq1(∆4)Lq2(∆4) (9.8) Donde los subindices i , i1 e i2 hacen referencia a la variable m0 , los subindices j , j1 y j2 hacen referencia a la variable CD0 , los subindices k , k1 y k2 hacen referencia a la variable c y los subindices q , q1 y q2 hacen referencia a CD2 . La forma de proceder para simplificar el sistema de ecuaciones es tratar cada variable por separado, como en los capítulos anteriores, multiplica por Ll(∆1) y se toma esperanza respecto a ∆1 , se reordena términos y se repite la misma operación con Lm(∆2) , Ln(∆3) y Lo(∆4) . Dando lugar a un sistema de (P+1)4ecuaciones diferenciales. 9.1 Cálculo de los coeficientes del GPC 61 ˙ hlmno(t) = −Aδ0lδ0o(¯ CD0δ0m+δCD0δ1m)(¯cδ0n+δcδ1n) −B P ∑ i1=0 P ∑ i2=0 P ∑ j1=0 P ∑ j2=0 P ∑ k1=0 P ∑ k2=0 P ∑ q1=0 P ∑ q2=0 hi1j1k1(t)hi2j2k2(t)Ci1i2l1Cj1j2m(¯cX0 k1k2n+δcX1 k1k2n) ×(¯ CD2X0 q1q2o+δCD2X1 q1q2o) (9.9) Para poder obtener los coeficientes hlmn0 queda establecer las condiciones iniciales, estas son las siguientes: h0mno(0) = ¯m0,h1mno(0) = δmhlmno(0) = 0,para l =2,...,P;m=0,...,P;n=0,...,P;o=0,...,P (9.10) Para este caso, se vuelve a usar el valor de P=3, dando lugar a 254 coeficientes aleatorios. Estos coeficientes siguen la misma tendencia de ir disminuyendo de magnitud. En concreto, a medida que aumenta la suma de los subindices de los coeficientes en una unidad con respecto al anterior, el orden de magnitud cae en dos unidades. En la figura 9.2 se representa los cuatro grupos de coeficientes con mayor orden de magnitud, que son los mas representativos. Es decir, la primera gráfica representan aquellos coeficientes que sus subindices suman 0, en la segunda aquellos cuyos subindices suman 1, en la tercera los que suman 2 y en la cuarta, los que suman 3. Figura 9.2 Coeficientes para distribución uniforme con incertidumbre en m0,CD0cyCD2. En la figura 9.3 se muestra el tiempo necesario para resolver el sistema de ecuaciones en función de P, para el caso en que P vale 3, se necesita alrededor de unos 100 minutos. Este tiempo ya es considerablemente alto, sin embargo, proporciona una precisión bastante adecuada. En cuanto al espacio que ocupa en memoria, con P=3, la variable que almacena los coeficientes aleatorios 68 Capítulo 10. Conclusiones En cuanto a posibles trabajos o ampliaciones varias lineas de trabajo pueden ser: • En primer lugar, se puede estudiar el efecto de otras variables sobre la masa. Es decir,extender el métodos a otras variables que no se han tenido en cuenta en este trabajo, como puede ser la velocidad o la densidad del aire. • En segundo lugar, se puede aplicar este mismo método a otros segmentos de vuelo, como puede ser descenso, ascenso o viraje. Estudiando el efecto sobre la masa de las mismas variables consideradas en este trabajo o otras que resulten de interés. • En tercer lugar, estudiar el efecto del viento, ya que el viento este afecta durante toda la fase de crucero. Al ser esta la fase de vuelo mas larga, puede dar lugar a diferencias importantes en el consumo de combustible. Bibliografía [1] Rafael Vazquez. http://www.complexworld.eu/the-ubiquity-of-uncertainty-in-atm/ , the ubiquity of uncertainty in atm. [2] R. Vazquez and D. Rivas, “Propagation of initial mass uncertainty in aircraft cruise flight,” Journal of Guidance, Control, and Dynamics, vol. 36, no. 2, pp. 415–429, 2013. [3] M. Á. Zapata Habas, “Análisis de trayectorias de crucero de aviones comerciales sujetas a incertidumbre en los datos,” 2015. [4] D. Xiu and G. E. Karniadakis, “The wiener–askey polynomial chaos for stochastic differential equations,” SIAM journal on scientific computing, vol. 24, no. 2, pp. 619–644, 2002. [5] J. E. C. Mendoza and G. C. G. Parra, “Comparación de caos polinomial y monte carlo para ecuaciones diferenciales ordinarias aleatorias,” Ciencia e Ingeniería, vol. 33, no. 1, pp. 9–20, 2012. 69