Comparativa de los métodos GUM clásico y el método de simulación Montecarlo

Desarrollo del método de Montecarlo en la respuesta a escalón de sensores de presión

Héctor Fuentes González y Miguel Antonio Marazuela Zapata, Instituto Nacional de Técnica Aeroespacial ‘Esteban Terradas’03/01/2018

A continuación, se desarrolla el cálculo de la incertidumbre de medida de un sistema de caracterización dinámica de sensores de presión ante entradas de tipo escalón, por dos técnicas diferentes: el método GUM clásico y el método de simulación de Montecarlo. Se establece una comparativa de la aplicación de cada uno de los métodos, aplicados al cálculo de la incertidumbre del tiempo de respuesta de un transmisor. Finalmente se extraen conclusiones acerca de cuál es el método idóneo para esta aplicación en concreto y se detalla la importancia de estimar adecuadamente la variabilidad de la sensibilidad de los parámetros de amortiguamiento respecto a las variables de salida y observar una adecuada distribución de probabilidad de la magnitud de salida que permita asegurar la idoneidad del método.

Este trabajo se presentó en el Congreso de Metrología de 2017

1. Introducción

En la calibración estática (régimen permanente) de sensores, los parámetros de interés son la sensibilidad y la linealidad. Sin embargo, en las aplicaciones dinámicas de estos sensores, el objetivo es determinar su comportamiento en un corto periodo de tiempo (régimen transitorio), en el que se aplica una función de entrada.

La manera en la que una magnitud varía con el tiempo constituye una información muy importante en muchos procesos de medida y control en condiciones dinámicas. En estos casos, tan importante como medir el valor exacto de la magnitud en cuestión es determinar cómo las señales de los sensores que miden esa magnitud reproducen el comportamiento en el tiempo de la misma.

El Laboratorio de Presión y Masa de INTA dispone de métodos de ensayo de caracterización dinámica de sensores de presión, que permiten la determinación de tiempos de respuesta, de subida, de establecimiento, tiempos muertos, tiempo de pico, sobreoscilación, factor de amortiguamiento, etc.; ante entradas de presión de tipo escalón. Tras la integración de centros tecnológicos que recientemente ha sufrido el Ministerio de Defensa, se está ampliando la aplicabilidad del método a otras funciones de entrada: entradas de tipo pulso y rampa, principalmente, que ya fueron aplicadas con éxito en el seno de la OTAN en la calibración de sensores piezoeléctricos. En concreto, en los años 90, España y el Laboratorio de Metrología Mecánica del ITM obtuvieron la homologación de la Alianza de dos de sus sensores como patrones internacionales.

En los siguientes apartados, tras describir muy brevemente el método, se desarrolla el cálculo de la incertidumbre de medida de uno de los parámetros de interés, el tiempo de respuesta, aplicando el método analítico GUM clásico, y posteriormente el método de simulación numérica de Montecarlo, de forma complementaria. Del análisis de los resultados obtenidos se obtienen las conclusiones mostradas en el apartado 4.

2. Desarrollo del Método de Medida de Respuesta a Escalón

2.1 El sistema físico de medida

La implementación física del sistema de ensayo está basada en la utilización de un dispositivo de rápida apertura (‘fast acting device’). Actualmente, el procedimiento se ve imitado por las presiones máximas que soportan las electroválvulas disponibles para materializar el escalón. De una manera muy esquemática, el sistema es el siguiente:

Figura 1: Equipo completo de ensayo. Imagen del montaje en el laboratorio.

El sistema originalmente diseñado ha ido implementando continuas mejoras, sobre todo en lo que a la adquisición de datos y tratamiento de resultados se refiere. También en el uso de sensores de presión piezoeléctricos con frecuencias de respuesta de hasta 3,5 kHz, para monitorizar la perturbación de entrada y sus desviaciones respecto a las funciones ideales. Para obtener resultados fiables y consistentes, y para tener una idea de la repetibilidad del método, la prueba se repetirá n veces en condiciones compatibles.

2.2. La función modelo del sistema. Magnitudes de entrada y magnitud de salida

El objetivo del ensayo es determinar el tiempo de respuesta o de establecimiento, t_s, (al 63% o más comúnmente, al 95%). Otro parámetro de habitual interés industrial es la sobreoscilación M_p y el tiempo en el que tiene lugar, t_p, o tiempo de pico. Estos parámetros deben calcularse sobre datos obtenidos experimentalmente obteniendo la curva de respuesta temporal, c(t) a la excitación de tipo escalón aplicada. Seguidamente se modeliza la curva para el tratamiento de los parámetros y el cálculo de incertidumbres, asimilando el sistema a uno de segundo orden subamortiguado, que se rige por las siguientes ecuaciones.

Aplicando la transformada de Laplace:

El comportamiento dinámico del sistema de segundo orden se puede escribir en términos de dos parámetros ζ (factor de amortiguamiento) y ω_n (frecuencia natural no amortiguada). Para sistemas subamortiguados la respuesta temporal responde a:

Su derivada respecto del tiempo, llamada función de sensibilidad o ecuación de pendientes será:

Los términos coseno se cancelan, y la ecuación puede simplificarse como:

Esta función de sensibilidad ꝺc(t)/ꝺt es fundamental para el desarrollo del cálculo analítico de incertidumbres de la respuesta temporal del sistema, según el método GUM descrito en el apartado 2.3.

Como se explica más adelante, de la dificultad de obtener estas derivadas, se justificará la elección de un desarrollo de incertidumbre basado en Montecarlo. El siguiente paso es determinar el valor de ζ y ω_n que ajusten mejor la curva modelo con la experimental. Para ello se pueden emplear métodos de ajuste por mínimos cuadrados o métodos iterativos de aproximación, utilizando como variable objetivo el valor mínimo de la máxima diferencia entre la curva observada experimentalmente y la modelizada, las variables son ζ y ω_n y las restricciones: 0,5 ≤ ζ ≤ 0,999 y 5 ≤ ω_n ≤ 6,5. Una vez calculados ζ y ω_n obtenemos la función modelo de la respuesta temporal ante una entrada de escalón de presión del instrumento bajo prueba. El ejemplo concreto mostrado en la figura siguiente está basado en un ensayo real realizado por el Laboratorio de Presión y Masa, sobre un transmisor de presión absoluta de 600 kPa de fondo de escala y salida 4-20 mA DC, ante una entrada escalón de 10% del span de presión inicial al 90% del span de presión final.

Figura 2: Representación gráfica de la señal obtenida experimentalmente y el escalón generado.

2.3. Cálculo de incertidumbres basado en el método GUM clásico

Tenemos dos magnitudes de salida, por un lado, el tiempo de respuesta, t_s, y por otro lado la sobreoscilación, M_p. Iniciamos el desarrollo a partir de la función modelo y su ecuación de sensibilidad o ecuación de pendientes:

Función modelo:

Ecuación de sensibilidad:

El tiempo de respuesta obtenido al 0,95% de la amplitud del escalón determinado es:

El tiempo de pico y la sobreoscilación son:

El balance de incertidumbres obtenido en la aplicación del método GUM clásico es el mostrado en la siguiente tabla.

Tabla 1: Tabla balance de incertidumbres del tiempo de respuesta, t_s.

No obstante, para la obtención de estos resultados, se ha considerado:

a) por un lado que se cumple el Teorema del Límite Central. Se cuenta con un buen número de distribuciones rectangulares y una normal tipo A que nos invitan a pensar que, por la Ley de Propagación de Incertidumbres, tendremos como resultado de la convolución una distribución muy parecida a una normal. Sin embargo, dada la complejidad de la función modelo, esta hipótesis podría ser incierta si una de las distribuciones uniformes tiene un carácter dominante, lo cual, a priori, se desconoce;

b) por otro lado, que los valores de ζ y ω_n que optimizan el ajuste de la respuesta temporal registrada experimentalmente con la respuesta modelizada a partir de la ecuación del sistema están carentes de un cierto valor de incertidumbre.

En este sentido, se ha realizado un análisis de sensibilidad sobre un determinado rango de valores de variables ζ, para un valor de ω_n constante, y viceversa, con un ζ constante variando ω_n en rangos acotados en la zona de subamortiguamiento. Analizando el efecto de dichas variaciones sobre el ajuste de ambas respuestas temporales, experimental y modelo, se determinan los valores de sensibilidad mostrados en la tabla siguiente.

Tabla 2: Análisis de sensibilidad de los parámetros de Tabla amortiguamiento sobre las variables de salida.

De estos surge la duda de que, al no considerar la variabilidad de estos parámetros, se está minusvalorando la incertidumbre de las variables de salida.

Resultaría muy complejo calcular las derivadas parciales de la función modelo respecto de ζ y ω_n para calcular los coeficientes de sensibilidad e incluir estos dos parámetros en el balance. Por este motivo y por determinar la distribución de probabilidad real que sigue la variable de salida, se opta por abordar el problema por el método de cálculo de incertidumbres basado en simulaciones de Montecarlo, según lo descrito en el documento JCGM 101:2008 [7], traducido al español por el Centro Español de Metrología y disponible en www.cem.es.

2.4. Cálculo de incertidumbres basado en el método de simulación de Montecarlo

Para la aplicación del método de simulación de Montecarlo se define la función modelos y las magnitudes y distribuciones de densidad de probabilidad de entrada y salida mostradas en la tabla 3.

En la propia función ya tenemos identificadas las dos variables de entrada cuya problemática sobre la consideración de su variación se planteó en 2.3, ζ y ω_n.

Las mismas estimaciones se han utilizado tanto en el desarrollo del cálculo de incertidumbre por el método GUM clásico como por el basado en las simulaciones de Montecarlo.

Tabla 3: Parametrización del desarrollo de MCM.

3. Análisis de resultados

Se obtienen los siguientes resultados para el caso concreto mostrado y una probabilidad de cobertura del 95,45%.

Se observa el carácter gaussiano de la FDP que sigue la variable de salida. Se ha superpuesto una FDP normal teórica que ilustra la similitud entre ambas. Se ha ejecutado el método varias veces para comprobar su consistencia, obteniendo resultados satisfactorios.

Figura 3: Resultados obtenidos de la implementación del MCM para la determinación del tiempo de respuesta y su incertidumbre.

Las desviaciones típicas encontradas en las pruebas de repetibilidad y reproducibilidad fueron satisfactorias, encontrándose los siguientes resultados:

Repetibilidad (n = 4): < 5 ms

Reproducibilidad: e_n< 0,29

Figura 4: Comprobación de la consistencia de los resultados obtenidos repitiendo la ejecución de MCM.

La incertidumbre obtenida para 2s (95,45%) es de 21 ms, frente a los 14 ms obtenidos en el método GUM. Esto es debido a la consideración en el cálculo por MCM de las incertidumbres típicas de ? y ?_n, que no se habían tenido en cuenta en el método GUM.

Con la herramienta desarrollada y la técnica de cálculo de incertidumbres de MCM podemos comprobar cómo afecta a los resultados obtenidos y a su FDP resultante, la no consideración de los parámetros de amortiguamiento, introduciendo valores cero en las incertidumbres definidas para ? y ?_n. El resultado sería el siguiente:

Se puede apreciar que los resultados obtenidos por ambos métodos, 14 ms (GUM), y 15 ms (MCM), son similares. Sin embargo, con la implementación de MCM comprobamos lo que ya se intuía al calcular la incertidumbre por GUM: que la FDP de salida no sigue la forma de una distribución normal.
Se observa que la consideración de u(?) y u(?_n) deforma la FDP de salida asemejándola más a una distribución rectangular o una especie de trapezoidal, por lo que no se podría afirmar que con el cálculo de k y los grados efectivos de libertad correspondientes se pudiera calcular una incertidumbre expendida que cubra el 95,45% de probabilidad.

Figura 5: Simulación por MCM del sistema sin tener en cuenta u(ζ) y u(ω_n).

4. Conclusiones

El Laboratorio de Presión y Masa de INTA ha desarrollado y puesto en práctica en los últimos años, métodos de ensayo que permiten la determinación de propiedades dinámicas de instrumentos de medida de presión. Estos métodos van experimentado mejoras sustanciales respecto a las sistemáticas diseñadas originalmente, permitiendo un alto grado de automatización y una rigurosa presentación de los resultados y su incertidumbre.

Las incertidumbres obtenidas son adecuadas para las características del ensayo. Se han calculado por dos métodos diferentes:

el método GUM, analítico, utilizado en la práctica totalidad de sistemas en el ámbito de la calibración. Si bien los resultados fueron, a priori satisfactorios, existen motivos para considerar evaluar la incertidumbre por otro método;
el método de Montecarlo (MCM), estadístico, preferible en aquellos casos donde, bien es complicado aplicar la GUM (como éste) o bien no se dan las condiciones necesarias para cumplir el Teorema del Límite Central (también como es este caso).

Se ha desarrollado un software de cálculo que implementa el cálculo de la incertidumbre del tiempo de respuesta de sensores de presión por el Método de simulación de Montecarlo. De su ejecución y análisis de resultados obtenidos, se concluye que el método que mejor se adapta a la evaluación de incertidumbres en este tipo de ensayos es el MCM.

El hecho de haber desarrollado ambos métodos ha permitido evidenciar que la no consideración de la variabilidad de todos los parámetros característicos de sistema arroja resultados incompletamente satisfactorios, al subestimar la incertidumbre final y no poderse asumir la normalidad necesaria para dar por válido el desarrollo por método GUM.

5. Referencias

[1] H. Fuentes. Desarrollo y aplicación de un procedimiento de caracterización de la respuesta dinámica de sensores de presión y su cálculo de incertidumbres por el método de Montecarlo. Universidad de León. Diciembre de 2015

[2] A.G. Gómez. Calibración de piezoeléctricos por comparación mediante sistema de generación de alta presión dinámico. Universidad Politécnica. Mayo de 2010.

[3] J.P. Damion. Means of Dynamic Calibration for Pressure Transducers. Metrologia, Vol. 30-743-746. 1994.

[4] Eichstädt, S. Analysis of Dynamic Measurements. Evaluation of Dynamic measurement uncertainty. Tesis Doctoral. Berlín 2012.

[5] JCGM (2008). Evaluation of Measurement Data - Supplement 1 to the Guide to the expression of uncertainty of measurement – Propagation of distributions using a Monte Carlo method. Joint Committee for Guides in Metrology, BIPM, JCGM 101, 2008.

Agradecimientos

A mis compañeros del Laboratorio de Presión y Masa por su esfuerzo y profesionalidad, del Laboratorio de Tiempo y Frecuencia por su colaboración, especialmente a Javier Cebollero, y a Ana Gloria Gómez y Jorge Barquero por su valiosa información de experiencias previas.