Estimación de límites de tolerancia. Aplicación al análisis de incertidumbre de códigos termohidráulicos

UNIVERSIDAD POLITÉCNICA DE VALENCIA Máster en Ingeniería de Análisis de Datos, Mejora de Procesos y Toma de Decisiones ESTIMACIÓN DE LÍMITES DE TOLERANCIA. APLICACIÓN AL ANÁLISIS DE INCERTIDUMBRE DE CÓDIGOS TERMOHIDRÁULICOS. PROYECTO FINAL DE MÁSTER DIRECTORA Ana Isabel Sánchez Galdón AUTOR Jose Miguel Prades Marco Valencia, Septiembre de 2016 ÍNDICE 1 INTRODUCCIÓN ..................................................................................................................... 3 2 INTERVALOS DE TOLERANCIA................................................................................................ 6 3 MÉTODO DE WILKS ............................................................................................................... 7 4 METAMODELOS ................................................................................................................... 10 4.1 VECINO MÁS PRÓXIMO ............................................................................................... 13 4.2 MÁQUINA DE SOPORTE VECTORIAL (SVM) ................................................................. 15 4.3 MEDIDAS DE EVALUACIÓN DEL METAMODELO ......................................................... 18 4.4 TRATAMIENTO DEL ERROR .......................................................................................... 19 5 METODOLOGÍA .................................................................................................................... 22 6 CASO DE APLICACIÓN .......................................................................................................... 25 6.1 DESCRIPCIÓN ............................................................................................................... 27 6.2 RESULTADOS MÉTODO DE WILKS ............................................................................... 28 6.3 RESULTADOS OBTENIDOS METAMODELOS ................................................................ 33 6.3.1 RESULTADOS OBTENIDOS CON VECINO MÁS PRÓXIMO .................................... 33 6.3.2 RESULTADOS OBTENIDOS CON MÁQUINA DE SOPORTE VECTORIAL ................. 38 7 CONCLUSIONES ................................................................................................................... 42 8 REFERENCIAS ....................................................................................................................... 43 9 ANEXO ................................................................................................................................. 45 3 1 INTRODUCCIÓN La principal preocupación a la hora de diseñar una central nuclear, es asegurar que los combustibles empleados en la obtención de energía permanezcan confinados dentro del sistema, evitando la liberación de radioactividad la cual podría poner en peligro no sólo a los trabajadores de la central, sino también al entorno donde se ubica y a las poblaciones cercanas. Es por ello que en el ámbito de la industria nuclear se ha establecido el concepto de “Defensa en profundidad”, la cual hace referencia a la implantación de numerosas medidas para asegurar que se adoptan todas las precauciones necesarias para mantener la seguridad. Las bases de esta filosofía son: 1. Un buen diseño, construcción y vigilancia de la planta para prever accidentes. 2. Utilización de diversos sistemas de seguridad y dispositivos que minimicen los daños en caso de producirse un accidente. 3. Adopción de medidas de seguridad adicionales, como sistemas de refrigeración de emergencia del núcleo o fuentes de energía independientes del reactor, para operar los sistemas de seguridad. Los tres puntos anteriores se pueden resumir en hacer uso de la redundancia, la diversidad y grandes márgenes a la hora de realizar los análisis deterministas de seguridad. Estos análisis utilizan simulaciones realizadas mediante programas informáticos (por ejemplo TRACE, RELAP, ..) para estudiar posibles accidentes en la central, escenarios base de diseño, y establecen unos criterios reguladores de aceptación, que son restricciones de seguridad que deben cumplirse. Las hipótesis extremadamente pesimistas suelen ser más simples de modelizar ya que no necesitan tener un análisis de incertidumbre de sus resultados. Mientras que las hipótesis más realistas, se basan en la utilización de mejores estimaciones y datos junto con la evaluación de las incertidumbres para obtener unas mejores predicciones. Esta metodología se conoce como BEPU (Best Estimate Plus Uncertainty). Las aproximaciones BEPU para el análisis de un accidente base de diseño específico asumen que la incertidumbre en las variables de salida, es decir, las figuras de mérito (FOMs) involucradas en los criterios de aceptación del análisis, derivan de las incertidumbres asociadas a las variables de entrada al modelo (condiciones iniciales y de contorno). Estas FOMs son habitualmente valores extremos (máximo y mínimo) de las variables de seguridad durante el transitorio, por ejemplo, la máxima temperatura de vaina. La metodología BEPU cuando se desarrolla para el análisis de accidentes base de diseño se suele denominar evaluación de márgenes debido a que el principal objetivo es comparar la mejor 4 respuesta estimada más su incertidumbre con una figura de mérito predeterminada con el objetivo de establecer el margen de seguridad. Existen diferentes aproximaciones BEPU en la literatura (Martin & O’Dell, 2005) pudiendo ser resumidas las diferentes etapas que constituyen dicha aproximación como: 1. Selección de escenario accidental objeto de estudio. 2. Selección del criterio de seguridad vinculado al escenario accidental objeto de estudio y la FOM involucrada en el criterio de aceptación. 3. Identificación y clasificación de los fenómenos físicos relevantes en base al criterio de seguridad. 4. Selección de los parámetros termohidráulicos significativos. 5. Identificación de los sistemas de seguridad relevantes involucrados en el escenario accidental estableciendo supuestos conservadores respecto a la disponibilidad de dichos sistemas. 6. Desarrollo del modelo de simulación del accidente utilizando códigos termohidráulicos como, por ejemplo, el código TRACE o RELAP. 7. Caracterización de la función de densidad (f.d.d.) para las diferentes parámetros termohidráulicos seleccionados. 8. Ejecutar N simulaciones del código termohidráulico, seleccionando los valores de los parámetros termohidráulicos mediante m.a.s., para obtener la figura de mérito de cada simulación. El número de simulaciones (N) dependerá de diferentes aspectos tales como el número de outputs de interés, el criterio de aceptación seleccionado, etc. 9. Procesar los resultados obtenidos de las N simulaciones para obtener la distribución de probabilidad de la FOM u otro estadístico de interés como, por ejemplo, un percentil, intervalo de tolerancia, etc. 10. Verificar el cumplimiento de los criterios de aceptación para cada FOM en función del método y criterio de aceptación adoptado. Muchos cálculos se desarrollan para estimar la distribución de probabilidad de la figura de mérito o de algún parámetro característico de la distribución, como por ejemplo, un intervalo de tolerancia. Las guías de seguridad recomiendan, en concordancia con la práctica regulatoria, que el valor que se debe comparar con el criterio de aceptación es el límite superior del intervalo de tolerancia del 95/95. 5 La mayoría de aproximaciones BEPU hacen uso, tal y como recomiendan las guías rgeuladoras, del método no paramétrico de Wilks, basado en determinar el número de cálculos mínimo necesarios para verificar los criterios de aceptación de los niveles estándar de tolerancia (standard tolerance levels, STL), que las autoridades reguladoras (Consejo de Seguridad Nuclear, CSN) han fijado, tal y como se ha mencionado anteriormente, en 95/95. Por ejemplo, se encuentra extendido el uso del estadístico de orden 1 que garantiza con un tamaño de muestra de 59 la obtención del intervalo de tolerancia unilateral 95/95. La principal ventaja de usar el estadístico de orden 1 del método de Wilks es el reducido número de simulaciones necesarias para la obtención del intervalo de tolerancia dado que se intenta hacer el menor número de simulaciones como consecuencia del elevado coste computacional. Sin embargo, el valor obtenido es excesivamente conservador. Uno de los retos que se plantean, en el contexto de la metodología BEPU, es el uso de resultados conservadores, pero más realistas y precisos con un coste computacional aceptable. Por ello, se ha propuesto en la literatura diferentes alternativas con el objetivo de obtener resultados más realistas como son el uso de estadísticos de orden superior o el uso de metamodelos, también llamados modelos sustitutivos o surrogados, que sustituyan el código termohidráulico. En este contexto, el presente trabajo se centra en la comparación de los resultados obtenidos en la estimación de intervalos de tolerancia unilaterales 95/95 utilizando el método de Wilks y el uso de metamodelos. Concretamente, se estudia el caso de un accidente base de diseño con pérdida de refrigerante debido a una rotura grande (Large-Break Loss of Coolant Accident, LBLOCA) en una tubería de la rama fría del sistema de refrigeración de un reactor de agua presurizada (Pressurized Water Reactor, PWR). La variable de interés analizada es la temperatura máxima de vaina (MPCT). El documento se estructura como se detalla a continuación. En la sección 2 se describe el concepto de intervalo de tolerancia. En la sección 3 se presenta el Método de Wilks, método no paramétrico que permite la estimación de intervalos de tolerancia y cuyo uso se encuentra extendido en el ámbito nuclear. En la sección 4 se describen los dos metamodelos utilizados en el presente trabajo (vecino más próximo y máquinas de soporte vectorial) y el tratamiento del error asociado al uso de un metamodelo. En la sección 5 se presenta la metodología utilizada y en la sección 6 el caso de aplicación desarrollado. Finalmente, en la sección 7 se resumen las conclusiones obtenidas. 6 2 INTERVALOS DE TOLERANCIA En el contexto de la metodología BEPU y con el objetivo de identificar unos límites de seguridad más realistas, la industria nuclear internacional ha adoptado la aproximación basada en intervalos de tolerancia. Un intervalo de tolerancia (Krishnamoorthy, 2009) es un intervalo, basado en una muestra aleatoria, que se espera contenga una proporción especificada de la población muestreada. Concretamente, en el análisis de seguridad de códigos termohidráulicos resulta de interés el límite superior de tolerancia, el cual se evalúa sujeto a la condición de que al menos un p% de la población, p.e. el 95%, estará por debajo del límite con un cierto nivel de confianza , p.e. el 95%. Así, sea 𝑋 una variable aleatoria continua con función de distribución (cdf) 𝐹(𝑥)=𝑃(𝑋≤𝑥). Para un 𝛾(0<𝑝<1) especificado la inversa de la cdf se define por: 𝐹−1(𝛾)=𝑖𝑛𝑓{𝑥:𝐹(𝑥)≥𝛾} (1) siendo 𝐹−1(𝛾) el valor de 𝑥 para el cual 𝐹(𝑥)=𝑃(𝑋≤𝑥)=𝛾 (2) Sea 𝑿=(𝑋1,𝑋2,…,𝑋𝑛) una muestra aleatoria de la población. Para definir un intervalo de tolerancia necesitamos especificar un nivel de cobertura, 𝛾, y un nivel de confianza, 𝛽. En las aplicaciones prácticas 𝛾 y 𝛽 toman valores en el conjunto {0.90,0.95,0.99}. El intervalo se construirá usando una muestra aleatoria 𝑿 y se requiere que contenga un 𝛾% de la población o superior con un nivel de confianza 𝛽. Formalmente, un intervalo de tolerancia unilateral (𝛾,𝛽) de la forma (−∞,𝑈(𝑿)) satisface la condición 𝑃{𝑃(𝑋≤𝑈(𝑿)|𝑿)≥𝛾}=𝛽 (3) Esto es 𝑈(𝑿) se determina de modo que al menos un 𝛾% de la población sea menor o igual a 𝑈(𝑿) con un nivel de confianza 𝛽. El intervalo (−∞,𝑈(𝑿)) se denomina intervalo de tolerancia unilateral y 𝑈(𝑿) límite superior de tolerancia unilateral. La ecuación (3) puede ser escrita como 𝑃{𝑞𝛾<𝑈(𝑿)}=𝛽 (4) Siendo 𝑞𝛾 el percentil 𝛾. 7 En la literatura se han propuesto diferentes métodos para el cálculo de intervalos de tolerancia. Autores como Glaeser, Guba, Makai, Pál, o Wallis entre otros (Gaeser, H. 2000; Guba, et al., 2003; Nutt & Wallis, 2004; Wallis, G.B. 2005), centran sus estudios en estimar los límites de tolerancia de variables aleatorias cuya función de probabilidad se desconoce a partir de la metodología de estadísticos ordenados de muestras aleatorias. Se trata de un método no paramétrico introducido por S.S. Wilks (Wilks, S.S. 1967). Métodos exactos para evaluar intervalos de tolerancia asumiendo normalidad de la variable de interés pueden consultarse en (Owen, 1964). Asimismo, los factores asociados al cálculo de intervalos de tolerancia para diferentes tamaños de muestra y diferentes softwares para evaluar los factores de tolerancia se encuentran también disponibles. No obstante, existen otras aproximaciones tanto paramétricas como no paramétricas que permiten estimar un intervalo de tolerancia. 3 MÉTODO DE WILKS En el análisis de códigos termohidráulicos se encuentra extendido el uso del método de Wilks (orden 1) para la estimación de intervalos de tolerancia principalmente por el número reducido de simulaciones necesarias (cada simulación tiene un elevado coste computacional) y por su carácter conservador. En el método de Wilks, sea w el límite superior de tolerancia, entonces 𝑃(−∞≤𝑋𝑖≤𝑤)≥𝛽 (5) la estimación de w puede escribirse como w=w(X). Sobre todas las posibles muestras se puede requerir que 𝑃(−∞≤𝑋𝑖≤𝑤)≥𝛽)≥𝛾 (6) La expresión anterior define un intervalo de tolerancia /. La interpretación de la ecuación anterior es que al menos 100x % de la población es cubierta por el intervalo de tolerancia con una probabilidad . El límite superior w del intervalo de tolerancia dado por la Ecuación anterior se interpreta como un límite de tolerancia conservador, es decir, el límite indica el valor extremo que cubre al menos 100 de la población con una confianza . El método de Wilks (Wilks, S.S, 1967; UPV-CSN, 2006) permite estimar los intervalos de tolerancia de una determinada variable aleatoria, a partir de una muestra cuyo tamaño mínimo es función de la cobertura () /confianza () deseado. Es por tanto interesante que el tamaño de muestreo sea pequeño. Cuando los elementos de la muestra aleatoria (de tamaño N) de la 8 variable se ordenan de menor a mayor, el elemento que ocupa la posición r-ésima es el elemento r-ésimo menor y el (N-r+1) mayor de la muestra, al que se denomina estadístico de orden r. Los estadísticos de orden 1 y N son el mínimo y el máximo de la muestra. Los intervalos de tolerancia que proporcionan los estadísticos de orden no dependen de la función de distribución de la variable. La teoría de los estadísticos de orden indica el tamaño de muestra necesario para garantizar un cierto nivel /. Así, un límite superior de tolerancia de la variable aleatoria X para un nivel 95/95 es también un valor aleatorio obtenido de los valores de la muestra, el cual es mayor que el percentil 95 de X con una confianza del 0.95. La NRC ha reconocido como adecuado un nivel 95/95 para la estimación de la MPCT que pueda compararse con el límite de seguridad (Martin, R.P, & O’Dell, L.D. 2005). El punto de partida de la configuración del problema es que se obtiene una muestra de tamaño n de los parámetros de entrada según su distribución de probabilidad correspondiente. Esta muestra se utiliza como entrada al código de simulación y se obtienen un número n de salidas de la variable de interés. La distribución de probabilidad de la f(y) de salida es una función desconocida. Los límites de tolerancia (superior, U, e inferior, L) se obtienen utilizando el método de Wilks como: 𝑃(∫𝑓(𝑦)𝑑𝑦>𝛾 𝑈 𝐿)=𝛽 (7) Donde L o U se seleccionan, respectivamente, como –∞ o +∞ para intervalos de tolerancia unilaterales. La expresión (7) define un intervalo de tolerancia (L,U) con un nivel de cobertura/confianza /. La interpretación de la ecuación anterior es que al menos 100 % de la población es cubierta por el intervalo de tolerancia con una probabilidad . En (Wilks, 1967) se presenta la formulación detallada para intervalos de tolerancia unilaterales y bilaterales. En el caso concreto de intervalos de tolerancia unilaterales considerando el estadístico de orden r se verifica: 𝛽≤1 − ∑(𝑛 𝑘)𝛾𝑘(1−𝛾)𝑛−𝑘 𝑛 𝑘=𝑛−𝑟+1 (8) Si se utiliza el estadístico de primer orden se obtiene la siguiente expresión 1−𝑛≥𝛽 (9) Por lo tanto, si se ordena la muestra de la salida de menor a mayor, el valor máximo de la muestra infiere el percentil  de la población de salida con una confianza de β. Por ejemplo, de acuerdo con las prácticas de la reglamentación actual, para un nivel estándar de tolerancia (STL) 9 95/95, utilizando el estadístico de primer orden, se debería seleccionar una muestra de 59 ejecuciones del código. En el análisis de incertidumbre de transitorios se encuentra extendido el uso del método de Wilks con el estadístico de orden 1 para un intervalo de tolerancia/confianza 95/95. Sin embargo, se trata de una estimación muy conservadora. Con el objetivo de obtener aproximaciones menos conservadoras es posible utilizar estadísticos de orden superior a costa de aumentar el número de ejecuciones del código termohidráulico. En la Tabla 3.1 se presentan los tamaños de muestra n necesarios para la estimación de un intervalo de tolerancia unilateral 95/95 utilizando estadísticos de orden superior obtenidos a partir de la Ec. (8). Tabla 3.1: Mínimo número de simulaciones método de Wilks intervalo de tolerancia unilateral 95/95 Orden Wilks Mínimo número de simulaciones (n) Orden Wilks Mínimo número de simulaciones (n) 1 59 4 153 2 93 5 181 3 124 16 Por otro lado, se conoce con el nombre de atributo a la variable predictora y característica a un atributo transformado que es usado para definir el hiperplano. La elección de la representación más adecuada del universo estudiado, se realiza mediante un proceso denominado selección de características. Support Vector Regression (SVR) es una particularización de la máquina de soporte vectorial. SVR utiliza los mismos principios que la SVM para la clasificación que caracteriza el algoritmo de máximo margen: que consiste en realizar un mapeo de los datos de entrenamiento x ∈ X, a un espacio de mayor dimensión F a través de un mapeo no lineal φ: X → F, donde sí se puede llevar a cabo una regresión lineal. En SVR, el objetivo es encontrar 𝑦(𝒙) que tiene las mayores desviaciones  (at most  deviations) de cada una de las salidas de los inputs de entrenamiento. El modelo SVR viene dado por: 𝑦(𝒙)=∑ (𝑎𝑖−𝑎𝑖∗)𝐾(𝒙𝑖𝒙)+𝑏 (17) donde 𝐾(𝒙𝑖𝒙) es la función núcleo (kernel), 𝒙𝑖 son los diferentes puntos del diseño de experimentos original y 𝒙 es el punto del espacio en el cual es modelo surrogado es evaluado. Los parámetros 𝑎𝑖,𝑎𝑖∗ 𝑦 𝑏 se obtienen durante el proceso de ajuste. La metodología base de SVR se puede resumir de la siguiente forma. En SVR, la estimación de 𝑓󰆹 se selecciona entre la familia de funciones 𝑓:𝑥∈𝑅𝑑→〈𝒘,∅(𝒙)〉𝐻+b (18) donde ∅ es una función de Rd en un espacio de Hilbert, w  H y b  R son parámetros obtenidos a partir del conjunto de datos de entrenamiento. En la literatura se han desarrollado diferentes estrategias para el aprendizaje de los parámetros w y b. Los coeficientes w y b se obtienen minimizando la siguiente ecuación estructural: 𝑅𝑆𝑉𝑀(𝐶)=𝐶∑𝑐(𝑓󰆹(𝑋𝑖),𝑦𝑖 )+1 2 𝑛 𝑖=1 ‖𝑤‖2 (19) donde 𝑐(𝑓󰆹(𝑋),𝑦)={|𝑓󰆹(𝑋)−𝑦|− 𝜀 𝑠𝑖|𝑓󰆹(𝑋)−𝑦|≥𝜀 0 𝑜𝑡𝑟𝑜𝑠 (20) En (19) el primer término es el riesgo empírico donde n representa el número total de observaciones. El riesgo empírico se mide por la función pérdida -insensitiva en (20). En la aproximación original se utiliza una función de pérdida  como un criterio de calidad para la 17 regresión. El segundo término en (19) es el término de regularización. El parámetro C, el cual es definido por el usuario, es el parámetro de regularización el cual determina la frontera entre el riesgo empírico y la categoría. En la función pérdida  -insensitiva en (20),  es el tamaño del tubo que controla el rango de penalización. Las desviaciones mayores que  son penalizadas, al igual que C el valor de  se selecciona por el usuario. Tras resolver el problema de minimización, (19) se obtiene la función final de predicción como: 𝑓(𝑋,𝛼𝑖,𝛼𝑖∗)=∑ (𝑎𝑖−𝑎𝑖∗)𝐾(𝑋,𝑋𝑖)+𝑏 𝑛 𝑖=1 (21) donde 𝐾(𝒙𝑖𝒙) es la función núcleo (kernel) y los parámetros 𝑎𝑖,𝑦 𝑎𝑖∗ se obtienen durante el proceso de ajuste. En la práctica se utilizan una serie de funciones Kernel entre las más utilizada en la práctica se encuentran la polinomial, base radial gaussiana y base radial exponencial. Se define el kernel polinomial como: 𝐾(𝑋,𝑋′)=(𝑋 𝑋′+𝑐)𝑑 𝑐>0 (22) Donde d es el grado del polinomio. La función base radial en su forma gaussiana es uno de los kernels más ampliamente utilizados. Dicho Kernel viene dado por: 𝐾(𝑋,𝑋′)=𝑒𝑥𝑝(−𝛾‖𝑋 −𝑋′‖2) (23) donde el parámetro 𝛾 controla la flexibilidad del kernel de forma similar al parámetro d en el kernel polinomial. Las funciones de base radial exponencial son de la forma: 𝐾(𝑋,𝑋′)=𝑒𝑥𝑝(−‖𝑋 −𝑋′‖‖𝑋 −𝑋′‖ 2𝜎 ) (24) 18 La figura 4.2-1 muestra un ejemplo de la función de regresión lineal unidimensional con - épsilon intensiva - banda. Las variables miden el costo de los errores en los puntos de entrenamiento. Estos son cero para todos los puntos que están dentro de la banda. Figura 4.2-1 Regresión lineal unidimensional con los márgenes épsilon 4.3 MEDIDAS DE EVALUACIÓN DEL METAMODELO Con el objetivo de evaluar el funcionamiento de los diferentes metamodelos se utilizan medidas de calidad. Estas medidas miden aspectos tales como la precisión, la robustez y la eficiencia. Otras medidas adicionales que pueden ser de interés son las relacionadas con la selección de parámetros y la interpretación de la estructura del modelo. La precisión es un importante indicador del funcionamiento de un metamodelo. Las métricas de precisión deben reflejar la salida del metamodelo y su desviación desde la salida de la simulación. La calidad de la aproximación ha sido evaluada a partir del cálculo de las métricas que se describen a continuación. El error cuadrático medio (MSE) que proporciona una evaluación de la precisión de la Estimación el cual es evaluado como: 𝑀𝑆𝐸=1 𝑁∑(𝑦𝑖  𝑁 𝑖=1 −𝑦𝑖)2 (25) Donde N es el tamaño de muestra utilizado para validación. El coeficiente de determinación 𝑅2 19 𝑅2=1− 𝑀𝑆𝐸 𝑉𝑎𝑟(𝑦) (26) siendo 𝑉𝑎𝑟(𝑦) la varianza de la muestra de validación. El Máximo Error Absoluto (MAE) el cual refleja la presencia de predicciones de mala calidad en áreas locales y se calcula como: 𝑀𝐴𝐸=max|𝑦𝑖−𝑦𝑖 | (27) El Error porcentual absoluto medio (MAPE, Mean absolute percentage error) que expresa la exactitud como un porcentaje del error y es evaluado como: 𝑀𝐴𝑃𝐸=1 𝑁∑ |(𝑦𝑖−𝑦𝑖 )/𝑦𝑖|∙100 𝑁 𝑖=1 (28) 4.4 TRATAMIENTO DEL ERROR A la hora de obtener unos buenos resultados en los análisis de un sistema, en el ámbito de la ingeniería, existen múltiples fuentes de errores a las que hay que hacer frente. Realizando una predicción conservadora se garantiza que se tienen en cuenta las incertidumbres y los errores en los análisis del sistema, a través del uso de cálculos o aproximaciones que tienden a estimar con seguridad la respuesta del sistema. Tradicionalmente, los factores de seguridad se han utilizado ampliamente para este propósito, a pesar de que su eficacia es cuestionable. Cuando se utilizan metamodelos para predecir los valores de las variables a estimar, se sabe muy poco en la práctica como para proporcionar estimaciones conservadoras, y de qué manera las estrategias conservadoras afectan en el proceso de diseño. La mayoría de los metamodelos están diseñados de modo que existe una posibilidad del 50% de que el valor de la predicción obtenida sea mayor que el valor real. Es por ello que se buscan alternativas para modificar los metamodelos tradicionales de manera que este porcentaje se incremente, obteniéndose un valor más conservador pero ocasionando el menor impacto en la precisión. Dado que los metamodelos son conservadores tienden a sobreestimar la respuesta real, por ello existe un compromiso entre la exactitud y el conservadurismo. Existe pues un equilibrio entre el riesgo de sobredimensionamiento y el riesgo de un diseño débil. La estrategia conservadora más clásico consiste en influir en la predicción de la respuesta mediante una constante multiplicativa o aditiva. Tales enfoques son llamados empírica porque la elección de la constante es de alguna manera arbitraria y se basa en el conocimiento previo del problema de ingeniería considerado. Sin embargo, es difícil predecir como de eficiente es su aplicación a los metamodelos, y cómo es posible fijar dichas cantidades. 20 Alternativamente, es posible utilizar el conocimiento estadístico del metamodelo para construir intervalos de confianza unilateral en la predicción. A la hora de definir los predictores conservadores se supone que un estimador conservador es un estimador que no subestima el valor real. Por lo tanto, un estimador conservador se puede obtener por simple adición de un margen de seguridad positiva al estimador insesgado: 𝑦𝑆𝑀(𝑥)=𝑦(𝑥)+𝑆𝑚 (29) Alternativamente, es posible definir un estimador conservador utilizando factores de seguridad: 𝑦𝑆𝐹(𝑥)=𝑦(𝑥)∗𝑆𝐹 (30) En el contexto de la obtención de los modelos sustituto, los márgenes de seguridad son más convenientes. De hecho, cuando se utiliza un factor de seguridad, el nivel de seguridad depende del valor de respuesta, mientras que en la mayoría de metamodelos asumen que el error es independiente de la media de la respuesta. Por lo tanto, en este trabajo se utiliza únicamente los márgenes de seguridad para definir un estimador conservador. El valor del margen repercute directamente en el nivel de seguridad y en el error en la estimación. Por lo tanto, hay que encontrar el margen que garantiza un nivel específico de seguridad con el menor impacto en la exactitud de las estimaciones. El margen de seguridad puede ser constante, o depender de la ubicación. A continuación, se proporcionan varias maneras de diseñar el margen, usando métodos paramétricas y no paramétricas, de punto a punto y de márgenes constantes. Las estimaciones conservadoras son parciales, y un alto nivel de conservadurismo sólo puede conseguirse a costa de la precisión. Por lo tanto, la calidad de un método sólo se puede medir como un compromiso entre el conservadurismo y la precisión. Con el fin de evaluar un desempeño global de los métodos, se propone definir los índices de precisión y conservadurismo. La medida más ampliamente utilizada para comprobar la exactitud de un modelo sustitutivo es la raíz del error cuadrático medio (eRMS), que se define como: 𝑒𝑅𝑀𝑆 =(∫𝑒𝑖2𝑑𝑥 𝐷)1 2 (31) 21 Donde 𝑒𝑖=𝑦(𝑥)−𝑦(𝑥) sabiendo que 𝑦(𝑥) es el valor de la predicción e 𝑦(𝑥) es el valor real de la variable respuesta para el caso i. Margen de Seguridad constante mediante técnicas de validación cruzada: En este caso, se tiene en cuenta el diseño de los metamodelos conservadores, cuando se aplica el mismo margen de seguridad por todas partes en el dominio de diseño. De ahora en adelante tales estimadores se conocerán como CSM (Margen Constante de Seguridad). En términos de la función de distribución acumulativa (CDF) de los errores, Fe, el margen de seguridad Sm para un valor de conservadurismo, %c, se da como: 𝑆𝑚=𝐹𝑒−1∗(%𝑐 100) (32) El objetivo es diseñar el margen de seguridad de tal manera que se asegura la ecuación anterior. La distribución de error real es, en la práctica, desconocida. Se propone estimar empíricamente mediante técnicas de validación cruzada para obtener el margen. La validación cruzada (CV) es un proceso de estimación de errores mediante la construcción de modelos sustitutivos sin algunos de los puntos, y el cálculo de los errores en estos puntos quedan fuera. El proceso se repite con diferentes conjuntos de puntos que se han quedado fuera con el fin de obtener estimaciones estadísticamente significativas de los errores. El proceso prosigue al dividir el conjunto de puntos de datos n en subconjuntos. El modelo sustitutivo se ajusta a todos los subgrupos, excepto uno, y el error se comprueba en el subconjunto que quedó fuera. Este proceso se repite para todos los subconjuntos para producir un vector de errores de validación cruzada, eXV. La CDF (distribución uniforme continua) empírica FXV, definida por los valores n de eXV, son una aproximación de la verdadera distribución Fe. Con el fin de diseñar el margen, reemplazamos Fe en la ecuación NN por FXV: 𝑆𝑚=𝐹𝑋𝑉−1∗(%𝑐 100) (33) 22 5 METODOLOGÍA En esta sección se presenta la metodología de análisis de incertidumbres coherente con la metodología de análisis “Best Estimate Plus Uncertainty, BEPU” descrita en la sección 1, que hace uso de códigos de sistemas realistas y tienen en cuenta el efecto de las incertidumbres, aplicando en algunos casos herramientas estadísticas. La metodología consiste en 4 etapas las cuales se describen brevemente a continuación. Primera etapa. Definición del escenario de estudio. El objetivo es la selección del escenario accidental de interés, escogiéndose el iniciador y la(s) secuencia(s). En el caso de aplicación que se presenta en el presente trabajo se considera la secuencia de éxito de un LBLOCA Segunda etapa. Identificación y caracterización de las variables influyentes. El objetivo es la elección a priori de las variables “input” que se consideran más influyentes en los resultados del accidente. La elección correcta de “inputs” se puede derivar de un estudio completo de las fenomenologías que resultan fundamentales en el caso de estudio. En este apartado se puede aprovechar los resultados de “peer reviews” previos sobre estudio fenomenológico del escenario accidental bajo análisis, por ejemplo un LBLOCA en el caso de aplicación. Adicionalmente resulta necesario la caracterización estadística de las variables “input” del modelo sujetas a incertidumbre y cuyo efecto en los resultados de salida “outputs” se pretende estudiar. Para ello, se establecen funciones de densidad de probabilidad (pdf) que las representen. Tercera etapa: Análisis de sensibilidad (screening) Cuando el número de variables inputs es elevado se puede realizar un screening previo que permita identificar dentro de todo el conjunto de variables “input” aquellas que sean significativamente importantes en el comportamiento final de las variables “output”. Cuando se considere despreciable el efecto de posibles interacciones entre variables es posible utilizar diseños de experimentos, como el método Placket-Burman, que con muy pocas simulaciones consiguen fácilmente identificar los efectos principales. Sólo los “inputs” seleccionados en esta etapa se considerarán en las siguientes etapas del análisis. 23 Cuarta etapa. Análisis de incertidumbre El análisis de incertidumbre se puede abordar mediante dos métodos distintos: obtención de límites de tolerancia a partir de estadísticos de orden y el uso de metamodelos. En el primer caso, el objetivo es calcular los intervalos de tolerancia para las variables aleatorias de salida “outputs”. Respecto al uso de metamodelos el objetivo se centra en establecer modelos de predicción de las variables de interés en la estimación de las variables de salida de interés y la propagación de la incertidumbre a través de las diferentes ecuaciones de predicción. En este caso En la Figura 5-1 se presenta el diagrama de flujo de la metodología propuesta. El presente trabajo se centra en la cuarta etapa de la metodología. 24 Figura 5-1 Etapas de la metodología propuesta para el análisis de incertidumbres En la Figura 5-2 se muestra el proceso de sustitución del código termohidráulico por un metamodelo. Figura 5-2 Metamodelo como alternativa al código termohidráulico. COMPARACIÓN Definición del escenario de estudio Identificación y caracterización de variables influyentes Muestreo variables influyentes Código termohidráulico Resultados simulaciones Sensibilidad (screening) Análisis incertidumbre Estadísticos de orden Metamodelos METAMODELO MPCT Error Predicción MPCTIntervalo de Tolerancia 95/95 X11 ………. X1m PCT1 ……….. Xn1 ………… Xnm PCTn 6 CASO DE APLICACIÓN Este trabajo se centra en uno de los Accidentes de Diseño Básicos más comunes (Design Basis Accidents) en los reactores de agua a presión (Power Water Reactors, PWRs), la pérdida de refrigerante producida por una rotura grande (Large-Break Loss of Coolant Accident, LBLOCA) en la tubería de la rama fría de refrigeración. A consecuencia del transitorio, se produce una rápida despresurización del sistema primario y consecuentemente se activa la correspondiente señal de SCRAM que indica que hay baja presión, continuando con los acumuladores de inyección y después la inyección del refrigerante frío de emergencia desde el Sistema de Inyección de Baja Presión (Low Pressure Injection System, LPIS), para evitar la no recuperación del núcleo. Con el SCRAM se consigue una reducción de la energía térmica y con el sistema de inyección de agua mediante los acumuladores y el Sistema de Inyección de Baja Presión se reduce la temperatura del núcleo. El calor disipado a través del sistema secundario no es considerado debido a la rápida despresurización del sistema primario. La planta seleccionada ha sido una típica de 4 lazos PWRWestinghouse, cuya referencia es Zion Nuclear Power Plant (NPP), y el código informático para modelizar el sistema termohidráulico es el TRACE V5.0 Patch 4. Los sistemas de seguridad que actúan cuando se produce el transitorio accidental son los sistemas de parada de emergencia, acumuladores y sistemas de inyección de baja presión. El Esquema 6-1 muestra el sistema primario modelado para TRACE V5.0 Patch 4 mediante el paquete de SNAP, que incluye un componente tridimensional de tipo recipiente, que representa la vasija de presión del reactor incluyendo el núcleo. También incluye la parte principal de los cuatro circuitos de refrigeración (tubos, 4 SGS, PRZ y 4 bombas). Además, también se han modelado los sistemas de seguridad necesarios en el transitorio de apoyo al sistema primario (4 acumuladores, 4 inyecciones de baja presión, todo en piernas frías). La gran ruptura se simula como una ruptura de doble guillotina por tres válvulas y dos componentes de rotura, en la rama fría. El Esquema 6-2 describe el sistema secundario, que incluye generadores de vapor (4 SG) asociados cada uno con los circuitos de refrigeración correspondientes en el sistema primario. 32 Figura 6.7. Distribución de los 59 valores MPCT Figura 6.8. Histograma MPCT (59) En la figura 6.7, se pude ver el valor que toman cada una de las 59 simulaciones para la variable de salida MPCT, y con una línea a trazos roja el IT 95/95. Mientras que en la figura 6.8 se puede ver el histograma para la variable MPCT y al igual que en el primer gráfico se ha trazado una línea discontinua roja para indicar el IT 95/95. El valor obtenido para el límite superior del intervalo de tolerancia 95/95 es de 1294.7 K. En la Tabla 6.1 se resumen los resultados obtenidos mediante la aplicación del método de Wilks con estadísticos de orden superior. Orden del estadístico Tamaño de muestra Límite superior del intervalo de tolerancia 95/95 1 59 1294.7 2 93 1286.9 3 124 1268.1 4 153 1266.5 Tabla 6.1. Límites superior del intervalo de tolerancia para estadísticos de orden 1, 2, 3 y 4. En la tabla 6.1 se muestra como el aumento en el orden del estadístico aumenta la proximidad de la estimación puntual al valor de referencia obtenido con las 1000 simulaciones (1259.5 K) manteniéndose en todos los casos por debajo del valor de seguridad de 1477 K. 33 6.3 RESULTADOS OBTENIDOS METAMODELOS Una alternativa al uso de estadísticos de orden en el análisis de incertidumbre es el uso de metamodelos. Con el objetivo de comparar el análisis de incertidumbre bajo los dos enfoques (estadísticos de orden superior vs metamodelos) se han ajustado diferentes metamodelos utilizando los tamaños de muestras correspondientes a los estadísticos de orden 2, 3 y 4 Los metamodelos seleccionados para realizar la comparativa son, como se ha indicado anteriormente, el método del vecino más próximo (KNN) y la máquina de soporte vectorial (SVR). Las variables de entrada utilizadas para entrenar los dos metamodelos corresponden a las 45 variables de entrada identificadas como importantes siendo el output analizado la MPCT. Las simulaciones utilizadas para entrenar los metamodelos corresponden a m.a.s. de tamaños 93, 124 y 153. La función de densidad de la FOM seleccionada como variable de interés, en este caso la MPCT, se obtiene por muestreo de Monte Carlo de los parámetros de entrada propagando la incertidumbre utilizando el metamodelo. El resultado final de la propagación de incertidumbre se ha realizado con un tamaño de muestra de 1000. El uso del metamodelo impone, como se ha comentado anteriormente, una nueva fuente de incertidumbre. Tanto la bondad del ajuste del modelo a los datos como la precisión de los nuevos outputs son dos aspectos importantes a considerar para el modelo surrogado. 6.3.1 RESULTADOS OBTENIDOS CON VECINO MÁS PRÓXIMO El algoritmo KNN se ha implementado utilizando la librería de R FNN (Beygelzimer et a. 2015). Los modelos han sido entrenados con tamaños de muestra de 153, 124 y 93 con el objetivo de comparar los resultados obtenidos con el metamodelo y el método de Wilks. El valor de K, es decir, el número de inputs más cercanos que se consideran para realizar la estimación de la variable respuesta ha sido en los tres casos igual a 3 y la métrica de distancia la distancia euclídea. La validación del modelo se ha realizado, puesto que se ha considerado utilizar el total de datos como muestra de entrenamiento, utilizando “Leave one out crossvalidation”. En la figura 6.9 se muestran las predicciones obtenidas vs a los valores reales de la variable respuesta, MPCT, con un tamaño de muestra de 153. 34 Figura 6.9. Predicciones de MPCT frente a valores reales (153 muestras) En las Figuras 6.10 y 6.11 se presentan la distribución y el histograma de los errores del modelo, respectivamente. En la misma se observa un valor aislado con un error de -10 K y varios puntos con errores superiores en valor absoluto de 3 K. Figura 6.10. Residuos obtenidos con el modelo de Vecino Más Próximo (153 muestras) 35 Figura 6.11. Histograma de los residuos (153 muestras) Las métricas de calidad del ajuste se muestran en la Tabla 6.2. Modelo MAE MSE RMSE MAPE VMP n=153 1.299782135 3.195236020 1.787522313 0.001076358 Tabla 6.2. Medidas de calidad del ajuste VMP (n=153) A partir del modelo ajustado se ha obtenido 10000 predicciones de la PCT. Con el objetivo de obtener un metamodelo conservador se ha utilizado un margen de seguridad, tal como se ha descrito en la sección 4.4. utilizando las ecuaciones (26) y (29), para un valor de conservadurismo del 95%. Utilizando las 10000 predicciones obtenidas con el metamodelo conservador se obtiene un valor del límite superior de tolerancia 95/95 de 1263.2 K. En las Figuras 6.12 y 6.13 se muestran, respectivamente, el histograma de los valores obtenidos con el metamodelo conservador y el valor del límite de tolerancia superior y la función de distribución obyenida con la 1000 simulaciones obtenidas del código termohidráulico TRACE y las 10000 simulaciones utilizando el metamodelo VMP. 36 Figura 6.12. Histograma de la MPCT obtenida con el metamodelo conservador y límite de tolerancia superior 95/95 (153 muestras) El mismo procedimiento seguido para la muestra de tamaño 153 se ha seguido para tamaños de muestra de 124 y 93. No se ha realizado el ajuste con tamaño de muestra de 59 dado el elevado número de variables de entrada. En las Figura 6.13 se muestra las predicciones obtenidas de la MPCT frente a los valores reales para los nuevos tamaños de muestra mientras en la Figura 6.12 se muestra la distribución de los residuos. Figura 6.13. Predicciones de la MPCT frente valores reales con tamaños de muestra n=124 (izquierda) y n=93 (derecha). 95-95,2871 Límites LST: 1261,28 1100 1140 1180 1220 1260 1300 1340 0 200 400 600 800 1000 frecuencia Límite tolerancia superior 95/95 37 Figura 6.14. Residuos obtenidos con VMP y tamaños de muestra n=124 (izquierda) y n=93 derecha). Las métricas de calidad del ajuste para un tamaño de muestra de 124 y 93 se presentan en la Tabla 6.3 Modelo MAE MSE RMSE MAPE VMP n=124 1.531273469 7.814239907 2.795396199 0.001272935 VMP n=93 1.684506398 11.806721863 3.436091073 0.001376384 Tabla 6.3. Medidas de calidad del ajuste (n=124, n=93) Los valores obtenidos del límite superior de tolerancia 95/95 para n=124 y n=93 son respectivamente de 1264.26 K y 1262.5 K, respectivamente. Si comparamos dichos valores con el de referencia observamos que, en todos los casos, el intervalo de toelrancia contiene el valor del percentil de referencia. Con el objetivo de verificar el comportamiento observado para el modelo VMP se han ajustado 50 modelos para cada uno de los diferentes tamaños de muestra y obtenido el límite superior de tolerancia utilizando 10000 predicciones. La Figura 6.15 muestra el diagrama Box-Whisker de los límites superior de tolerancia 95/95 obtenidos para los diferentes tamaños de muestra. 38 Figura 6.15. Box-Whisker límites superior de tolerancia 95/95 de MPCT para diferentes tamaños de muestra. Si se comparan los resultados con los obtenidos con el método de Wilks para los diferentes tamaños de muestra se observa que se obtienen, en todos los casos, límites conservadores pero más realistas utilizando el metamodelo VMP. 6.3.2 RESULTADOS OBTENIDOS CON MÁQUINA DE SOPORTE VECTORIAL El algoritmo SVR se ha implementado utilizando la librería de R “e1071” (Meyer et al., 2015). Los modelos han sido entrenados con tamaños de muestra de 153, 124 y 93 con el objetivo, al igual que en el caso del apartado anterior de comparar los resultados obtenidos con el metamodelo SVR y el método de Wilks. En el caso de SVR, resulta necesario la selección de un conjunto de parámetros como son la función pérdida , el parámetro de regularización C y el parámetro del kernel Gaussiano, . La Tabla 6.4 muestra los parámetros seleccionados en el caso de aplicación. Parámetro Tipo Núcleo Base radial  0.02 Coste 1  0.1 Tabla 6.4 Parámetros SVR 39 En la figura 6.16 se muestran las predicciones obtenidas de la variable respuesta MPCT, a partir de las 153 muestras. Se aprecia como los valores predichos no se acaban de aproximan a los valores reales. Por este motivo los puntos en la gráfica no se ajustan perfectamente a la recta, se pueden observar dos zonas distintas respecto al comportamiento del error. La primera de ellas, que corresponde a valores por debajo de 1200 se sobreestiman las predicciones de MPCT, es decir, los valores predichos son siempre superiores a los valores reales. Mientras que en el segunda zona ocurre exactamente lo contrario, se subestiman las predicciones de MPCT, es decir, los valores predichos adopta siempre valores inferiores a los reales. Figura 6.16. Predicciones de MPCT frente a valores reales (153 muestras) En la Figura 6.17 se puede observar la distribución de los errores del modelo anterior. Se aprecian como la mayoría de los residuos toman valor 3 ó - 3, lo que corresponde con los dos tramos de valores predichos para la variable respuesta que se ha mencionado anteriormente. 1150 1200 1250 1140 1160 1180 1200 1220 1240 1260 1280 y y ^ 40 Figura 6.17. Residuos obtenidos con el modelo de SVM (153 muestras) En la Tabla 6.5 se muestran el valor de las diferentes métricas para el modelo SVR y tamaño de muestra 153. Modelo MAE MSE RMSE MAPE SVM (e1071) 2.644030353 7.132441896 2.670663194 0.002188595 Tabla 6.5. Métricas de calidad modelo SVR (n=153) Siguiendo el mismo procedimiento que en el caso de KNNse ha obtenido el límite superior del intervalo de tolerancia 95/95 obteniéndose un valor de 1259.7 K. Al igual que en el caso anterior se ha ajustado el modelo con tamaños de muestra de 93 y 125. En las Figuras 6.18, 6.19 y Tabla 6.6 se muestran el ajuste obtenido y las métricas de calidad para los dos tamaños de muestra. 41 Figura 6.18. Predicciones de la MPCT obtenidas con SVR frente valores reales con tamaños de muestra n=124 (izquierda) y n=93 (derecha). Figura 6.19. Residuos obtenidos con SVR y tamaños de muestra n=124 (izquierda) y n=93 derecha). Modelo MAE MSE RMSE MAPE VMP n=124 2.631038866 7.233376236 2.689493676 0.002187407 VMP n=93 2.692395600 7.389489383 2.718361526 0.002231905 Tabla 6.5. Medidas de calidad del ajuste SVR (n=124, n=93) Los valores obtenidos del límite superior de tolerancia 95/95 para n=124 y n=93 son respectivamente de 1259.7 K y 1258.1 K, respectivamente. Si comparamos dichos valores con el 1150 1200 1250 1150 1200 1250 y y ^ 1160 1180 1200 1220 1240 1260 1280 1160 1180 1200 1220 1240 1260 1280 y y ^