ER2.10. Análisis de sensibilidad e incertidumbre de un modelo masa-resorte

Introducción (planteamiento)

Dado el siguiente modelo lineal, implementar en MATLAB los métodos OAT, fuerza bruta y Saltelli para el análisis de sensibilidad con 4 factores (3 parámetros y 1 condición inicial) con respecto a dos salidas escalares (error cuadrático de ajuste a la salida nominal y el tiempo que tarda en llegar al valor máximo, también llamado tiempo de pico) y una salida vectorial (respuesta temporal con el modelo nominal). Probar con porcentajes de incertidumbre del 20%.

$\begin{cases} \dot{x}_1=x_2\\ \dot{x}_2=-\frac{k}{m}x_1-\frac{f}{m}x_2 + \frac{u(t)}{m}\\\end{cases}\\  x_1(0)=x_{10}, x_2(0)=0\\ k=0.15, f=0.084, m=0.54, x_{10}=0.20$

Método (plan de solución)

Pasos: (i) utilizar el GSUA Toolbox de MATLAB para obtener los índices de sensibilidad con un número diferente de muestras, utilizando como función escalar la función de coste de ajuste a la salida nominal (con los valores dados arriba); (ii) utilizar el GSUA Toolbox de MATLAB para obtener los índices de sensibilidad utilizando el tiempo de pico como función escalar (cambiar las condiciones iniciales cerca a cero y poner la entrada en 0.1 N); (iii) analizar los resultados.

Resultados (solución)

Los archivos para la ejecución de las pruebas con MATLAB de este ejercicio con el GSUA se encuentran dentro de los ejemplos del GSUA Toolbox.

Resultados con el método OAT (N = 1000):

Resultados del método de fuerza bruta:

Resultados del método de Saltelli (N = 80000):

Cálculo para una salida escalar igual al tiempo de pico de la respuesta temporal:

Cálculo para una salida escalar igual al valor máximo de la respuesta temporal:

Cálculo para una salida escalar igual al valor final (en estado estacionario) de la respuesta temporal:

Cálculo para una salida escalar igual al cuasiperíodo de la respuesta temporal:

Discusión y verificación

En relación con el gráfico de incertidumbre (sólo se muestra para el método de Saltelli, pues en ese caso se utilizó un número grande de 80000 muestras), los resultados muestran que el sistema siempre se estabiliza (aunque en diferentes valores) y presenta un comportamiento subamortiguado. Para una incertidumbre del los factores del 20%, la incertidumbre con relación a la respuesta temporal nominal supera el 20% en la respuesta transitoria, pero no la estacionaria. El valor del tiempo de pico, sobreimpulso máximo y valor final cambian.

Los gráficos de sensibilidad vectorial muestran el cambio de los índices de sensibilidad en el tiempo y sirven para determinar la influencia mayor o menos de los factores para ciertas características, como el sobreimpulso máximo (o el valor máximo) o el valor final. Los métodos de fuerza bruta y Saltelli dan una suma de los valores de los índices de sensibilidad de primer orden (S1) mayores que 1, lo cual indica que se debe tomar un número más grande muestras. Los diferentes gráficos permiten deducir lo siguiente:

• $x_1(0)$ es importante al inicio de la respuesta temporal y $f$ al final (la franja roja se hace mayor)
• En todos los puntos de inflexión $(y=0)$ los factores $\{k,m \}$ son fundamentales, dado que la elasticidad y la masa determinan la fuerza e inercia del modelo para retornar al punto de equilibrio
• En el primer pico negativo alrededor de $t=6.1~\mathrm{seg}$ la variable que más incide es $x_1(0)$ seguido por $f$, debido a que la condición inicial y la fricción son determinantes para alcanzar dicho pico
• Con el tiempo, por ejemplo en $t=30.3~\mathrm{seg}$, los picos positivos y negativos dependen cada vez menos de la condición inicial y de la fricción (los picos en morado y rojo se han cada vez más pequeños) y más de la elasticidad y la masa (franjas azul y amarilla).

Con respecto a los índices de sensibilidad escalar de primer orden (S1s), tomando como salida escalar la función de coste de ajuste de cualquier respuesta temporal a la salida nominal, con cualquier método y cualquier nivel de incertidumbre, el orden de importancia de los factores es: $\{ k,m,f,x_1(0)\}$, siendo la sensibilidad de $\{k\}$ muy similar y superior a la de los otros factores. Este es el principal resultado. Por lo tanto, es necesario medir de una manera más exacta el factor $\{k\}$. De otro lado, se observa que el modelo es aditivo (la suma de los índices es de más del 90%), por lo cual hay poca interacción entre los factores, lo cual se valida con los resultados muy similares del método OAT; es decir, en este modelo el método local OAT es suficiente. Comparando los índices de sensibilidad de primer orden con los totales se observa un mayor cambio en el índice para $m$, lo cual indica que este factor es el más decisivo en las interacciones.

La suma de los índices de sensibilidad vectorial de primer orden mayor que 1 se debe también al hecho que la varianza en cada instante del tiempo de las respuestas temporales generadas por el muestreo de Montecarlo se toma con respecto a la media y no a la salida nominal (como en el cálculo de los índices escalares). Al cambiar la varianza calculada con MATLAB con la obtenida con referencia a la salida nominal, los índices de primer orden serán menores que 1, pero la interpretación será diferente.

Todas las pruebas se realizan con un número determinado de muestras $N$, pero para ser más precisos se deberían hacer cálculos con distintos valores para observar la convergencia de los índices de sensibilidad, o, mejor aún, calcular los índices un número grande de veces y calcular los intervalos de confianza de cada índices de sensibilidad (no se hace esto aquí por ser un ejercicio meramente académico).

Comparando los métodos entre sí, se observa que en el método de fuerza bruta con un número de muestras de N = 100 se obtienen resultados aceptables, por lo que este método debería ejecutarse siempre en primer lugar para tener una idea global del comportamiento de los índices de sensibilidad. Con 1000 muestras, el tiempo de cálculo por el método de fuerza bruta tiene un incremento de 1 minuto a más de 2 horas. Con el método de Saltelli y N = 80000 se llega a un resultado similar en menos de 16 minutos, y con N = 10000 y 1 minuto y 30 segundos se obtienen valores aceptables de los índices de sensibilidad. Por lo tanto, el uso del método de Saltelli con N = 10000 parece un enfoque adecuado y no demasiado costoso computacionalmente.

El análisis de sensibilidad puede usarse para determinar la dependencia de una característica escalar de una respuesta temporal en función de los parámetros. En este ejercicio se probó con cuatro características escalares: tiempo de pico, valor máximo, valor final y cuasiperíodo. 

• En el caso del tiempo de pico los factores que más influyen son $k$ y $m$, lo cual tiene sentido, dado que la velocidad de movimiento depende de la constante elástica (a mayor valor, mayor velocidad) y de la masa (a mayor inercia mayor será el tiempo antes de que se frene la masa). Los valores negativos son comunes en el método de Saltelli cuando un índice es cercano a cero. Este resultado es coherente con el gráfico de sensibilidad de la forma analítica (cálculos detallados): 

$t_p=\frac{\pi}{\sqrt{\frac{k}{m}-\frac{f^2}{4m^2}}}$

• En el caso del valor máximo el factor determinante es $k$ y luego $f$, algo que también se observa en los gráfico de los índices vectoriales. 
• Para el valor final el factor determinante es únicamente $k$, lo cual se corrobora con el gráfico de sensibilidad vectorial. Teniendo en cuenta las características temporales de un sistema de segundo orden (sección 3.6.2 del libro), se deduce que la ganancia es inversamente proporcional a $k$ y depende sólo de $k$: a mayor $k$ menor elongación del resorte.
• El cuasiperíodo depende principalmente de $k$ y luego de $m$. La fórmula analítica es $P^*=2t_p$, similar al primer caso.

Finalmente, es importante tratar de comprender el comportamiento anterior de los índices de sensibilidad a la luz del fenómeno físico y la forma del modelo. Por ejemplo, ¿por qué es la incertidumbre de la constante elástica la de mayor peso?, ¿por qué la incertidumbre de la condición inicial es la de menor importancia?, ¿por qué hay una interacción entre algunos factores en particular? Puede ser que se obtenga una respuesta a dichas preguntas o es precisamente el análisis de sensibilidad la herramienta para responderlas, o, al menos, para abordar el proceso respectivo de análisis. En el ejemplo, el modelo es lineal y existe una relación $k/m$ entre los factores más sensibles y multiplicados por la variable de posición, que es diferente de cero en el punto de inicio. Aunque los factores tienen valores diferentes, el porcentaje de incertidumbre de cada uno es igual. La condición inicial $x_1(0)$ está más aislada de los otros factores y $f$ es un valor mucho más pequeño que la masa.