Parameter estimation of dynamic systems: application to civil engineering structures

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 Tesis Doctoral Ingeniería de Caminos, Canales y Puertos Parameter estimation of dynamic systems: application to civil engineering structures Autor: Javier Naranjo Pérez Directores: Javier Fernando Jiménez Alonsoy Andrés Sáez Pérez Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Tesis Doctoral Programa de Doctorado en Ingeniería Mecánica y de Organización Industrial Parameter estimation of dynamic systems: application to civil engineering structures Estimación de parámetros de sistemas dinámicos: aplicación a estructuras de ingeniería civil Autor: Javier Naranjo Pérez Tesis para la obtención del grado de Doctor por la Universidad de Sevilla Director: Javier Fernando Jiménez Alonso Profesor Ayudante Doctor de la Universidad Politécnica de Madrid Codirector: Andrés Sáez Pérez Catedrático de Universidad de la Universidad de Sevilla Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2020 Tesis Doctoral: Parameter estimation of dynamic systems: application to civil engineering structures Autor: Javier Naranjo Pérez Director: Javier Fernando Jiménez Alonso Codirector: Andrés Sáez Pérez El tribunal nombrado para juzgar la Tesis arriba indicada, compuesto por los siguientes doctores: Presidente: Vocales: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Prefacio Esta Tesis Doctoral se ha realizado en el Departamento de Mecánica de Medios Continuos y Análisis Estructural de la Universidad de Sevilla. Este trabajo ha sido supervisado por los profesores Javier Fernando Jiménez Alonso y Andrés Sáez Pérez. La tesis se basa en los seis artículos científicos indicados a continuación. Partes de los artículos se emplean en el resumen ampliado de la tesis. Los artículos incluidos como apéndices han sido preparados en colaboración con los autores. La lista de los seis artículos científicos es la siguientes: A: Jiménez-Alonso, J. F., Naranjo-Pérez, J., Pavic, A., Sáez, A. Maximum likelihood finite-element model updating of civil engineering structures using nature-inspired computational algorithms. Structural Engineering International. 2020. (Accepted). B: Naranjo-Pérez, J., Sáez, A., Jiménez-Alonso, J.F., Pachón, P., Compán, V. (2019). A hybrid UKF-MGA algorithm for finite element model updating of historical constructions. IABSE Symposium, Guimaraes 2019: Towards a Resilient Built Environment Risk and Asset Management – Report pp. 29-36. C: Naranjo-Pérez, J., Jiménez-Alonso, J. F., Pavic, A., Sáez, A. Finite-element model updating of civil engineering structures using a hybrid UKF-HS algorithm. Structure and Infraestructure Engineering. 2020. D: Naranjo-Pérez, J., Jiménez-Alonso, J. F., Sáez, A. Parameter identification of the dynamic Winkler soil-structure interaction model using a hybrid UKF-MHS algorithm. Advances in Structural Engineering. 2020. E: Infantes M, Naranjo-Pérez J, Jiménez-Alonso JF, Sáez A (2019). Determining the Best Pareto-solution in a Multi-Objective Approach for Model Updating. IABSE Symposium, Guimaraes 2019: Towards a Resilient Built Environment Risk and Asset Management – Report pp. 523-530. F: Naranjo-Pérez, J., Infantes, M., Jiménez-Alonso, J.F., A. Sáez. A collaborative machine learning-optimization algorithm to improve the finite element model updating of civil engineering structures. Submitted for publication. I A mi familia Índice Resumen VII Abstract IX Índice de Figuras XV Índice de Tablas XVII 1 Introducción 1 1.1 Antecedentes 1 1.2 Estado del arte 2 1.3 Motivación 5 1.4 Objetivos de la tesis 6 1.5 Publicaciones relacionadas con la tesis 6 1.5.1 Artículos en revistas internacionales 6 1.5.2 Artículos en congresos 7 1.6 Organización del texto 7 2 Estimación de parámetros a través de problema inverso 9 2.1 Enfoques del problema inverso 9 2.2 Filtro de Kalman 13 2.3 Estimador de máxima verosimilitud 14 2.4 Algoritmos de optimización 15 2.4.1 Algoritmos genéticos 18 2.4.2 Algoritmo Harmony Search 20 2.4.3 Algoritmo de enjambre de partículas 20 2.5 Comparación de algoritmos metaheurísticos para estimación de parámetros 23 2.5.1 Actualización de modelos de EF 23 2.5.2 Ejemplo de aplicación 24 2.5.3 Comparación de los tres algoritmos con un único objetivo 26 XI XII Índice 2.5.4 Comparación de los tres algoritmos multiobjetivo 27 2.5.5 Comparación de los planteamientos de único objetivo y multiobjetivo 29 2.6 Conclusiones 29 3 Algoritmo híbrido UKF-HS para estimación de parámetros 31 3.1 Fuentes de error 31 3.2 Unscented Kalman Filter 32 3.2.1 Actualización del modelo de EF de una cúpula de fábrica 34 3.3 Implementación del algoritmo híbrido UKF-MHS 37 3.3.1 Validación del algoritmo propuesto 38 Actualización utilizando el algoritmo UKF-HS 38 Actualización utilizando el algoritmo UKF-MHS 40 3.4 Aplicación para estimar los parámetros del modelo dinámico de Winkler 41 3.5 Conclusiones 45 4 Algoritmo colaborativo de inteligencia artificial para estimación de parámetros 47 4.1 El problema de la toma de decisiones en optimización multiobjetivo 47 4.2 Descripción del algoritmo colaborativo 50 4.2.1 Análisis de componentes principales 51 4.2.2 Redes neuronales artificiales 51 4.2.3 Algoritmo colaborativo para estimación de parámetros y toma de decisiones 53 4.3 Ejemplo de aplicación 54 4.3.1 Actualización del modelo de EF aplicando el algoritmo colaborativo 56 4.3.2 Comparación de resultados y tiempos de ejecución con un algoritmo convencional 59 4.4 Conclusiones 60 5 Conclusiones y trabajos futuros 63 5.1 Conclusiones 63 5.2 Líneas de investigación futuras 65 Bibliografía 67 Apéndice A Artículo A 77 Apéndice B Articulo B 105 Apéndice C Articulo C 115 Apéndice D Articulo D 151 Índice XIII Apéndice E Articulo E 179 Apéndice F Articulo F 189 Apéndice G Articulo G 191 Índice de Figuras 1.1 Diagrama general del proceso de estimación de parámetros 2 2.1 Comparación de enfoques: a) un único objetivo y b) multiobjetivo 16 2.2 Diagrama de flujo del algoritmo GA 19 2.3 Diagrama de flujo del algoritmo HS 21 2.4 Diagrama de flujo del algoritmo PS 22 2.5 Pasarela de laboratorio de la Universidad de Exeter y configuración del ensayo de vibración forzada. [Tomada del Artículo del Apéndice A] 25 2.6 Modelo de EF de la pasarela de laboratorio de la Universidad de Exeter, configuración del ensayo de vibración forzada y parámetros de actualización. [Tomada del Artículo del Apéndice C] 26 2.7 Análisis de sensibilidad de los parámetros del modelo de EF: matriz de sensibilidad representada en un diagrama de barras 27 2.8 Comparación de los valores residuales en el espacio funcional multiobjetivo, siendo f1la componente de la función multiobjetivo asociada a los residuos de las frecuencias naturales y f2la componente asociada a los residuos de los modos de vibración, para los tres algoritmos comparados, a) enfoque de un único objetivo y b) enfoque multiobjetivo 28 2.9 Comparación de los valores residuales en el espacio funcional multiobjetivo, siendo f1la componente de la función multiobjetivo asociada a los residuos de las frecuencias naturales y f2la componente asociada a los residuos de los modos de vibración, para los tres algoritmos comparados siguiendo ambos enfoques 30 3.1 a) Vista interior de la capilla de la Residencia Würzburg, b) estructura de la cúpula y c) modelo de EF de la estructura. [Tomada del Artículo del Apéndice B] 35 XV XVI Índice de Figuras 3.2 Ilustración de los tres modos de vibración de la cúpula 36 3.3 Diagrama de flujo del algoritmo híbrido UKF-MHS propuesto en este capítulo para resolver problemas de estimación de parámetros 39 3.4 Pasarela para la estimación de los parámetros del modelo dinámico de Winkler: a) vista frontal, b) plano de la sección y c) plano de alzado 43 3.5 Modelo de interacción suelo-estructura de cimentaciones superficiales y profundas y modelo de EF de la pasarela. [Tomada del Artículo del Apéndice D] 44 4.1 Representación de compromiso en una optimización de dos objetivos. Tomada del Artículo E 48 4.2 Comparación gráfica de los knee points obtenidos por los diferentes métodos 50 4.3 Arquitectura del MLP con una capa de entrada, jcapas ocultas y la capa de salida. La capa de entrada tiene mneuronas (parámetros de entrada) y la capa de salida tiene nneuronas (parámetros de salida) 52 4.4 Esquema del algoritmo colaborativo de inteligencia artificial y optimización multiobjetivo 55 4.5 Pasarela en estudio: a) ilustración de la pasarela, b) modelo de EF y c) configuración del ensayo de vibración ambiental 56 4.6 Representación gráfica de la matriz de sensibilidad para la elección de los parámetros a actualizar 57 4.7 Frente de Pareto del problema de estimación de parámetros de la pasarela: a) Original poco poblado, b) frente de Pareto procesado en el espacio de componentes principales y aproximado por la red neuronal artificial 59 4.8 Comparación entre las soluciones del algoritmo propuesto y el convencional para a) 150 iteraciones, b) 200 iteraciones, c) 250 iteraciones y d) 300 iteraciones 60 Índice de Tablas 2.1 Propiedades modales numéricas y experimentales de la pasarela del laboratorio de la Universidad de Exeter 25 2.2 Valores medios de los residuos en términos de los factores de peso para la actualización del modelo de EF basado en un único objetivo para los tres algoritmos metaheurísticos 28 3.1 Propiedades modales numéricas y experimentales de la cúpula de la capilla de la Residencia Würzburg antes de la actualización. fnum yfexp son las frecuencias naturales numéricas y experimentales, respectivamente 35 3.2 Propiedades modales numéricas y experimentales de la cúpula de la capilla de la Residencia Würzburg después de la actualización usando el algoritmo híbrido. fnum,act.son las frecuencias naturales numéricas actualizadas 36 3.3 Propiedades modales numéricas y experimentales de la cúpula de la capilla de la Residencia Würzburg después de la actualización usando el GA 37 3.4 Comparación de los tres algoritmos para la actualización del modelo de EF de la pasarela de laboratorio siguiendo el enfoque de un único objetivo. Iter es el número de iteraciones, Pob es el tamaño de la población, Eval es el número de evaluaciones de la función objetivo y fel valor de la función objetivo 40 3.5 Comparación de los tres algoritmos para la actualización del modelo de EF de la pasarela de laboratorio siguiendo el enfoque multiobjetivo. Iter es el número de iteraciones, Pob es el tamaño de la población, N.Pob es el tamaño de la población que se crea en cada iteración, Eval es el número de evaluaciones de la función objetivo, Dist es la distancia euclídea al origen y Sum es la suma de f1yf242 XVII XVIII Índice de Tablas 3.6 Propiedades modales numéricas y experimentales de la pasarela integral antes de la actualización 44 3.7 Propiedades modales numéricas y experimentales de pasarela integral después de la actualización 45 4.1 Propiedades modales numéricas y experimentales de la pasarela 57 4.2 Valores de la función multiobjetivo para el modelo mejor actualizado, distancia entre las mejores soluciones obtenidas por ambos algoritmos y tiempo total de simulación 61 1 Introducción 1.1 Antecedentes Los sistemas dinámicos se encuentran presentes en la práctica totalidad de ámbitos de la vida real, desde sistemas mecánicos y eléctricos hasta los sistemas biológicos. Conocer el comportamiento de dichos sistemas es crucial y los modelos, ya sean analíticos o numéricos, representan uno de los modos más frecuentes de simularlos. No obstante, los sistemas presentan propiedades desconocidas o inciertas que deben ser asumidas cuando se diseña el modelo. Por tanto, los modelos solo quedan bien definidos cuando los parámetros que lo caracterizan son calibrados, es decir, adoptan un valor que presenta la menor incertidumbre posible. Para definir los parámetros hay que resolver un problema de estimación de parámetros a partir del conocimiento de las entradas y salidas del sistema dinámico. Tal y como se ilustra en la Figura 1.1, la idea principal de estos problemas es minimizar el error (definido previamente) entre las salidas que el modelo ofrece, ˆy , y las medidas (salidas del sistema dinámico) registradas experimentalmente, z , mediante la modificación de los parámetros del modelo. Las fuentes de error, como el ruido en las medidas, pueden incluirse en el análisis. El constante avance que desarrolló la ingeniería informática en el último tercio del siglo pasado trajo consigo la aparición de nuevos algoritmos que permitían estimar los parámetros de modelos cada vez más complejos, siendo la mayor parte aplicados a problemas de ingeniería aeronáutica donde los parámetros de la aeronave se estiman a partir de los datos registrados durante el vuelo [61, 62]. En el ámbito de la ingeniería civil, la estimación de parámetros se ha aplicado de manera más reciente al caso particular de la actualización de modelos de elementos finitos (EF). Para modelizar las estructuras civiles se recurre a modelos numéricos, como los modelos de elementos finitos, que se definen a partir de parámetros mecánicos y geométricos que pueden presentar incertidumbre. Es en este aspecto donde el problema de estimación de parámetros permite definir los parámetros del modelo a partir del comportamiento realmente observado de la estructura. En el caso de estructuras se recurre habitualmente a los datos obtenidos en ensayos de vibración, como las propiedades modales, las funciones de respuesta en frecuencia o las historias temporales de aceleración [81, 102]. 1 8Capítulo 1. Introducción El Capítulo 5 cierra la tesis incluyendo las principales conclusiones y con una discusión de los trabajos futuros. El Apéndice A contiene al artículo “Jiménez-Alonso, J. F., Naranjo-Pérez, J., Pavic, A., Sáez, A. (2020). Maximum likelihood finite-element model updating of civil engineering structures using nature-inspired computational algorithms”. El Apéndice B contiene al artículo “J. Naranjo-Pérez, J.F. Jiménez-Alonso, P. Pachón, V. Compán, A. Sáez (2019). A hybrid UKF-MGA algorithm for finite element model updating of historical constructions”. El Apéndice C contiene al artículo “J. Naranjo-Pérez, J.F. Jiménez-Alonso, A. Pavic, A. Sáez (2020). Finite element model updating of civil engineering structures using a hybrid UKF-HS algorithm”. El Apéndice D contiene al artículo “J. Naranjo-Pérez, J.F. Jiménez-Alonso, A. Sáez (2020). Parameter identification of the dynamic Winkler soil-structure interaction model using a hybrid UKF-MHS algorithm”. El Apéndice E contiene al artículo “M. Infantes, J. Naranjo-Pérez, J.F. Jiménez-Alonso, A. Sáez (2019). Determining the best Pareto-solution in a multi-objective approach for model updating”. El Apéndice F contiene al artículo “J. Naranjo-Pérez, M. Infantes, J.F. Jiménez-Alonso, A. Sáez (2020). A collaborative machine learning-optimization algorithm to improve the finite element model updating of civil engineering structures”. El Apéndice G contiene al artículo “J. Naranjo-Pérez, J.F. Jiménez-Alonso, I. Díaz M., G. Quaranta, A. Sáez (2020). Motion-Based Design of Passive Damping Systems to Reduce Wind-Induced Vibrations of Stay Cables under Uncertainty Conditions. 2 Estimación de parámetros a través de problema inverso Cuando el problema de estimación de los parámetros de un modelo se lleva a cabo a partir de medidas o datos observados se le conoce como problema inverso. Este tipo de problemas puede plantearse según tres métodos: métodos basados en teoría de conjuntos difusos, el método Bayesiano (estimadores de intervalo) o usando estimadores puntuales. Los tres se resuelven de manera iterativa recurriendo a algoritmos de optimización. El presente capítulo pretende arrojar luz sobre los algoritmos que presentan mejor rendimiento para resolver problemas de ingeniería civil. En este sentido, el capítulo se estructura de la siguiente manera. En primer lugar se describen los enfoques del problema inverso, haciendo hincapié en los dos estimadores que se usan en esta tesis, el estimador de máxima verosimilitud y el filtro de Kalman. Seguidamente, se exponen los distintos tipos de algoritmos de optimización basados en métodos metaheurísticos, describiendo con mayor detalle los tres que se emplean en los siguientes capítulos. Finalmente se realiza un estudio comparativo de los tres algoritmos descritos y de los dos planteamientos que puede adoptar el problema, único objetivo o multiobjetivo, mediante la aplicación a un problema de actualización de modelo de EF de una pasarela de laboratorio. Partes de este capítulo han sido publicado en el Artículo del Apéndice A. 2.1 Enfoques del problema inverso Un problema inverso se define como aquel donde los parámetros de un modelo o sistema se derivan a partir de los datos observados o medidas. La variación de los valores de los parámetros provoca un cambio en la salida del sistema, por lo que de esta manera, el problema consiste en obtener aquellos valores de los parámetros que proporcionen valores de la salida del sistema más cercanos a los datos observados. 9 10 Capítulo 2. Estimación de parámetros a través de problema inverso Existen tres principales enfoques a la hora de resolver un problema de estimación de parámetros resolviendo un problema inverso: aplicando teoría de conjuntos difusos, aplicando métodos Bayesianos o aplicando estimadores puntuales. En el enfoque de la teoría difusa, la incertidumbre en el valor de los parámetros se trata de manera no probabilística, lo que es especialmente apropiado para modelizar incertidumbre sistemática (no aleatoria). La mayoría de los métodos no probabilísticos se basan en el análisis por intervalos, donde la incertidumbre de los parámetros se representa por un determinado rango de valores acotado. La teoría difusa fue propuesta en 1965 por Zadeh [ 116 ] con el objetivo de extender los operadores lógicos clásicos a un marco conceptual donde tratar la incertidumbre no probabilística. La idea de este planteamiento es dar un grado de pertenencia haciendo que el rango de valores de pertenencia de un elemento a un conjunto pueda variar en el intervalo [0,1] en lugar de seguir la lógica clásica, es decir, tomar el valor 0 o 1. Es un método muy adecuado para tratar la incertidumbre en situaciones donde, además de los límites de los intervalos, se dispone de valores de referencia de los parámetros a estimar. La teoría difusa ha sido aplicada a varios problemas de ingeniería estructural como, por ejemplo, dinámica estructural [ 1 ], análisis de fiabilidad [ 97 ], actualización de modelos de EF [59] o diseño estructural [77]. Los métodos Bayesianos, que permiten conocer la función de distribución de los parámetros, están basados en el teorema de Bayes, que puede formularse para la identificación de parámetros como [76]: P(θ|D) = P(D|θ)(P(θ) P(D)(2.1) donde θ representa el vector de parámetros a identificar y D es el vector que recoge los datos observados. Por tanto, cada término de la ecuación se define del siguiente modo: P(θ|D) es la función de densidad de probabilidad posterior de los parámetros condicionada por el valor de los datos observados, P(D|θ) es la función de verosimilitud de los datos condicionada por el valor de los parámetros, P(θ) es la función de densidad de probabilidad anterior, es decir, en ausencia de los datos observados y P(D) es un factor de normalización. La esperanza posterior condicionada de una función f(θ) se calcula como: E[f(θ)|D] = Rf(θ)P(D|θ)P(θ)dθ P(D).(2.2) Resolver la Ecuación 2.2 analíticamente puede no ser posible para sistemas complejos ya que la integral depende de la función de densidad de probabilidad posterior. Para salvar esta limitación, existen varias alternativas que permiten obtener una aproximación de la solución de la ecuación. La primera alternativa consiste en resolverla mediante evaluaciones numéricas, utilizando por ejemplo el método de aproximación asintótica, el cual considera una aproximación gaussiana de la función de densidad de probabilidad posterior. Sin embargo, si el tamaño de la función de densidad es grande este método no es preciso [11]. La segunda alternativa la constituyen las aproximaciones analíticas, como las aproximaciones de máxima verosimilitud, máximo a posteriori y de Laplace. El inconveniente de 2.1 Enfoques del problema inverso 11 estos métodos es su carácter local, por lo que la convergencia puede no ser a un mínimo global. Finalmente, las técnicas de muestreo, como los métodos de Monte Carlo basados en cadenas de Markov, son las más empleadas ya que dan como resultado una mejor aproximación para funciones complejas [ 9 ]. A continuación se esbozarán solo brevemente las técnicas de muestreo más comunes para la identificación de parámetros en ingeniería civil. Los métodos de Monte Carlo tratan de estimar la distribución de probabilidad mediante repeticiones aleatorias de muestreo [ 94 ]. El proceso comienza por crear de manera aleatoria muestras de los parámetros a partir de la función de distribución posterior. Esta nueva muestra es aceptada o rechazada al azar y se repite iterativamente hasta que el número de muestras ha alcanzado un valor prefijado. El valor medio de la función de distribución se calcula a partir del conjunto de muestras aceptadas. El método de Monte Carlo basado en cadenas de Markov es una técnica en la cual la función de distribución de probabilidad se identifica mediante una cadena de Markov. En este sentido, se asume que la evolución del sistema se caracteriza a través de un proceso estocástico de variables aleatorias (cuyo dominio forman el espacio estado), Ki , donde cada variable es un vector que contiene los parámetros a identificar. Como las variables aleatorias forman una cadena de Markov, la probabilidad de que el sistema esté en un estado Kn en un tiempo n depende únicamente de que el sistema estuviera en un estado Xn−1 en el tiempo n−1 . El estado actual se acepta o rechaza aplicando el algoritmo Metropolis. Para llevar a cabo una simulación de Monte Carlo basada en cadenas de Markov existen diferentes algoritmos, como el algoritmo Metropolis-Hastings, el algoritmo Slice Sampling o el algoritmo Metropolis-Hastings adaptativo [76]. Los métodos Bayesianos requieren, en general, de un tiempo de computación excesivamente alto para calcular la función de distribución de los parámetros. Además, el cálculo de la media de los parámetros no es muy preciso, por lo que debe calcularse previamente usando estimadores para asegurar que la función de distribución sea adecuada al problema. Este segundo enfoque mediante estimadores puntuales es en el que se basa la presente tesis. Uno de los aspectos fundamentales de los estimadores puntuales es el método para definir el error que se considere [ 92 ]. El primer método es el de error de ecuación (equation error), donde el error es la diferencia entre la medida y la salida obtenida por la ecuación del modelo con los valores de los parámetros estimados. La función de coste es lineal en los parámetros. El segundo método es el de error en la salida (output error), donde el error es la diferencia entre la salida del modelo y la medida del sistema real. Para la aplicación de este método, la entrada al modelo y al sistema debe ser la misma. Finalmente, el método de error de parámetros (parameter error) se define como la diferencia entre los valores estimados de los parámetros y sus valores reales. De este modo, el error de parámetro sólo puede ser obtenido si el valor real es conocido, lo que limita su aplicación para problemas reales. Entre los estimadores más comunes se encuentran el estimador de mínimos cuadrados, el estimador de máxima verosimilitud y los filtros de Kalman. 12 Capítulo 2. Estimación de parámetros a través de problema inverso El estimador de mínimos cuadrados se basa en un modelo de ecuación de la medida que toma la siguiente forma: z=y+v=Xθ θ θ+v(2.3) siendo z el vector que recoge las medidas, y las salidas del modelo, v un vector que representa el ruido/error en las medidas, que se asume como gaussiano y con media cero, θ θ θ el vector de parámetros y X una matriz que contiene las variables explicativas del problema, las cuales representan el comportamiento del modelo y son perfectamente conocidas. El objetivo del estimador de mínimos cuadrados es encontrar el valor de los parámetros que minimice la función de coste, J , definida como la suma de los errores al cuadrado entre las medidas y las predicciones del modelo: J∼ = M ∑ j=1 v2 j= (z−Xθ θ θ)T(z−Xθ θ θ)(2.4) donde M es el número de medidas. La minimización con respecto a θ θ θ de la función J conduce a una expresión analítica de la estimación de mínimos cuadrados de θ θ θ : ˆ θ θ θ= (XTX)−1XTz . Por consiguiente, como el vector de medidas z y la matriz H son conocidos, las estimaciones se obtienen de manera directa. El estimador de mínimos cuadrados presenta tres propiedades principales. La primera es que la estimación ˆ θ θ θ es una función lineal del vector de medidas z . La segunda es que el error en el estimador, θ θ θ−ˆ θ θ θ , es una función lineal de los errores en las medidas v . La tercera es que si la esperanza matemática del error en las medidas v es cero, la estimación ˆ θ θ θ es no sesgada ya que la esperanza del error en la estimación , E[θ θ θ−ˆ θ θ θ] , también es nula. Por tanto, si z=y , ˆ θ θ θ es una cantidad determinista igual a θ θ θ . No obstante, en presencia de ruido en la medida la estimación ˆ θ θ θ es sesgada y el estimador de mínimos cuadrados es empleado como punto de partida para otros métodos de estimación como el estimador de mínimos cuadrados generalizado o el método de error en la salida basado en el estimador de máxima verosimilitud. El método de mínimos cuadrados generalizado, también llamado mínimos cuadrados ponderado, consiste en incluir una matriz de peso simétrica y definida positiva en la definición de la función J que sirve para controlar la influencia de medidas específicas en las estimaciones de los parámetros. El estimador de mínimos cuadrados puede aplicarse para sistemas no lineales asumiendo que las salidas del modelo (también no lineal) adopta la forma y=h(θ θ θ) , donde h es una función no lineal conocida. Para este caso, la estimación ˆ θ θ θ se obtiene a través de un proceso iterativo. El estimador de máxima verosimilitud es el caso general del estimador de mínimos cuadrados ya que la estimación es idéntica cuando se trata de un sistema lineal con ruido gaussiano. La idea principal es definir una función de las medidas y de los parámetros, llamada función de verosimilitud, y maximizarla con respecto a los parámetros. La función de verosimilitud representa la densidad de probabilidad de las medidas dados los parámetros. El estimador de máxima verosimilitud puede aplicarse para los métodos de error de ecuación (se considera ruido únicamente en el proceso de estimación) y de error en la salida (se considera ruido únicamente en las medidas). No obstante, para poder considerar ambos ruidos en el problema es necesario recurrir a otras técnicas. El filtro de Kalman [ 53 ] es un método de estimación de parámetros que permite tener 2.2 Filtro de Kalman 13 en cuenta ruido en el proceso de estimación y en las medidas. En el caso de sistemas no lineales, la estimación del filtro de Kalman no es óptima y hay que recurrir a variaciones del filtro de Kalman, como el filtro de Kalman extendido (EKF de Extended Kalman Filter) [ 49 ]. No obstante, su aplicación para estimación de parámetros revela que los resultados son muy sensibles a los valores iniciales de las matrices de covarianza del ruido en las medidas [92]. La otra alternativa consiste en aplicar los denominados métodos de error de filtro, los cuales implementan el filtro de Kalman dentro de la estructura del método de error en la salida basado en el estimador de máxima verosimilitud, de modo que permite manejar ruidos en las medidas y en el proceso de estimación. Esta idea justifica el planteamiento del algoritmo híbrido que se propone en esta tesis. En concreto, se emplea una variante del filtro de Kalman para sistemas no lineales y el estimador de máxima verosimilitud por lo que serán explicados con mayor detalle en los siguientes apartados. 2.2 Filtro de Kalman El filtro de Kalman es un algoritmo que permite estimar las medias y las covarianzas de los parámetros llevando a cabo un proceso iterativo en el que se reduce la covarianza del error en el proceso de estimación hasta converger en una solución en los parámetros no sesgada. La representación matemática de un problema de identificación de parámetros se basa en el modelo de espacio de estados donde la Ecuación 2.5 representa el proceso de transición estacionario en la estimación y la Ecuación 2.6 representa la ecuación del modelo: θk=θk−1+wk−1(2.5) yk=h(θ)+vk(2.6) siendo h una función de modelización, w el ruido en el proceso de estimación y v el ruido en la modelización. Ambos términos de ruido se consideran ruido blanco gaussiano con media cero y matrices de covarianza Q y R , respectivamente. La matriz R presenta dos contribuciones, el ruido en la modelización y el ruido en las medidas. El filtro de Kalman se basa en un proceso iterativo en el que cada iteración ejecuta dos pasos: predicción y corrección. El primer paso se basa en las estimaciones del paso anterior (Ecuación 2.5) para predecir las estimaciones del paso actual. El segundo paso tiene en cuenta las medidas para corregir las estimaciones de manera que se reduzca la covarianza del error. Para problemas no lineales se desarrolló el EKF que tiene mejores prestaciones para sistemas ligeramente no lineales. El EKF requiere la linelización de la Ecuación 2.6 (la Ecuación 2.5 ya lo es), por lo que es necesario calcular el Jacobiano de la función de modelización h . Sin embargo, cuando el sistema presenta una fuerte no linealidad, el rendimiento del EKF disminuye. En el caso de que la función h se defina mediante un modelo numérico no hay una expresión analítica para calcular el Jacobiano y hay que aproximarlo numéricamente. 14 Capítulo 2. Estimación de parámetros a través de problema inverso Más recientemente se desarrolló el Unscented Kalman Filter (UKF) [ 52 ]. Se trata de un algoritmo que no requiere ninguna operación diferencial ya que no hay que linealizar ninguna función [ 105 , 84 ]. La premisa del UKF es que un conjunto de puntos pueden capturar y propagar la media y covarianza de los parámetros a través de la función h para obtener la media y covarianza posterior, que se usan en la siguiente iteración para predecir la media y covarianza anterior. Varios estudios han demostrado que el algoritmo UKF supera al EKF. Mariani y Ghisi los compararon para la calibración de la ley constitutiva de un sistema dinámico de un grado de libertad [ 73 ]. Más recientemente, Nguyen y Nestorović compararon estos dos algoritmos para la estimación de parámetros geomecánicos en modelos del suelo para problemas de túneles [ 85 ]. Por tanto, en esta tesis se considera el algoritmo UKF para su implementación en problemas de estimación de parámetros. 2.3 Estimador de máxima verosimilitud El estimador de máxima verosimilitud se basa en la definición de una función, llamada función de verosimilitud, que depende de las medidas experimentales del sistema y de los parámetros del modelo. El objetivo de este estimador es maximizar dicha función, es decir, buscar aquellos parámetros que hacen que de manera más probable la función conduzca a las medidas. Como se dijo anteriormente, la función de verosimilitud es la densidad de probabilidad de las medidas y se calcula para las medidas estadísticamente independientes donde el error se asume que está normalmente distribuido. La densidad de probabilidad conjunta se halla como el producto de las probabilidades de las medidas individuales. La función de verosimilitud puede escribirse como el logaritmo de la función de densidad de probabilidad ya que esta función y su logaritmo siguen teniendo el máximo en el mismo argumento. Esta función depende del error entre las medidas y las salidas del modelo, asumido como normalmente distribuido, y está definida de manera que la minimización del error conlleve la maximización de la probabilidad, es decir, de la función de verosimilitud. Por tanto, el problema de maximización se reduce a un problema de minimización. A diferencia del estimador de mínimos cuadrados, la estimación no es directa y es necesario realizar un proceso iterativo. El proceso comienza por realizar una primera estimación de los parámetros y calcular a partir de ellos la salida del modelo. Teniendo en cuenta las medidas y la salida del modelo se evalúa el error y la función objetivo. En la siguiente iteración, los parámetros se actualizan a partir de los parámetros del paso anterior y del valor de la función objetivo. Este proceso continua hasta que se logra la convergencia. Esta metodología puede ser aplicada a sistemas no lineales. Además, el estimador de máxima verosimilitud presenta varias ventajas matemáticas. En primer lugar, la estimación es consistente. Segundo, en diversas condiciones el estimador tiene una distribución asintótica normal. Finalmente, el estimador de máxima verosimilitud es eficiente [92]. Estos motivos justifican su elección para la resolución de problemas de estimación de parámetros. El proceso iterativo puede definirse como un problema de optimización, donde el objetivo es encontrar los valores de los parámetros que optimicen la función, expresado 2.4 Algoritmos de optimización 15 matemáticamente como: min f(θ θ θ)(2.7) sujeto a: θ θ θl≤θ θ θ≤θ θ θu(2.8) donde f es la función objetivo a minimizar, θ θ θ es el vector de parámetros a estimar y θ θ θl y θ θ θuson los límites inferiores y superiores de los parámetros. Para la ejecución del proceso iterativo se recurre a algoritmos de optimización que solucionen el problema de manera automática [ 86 ]. En el siguiente apartado se describirán los principales algoritmos de optimización y su aplicación a los distintos enfoques del problema de estimación de parámetros. 2.4 Algoritmos de optimización Cuando se elige el algoritmo de optimización hay que considerar el carácter local o global de la solución que obtiene. Un óptimo local es una solución óptima en un entorno próximo, mientras que un óptimo global es la solución óptima de todas las posibles soluciones. Por consiguiente, según sea la función objetivo (o de coste), puede resultar más apropiado elegir un tipo de algoritmo u otro. Los algoritmos locales convergen a óptimos locales ya que se basan en el gradiente y son los más adecuados cuando la función objetivo tiene forma cóncava/convexa para las que el óptimo local es la solución óptima del problema. Por el contrario, los algoritmos globales permiten hallar óptimos globales y su aplicación es adecuada cuando la función objetivo presenta múltiples mínimos o máximos locales [86]. A la hora de definir la función objetivo que hay que optimizar hay que tener en cuenta los errores residuales entre las salidas del modelo y las medidas. La función objetivo se define a partir de uno o varios residuos según el número de criterios de optimización, asumiendo en este último caso unos factores de peso que representan la influencia de cada residuo (Figura 2.1). A este tipo de planteamiento se le conoce como optimización de un único objetivo. El objetivo de este tipo de planteamiento es sencillo ya que se trata de hallar el valor de los parámetros que encuentre el óptimo de dicha función. En múltiples casos, no es adecuado o posible unir los residuos en una misma función objetivo. En estos casos se recurre a la optimización multiobjetivo, donde cada función objetivo está definida por un residuo. Para este planteamiento se recurre al concepto de óptimo de Pareto formulado en el siglo XIX. El objetivo de este enfoque es encontrar una solución Pareto-óptima cuando no existe otra solución que mejore un objetivo sin sacrificar al menos uno de los otros. Al conjunto de soluciones que cumplen este criterio se les denomina frente de Pareto o soluciones no dominadas. El resultado de este tipo de problemas es una curva formada por un conjunto de puntos donde cada uno de ellos es una posible solución óptima al problema. Por tanto, para finalizar el problema es necesario elegir el mejor punto mediante técnicas de toma de decisiones (Figura 2.1). Es preciso remarcar en este punto que el carácter local o global del algoritmo no está directamente relacionado con el enfoque de un único objetivo o multiobjetivo, ya que un 16 Capítulo 2. Estimación de parámetros a través de problema inverso a) Enfoque único objetivo Función objetivo: suma ponderada residuos Optimización Análisis sensibilidad pesos Elección mejor solución Modificar pesos b) Enfoque multiobjetivo Función multiobjetivo Optimización Frente de Pareto Toma de decisiones Figura 2.1 Comparación de enfoques: a) un único objetivo y b) multiobjetivo. único objetivo puede presentar un único óptimo local (también global) o varios óptimos locales. En el caso de que la relación entre los parámetros y la función objetivo sea no lineal o el espacio de solución sea muy amplio, las técnicas metaheurísticas representan la mejor opción para resolver el problema frente a los algoritmos deterministas convencionales [ 74 ]. Son métodos que enfocan el problema desde un punto de vista probabilístico, como el método de Monte Carlo. Los métodos metaheurísticos se basan en criterios que permiten discernir qué alternativa es la más efectiva para logar un objetivo, es decir, qué candidato a solución es el más efectivo para seguir considerándolo en la generación de nuevos candidatos. Los algoritmos metaheurísticos se pueden clasificar según su estrategia de búsqueda por puntos singulares o basada en poblaciones. Los primeros generan una única solución en cada iteración. Los segundos mantienen un conjunto de candidatos que van evolucionando a raíz de la exploración en el espacio de búsqueda. Las principales ventajas de este tipo de algoritmos es su independencia del valor inicial de los candidatos (que son estimados aleatoriamente dentro de un rango definido a partir del conocimiento físico del problema) y su habilidad para encontrar el óptimo global. No obstante, presentan un claro inconveniente: el alto tiempo de simulación requerido para resolver el problema de optimización. Aunque la mayoría de los algoritmos metaheurísticos son algoritmos globales, existen también algoritmos locales como es el caso del algoritmo Hill Climbing (escalada simple o escalada de colinas) [ 39 ] que comienza por una solución arbitraria y busca encontrar una mejor solución variando incrementalmente un único elemento (parámetro) de la solución. El proceso continua hasta que no se encuentra ninguna solución mejor. Hay una amplia variedad de algoritmos metaheurísticos globales que han sido usados para aplicaciones de ingeniería civil [ 107 , 114 ]. Moradi et al. [ 80 ] exploraron la eficiencia del algoritmo de colonia de abejas para la resolución de un problema de actualización 2.4 Algoritmos de optimización 17 del modelo de EF de un sistema de tuberías. Las redes neuronales artificiales también han demostrado su eficiencia para resolver problemas de optimización no lineales [ 63 , 88 , 42 , 40 , 68 ]. El algoritmo de colonia de hormigas ha sido igualmente usado para problemas de optimización [ 25 ]. Los algoritmos genéticos han sido ampliamente usados para aplicaciones de ingeniería, como es el caso de problemas de actualización de modelos o para detección de daño [ 51 , 67 , 50 ]. El algoritmo Harmony Search, ideado en 2001 [ 38 ], sigue un proceso muy similar a los algoritmos genéticos aunque más simplificado ya que no se basa en poblaciones. Aunque existan múltiples aplicaciones para estimación de parámetros en ingeniería estructural, su uso en problemas de actualización de modelos de EF es escaso [ 71 ]. El algoritmo de enjambre de partículas, a diferencia de los algoritmos genéticos, mantiene la misma población la cual va evolucionando durante el proceso iterativo. En [ 96 ] se aplica este algoritmo para la actualización del modelo de EF de vigas con el objetivo de identificar el daño. El algoritmo de evolución diferencial puede considerarse como una versión más desarrollada de los algoritmos genéticos donde las ecuaciones de actualización de la población son explícitas. Quaranta et al. [ 91 ] emplearon este algoritmo para el diseño óptimo (en términos económicos) de vigas de hormigón pretensadas. El algoritmo Cuckoo search es uno de los más recientes. Fue desarrollado por Yang y Deb en 2009 [ 115 ] y está basado en el parasitismo de algunas especies de cuco que dejan sus huevos en nidos de otras especies. Gandomi et al. [ 33 ] emplearon este algoritmo para el diseño óptimo de estructuras en celosía. El algoritmo colliding bodies (cuerpos en colisión) se basa en las leyes de conservación del momento y la energía de una colisión unidimensional [ 54 ]. Cada cuerpo posee una masa, posición y velocidad y la colisión con otros cuerpos actualiza la posición de acuerdo a unas leyes de colisión. El algoritmo de búsqueda gravitacional está basado en las leyes de la gravitación de Newton. Este algoritmo ha sido aplicado a una amplia variedad de problemas de optimización, como el diseño óptimo de la forma de presas de gravedad de hormigón [ 57 ] o el diseño basado en la fiabilidad de structuras de hormigón reforzado considerando interacción suelo-estructura [58]. Los resultados obtenidos por el algoritmo de enjambre de partículas y el algoritmo genético fueron comparados por Perera et al. [ 89 ] para la actualización del modelo de EF de un pórtico de una planta de hormigón armado. Igualmente, Levin y Lieven compararon los resultados entre los algoritmos genéticos y el algoritmo simulated annealing para la actualización de un modelo de EF considerando datos simulados [ 64 ]. Estos tres algoritmos, enjambre de partículas, genéticos y simulated annealing, fueron comparados en un mismo estudio por Marwala, concluyendo que el algoritmo de enjambre de partículas era superior a los otros dos [75]. Para vencer el problema del alto tiempo de simulación requerido, la tendencia actual es usar algoritmos híbridos aprovechando la sinergia del algoritmo resultante. Algunos ejemplos de esta tendencia es la combinación de algoritmo genético y simulated annealing propuesto por Feng et al. [ 27 ], la combinación del algoritmo simulated annealing con el UKF [ 4 ] o la combinación del algoritmo de enjambre de partículas con un surrogate model [90]. Dentro del amplio abanico de algoritmos de optimización existentes, los algoritmos genéticos constituyen la opción más empleada para resolver problemas de actualización de modelos de EF en problemas prácticos de ingeniería civil. El algoritmo de enjambre de 24 Capítulo 2. Estimación de parámetros a través de problema inverso rf j=fnum,j−fexp,j fexp,j (2.10) (rm j)2=(1−pMACj)2 MACj (2.11) donde fnum,j y fexp,j son la frecuencia natural numérica y experimental j , respectivamente, y MACj es un parámetro para medir la correlación entre dos modos de vibración que se calcula como [3]: MACj=(φnum,j·φexp,j)2 (φT num,j·φnum,j)·(φT exp,j·φexp,j)(2.12) siendo φnum,j y φexp,j los vectores con las coordenadas modales del modo de vibración j numérico y experimental, respectivamente. En segundo lugar, el enfoque multiobjetivo se formula como [51]: min f(θ θ θ) = min f1=1 2"n ∑ j (rf j)2#1/2 ,f2=1 2"n ∑ j (rm j)2#1/2  T (2.13) siendo f1 y f2 las dos componentes de la función multiobjetivo, estando la primera de ellas definida en términos de los residuos de las frecuencias naturales y la segunda a partir de los residuos de los modos de vibración. La solución de este planteamiento es el frente de Pareto del problema, donde cada punto representa una posible solución, es decir, un vector de parámetros del modelo de EF que permite reducir las diferencias entre las propiedades modales numéricas y experimentales. 2.5.2 Ejemplo de aplicación La estructura de laboratorio consiste en una pasarela porticada de un solo vano de 15 m de longitud situada en el laboratorio de la sección de Vibraciones de la Universidad de Exeter [ 45 ]. Dos vigas de acero separadas 2.5m en la dirección transversal a la pasarela la recorren longitudinalmente. Estas vigas se encuentran sujetas por diafragmas de acero colocados cada 1.25 m cuya sección transversal es una placa rectangular de 200 ×12 mm . Toda la pasarela se apoya en cuatro pilares de acero situados en los extremos de las dos vigas laterales que están anclados al suelo. El tablero de la pasarela está formado por un panel compuesto atornillado a las vigas laterales y a los diafragmas. En la Figura 2.5 se ilustra la pasarela en estudio. El modelo inicial de EF de la estructura se construye en el software ANSYS [ 78 ]. Para las dos vigas laterales, los diafragmas y los paneles compuestos se consideran elementos tipo placa de cuatro nodos con seis grados de libertad por nodo (elemento Shell181). Los tornillos que unen los paneles a los elementos metálicos se modelizan mediante elementos tipo viga de dos nodos y seis grados de libertad por nodo (elemento Beam188). Finalmente, los apoyos se modelizan a través de elementos tipo muelle en la dirección longitudinal 2.5 Comparación de algoritmos metaheurísticos para estimación de parámetros 25 Figura 2.5 Pasarela de laboratorio de la Universidad de Exeter y configuración del ensayo de vibración forzada. [Tomada del Artículo del Apéndice A]. Tabla 2.1 Propiedades modales numéricas y experimentales de la pasarela del laboratorio de la Universidad de Exeter. Modo fnum [Hz]fexp [Hz]∆f[%]MAC [-] 1 3.638 3.810 -4.509 0.999 2 5.329 5.144 3.600 0.994 3 10.185 8.485 20.033 0.990 4 11.310 12.366 -8.540 0.877 5 17.364 18.605 -6.670 0.985 6 20.238 20.459 -1.080 0.993 7 21.105 22.980 -8.159 0.910 y transversal (elemento COMBIN14). El modelo de EF (ver Figura 2.6) consiste en una malla de 31903 elementos. El análisis modal numérico realizado al modelo de EF de la pasarela da como resultado las siete frecuencias numéricas y sus correspondientes modos de vibración que se muestran en la Tabla 2.1. Las propiedades modales experimentales se obtienen a partir del análisis modal experimental de las señales de aceleración registradas durante un ensayo de vibración forzada, llevado a cabo en la pasarela cuya configuración se representa en la Figura 2.6 (consúltese el Artículo del Apéndice A para mayor información). En la Tabla 2.1 se indican los valores de las frecuencias naturales experimentales. Para analizar la correlación entre las propiedades modales numéricas y experimentales se indican en la Tabla 2.1 los valores de las diferencias relativas entre frecuencias y del parámetro MAC. A pesar de los buenos valores del MAC, las diferencias relativas de cuatro de las siete frecuencias son todavía grandes [ 117 ]. Por tanto, para mejorar el modelo de EF de manera que se ajuste mejor al comportamiento real de la pasarela es necesario actualizarlo. Para elegir los parámetros a actualizar, que representan las variables de diseño del problema, se realiza un análisis de sensibilidad de varios parámetros para conocer la influencia de cada uno de ellos en una característica de referencia de la estructura. Siguiendo la formulación de Fox y Kapoor se adopta la energía de deformación modal como esta característica de referencia [ 29 ]. De esta manera, los parámetros seleccionados son aquellos que tienen una mayor influencia en la energía de deformación modal para cada modo de 26 Capítulo 2. Estimación de parámetros a través de problema inverso Figura 2.6 Modelo de EF de la pasarela de laboratorio de la Universidad de Exeter, configuración del ensayo de vibración forzada y parámetros de actualización. [Tomada del Artículo del Apéndice C]. vibración. En la Figura 2.7 se muestra la matriz de sensibilidad en forma de diagrama de barras. De acuerdo a dicha matriz, las parámetros físicos más relevantes son los diez siguientes: el módulo de Young del acero de las vigas longitudinales en 6 secciones ( θ1−θ6 ), el módulo de Young del poliuretano del panel compuesto ( θ7 ), el módulo de Young del acero de los tornillos ( θ8 ), la rigidez longitudinal equivalente de los apoyos ( θ9 ) y la rigidez transversal equivalente de cada apoyo ( θ10 ). En la Figura 2.6 se ilustran los parámetros de actualización. En el Artículo del Apéndice A se puede encontrar los límites inferiores y superiores impuestos a cada parámetro para garantizar un significado físico adecuado de los parámetros. 2.5.3 Comparación de los tres algoritmos con un único objetivo En primer lugar se soluciona el problema de actualización del modelo de EF siguiendo el planteamiento de un único objetivo mediante la implementación de los tres algoritmos descritos anteriormente, GA, HS y PS. Se consideran dos criterios de comparación para validar el rendimiento de los algoritmos: el coste computacional (tiempo de simulación) y la exactitud del ajuste. Para todos los algoritmos las variables siguientes se consideran comunes: tamaño de la población (100), número máximo de iteraciones (50) y tolerancia de la función objetivo ( 10−5 ). Los parámetros intrínsecos de cada algoritmo se indican con detalle en el Artículo del Apéndice A. Se consideran once casos distintos en función de los valores de los pesos asociados a las frecuencias naturales, ∑wf j , y a los modos de vibración, ∑wm j , para considerar su impacto en los residuos [ 104 ]. Además, cada caso se ejecuta diez veces para tener en cuenta la aleatoriedad con la que la población inicial se genera. En la Tabla 2.2 se muestra la solución al problema para los once casos y los tres 2.5 Comparación de algoritmos metaheurísticos para estimación de parámetros 27 Figura 2.7 Análisis de sensibilidad de los parámetros del modelo de EF: matriz de sensibilidad representada en un diagrama de barras. algoritmos metaheurísticos. Además, en la Figura 2.8a se muestran estas soluciones en el espacio funcional multiobjetivo. El coste computacional de cada algoritmo se calcula como la media de los tiempo de simulación de los once casos, dando como resultado los valores siguientes. Para el GA: t=74072 s, para el HS, t=34380 s y para el PS, t=43857 s. De acuerdo con los valores obtenidos y analizando la exactitud en la solución y el tiempo de simulación se pueden extraer dos conclusiones principales. La primera es que la precisión de los tres algoritmos es similar. La segunda es que el algoritmo HS es más eficiente, ya que permite reducir el tiempo de computación requerido para ejecutar la actualización del modelo de EF sin comprometer la precisión en el ajuste. 2.5.4 Comparación de los tres algoritmos multiobjetivo Los dos mismos criterios se establecen para la comparación de la solución dada por los tres algoritmos cuando se aplican a la actualización del modelo de EF siguiendo el enfoque multiobjetivo. Los parámetros relativos a los tres algoritmos se asumen iguales al caso anterior. De igual manera, se ejecuta diez veces cada algoritmo para tener en cuenta la generación aleatoria de la población inicial. En la Figura 2.8b se representa el frente de Pareto obtenido por los tres algoritmos. 28 Capítulo 2. Estimación de parámetros a través de problema inverso Tabla 2.2 Valores medios de los residuos en términos de los factores de peso para la actualización del modelo de EF basado en un único objetivo para los tres algoritmos metaheurísticos. Caso ∑wf j∑wm j GA (·10−3) PS (·10−3) HS (·10−3) ∑(rf j)2∑(rm j)2∑(rf j)2∑(rm j)2∑(rf j)2∑(rm j)2 1 1 0 3,63 2,69 3,56 2,68 3,55 2,69 2 0,9 0,1 3,57 2,66 3,56 2,68 3,55 2,69 3 0,8 0,2 3,59 2,62 3,58 2,51 3,58 2,66 4 0,7 0,3 3,58 2,53 3,57 2,52 3,61 2,52 5 0,6 0,4 3,57 2,56 3,57 2,5 3,62 2,5 6 0,5 0,5 3,58 2,51 3,58 2,51 3,58 2,57 7 0,4 0,6 3,58 2,52 3,58 2,49 3,62 2,51 8 0,3 0,7 3,59 2,48 3,58 2,49 3,63 2,51 9 0,2 0,8 3,58 2,51 3,58 2,49 3,62 2,49 10 0,1 0,9 3,58 2,53 3,59 2,47 3,7 2,47 11 0 1 3,59 2,5 3,59 2,49 3,64 2,46 a) 45 6 7 ·10−3 2.2 2.3 2.4 2.5 2.6 2.7·10−3 f1 f2 GA PS HS b) 4567 ·10−3 2.2 2.3 2.4 2.5 2.6 2.7·10−3 f1 f2 MGA MPS MHS Figura 2.8 Comparación de los valores residuales en el espacio funcional multiobjetivo, siendo f1 la componente de la función multiobjetivo asociada a los residuos de las frecuencias naturales y f2 la componente asociada a los residuos de los modos de vibración, para los tres algoritmos comparados, a) enfoque de un único objetivo y b) enfoque multiobjetivo. Los tiempos que ha tardado cada simulación son: t=307057 s para el GA, t=121676 s para el MHS y t=240663 s para el PS. De acuerdo a los resultados para las diez simulaciones pueden sacarse tres conclusiones. La primera es que la precisión de los algoritmos PS y MHS es mayor que la del GA ya que 2.6 Conclusiones 29 el frente de Pareto de ambos permite reducir el valor total de los residuos. La segunda es que la precisión entre el PS y el MHS es similar. La última es que el tiempo de simulación del MHS es significativamente menor al de los otros dos algoritmos. Por tanto, nuevamente el MHS es el algoritmo más eficiente para resolver la actualización del modelo de EF siguiendo un planteamiento multiobjetivo. 2.5.5 Comparación de los planteamientos de único objetivo y multiobjetivo La comparación entre ambos enfoques se efectúa en base a tres criterios: la capacidad de rastreo, el coste computacional y el ajuste de la solución. En la Figura 2.9 se representa la solución obtenida por los tres algoritmos para los dos enfoques. Puede apreciarse que a pesar de que las soluciones obtenidas siguiendo el planteamiento de un único objetivo pertenecen al frente de Pareto, se centran en una sola rama de la curva. De este modo, la capacidad de rastreo del enfoque de un único objetivo es menor. Por otro lado, el conjunto de soluciones del enfoque multiobjetivo presentan un mejor ajuste ya que la distancia al origen de coordenadas es menor que la de los puntos del enfoque de un objetivo, que se encuentran más lejos del punto de la curva que mejor equilibra ambos residuos (llamado knee point) [23]. Finalmente, para comparar los tiempos de simulación es preciso remarcar que bajo el enfoque de un único objetivo se considera la suma de los once casos ya que el estudio de sensibilidad de los pesos es fundamental y necesario para cada problema. De este modo, se obtiene que el tiempo de simulación del algoritmo MHS es menor que el del algoritmo HS (121676 s frente a 378180 s). Por todos estos motivos se puede justificar que el enfoque multiobjetivo del algoritmo HS (MHS) representa la mejor opción para llevar a cabo la actualización del modelo de EF de estructuras civiles. 2.6 Conclusiones En este capítulo se han descrito los dos enfoques más generales que pueden seguirse para estimar parámetros resolviendo un problema inverso. El estimador de máxima verosimilitud es el empleado en esta tesis debido a su buena relación coste computacional/precisión en la solución. La estimación de los parámetros se realiza siguiendo un proceso iterativo que puede plantearse de dos maneras: considerando una única función objetivo o una función multiobjetivo. Existen diversos algoritmos de optimización para resolver el proceso iterativo, siendo los métodos metaheurísticos los más adecuados, ya que no tienen dependencia del valor inicial y logran encontrar el óptimo global. En este capítulo se ha llevado a cabo un análisis comparativo para estimar qué algoritmo es el que mejor rendimiento ofrece. Además, también se han analizado y comparado los dos enfoques del problema, el de un único objetivo y el de multiobjetivo. Como resultado de esta comparación se han obtenido dos conclusiones: la primera es que el algoritmo HS es el más eficiente en los dos enfoques y la segunda es que el enfoque multiobjetivo permite obtener una mejor solución al problema. 30 Capítulo 2. Estimación de parámetros a través de problema inverso 3.544.5 5 5.5 6 6.577.5 ·10−3 2.2 2.3 2.4 2.5 2.6 2.7·10−3 f1 f2 GA PS HS MGA MPS MHS Figura 2.9 Comparación de los valores residuales en el espacio funcional multiobjetivo, siendo f1 la componente de la función multiobjetivo asociada a los residuos de las frecuencias naturales y f2 la componente asociada a los residuos de los modos de vibración, para los tres algoritmos comparados siguiendo ambos enfoques. Estas dos conclusiones son la base para el desarrollo teórico de los algoritmos que se proponen en esta tesis en los Capítulos 3 y 4. 3 Algoritmo híbrido UKF-HS para estimación de parámetros En este capítulo se aborda la consideración de los errores en el proceso de estimación y en las medidas que aparecen en el problema de actualización de parámetros de modelos. En primer lugar se esbozan las distintas fuentes de error que aparecen al construir un modelo numérico y estimar sus parámetros. A continuación se explica la teoría que encierra el filtro de Kalman empleado y su formulación, así como una primera aplicación. En este ejemplo se implementa un algoritmo híbrido que combina el algoritmo GA y el filtro de Kalman explicado anteriormente. Seguidamente y atendiendo a las conclusiones del Capítulo 2, se describe la implementación de un algoritmo híbrido basado en el MHS que resuelva el problema de optimización. Para tal fin, el algoritmo MHS se combina con el filtro de Kalman explicado. La ventaja de la implementación de este algoritmo híbrido es doble. En primer lugar el tiempo de computación es reducido y, en segundo lugar, el error en la medida y en el proceso de estimación no se ignoran en la resolución del problema. La validación del algoritmo se realiza resolviendo un problema de actualización del modelo de EF de una pasarela de laboratorio. Los resultados se comparan con los obtenidos por dos algoritmos convencionales. Finalmente, se describe una aplicación adicional del algoritmo híbrido propuesto basado en el MHS y en el filtro de Kalman, para la estimación de los parámetros del modelo dinámico de Winkler en un problema de interacción suelo-estructura de una pasarela. Partes de este capítulo se han publicado en los Artículos B, C y D. 3.1 Fuentes de error Cuando se construye el modelo numérico/matemático de un sistema entran en juego numerosos tipos de error, algunos de los cuales se corrigen mediante la resolución del problema de estimación de parámetros (incertidumbre de los parámetros). En cambio, hay otros a los que hay que prestar atención. Dentro de este segundo grupo se encuentran las dos fuentes de error que se indican a continuación [ 82 ]. En primer lugar, los errores debido 31 32 Capítulo 3. Algoritmo híbrido UKF-HS para estimación de parámetros a las hipótesis asumidas para caracterizar al sistema. Estos errores generalmente surgen de simplificaciones en la modelización (por ejemplo, una placa modelizada como una viga) o asumir comportamiento lineal cuando realmente es no lineal. En segundo lugar, los errores generados por la distretización de la que hacen uso los métodos numéricos, aunque suelen ignorarse al ser de un orden muy inferior a los errores en la modelización [99]. Sin embargo, cuando se formula un problema de estimación de parámetros no solo aparecen los errores citados arriba asociados al modelo numérico sino que también existe incertidumbre en los datos experimentales. Estos errores pueden ser aleatorios o sistemáticos debido a imperfecciones en los equipos de medida o en el procesamiento de las señales. Finalmente, cuando se realiza la actualización de parámetros del modelo aparece un error, denominado error en el proceso de estimación. Este error se utiliza para establecer las estimaciones de una iteración basadas en la bondad de las estimaciones del paso anterior. Para tener en cuenta las fuentes de error, tanto del modelo como de las medidas y proceso de estimación, el algoritmo híbrido que se propone en este capítulo hace uso de las ventajas que presenta el algoritmo UKF. En el siguiente apartado se explica la formulación que sigue este algoritmo. 3.2 Unscented Kalman Filter Debido a las ventajas sobre el EKF (Extended Kalman Filter), el algoritmo de minimización local usado para acelerar al algoritmo harmony search multiobjetivo (MHS) es el UKF (Unscented Kalman Filter). Se trata de un algoritmo que no precisa de realizar ninguna operación diferencial ya que no es necesario calcular el Jacobiano ni el Hessiano. El UKF se desarrolló para resolver problemas de estimación de parámetros de la mano de Wan y van der Merwe [109, 108], asumiendo que puede ser formulado como: θ θ θk=θ θ θk−1+wk−1(3.1) zk=h(θ θ θk)+vk(3.2) con θ θ θ siendo el vector de parámetros a estimar, h la función de modelización no lineal, z z z el vector que agrupa las salidas del sistema dinámico y w w w y v v v los ruidos estadísticos del proceso de estimación y de las observaciones, respectivamente. Para la definición de estos dos últimos términos se asume que no están correlacionados y que se trata de ruido blanco con media cero y matrices de covarianza Q y R , respectivamente. La matriz R se calcula como la suma de dos componentes, el ruido en la medida y el ruido en la modelización [101]. La estimación se realiza por medio de 2nd+1 , siendo nd el número de parámetros a estimar, puntos de muestreo deterministas o puntos sigma que toman la anterior media y covarianza real de las cantidades a estimar. La propagación de estos puntos a través de la función no lineal h permite derivar la media y la covarianza posterior. Los puntos sigma se calculan mediante la factorización raíz cuadrática de la matriz de covarianza posterior, P . Por ello, esta matriz debe ser semidefinida positiva en cada paso, lo que supone el principal esfuerzo computacional del algoritmo. La factorización de Cholesky 3.2 Unscented Kalman Filter 33 puede emplearse de manera eficiente para realizar este proceso. No obstante, la matriz P ha de actualizarse en cada paso y la aparición de errores numéricos puede dar lugar a una matriz no semidefinida positiva. El algoritmo square-root UKF propuesto por van der Merwe y Wan [ 106 ] supera esta limitación ya que se evita factorizar en cada paso. En su lugar, la factorización de Cholesky ( A=√P=chol(P) , siendo P=AAT ) se propaga directamente asegurándose que la matriz de covarianza es semidefinida positiva. El algoritmo consiste en los dos siguientes pasos. El paso de predicción se lleva a cabo considerando los resultados anteriores para calcular los puntos sigma, predecir las estimaciones de la matriz de covarianza del error en la estimación, Aθ θ θ , y calcular las salidas del modelo zk|k−1. El paso de corrección o actualización estima la media y covarianza del error en la estimación posterior en términos de la matriz de ganancia de Kalman, K , las medidas, zobs, y las salidas del modelo ponderadas, ˆ z. La formulación del algoritmo UKF es la siguiente (en el Artículo del Apéndice C se encuentran más detalles de los fundamentos matemáticos de la formulación): ˆ θ θ θ=θ θ θprior (3.3) Aθ 0=chol(Pθ 0)(3.4) Bucle principal: parak =1 : NUKF (número de iteraciones) Predicción Calcula los 2nd+1puntos sigma :(χk−1)i(3.5) (χk|k−1)i= (χk−1)i(3.6) θ θ θk|k−1= 2nd ∑ 0 Wi·(χk|k−1)i(3.7) Aθ k|k−1=γ−0.5Aθ k−1|k−1(3.8) (zk|k−1)i=h((χk−1)i)(3.9) ˆzk|k−1= 2nd ∑ 0 Wi·(zk|k−1)i(3.10) Corrección Sz k|k−1=qr([qW1:2nd·[(zk|k−1)1:2nd−ˆ zk|k−1]√R]) (3.11) Sz k|k−1=cholupdate(Sz k|k−1,(zk|k−1)0−ˆ zk|k−1,sgn(W0)) (3.12) 40 Capítulo 3. Algoritmo híbrido UKF-HS para estimación de parámetros Tabla 3.4 Comparación de los tres algoritmos para la actualización del modelo de EF de la pasarela de laboratorio siguiendo el enfoque de un único objetivo. Iter es el número de iteraciones, Pob es el tamaño de la población, Eval es el número de evaluaciones de la función objetivo y fel valor de la función objetivo. Iter Pob NUKF Pθ 0Rj j Eval Tiempo f [s] GA 50 20 - - - 1020 9790 0.0036 HS 500 20 - - - 520 6640 0.0033 UKF-HS 2 5 1 diag((θ θ θu−θ θ θl 2000 )2)0.001 147 2100 0.0037 2 5 3 diag((θ θ θu−θ θ θl 500 )2)0.001 441 5680 0.0034 2 5 3 diag((θ θ θu−θ θ θl 2000 )2)0.001 441 5710 0.0033 2 5 3 diag((θ θ θu−θ θ θl 4000 )2)0.001 441 5690 0.0036 2 5 3 diag((θ θ θu−θ θ θl 2000 )2)0.1441 5680 0.0060 2 5 3 diag((θ θ θu−θ θ θl 2000 )2)0.0001 441 5640 0.0040 2 5 5 diag((θ θ θu−θ θ θl 2000 )2)0.001 735 9440 0.0033 siguientes valores: NUKF =3 , Rj j =0.001 y Pθ 0=diag((θ θ θu−θ θ θl 2000 )2) , ya que mejoran la precisión sin sacrificar el tiempo de computación. Los resultados se comparan con los proporcionados por dos algoritmos diferentes, el HS y el GA. Para el primero se adoptan los mismos valores de los parámetros hmcr , par y bw. En el segundo caso, se asume que la relación de cruce es de 0.8 y la de mutación 0.4 [ 74 ]. Esta comparación se establece en base a dos criterios, la velocidad de convergencia, analizada a través del tiempo de simulación y el número de evaluaciones de la función objetivo, y la precisión del ajuste, analizada a través del valor de la función objetivo. Los resultados obtenidos mediante estos dos algoritmos se incluyen en la Tabla 3.4. A la vista de los resultados, el algoritmo híbrido UKF-HS demuestra ser el más eficiente ya que presenta la misma precisión en el resultado que el HS pero siendo el tiempo de simulación y el número de evaluaciones de la función objetivo menor. Además, es preciso incidir en que el GA es el menos eficiente de los tres comparados para resolver este problema de actualización. Actualización utilizando el algoritmo UKF-MHS La actualización del modelo de EF de la pasarela se lleva a cabo también aplicando el algoritmo híbrido basado en un enfoque multiobjetivo del problema de optimización. Los parámetros del algoritmo híbrido son los mismos que los adoptados para la resolución del problema de un único objetivo. La selección del punto óptimo de entre todos aquellos que forman el frente de Pareto se hace aplicando el método de intersección normal al 3.4 Aplicación para estimar los parámetros del modelo dinámico de Winkler 41 contorno [ 23 ]. Con el fin de comparar los resultados con otros algoritmos, la actualización del modelo de EF se realiza también aplicando los algoritmos MHS y el GA multiobjetivo (MGA), adoptando los mismos valores de los parámetros de los dos algoritmos que en el caso anterior. En primer lugar, se realiza un análisis de sensibilidad de los parámetros del algoritmo híbrido. Para ello, se consideran los tres mismos parámetros, NUKF , Rj j y Pθ 0 . Los resultados de este análisis se muestran en la Tabla 3.5. En relación al parámetro NUKF se obtiene una relación directa entre el valor de este parámetro y la precisión del ajuste y una relación inversa entre este parámetro y el tiempo de simulación. Respecto al parámetro Rj j , un valor intermedio, Rj j =0.001 , mejora la precisión del ajuste sin comprometer el tiempo de simulación. Finalmente, el parámetro Pθ 0 dota de una mayor precisión en el ajuste cuando su valor disminuye aunque su influencia no es significativa. Los valores que se adoptan para los tres parámetros son los siguientes: NUKF =3 , Rj j =0.001 y Pθ 0=diag((θ θ θu−θ θ θl 2000 )2). Finalmente, el rendimiento del algoritmo híbrido aplicado a la actualización del modelo de EF se compara con el de los dos algoritmos convencionales mencionados anteriormente, el MHS y el GA. Como criterios de comparación se tienen en consideración el número de evaluaciones de la función objetivo y el tiempo de simulación para analizar la velocidad de convergencia y, por otro lado, la precisión en el ajuste evaluada en términos de la distancia entre el punto seleccionado del frente de Pareto y el origen y de la suma de los dos términos de la función objetivo. En la Tabla 3.5 se detallan los resultados obtenidos para la actualización del modelo de EF implementando los tres algoritmos. Atendiendo a los resultados puede comprobarse que el algoritmo UKF-MHS es más eficiente, ya que tiene un tiempo de simulación menor alcanzando la misma precisión en el ajuste. Por otro lado, es nuevamente destacable que el algoritmo MGA es el que proporciona un peor ajuste, además, con un tiempo de simulación mayor. Mediante estos dos problemas de optimización se ha validado el algoritmo híbrido propuesto tanto en la versión para un único objetivo como la de multiobjetivo. 3.4 Aplicación para estimar los parámetros del modelo dinámico de Winkler Además de estimar los parámetros del modelo de EF de una estructura, el algoritmo híbrido propuesto en este capítulo puede implementarse para resolver cualquier problema donde haya que estimar parámetros de un sistema donde hay involucradas mediciones. En este apartado, el algoritmo se aplica para estimar los parámetros del modelo dinámico de Winkler para la interacción suelo-estructura de una pasarela integral de estructura mixta de acero y hormigón situada en la provincia de Sevilla (ver Figura 3.4). El modelo de Winkler es un modelo simplificado para simular la interacción sueloestructura. Su formulación depende del tipo de cimentación que tenga la estructura, que grosso modo, pueden clasificarse en cimentaciones directas y cimentaciones profundas. El primer caso se caracteriza por que las cimentaciones son, generalmente, estructuras masivas rígidas que se modelizan a través de elementos muelle-amortiguador equivalentes 42 Capítulo 3. Algoritmo híbrido UKF-HS para estimación de parámetros Tabla 3.5 Comparación de los tres algoritmos para la actualización del modelo de EF de la pasarela de laboratorio siguiendo el enfoque multiobjetivo. Iter es el número de iteraciones, Pob es el tamaño de la población, N.Pob es el tamaño de la población que se crea en cada iteración, Eval es el número de evaluaciones de la función objetivo, Dist es la distancia euclídea al origen y Sum es la suma de f1yf2. Iter Pob N.Pob NUKF Pθ 0Rj j Eval Tiempo [s] Dist Sum (·10−5) (·10−4) (·10−4) MGA 25 50 - - - - 1450 0.18 50.7 68.6 45 50 - - - - 2570 0.33 50.2 66.9 150 10 - - - - 16600 2.12 52.0 68.0 MHS 50 50 25 - - - 1300 0.17 45.4 62.1 100 100 25 - - - 2600 0.33 45.1 61.8 300 200 50 - - - 15200 1.98 43.7 60.6 UKF-MHS 3 10 3 3 diag((θ θ θu−θ θ θl 2000 )2)0.001 1197 0.16 45.2 62.1 3 10 3 3 diag((θ θ θu−θ θ θl 500 )2)0.001 1197 0.15 43.0 63.5 3 10 3 3 diag((θ θ θu−θ θ θl 4000 )2)0.001 1197 0.15 44.6 61.4 3 10 3 3 diag((θ θ θu−θ θ θl 2000 )2)0.11197 0.15 52.6 71.3 3 10 3 3 diag((θ θ θu−θ θ θl 2000 )2)0.0001 1197 0.15 145.0 177.6 8 15 3 3 diag((θ θ θu−θ θ θl 2000 )2)0.001 2457 0.32 44.3 61.1 20 40 10 3 diag((θ θ θu−θ θ θl 2000 )2)0.001 15120 1.96 43.7 60.4 3 8 3 1 diag((θ θ θu−θ θ θl 2000 )2)0.001 357 0.05 50.6 68.2 3 8 3 3 diag((θ θ θu−θ θ θl 2000 )2)0.001 1071 0.14 49.1 66.5 3 8 3 5 diag((θ θ θu−θ θ θl 2000 )2)0.001 1785 0.23 45.2 62.0 [ 35 ] (Figura 3.5). No obstante, para niveles bajos de vibración de la cimentación la componente de amortiguamiento puede ser excluida sin perjuicio de la validez en la simulación. La rigidez equivalente en los seis grados de libertad, tres de traslación y tres de giro, se calculan a partir de los datos geométricos de la cimentación y de los parámetros mecánicos del suelo ya sean cimentaciones embebidas en el suelo o superficiales según la formulación propuesta por Gazetas [ 35 ]. En el Artículo del Apéndice D se pueden consultar las fórmulas que permiten obtener las rigideces del modelo de interacción suelo-estructura. Por otro lado, las cimentaciones profundas se caracterizan por ser estructuras flexibles y esbeltas cuyo comportamiento estructural depende principalmente de su longitud. En este caso es preciso tener en cuenta la componente de amortiguamiento para simular adecuadamente la interacción suelo-estructura ya que los efectos dinámicos tienen una mayor influencia (Figura 3.5). El modelo dinámico de Winkler es el modelo más simple debido a la buena relación entre complejidad y precisión en los resultados que otorga. La rigidez se calcula en función de las propiedades del suelo y de la cimentación. El amortiguamiento se calcula como la suma de dos mecanismos de disipación de energía, a saber, amortiguamiento histerético (debido al material) y amortiguamiento por radiación [ 36 ]. Se asume que la contribución del amortiguamiento histerético no se tiene en cuenta 3.4 Aplicación para estimar los parámetros del modelo dinámico de Winkler 43 a) b) c) Figura 3.4 Pasarela para la estimación de los parámetros del modelo dinámico de Winkler: a) vista frontal, b) plano de la sección y c) plano de alzado. en la formulación del problema ya que presenta dos órdenes de magnitud inferior al amortiguamiento de radiación en el rango de frecuencias de interés. En el Artículo del Apéndice D se encuentran las ecuaciones empleadas para calcular la componente de rigidez y de amortiguamiento del modelo dinámico de Winkler. La pasarela integral tiene un único vano de 12 m de longitud y 2.5m de ancho. La sección transversal es mixta con un cajón de acero de 0.4m de alto y un tablero de hormigón de 0.2m de alto y 2.5m de ancho. Los estribos de los extremos se apoyan sobre 4 micropilotes de 0.2m de diámetro y 15 m de longitud. El modelo de EF se implementa en el software ANSYS [ 78 ]. El tablero de hormigón, la sección cajón de acero, la barandilla y los micropilotes se modelizan mediante elementos tipo viga de dos nodos y seis grados de libertad por nodo (Beam188). La interacción suelo-estructura se simula mediante elementos tipo muelle-amortiguador unidimensionales (COMBIN14) colocados en cada nodo de los micropilotes y los estribos (Figura 3.5). El modelo de EF presenta un total de 2008 elementos. Para la caracterización de los elementos muelle-amortiguador se parte de los datos del terreno, que consiste en una primera capa de material de relleno de 5m de espesor y una segunda capa de arcilla de espesor indefinido y de los que se dispone el valor de la densidad, módulo de Young, coeficiente de Poisson y coeficiente de balasto. El análisis modal da como resultados las tres frecuencias naturales numéricas y sus amortiguamientos numéricos asociados, ξnum , que se muestran en la Tabla 3.6. Las propiedades modales experimentales se obtienen mediante el análisis modal operacional de las señales de 44 Capítulo 3. Algoritmo híbrido UKF-HS para estimación de parámetros Figura 3.5 Modelo de interacción suelo-estructura de cimentaciones superficiales y profundas y modelo de EF de la pasarela. [Tomada del Artículo del Apéndice D]. Tabla 3.6 Propiedades modales numéricas y experimentales de la pasarela integral antes de la actualización. Modo fnum [Hz]fexp [Hz]ξnum [%]ξexp [%]∆f[%]∆ξ[%]MAC [-] 1 15.322 15.602 3.103 4.650 1.795 33.268 0.975 2 32.179 35.971 1.121 1.049 10.542 6.863 0.692 3 38.614 40.017 0.700 4.846 3.506 85.555 0.907 aceleraciones medidas durante un ensayo de vibración ambiental. En el Artículo del Apéndice D se encuentran más detalles sobre la configuración del ensayo y el tratamiento de las señales. Las tres frecuencias naturales experimentales obtenidas se recogen en la Tabla 3.6. Además, también se han calculado los amortiguamientos modales experimentales, ξexp , [ 70 ]. Puede observarse que el modelo no simula el comportamiento real de la estructura ya que existen diferencias altas entre los valores numéricos y experimentales, por lo que se procede a actualizar el modelo [117]. El conjunto de parámetros seleccionado resulta de la combinación de aquellos que tienen una mayor influencia en la energía de deformación modal y los parámetros que 3.5 Conclusiones 45 Tabla 3.7 Propiedades modales numéricas y experimentales de pasarela integral después de la actualización. Modo fnum [Hz]fexp [Hz]ξnum [%]ξexp [%]∆f[%]∆ξ[%]MAC [-] 1 16.380 15.602 4.588 4.650 4.987 1.333 0.975 2 35.377 35.971 1.058 1.049 1.651 0.858 0.907 3 40.599 40.017 4.701 4.846 1.454 2.992 0.964 caracterizan el fenómenos de interacción suelo-estructura. Dentro de éstos últimos se hallan: el módulo de Young del relleno y de la arcilla y el coeficiente de balasto también de ambos materiales. Para la implementación del algoritmo híbrido se adoptan los siguientes valores. El número de armonías generadas aleatoriamente por el algoritmo MHS al inicio de la simulación es de 40, en cada iteración del MHS se crean 20 armonías nuevas y el número de iteraciones del MHS es de 20. Respecto a la implementación del algoritmo UKF se consideran 5 iteraciones, el vector de observaciones se construye a partir de las tres frecuencias naturales experimentales, los tres amortiguamientos experimentales y tres componentes unitarias relativas al residuo de los modos. Se adoptan los siguientes valores para los factores de escala: λ=0.0001 [4] y para γ=0.99 [106]. A diferencia de las otras aplicaciones, en este problema se consideran tres funciones objetivos que se definen a partir de los residuos de las frecuencias naturales, modos de vibración y amortiguamientos modales, respectivamente. La actualización del modelo de EF tiene de este modo dos ventajas. La primera es obtener un modelo de EF más preciso y la segunda es definir el modelo dinámico de Winkler a partir de los valores de los parámetros del terreno actualizados. Los parámetros del terreno obtenidos son: módulo de Young del relleno: 7.77·108N/m2 ; módulo de Young de la arcilla: 3.5·107N/m2 ; coeficiente de balasto del relleno: 4.16 · 107N/m3 y coeficiente de balasto de la arcilla: 1.27 ·108N/m3 . A partir de estos valores se estiman los parámetros del modelo dinámico de Winkler dando como resultado: rigidez en la capa del relleno: 4.16 ·106N/m ; amortiguamiento en la capa del relleno: 7.93 ·105sN/m ; rigidez en la capa de arcilla: 1.27 ·107N/m y amortiguamiento en la capa de arcilla: 1.68 ·105sN/m . La validación de estos resultados se realiza en base a la diferencias entre las propiedades modales experimentales y actualizadas. Como puede observarse en la Tabla 3.7, las diferencias relativas entre frecuencias naturales y amortiguamientos han disminuido y el valor del MAC ha aumentado. Además, al haber definido un rango de búsqueda para los parámetros, se garantiza que los valores encontrados mantienen un significado físico adecuado. 3.5 Conclusiones Los algoritmos de optimización usados para la estimación de parámetros de estructuras de ingeniería civil, generalmente, no tienen en consideración el error en el proceso de estimación ni en las medidas. En el presente capítulo se ha propuesto, desarrollado e 46 Capítulo 3. Algoritmo híbrido UKF-HS para estimación de parámetros implementado un algoritmo híbrido que permita resolver de manera eficiente el problema de estimación de parámetros bajo la hipótesis de errores asumidos como ruido blanco gaussiano en las medidas y en el proceso de estimación. Este algoritmo resulta de la combinación de dos técnicas, un algoritmo de búsqueda global, el algoritmo HS, y otro algoritmo de búsqueda local, el algoritmo UKF. Bajo este esquema, el carácter local del algoritmo UKF acelera la búsqueda de una solución global por parte del algoritmo HS. Para garantizar el buen funcionamiento del algoritmo en múltiples problemas, se ha desarrollado para dos casos distintos: el enfoque de optimización de un solo objetivo y el enfoque de optimización multiobjetivo. Esta propuesta se ha validado a través de la solución de un problema de actualización del modelo de EF de una pasarela donde se estiman los parámetros físicos del modelo de manera que las propiedades modales numéricas y experimentales de la pasarela sean lo más parecidas posible. El problema se ha resuelto asumiendo los dos enfoques citados anteriormente y los resultados se han comparado con los obtenidos por dos algoritmos metaheurísticos convencionales. La comparación ha demostrado que el algoritmo propuesto es más eficiente ya que logra resultados más precisos en un tiempo de simulación menor. Finalmente, se ha descrito una aplicación adicional del algoritmo para resolver un problema de estimación de parámetros. Este problema consiste en la estimación de los parámetros del modelo dinámico de Winkler para el modelo de una pasarela integral donde el fenómeno de interacción suelo-estructura tiene gran influencia. De esta forma se prueba que el algoritmo híbrido puede extenderse con éxito a otros tipos de problemas de ingeniería civil. 4 Algoritmo colaborativo de inteligencia artificial para estimación de parámetros El presente capítulo aborda la implementación de un nuevo algoritmo colaborativo que, además de reducir el tiempo de computación, permite elegir la mejor solución entre todo el conjunto de posibles soluciones óptimas que forman el frente de Pareto de una manera robusta. Partes del presente capítulo se pueden encontrar en los Artículos de los Apéndices E y F. La organización del capítulo es la siguiente. En primer lugar se expone el problema que presenta el uso de algoritmos de optimización multiobjetivo. Bajo este enfoque no se obtiene una solución única, sino un conjunto de soluciones y, por tanto, para finalizar el problema, hay que elegir el punto óptimo de entre todos ellos. Se describen los métodos más comunes encontrados en la literatura para resolver este problema. Seguidamente se describe el algoritmo que se propone en esta tesis para resolver problemas de estimación de parámetros empleando algoritmos multiobjetivo, reduciendo el coste computacional e implementando un método para resolver la toma de decisiones de manera robusta. Asimismo, se detallan las técnicas de inteligencia artificial y numéricas que se consideran en el algoritmo. Finalmente se expone una aplicación del algoritmo para su validación, aplicándolo a un problema de actualización del modelo de EF de una pasarela. 4.1 El problema de la toma de decisiones en optimización multiobjetivo La resolución de un problema de optimización multiobjetivo da como resultado el frente de Pareto, que está formado por la frontera de puntos para los que no existe otra solución tal que mejore un objetivo sin sacrificar el otro. Cada punto que forma esta curva es una posible solución al problema y la elección de la solución más adecuada generalmente no 47 48 Capítulo 4. Algoritmo colaborativo de inteligencia artificial para estimación de parámetros es inmediata ya que no existe un punto que minimice ambos objetivos simultáneamente (ver Figura 4.1). Para ayudar en la tarea de elegir la mejor solución se realiza un análisis teniendo en cuenta múltiples criterios que permite obtener la solución óptima de una manera sistemática. No obstante, como puede apreciarse en la Figura 4.1, la solución óptima es sensible a la discretización del frente de Pareto, ya que a menor número de puntos la probabilidad de que la solución óptima esté mal representada es mayor. feasible region Infeasible region GOOD ←− Objective 1 −→BAD GOOD ←− Objective 2 −→BAD Pareto-optimal front Dominated solutions Utopian solution Figura 4.1 Representación de compromiso en una optimización de dos objetivos. Tomada del Artículo E. A este tipo de análisis se les denomina problema o estrategia de toma de decisiones y pueden clasificarse según el momento en el que se aplica en [14]: • Métodos a priori: antes de realizar el proceso de optimización es necesario aportar suficiente información sobre la preferencia de los objetivos. • Métodos interactivos: consiste en expresar las preferencias en cada iteración del algoritmo de optimización. • Métodos a posteriori: Las preferencias acerca de los objetivos se introducen una vez realizado el proceso de optimización. Los métodos a posteriori constituyen la opción más empleada debido a dos razones. La primera de ellas es que la toma de decisión se realiza considerando toda la información de compromiso entre los objetivos. La segunda razón es que los algoritmos son más simples en comparación con aquellos empleados en los otros dos métodos ya que no es necesario modificarlos para incluir las preferencias. Las preferencias se formulan de manera que haya un compromiso entre los valores de los objetivos. En este sentido, la solución del frente de Pareto que para una pequeña mejora en un objetivo conlleve un gran deterioro en al menos otro objetivo es las más adecuada. Este punto se conoce como knee point [ 22 ]. 4.1 El problema de la toma de decisiones en optimización multiobjetivo 49 Por tanto, la definición del knee point es fundamental a la hora de resolver el problema de toma de decisiones. Existen en la literatura diversas definiciones del knee point, de las cuales se describirán a continuación las más empleadas. En el Artículo del Apéndice E se puede encontrar una definición más detallada de cada uno de los siguientes métodos. •Método de intersección normal al contorno [22]. La idea de este método es unir mediante una línea recta los dos puntos extremos del frente de Pareto normalizado para construir el contorno. El punto que se encuentre a una mayor distancia en la dirección normal al contorno es el knee point. •Método del ángulo de reflexión [15]. Este método define el knee point como aquel que tenga el máximo ángulo de reflexión con el punto inmediatamente próximo a la izquierda y a la derecha. •Método del ángulo de flexión [23]. Análogo al anterior, el knee point que proporciona este método es aquel que tiene el máximo ángulo de flexión positivo. •Método basado en la utilidad [15]. Este método requiere definir una función de utilidad lineal que es función del valor de los dos objetivos y de un parámetro, λ: U(f,λ) = λf1+(1−λ)f2(4.1) siendo f1 y f2 el valor de los dos objetivos de un punto del frente de Pareto. La aplicación del método consiste en los siguientes pasos: i) para cada solución del frente de Pareto se define un contador c que toma el valor nulo inicialmente, ii) se toma un vector con valores de λ∈[0,1] dividiendo el rango en un número dado, iii) para cada valor de λ se calcula el valor de U(f,λ) para todos los puntos del frente de Pareto, iv) para cada valor de λ a aquella solución del frente de Pareto que tenga el valor mínimo de U se le incrementa su contador en una unidad, v) se repiten los pasos ii)-iv) iterativamente. Para garantizar la precisión de este método es fundamental que el número de divisiones del rango realizado en el paso ii) no sea el mismo en cada iteración. La solución del frente de Pareto con mayor valor de ces el knee point. •Método de compromiso [23]. Finalmente, este método dota al usuario de un mayor control en la toma de decisiones ya que se especifican dos parámetros, α y β , para proporcionar la información sobre la preferencia en el balance de objetivos. En el Artículo del Apéndice E se realiza un análisis comparativo de los métodos anteriores para la obtención de la mejor solución de entre el conjunto de posibles soluciones que forman el frente de Pareto. Para tal fin, el modelo de EF de una pasarela real se actualiza en base a las propiedades modales experimentales empleando un algoritmo de optimización con dos funciones objetivos. En la Figura 4.2 se muestra el frente de Pareto resultante y la solución óptima según los cuatro primeros métodos anteriores, siendo NBI el método de 56 Capítulo 4. Algoritmo colaborativo de inteligencia artificial para estimación de parámetros Figura 4.5 Pasarela en estudio: a) ilustración de la pasarela, b) modelo de EF y c) configuración del ensayo de vibración ambiental. 4.3.1 Actualización del modelo de EF aplicando el algoritmo colaborativo Para resolver el problema de actualización del modelo de EF de la pasarela se eligen como variables de diseño los parámetros físicos del modelo que presentan mayor incertidumbre. La elección de estos parámetros a actualizar se hace de manera sistemática realizando un análisis de sensibilidad que permite conocer aquellos que tengan mayor influencia en alguna característica de referencia de la estructura. La energía de deformación modal, siguiendo la formulación de Fox y Kapoor [ 29 ], es el parámetro empleado para dilucidar cuáles son los parámetros físicos que más influyen en el comportamiento dinámico de la pasarela. En la Figura 4.6 se representa la matriz de sensibilidad, donde puede apreciarse aquellos parámetros que tienen una mayor influencia en la energía de deformación modal. En vista a los resultados, se toman los siguientes doce parámetros de actualización: rigidez de la interacción suelo-pilote de cada pilote (tres), rigideces en los estribos (longitudinal, transversal, giro sobre eje vertical y giro sobre eje transversal), módulo de Young del tablero y módulos de Young del acero de cada vano (cuatro). En el Artículo del Apéndice F 4.3 Ejemplo de aplicación 57 Tabla 4.1 Propiedades modales numéricas y experimentales de la pasarela. Modo fnum [Hz]fexp [Hz]∆f[%]MAC [-] 1. Lateral 1.704 1.872 -8.974 0.997 2. Lateral 3.560 3.506 1.540 0.969 3. Lateral 6.424 5.851 9.793 0.922 4. Longitudinal 6.610 6.914 -4.397 0.966 5. Lateral 8.868 8.217 7.923 0.958 6. Lateral 9.486 8.515 11.403 0.929 7. Vertical 9.483 9.189 3.199 0.947 8. Vertical 9.720 9.483 2.499 0.914 9. Vertical 11.788 11.240 4.875 0.875 10. Lateral 13.854 11.828 17.129 0.765 11. Vertical 14.068 12.812 9.803 0.816 se puede encontrar una descripción más detallada del análisis de sensibilidad llevado a cabo. Figura 4.6 Representación gráfica de la matriz de sensibilidad para la elección de los parámetros a actualizar. La estimación de parámetros empleando el algoritmo colaborativo comienza con la generación de un conjunto de posibles soluciones aleatorias. Para ello, 100 armonías iniciales se crean aleatoriamente usando el algoritmo MHS donde cada armonía es un vector de 12 componentes (los 12 parámetros físicos a actualizar). Para obtener el frente de Pareto poco poblado se adoptan los siguientes valores para los parámetros del algoritmo MHS: número máximo de iteraciones, 100; nuevas armonías en cada iteración, 25; hmcr =0.8 y 58 Capítulo 4. Algoritmo colaborativo de inteligencia artificial para estimación de parámetros par =0.4 [ 110 ]. La ejecución de este proceso iterativo tiene una duración de t1=11700 s , empleado un ordenador con un procesador de 3.6GHz y una memoria RAM de 8Gb . En la Figura 4.7a se muestra el frente de Pareto poco poblado obtenido mediante el algoritmo MHS. A continuación, se realiza el análisis de componentes principales del frente de Pareto. La matriz de covarianza de los datos se descompone aplicando el método de descomposición en valores singulares. En la Figura 4.7b se puede observar la representación en el espacio de componentes principales del frente de Pareto procesado. Como se indicó anteriormente, esta representación toma la forma de una curva convexa cuya principal ventaja es la simplicidad que toma la definición de knee point. La ejecución de este paso es inmediata asumiendo que el orden de magnitud del proceso completo será del orden de t1 . Por ello, el tiempo de ejecución de este paso se redondea a t2=1s. Al ser el frente de Pareto procesado una curva formada por pocos puntos, cabe la posibilidad de que no tenga representada la solución óptima del problema. La red neuronal artificial permite tener una función continua que aproxima al frente de Pareto procesado. La configuración de esta red neuronal consiste en una capa de entrada con 12 neuronas (parámetros a estimar), una capa de salida con 2 neuronas (número de funciones objetivo) y una capa oculta como se justificó anteriormente. El número de neuronas de esta capa oculta se calcula empleando la regla de Kermanshahi [56] neuronas =m+n 2+δ(4.2) siendo m=12 , n=2 y δ un factor de normalización que puede tomar el valor de 1 o 2. Asumiendo δ=1 se obtiene que el número de neuronas de la capa oculta de 8. El entrenamiento de la red para calcular los pesos de las conexiones entre neuronas se realiza usando el algoritmo de Levenberg-Marquardt, formulando la función de error mediante el error cuadrático medio [ 41 ]. La convergencia del algoritmo de entrenamiento se alcanza en la iteración 8 teniendo un tiempo de ejecución del orden de t3=1s . En la Figura 4.7b se representa el frente de Pareto aproximado por la red neuronal donde puede contemplarse la exactitud con la que la red aproxima los puntos. Finalmente, aprovechando las virtudes que el frente de Pareto aproximado presenta por sus propiedades convexas, la elección del knee point se realiza resolviendo un problema de minimización con restricciones usando el algoritmo de minimización local activeset. Las restricciones implementadas en el algoritmo representan los límites donde el frente de Pareto aproximado queda definido. El tiempo de simulación de este paso es de aproximadamente t4=2s . En el Artículo del Apéndice F se pueden consultar más detalles sobre la resolución del problema de actualización del modelo de EF empleando el algoritmo colaborativo propuesto. La solución obtenida es la solución óptima al problema de estimación de parámetros y los parámetros físicos asociados son aquellos que aseguran el mejor ajuste entre las propiedades modales numéricas y experimentales de la pasarela. 4.3 Ejemplo de aplicación 59 a) 8·10−29·10−20.1 0.11 0.12 5.5 6 6.5 7 7.5 8·10−2 f1 f2 b) −2·10−202·10−2 −4 −2 0 2 4 6·10−3 Primera componente principal Segunda componente principal Valor real Salida estimada Figura 4.7 Frente de Pareto del problema de estimación de parámetros de la pasarela: a) Original poco poblado, b) frente de Pareto procesado en el espacio de componentes principales y aproximado por la red neuronal artificial. 4.3.2 Comparación de resultados y tiempos de ejecución con un algoritmo convencional Para evaluar el rendimiento del algoritmo colaborativo propuesto, se comparan los resultados y el tiempo de simulación con los de la actualización del modelo de EF considerando el algoritmo MHS de manera clásica, es decir, generando como resultado un frente de Pareto con muchos puntos y procesándolo con alguno de los métodos de toma de decisiones. Se asume, por tanto, que para garantizar que el frente de Pareto está bien definido la población inicial del algoritmo MHS es de 200 armonías generadas aleatoriamente. Los valores de hmcr y par son los mismos indicados anteriormente. El problema de la toma de decisiones para obtener la solución óptima se resuelve empleando el método de intersección normal al contorno descrito en el apartado 4.1 [ 22 ]. Para enriquecer la comparación entre ambos métodos se lleva a cabo un estudio de sensibilidad para analizar la variación de la mejor solución en términos del número de iteraciones considerado. Con este propósito, la actualización del modelo de EF se resuelve para cuatro valores diferentes del número de iteraciones: 150, 200, 250 y 300. En las Figuras 4.8a-d se representa una comparación gráfica de la solución obtenida empleando el algoritmo colaborativo y las soluciones considerando el algoritmo MHS convencional. Además, en la Tabla 4.2 se muestra numéricamente la comparación entre las soluciones en términos de distancia entre los puntos obtenidos por el algoritmo colaborativo y el convencional y en términos del tiempo de simulación. El tiempo que tarda el algoritmo colaborativo se obtiene como la suma de t1,t2,t3yt4. Queda ilustrado el buen rendimiento del algoritmo colaborativo para abordar problemas de estimación de parámetros basados en optimización multiobjetivo debido a la clara convergencia que exhibe la solución con la solución obtenida mediante la aplicación del método convencional cuando el número de iteraciones aumenta. En particular, el algoritmo 60 Capítulo 4. Algoritmo colaborativo de inteligencia artificial para estimación de parámetros a) 8·10−20.1 0.12 6·10−2 7·10−2 8·10−2 9·10−2 0.1 0.11 f1 f2 Pareto poco poblado Algoritmo propuesto Pareto convencional Método NBI b) 9·10−20.1 0.11 0.12 0.13 5·10−2 6·10−2 7·10−2 8·10−2 f1 f2 Pareto poco poblado Algoritmo propuesto Pareto convencional Método NBI c) 9·10−20.1 0.11 0.12 6·10−2 7·10−2 8·10−2 f1 f2 Pareto poco poblado Algoritmo propuesto Pareto convencional Método NBI d) 9·10−20.1 0.11 0.12 6·10−2 7·10−2 8·10−2 f1 f2 Pareto poco poblado Algoritmo propuesto Pareto convencional Método NBI Figura 4.8 Comparación entre las soluciones del algoritmo propuesto y el convencional para a) 150 iteraciones, b) 200 iteraciones, c) 250 iteraciones y d) 300 iteraciones. propuesto permite reducir el tiempo de simulación en más de un 50%. De este modo, se ha validado el algoritmo propuesto. Además, se postula como una herramienta eficiente para resolver problemas de estimación de parámetros, reduciendo el tiempo de computación sin perder precisión en los resultados y abordando el problema de la toma de decisiones de manera robusta. 4.4 Conclusiones Los problemas de optimización multiobjetivo presentan dos principales inconvenientes: el primero es el alto coste computacional que requieren y el segundo es la necesidad de 4.4 Conclusiones 61 Tabla 4.2 Valores de la función multiobjetivo para el modelo mejor actualizado, distancia entre las mejores soluciones obtenidas por ambos algoritmos y tiempo total de simulación. Algoritmo Iteraciones f1(θ θ θ)f2(θ θ θ)distancia [−]t[s] Colaborativo 100 9.47 ·10−25.99 ·10−2- 11704 MHS 150 9.09 ·10−26.41 ·10−25.66 ·10−316336 MHS 200 9.78 ·10−25.78 ·10−23.74 ·10−323640 MHS 250 9.52 ·10−25.92 ·10−28.60 ·10−425194 MHS 300 9.48 ·10−25.94 ·10−25.10 ·10−429582 resolver un problema de toma de decisiones para elegir el mejor punto del frente de Pareto cuya precisión depende de la definición de un Pareto suficientemente poblado. A pesar de la existencia de varios métodos para resolver este problema de toma de decisiones, el objetivo de este capítulo es desarrollar e implementar un algoritmo colaborativo que afronta los dos inconvenientes anteriores de manera simultánea. Para ello se emplean técnicas de inteligencia artificial, estadísticas y de optimización. En concreto, el algoritmo propuesto se describe como una sucesión de pasos empleando el algoritmo MHS, el análisis de componentes principales, una red neuronal y finalmente un algoritmo de minimización local para obtener el punto óptimo. El rendimiento del algoritmo se ha validado mediante la aplicación del mismo para resolver un problema de actualización del modelo de EF de una pasarela y comparando los resultados con los obtenidos por un algoritmo convencional. Esta comparación reporta unas reducciones de tiempo de computación de más del 50% sin comprometer la precisión de la solución. 5 Conclusiones y trabajos futuros 5.1 Conclusiones En esta tesis se han propuesto, implementado y validado dos algoritmos de optimización para la estimación de parámetros, que de un modo u otro, presentan una mayor eficiencia que los algoritmos convencionales cuando se emplean para problemas de ingeniería civil. Para llevar a cabo estas propuestas se ha realizado, en primer lugar, un análisis comparativo para conocer el algoritmo de optimización que mejor rendimiento ofrece, así como cual de los dos enfoques posibles, la optimización de un único objetivo y la optimización multiobjetivo, presenta mejores resultados. Para la comparación entre algoritmos de optimización se ha optado por tres algoritmos metaheurísticos, dos de ellos ampliamente aplicados para la actualización de modelos de EF, el GA y el PS. El tercero es el algoritmo HS, del que existen muy pocos casos prácticos donde se aplique a la actualización de modelos de EF. Para realizar la comparación se ha resuelto un problema de actualización del modelo de EF de una pasarela de laboratorio. Los resultados de este análisis comparativo han reportado una mayor eficiencia del algoritmo HS, tanto en un único objetivo como en multiobjetivo, ya que obtiene una solución con la misma precisión que el PS pero con coste computacional menor. Es destacable también que el GA es el que presenta una peor eficiencia de los tres. Comparando los dos enfoques entre sí, se concluye que el enfoque multiobjetivo ofrece soluciones más precisas en términos de los valores de los residuos ya que los puntos de la optimización de un único objetivo se encuentran más lejos del knee point del frente de Pareto. Además, el tiempo de simulación es mayor para el caso de un único objetivo ya que es necesario llevar a cabo un análisis de sensibilidad de los pesos de los residuos. Estos dos motivos justifican la elección del algoritmo HS en optimización multiobjetivo (MHS) para la formulación de los dos algoritmos que se han propuesto en esta tesis. El primer algoritmo que se ha propuesto es un algoritmo híbrido UKF-MHS que aprovecha las virtudes que un algoritmo global (el MHS) y uno local (el UKF) presentan de manera independiente. El algoritmo UKF, al pertenecer a la familia de los filtros de Kalman, permite considerar el error tanto en las medidas (observaciones) como en el 63 64 Capítulo 5. Conclusiones y trabajos futuros proceso de estimación para converger a un valor esperado de los parámetros. Además, el carácter global del MHS asegura que el algoritmo no convergerá a un mínimo local de la función objetivo. De este modo, la ventajas de este algoritmo es doble: primero, el modelo actualizado es más preciso ya que al considerar los errores, las discrepancias entre el modelo numérico actualizado y la estructura real se reducen y, segundo, se aprovecha el carácter local-global para acelerar la búsqueda de soluciones sin caer en un mínimo local. La validación del algoritmo se ha realizado a través de la actualización del modelo de EF de una pasarela de laboratorio, comparando los resultados con los reportados por los algoritmos GA y HS. Como principal conclusión cabe indicar que el algoritmo propuesto presenta una mayor eficiencia ya que la solución elegida del frente de Pareto se encuentra a menor distancia del origen de coordenadas. El GA es el que presenta un peor rendimiento de los tres en términos de precisión y tiempo de simulación. Se ha presentado también una segunda aplicación del algoritmo para estimar los parámetros del modelo dinámico de Winkler en el estudio del fenómeno de interacción suelo-estructura de una pasarela integral. Los parámetros se han definido a partir de los valores de las propiedades mecánicas del terreno calibrados en base a las propiedades modales de la estructura. El enfoque multiobjetivo de la optimización da como resultado un conjunto de soluciones óptimas representado en el frente de Pareto. Elegir la adecuada exige resolver un problema de toma de decisiones donde la precisión de la solución depende del número de puntos que formen el Pareto. Para abordar esta cuestión, el segundo algoritmo propuesto en esta tesis se define como un algoritmo colaborativo, ya que las técnicas no se combinan sino que se aplican de manera sucesiva, eficiente y que resuelve el problema de toma de decisiones de manera sistemática sin necesidad de disponer de un elevado número de puntos en el Pareto. Este algoritmo colaborativo comienza por la aplicación del algoritmo MHS para obtener el frente de Pareto poco poblado. En segundo lugar, utiliza el análisis de componentes principales para representar el frente de Pareto en el espacio de componentes principales. Posteriormente, se emplea una red neuronal para simular el comportamiento de este nuevo frente de Pareto que presenta la característica de ser convexo. Esta red neuronal permite establecer una relación no lineal entre los parámetros a estimar y los valores de los residuos. La toma de decisiones se resuelve a través de un algoritmo de minimización local que encuentre el mínimo de la función convexa definida por la red neuronal. Este proceso otorga al algoritmo colaborativo dos ventajas principales. La primera es que el tiempo de computación se reduce significativamente ya que son precisas menos iteraciones para obtener el frente de Pareto poco poblado. La segunda hace referencia a la robustez del método para resolver el problema de toma de decisiones y hallar el knee point del problema. El rendimiento de este algoritmo se ha validado a través de la actualización del modelo de EF de una pasarela comparando los resultados con aquellos obtenidos por el algoritmo MHS convencional y resolviendo el problema de toma de decisiones mediante el método de intersección normal al contorno. Sin comprometer la precisión en la solución, el algoritmo colaborativo reduce el tiempo de simulación en un 50%. 5.2 Líneas de investigación futuras 65 5.2 Líneas de investigación futuras El alcance de esta tesis puede extenderse a través de nuevos trabajos de investigación. El punto más inmediato concierne al desarrollo e implementación de un algortimo que combine los dos algoritmos propuestos en esta tesis. En este contexto, este nuevo algoritmo se caracterizaría por calcular un frente de Pareto poco poblado utilizado el algoritmo híbrido UKF-MHS y a partir de él continuar con el proceso del algoritmo colaborativo para acabar teniendo la solución óptima del problema. Otra posible futura aplicación de los algoritmos propuestos es la actualización de modelos de EF de alta fidelidad de estructuras de ingeniería civil. Estos modelos de EF se caracterizan por representar con un alto nivel de detalle la estructura en consideración. De este modo, la incertidumbre es menor y el proceso de actualización del modelo tiene menor coste computacional. Las aplicaciones de los algoritmos propuestos pueden ir más allá de la actualización de modelos de EF, como por ejemplo, estimar los parámetros de los dispositivos de control pasivo de vibraciones aplicando el método de diseño basado en el movimiento (motionbased design). Este trabajo, incluido en el Apéndice G, ha sido abordado por los autores considerando como algoritmo de optimización el GA para calcular los parámetros de los amortiguadores pasivos que mitigan las vibraciones de los cables de puentes atirantados [ 83 ]. Como continuación a esta investigación se propone realizar el diseño de amortiguadores aplicando el algoritmo colaborativo ya que los análisis transitorios involucrados aumentan significativamente el tiempo de cálculo. 72 Bibliografía [67] T. Liu, Q. Zhang, T. Zordan, and B. Briseghella, Finite element model updating of canonica bridge using experimental modal data and genetic algorithm, Structural Engineering International 26 (2016), no. 1, 27–36. [68] Y. Lu and Z. Tu, A two-level neural network approach for dynamic FE model updating including damping, Journal of Sound and Vibration 275 (2004), no. 3-5, 931–952. [69] F. Magalhães and Á. Cunha, Explaining operational modal analysis with data from an arch bridge, Mechanical systems and signal processing 25 (2011), no. 5, 1431–1450. [70] F. Magalhães, A. Cunha, E. Caetano, and R. Brincker, Damping estimation using free decays and ambient vibration tests, Mechanical Systems and Signal Processing 24 (2010), no. 5, 1274–1290. [71] D. Manjarres, I. Landa-Torres, S. Gil-Lopez, J. Del Ser, M.N. Bilbao, S. Salcedo- Sanz, and Z.W. Geem, A survey on applications of the harmony search algorithm, Engineering Applications of Artificial Intelligence 26 (2013), no. 8, 1818–1831. [72] C. Mares, J.E. Mottershead, and M.I. Friswell, Stochastic model updating: part 1—theory and simulated example, Mechanical systems and signal processing 20 (2006), no. 7, 1674–1695. [73] S. Mariani and A. Ghisi, Unscented kalman filtering for nonlinear structural dynamics, Nonlinear Dynamics 49 (2007), no. 1-2, 131–150. [74] T. Marwala, Finite element model updating using computational intelligence techniques: applications to structural dynamics, Springer, London, 2010. [75] T. Marwala, Finite-element-model updating using computional intelligence techniques, ch. Finite-element-model Updating Using Particle-swarm Optimization, pp. 67–84, Springer, London, 2010. [76] T. Marwala, I. Boulkaibet, and S. Adhikari, Probabilistic finite element model updating using bayesian statistics: applications to aeronautical and mechanical engineering, John Wiley & Sons, 2016. [77] F. Massa, T. Tison, and B. Lallemand, A fuzzy procedure for the static design of imprecise structures, Computer methods in applied mechanics and engineering 195 (2006), no. 9-12, 925–941. [78] Ansys Mechanical 19.0, Ansys inc. 2019, Available online: http:// www.ansys.com/ . [79] S.V. Modak, T.K. Kundra, and B.C. Nakra, Comparative study of model updating methods using simulated experimental data, Computers & structures 80 (2002), no. 5-6, 437–447. [80] S. Moradi, L. Fatahi, and P. Razi, Finite element model updating using bees algorithm, Struct Multidisc Optim 42 (2010), no. 2, 283–291. Bibliografía 73 [81] J.E. Mottershead and M.I. Friswell, Model updating in structural dynamics: a survey, Journal of sound and vibration 167 (1993), no. 2, 347–375. [82] J.E. Mottershead, M. Link, and M.I. Friswell, The sensitivity method in finite element model updating: a tutorial, Mechanical systems and signal processing 25 (2011), no. 7, 2275–2296. [83] J. Naranjo-Pérez, J.F. Jiménez-Alonso, I.M. Díaz, G. Quaranta, and A. Sáez, Motionbased design of passive damping systems to reduce wind-induced vibrations of stay cables under uncertainty conditions, Applied Sciences 10 (2020), no. 5, 1740. [84] L.T. Nguyen and T. Nestorović, The extended kalman filter and the unscented kalman filter for material parameter identification with application in tunneling, PAMM 13 (2013), no. 1, 393–394. [85] L.T. Nguyen and T. Nestorović, Nonlinear kalman filters for model calibration of soil parameters for geomechanical modeling in mechanized tunneling, Journal of Computing in Civil Engineering 30 (2016), no. 2, 04015025. [86] J. Nocedal and S. Wright, Numerical optimization, Springer Science & Business Media, 2006. [87] P. Pachón, V. Compán, and J.F. Jiménez-Alonso, Ambient vibration testing, dynamic identification and model updating of a historical building. chapel of the würzburg residence (germany), International Operational Modal Analysis Conference (IOMAC) (Gijón, Spain), 2015. [88] Y.S. Park, S. Kim, N. Kim, and J.J. Lee, Finite element model updating considering boundary conditions using neural networks, Engineering Structures 150 (2017), 511–519. [89] R. Perera, S.E. Fang, and A. Ruiz, Application of particle swarm optimization and genetic algorithms to multiobjective damage identification inverse problems with modelling errors, Meccanica 45 (2010), no. 5, 723–734. [90] S. Qin, Y. Zhang, Y.L. Zhou, and J. Kang, Dynamic model updating for bridge structures using the kriging model and pso algorithm ensemble with higher vibration modes, Sensors 18 (2018), no. 6, 1879. [91] G. Quaranta, A. Fiore, and G.C. Marano, Optimum design of prestressed concrete beams using constrained differential evolution algorithm, Structural and Multidisciplinary Optimization 49 (2014), no. 3, 441–453. [92] J.R. Raol, G. Girija, and J. Singh, Modelling and parameter estimation of dynamic systems, vol. 65, Iet, 2004. [93] W.X. Ren and H.B. Chen, Finite element model updating in structural dynamics by using the response surface method, Engineering structures 32 (2010), no. 8, 2455–2465. 74 Bibliografía [94] C.P. Robert and G. Casella, Monte carlo statistical methods, London: Springer, 2004. [95] R.S. Rosenberg, Simulation of genetic populations with biochemical properties, Ph.D. thesis, University of Michigan, 1967. [96] M.M. Saada, M.H. Arafa, and A.O. Nassef, Finite element model updating approach to damage identification in beams using particle swarm optimization, Engineering optimization 45 (2013), no. 6, 677–696. [97] M. Savoia, Structural reliability analysis through fuzzy number approach, with application to stability, Computers & Structures 80 (2002), no. 12, 1087–1102. [98] D. Shan, Q. Li, I. Khan, and X. Zhou, A novel finite element model updating method based on substructure and response surface model, Engineering Structures 103 (2015), 147–156. [99] E. Simoen, G. De Roeck, and G. Lombaert, Dealing with uncertainty in model updating for damage assessment: A review, Mechanical Systems and Signal Processing 56 (2015), 123–149. [100] N. Srinivas and K. Deb, Multi-objective function optimisation using non-dominated sorting genetic algorithm, Evolutionary Comp 2(1994), no. 3, 221–248. [101] A. Tarantola, Inverse problem theory and methods for model parameter estimation, vol. 89, SIAM, 2005. [102] A. Teughels and G. De Roeck, Damage detection and parameter identification by finite element model updating, Archives of Computational Methods in Engineering 12 (2005), no. 2, 123–164. [103] A. Teughels, G. De Roeck, and J.A.K. Suykens, Global optimization by coupled local minimizers and its application to fe model updating, Computers & structures 81 (2003), no. 24-25, 2337–2351. [104] A. Teughels, J. Maeck, and G. De Roeck, Damage assessment by fe model updating using damage functions, Computers & structures 80 (2002), no. 25, 1869–1879. [105] R. Van Der Merwe, Sigma-point kalman filters for probabilistic inference in dynamic state-space models, Ph.D. thesis, OGI School of Science & Engineering at OHSU, 2004. [106] R. Van Der Merwe and E.A. Wan, The square-root unscented kalman filter for state and parameter-estimation, 2001 IEEE international conference on acoustics, speech, and signal processing. Proceedings (Cat. No. 01CH37221), vol. 6, IEEE, 2001, pp. 3461–3464. [107] A. Vasuki, Nature-inspired optimization algorithms, CRC Press, 2020. Bibliografía 75 [108] E.A. Wan and R. Van Der Merwe, The unscented kalman filter for nonlinear estimation, Proceedings of the IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium (Cat. No. 00EX373), IEEE, 2000, pp. 153– 158. [109] E.A. Wan, R. Van Der Merwe, and A.T. Nelson, Dual estimation and the unscented transformation, Advances in neural information processing systems, 2000, pp. 666– 672. [110] X. Wang, X.Z. Gao, and K. Zenger, An introduction to harmony search optimization method, Springer, 2015. [111] Y.Q. Wang, L. Zong, Y.J. Shi, and N. Yao, Damage detection and rehabilitation on a curvilinear steel box girder bridge by multistage model updating, Structure and Infrastructure Engineering 11 (2015), no. 11, 1420–1431. [112] S. Weng, H. Zhu, Y. Xia, J. Li, and W. Tian, A review on dynamic substructuring methods for model updating and damage detection of large-scale structures, Advances in Structural Engineering 23 (2020), no. 3, 584–600. [113] J. Wu, Q. Yan, S. Huang, C. Zou, J. Zhong, and W. Wang, Finite element model updating in bridge structures using kriging model and latin hypercube sampling method, Advances in Civil Engineering 2018 (2018). [114] X.S. Yang, Nature-inspired optimization algorithms, Elsevier, 2014. [115] X.S. Yang and S. Deb, Engineering optimisation by cuckoo search, Int. J. Mathematical Modelling and Numerical Optimisation 1(2010), no. 4, 330–343. [116] L.A. Zadeh, Fuzzy sets, Information and control 8(1965), no. 3, 338–353. [117] S. Živanović, A. Pavic, and P. Reynolds, Finite element modelling and updating of a lively footbridge: The complete process, Journal of Sound and Vibration 301 (2007), no. 1-2, 126–145. Apéndice A Artículo A Jiménez-Alonso, J.F, Naranjo-Pérez, J., Paviv, A., Sáez, A. Maximum likelihood finiteelement model updating of civil engineering structures using nature-inspired computational algorithms. Structural Engineering International. 2020. (Accepted). Revista: Structural Engineering International ISSN: 1683-0350 JCR (Journal Citation Reports) (2018): Factor de impacto: 0.608 •Civil Engineering: Q4 (120/132) •Construction and Building Technology: Q4 (54/63) 77 Maximum Likelihood Finite-Element Model Updating of Civil Engineering Structures Using Nature-Inspired Computational Algorithms Javier Fernando Jiménez-Alonsoa*, Javier Naranjo-Perezb, Aleksandar Pavicc and Andrés Sáezd a*Assistant professor. Department of Continuum Mechanics and Structures, E.T.S. Ingenieros de Caminos, Canales y Puertos. Universidad Politécnica de Madrid, Madrid (Spain). bPhD Candidate. Department of Continuum Mechanics and Structural Analysis, Universidad de Sevilla, Seville (Spain). cFull professor. Vibration Engineering, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter (U.K.). dFull professor. Department of Continuum Mechanics and Structural Analysis, Universidad de Sevilla, Seville (Spain). *Corresponding author: Assistant Professor: Javier Fernando Jiménez-Alonso. Department of Continuum Mechanics and Structures. E.T.S. Ingenieros de Caminos, Canales y Puertos. Universidad Politécnica de Madrid, Calle del Profesor Aranguren, 3, 28040 Madrid (Spain) Ph:+34 91 0674154. e-mail: [email protected]. Maximum Likelihood Finite-Element Model Updating of Civil Engineering Structures Using Nature-Inspired Computational Algorithms In finite-element model updating of civil engineering structures, the numerical finite-element model is calibrated in order to fit the experimentally measured dynamic properties (say mode shapes and natural frequencies) of the structure. Such updating process is usually performed under the maximum likelihood method in practical engineering applications. According to this approach, the updating problem is transformed into an optimization problem, where the values of the most relevant physical parameters of the structure are obtained via the minimization of an objective function. This function is usually defined in terms of the relative differences between the numerical and the experimental modal properties of the structure. To this aim, either (i) a single-objective or (ii) a multi-objective approach may be adopted. In the majority of applications, and due to the complexity of the problem, global optimizers are considered. Among these algorithms, natureinspired computational algorithms have been widely employed with satisfactory results. Nevertheless, such model updating approach presents two main limitations: (i) a clear dependence between the updated model and the objective function considered for the updating problem; and (ii) a high simulation cost in terms of computation time. In order to overcome these drawbacks, a detailed study has been performed herein both to assist in the selection of an efficient computational algorithm among several well-known nature-inspired algorithms and to further establish the most adequate (single- or multi-) objective function to tackle the updating problem. For this purpose, a laboratory footbridge has been considered as benchmark to conduct the updating process under different scenarios. Keywords: finite-element model updating, maximum likelihood method, singleobjective optimization, multi-objective optimization, nature-inspired computational algorithms, harmony search, civil engineering structures. Introduction Finite-element (FE) models are extensively used to simulate numerically the behaviour of civil engineering structures. However, these numerical models may not always reflect adequately the actual behaviour of the structure. To solve this issue, such models can be tuned based on the experimental results obtained from either field tests or continuous monitoring [1]. For practical engineering applications, the tuning of the numerical models is usually achieved via the modification of the most relevant physical parameters of the structure, in order to reduce the differences between its numerical and experimental behaviours. After the calibrating process, the resulting updated model can be used either to assess more accurately the behaviour of civil engineering structures [2, 3, 4] or to establish strategies for their maintenance or damage detection [5, 6]. Thus, given that one of the main objectives of FE model updating is to indirectly estimate the values of some relevant physical parameters of the structure, the updating process may be formulated as a parameter identification problem [7], where appropriate estimators need to be employed [8]. Generally speaking, estimators may be classified into two categories [9, 10]: (i) Point estimators, which return the expected value of each considered design parameter; and (ii) Interval estimators, which determine either an interval in which the value of each parameter lies or a probability density function for each parameter. Among the interval estimators, the Bayesian method [11] has prevailed, since the determination of the probabilistic density function of the design parameters is relevant when performing subsequent structural reliability analyses [12]. Meanwhile, among the point estimators, the maximum likelihood method (MLM) has been widely implemented due to its proved efficiency and accuracy when tackling model updating [8]. In this paper, we will focus on analysing the performance of FE model updating under the MLM, based on the advantages this approach exhibits for practical civil engineering applications, when compared to Bayesian methods, in particular in terms of: (i) ease of implementation, (ii) lower simulation times required; and (iii) straightforward use of the updated models for deterministic applications [2, 13]. In this manner, the parameter identification problem under the MLM, assuming a normal distribution of the adjustment errors, is equivalent to an ordinary least squares problem [14]. Therefore, the FE model updating aims at solving an optimization problem, whose objective is to determine the values of the design parameters that minimize the difference between the numerical (FE) modal properties of the structure and its experimental counterparts. The experimental data follow from the signal processing (applying either experimental or operational modal analysis) of a forced or an ambient vibration test [15]. At this stage, the formulation of the optimization problem requires to adopt two key decisions: (i) the definition of an adequate objective function; and (ii) the selection of an efficient computational optimization algorithm. On the one hand, the objective function is usually defined in terms of the residuals between the numerical and experimental modal properties of the structure. Since two types of modal properties are normally considered (natural frequencies and associated vibration modes), two types of residuals come into play. As both sets of modal properties have different nature, a new problem arises, namely, how to weigh the influence of each residual. Two approaches may be adopted, depending on whether a single-objective function (SOF) or a multi-objective function (MOF) [16] is considered. Whilst in the SOF approach a single objective is optimized with different weighting factors assigned to balance the influence of each residual, in the MOF approach the different terms of a multiobjective function are optimized (without the need to define weighting factors). Although both approaches have been used interchangeably for practical engineering applications [17, 18], the number of studies in which their performance has been compared is scarce. On the other hand, the selection of an appropriate optimization algorithm is a key ingredient for the success of the FE model updating. Due to the nonlinear relationship between the modal properties of the structure and the considered physical parameters, computational intelligence algorithms are usually considered to ensure the finding of the so-called Pareto front defining a set of possible solutions of the updating problem, is obtained. a) GA b) PS c) HS Fig. 2. Flowchart of the three considered NIC algorithms: a) genetic algorithms (GA); b) particle swarm (PS); and c) harmony search (HS). Similarly, the MOF approach of the PS algorithm is also based on the classification of the non-dominated particles [32]. Hence, once the initial particles’ population is created and the multi-objective function is evaluated, the non-dominated solutions are obtained and stored in the particles’ repository. A similar process, to define the Pareto front, is obtained after the mutation of the position. Thus, if the mutated particle dominates the non-mutated, the latter is updated. This step is applied iteratively to both the updated particle’s position and the best particle’s position. Then, the non-dominated solutions are added to the repository and finally, the domination is evaluated for the overall repository to dismiss the dominated particles. These steps are repeated iteratively until some convergence criterion is met. As results of this process, the Pareto font is obtained. Finally, the MOF approach of the HS algorithm has also been implemented [33] considering the main rules of the NSGA-II algorithm. In this manner, the HS algorithm under a MOF approach consists of the following steps: (i) an initial set of harmonies is generated (harmony memory matrix, 𝐇𝐌); (ii) the multi-objective function is evaluated for all of harmonies; (iii) the attributes of each individual are calculated using the crowded comparison operator; (iv) the harmonies are sorted based on the crowding distance; (v) a new set of harmonies is created based on the three mentioned design variables (𝐻𝑀𝐶𝑅, 𝑃𝐴𝑅 and 𝑏𝑤); (vi) the multi-objective function is evaluated for the new harmonies; (vii) the harmony memory matrix, 𝐇𝐌, is truncated using both the non-dominated sorting technique and the crowding distance; and (viii) the steps (v) to (vii) are repeated until some convergence criterion is met. Finally, as result of this process, the Pareto front is obtained. Initialize Population Evaluate Select Crossover Mutate Convergence? Stop Yes No Operations Update Population Convergence? Stop Yes No Position Velocity Operations Evaluate Initialize Population Initialize HM Convergence? Stop Yes No Operations Evaluate Select Improvise Random Update HM Preliminary discussion on the Performance of the Three Considered Nature-Inspired Algorithms As Fig. 2 illustrates, the three considered NIC algorithms under the SOF approach share three common characteristics: (i) the random initialization of a population; (ii) the implementation of several mathematical operations to modify this population; and (iii) the iterative application of this re-generation process until some convergence criteria is met. Similarly, they also share the use of the non-dominated sorting rules for their implementation under a MOF approach. However, some particular differences among these NIC algorithms make interesting comparing their performance when they are implemented for the FE model updating of civil engineering structures. Between the GA and the HS algorithm, the main difference lies in the different complexity of the operations that control the updating of the population: the operations associated with the HS algorithm (selection, improvement and randomization) are computationally easier than the operations associated with the GA (selection, crossover and mutation). Thus, it is expected that the computational time required to solve an optimization problem using the HS algorithm is lower than the one required employing GA. However, the simplicity of the operations can have a negative influence on both the accuracy of the solution and the sweeping capacity of the algorithm, reducing the performance of the HS algorithm when it is implemented for a practical engineering application. Additionally, among these two algorithms (GA and HS) and the PS algorithm, the main difference is related to the re-generation of the population. While the two mentioned algorithms (GA and HS) generate iteratively different populations (new populations) during the updating process; the PS algorithm maintains the initial population and it only updates iteratively its attributes. In this manner, the PS algorithm may reduce the simulation time required to solve an optimization problem. However, as on the one hand, the sweeping capacity of the algorithm is lower than the GA; and on the other hand, the computational operations needed to compute the position and the velocity of each particle are more complex than the re-generating operations of the HS algorithm, it is not possible to establish a priori whether the PS algorithm is the most efficient computational algorithm among the three considered NIC ones. For all these reason, it is not possible to preliminary conclude which will be the most efficient algorithm to solve a particular optimization problem, being necessary to perform a comparative study to analyse their performance. This will be done in the next section. Application Example: FE Model Updating of a Laboratory Footbridge A real case-study is next analysed to shed some light onto the two objectives of this paper: (i) establishing the most adequate approach (SOF or MOF) to formulate the updating problem under the MLM; and (ii) determining the most efficient NIC algorithm among the three selected above (GA, PS and HS). A comparison of the performance of the different approaches when applied to the FE model updating of a benchmark laboratory footbridge is conducted in this section. Preliminary FE Model and Numerical Modal Analysis of the Laboratory Footbridge The laboratory footbridge adopted for this study is a frame structure with a single span of 15 m. The structure is configured by two lateral steel beams separated transversally 2.5 m. These beams are braced by diaphragms separated longitudinally 1.25 m. The diaphragms consist of rectangular plates of 200x12 mm. The lateral beams are connected at both sides to steel columns. These columns are directly pinned to the ground. The deck of the structure is formed by composite SPS panels [34]. These composite SPS panels are bolted to both the longitudinal and transversal elements. An overall view of the structure and its main constitutive components is illustrated in Fig. 3. A more detailed description of this benchmark structure may be found in Ref [25]. The structure was modelled with three types of elements in the commercial FE package Ansys [35]. All the numerical simulations were performed using a laptop computer with a processor of 3.6 GHz and a RAM memory of 8 GB. The mesh consisted of 31903 elements (Fig. 5). The FE model has been realized according to the following scheme: (i) the two main lateral beams and the transversal plates have been modelled with a four node shell elements with six degrees of freedom at each node (SHELL181); (ii) the SPS panels have been modelled with the same type of element (SHELL181), considering its sandwich-type behaviour through the first-order shear-deformation theory; (iii) the bolts that configure the connection between the SPS panel and the steel structure have been modelled by 3-D beam elements (BEAM188); and (iv) the supports were modelled by means of longitudinal and lateral spring elements (COMBIN14) and assuming that the vertical displacement was constrained. The initial value for these spring parameters has been estimated from a simplified FE model of just the column that configure the support, resulting an equivalent stiffness of 𝑘𝑙𝑜𝑛=5.5∙107 N/m in the longitudinal direction and 𝑘𝑙𝑎𝑡=1.9∙107 N/m in the lateral direction. The mechanical properties of the constitutive materials adopt the following values: (i) for the steel [36], a density 𝛾𝑠=7850 kg/m3, a Young’s modulus 𝐸𝑠=2.1∙ 105 MPa and a Poisson’s ratio 𝜈𝑠=0.3; and (ii) for the polyurethane [37], a density 𝛾𝑝= 1100 kg/m3, a Young’s modulus 𝐸𝑝=750 MPa and a Poisson’s ratio 𝜈𝑝=0.5. This FE model has been used to perform a numerical modal analysis: The obtained first seven natural frequencies, 𝑓𝑛𝑢𝑚,𝑗(𝛉), are given in Table 1; whilst the associated numerical vibration modes, 𝝓𝑛𝑢𝑚,𝑗(𝛉), are depicted in Fig. 3 ( 𝑗 being the considered vibration mode). Laboratory Footbridge 𝑓𝑛𝑢𝑚,1(𝛉)=3.638 Hz SPS Panel Support Transversal Plate Longitudinal Beam Connection (Bolts) 2.5 m 𝑓𝑛𝑢𝑚,2(𝛉)=5.329 Hz 𝑓𝑛𝑢𝑚,3(𝛉)=10.185 Hz 𝑓𝑛𝑢𝑚,4(𝛉)=11.310 Hz 𝑓𝑛𝑢𝑚,5(𝛉)=17.364 Hz 𝑓𝑛𝑢𝑚,6(𝛉)=20.238 Hz 𝑓𝑛𝑢𝑚,7(𝛉)=21.105 Hz Fig. 3. Description of the structure and initial first seven numerical vibration modes, 𝝓𝑛𝑢𝑚,𝑗(𝛉), being 𝑗 the considered vibration mode. Experimental Identification of the Modal Properties of the Laboratory Footbridge Subsequently, a forced vibration test was conducted to identify experimentally the modal properties of the footbridge. A set of proof-mass actuators and rowing accelerometers were used for this purpose (Fig. 4). The actuators were simultaneously driven by uncorrelated random signals considering a Multi-input Multi-output (MIMO) configuration [15, 25]. At all the considered reference points (instrumented points), the Frequency Response Functions (FRFs) were determined considering an overlap of 50%. Later, probable mode locations were identified in the FRFs fitted curves via a complex mode indicator function. Finally, a global polynomial curve fitting method was used to identify experimentally the first seven natural frequencies and associated vibration modes based on the mentioned probable mode locations [15]. Both the forced vibration test and the corresponding experimental modal analysis are described with more detail in Ref [25]. Fig. 4. Forced vibration test of the laboratory footbridge [25]. The first seven experimental natural frequencies (Table 1), 𝑓𝑒𝑥𝑝,𝑗, and associated experimental vibration modes (Fig. 5), 𝝓𝑒𝑥𝑝,𝑗, were obtained as result of this experimental modal analysis (being 𝑗 the considered vibration mode). Correlation between the Preliminary Numerical and Experimental Modal Properties of the Laboratory Footbridge In order to analyse the correlation between the numerical and experimental modal properties of the laboratory footbridge, the relative differences, ∆𝑓𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗(𝛉), and the 𝑀𝐴𝐶(𝛉)𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗 ratios were determined. Table 1 shows the correlation between the numerical and experimental modal properties of the laboratory footbridge. Additionally, Fig. 5 illustrates a comparison between the numerical, 𝝓𝑛𝑢𝑚,𝑗(𝛉), and experimental, 𝝓𝑒𝑥𝑝,𝑗, vibration modes. Table 1. Correlation between the initial numerical (FE) and experimental modal properties of the laboratory footbridge (being 𝑓𝑛𝑢𝑚,𝑗(𝛉) the numerical natural frequencies; 𝑓𝑒𝑥𝑝,𝑗 the experimental natural frequencies; ∆𝑓𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗(𝛉), the relative differences between the numerical and experimental natural frequencies; and the 𝑀𝐴𝐶(𝛉)𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗 ratio of the 𝑗 considered vibration mode). Vibration Mode 𝑓𝑛𝑢𝑚,𝑗(𝛉) [Hz] 𝑓𝑒𝑥𝑝,𝑗 [Hz] ∆𝑓𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗(𝛉) [%] 𝑀𝐴𝐶(𝛉)𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗 [-] 1 (1st Bending) 3.638 3.810 -4.509 0.999 2 (1st Torsional) 5.329 5.144 3.600 0.994 3 (2nd Torsional) 10.185 8.485 20.033 0.990 4 (2nd Bending) 11.310 12.366 -8.540 0.877 5 (3rd Bending) 17.364 18.605 -6.670 0.985 6 (3rd Torsional) 20.238 20.459 -1.080 0.993 7 (4th Bending) 21.105 22.980 -8.159 0.910 Actuator Y X ZY ZX Forced Vibration Test Accelerometer 15 m 2.5 m Table 1 shows that, although the shapes of six of the seven identified vibration modes were in good agreement (with 𝑀𝐴𝐶(𝛉)𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗 ratios greater than 0.90), the relative differences, ∆𝑓𝑒𝑥𝑝,𝑗 𝑛𝑢𝑚,𝑗(𝛉), between the numerical and experimental natural frequencies of four of the seven vibration modes were too large to be acceptable. Therefore, it becomes necessary to improve the FE model of this structure in order to ensure a better correlation between its numerical and experimental modal properties and, consequently, to obtain a valid FE model that reflects the actual dynamic behaviour of this structure [5]. FE Model Updating of the Laboratory Footbridge In order to formulate the FE model updating problem of this laboratory footbridge under the MLM, the following elements must be defined: (i) the objective function; (ii) the design parameters; (iii) the search domain; and (iv) the optimization algorithm. Whilst the two approaches for the objective function (SOF/MOF) and the considered optimization algorithms (GA/PS/HS) were described previously, both the design parameters and their corresponding search domains will be established here. Hence, a set with the most relevant physical parameters of the model needs to be selected, 𝛉. To this end, the ratio between modal strain energy associated with each physical parameter and the overall modal strain energy of the structure was considered as selection criterion [38]. Following the sensitivity analysis carried out for this purpose, ten physical parameters of the structure were selected as design parameters (Fig. 7), 𝛉= [𝜃1,𝜃2,𝜃3,𝜃4,𝜃5,𝜃6,𝜃7,𝜃8,𝜃9,𝜃10]: (i) the Young’s modulus of the steel of the longitudinal beams in different sections, 𝜃1−𝜃6 [MPa]; (ii) the Young’s modulus of the polyurethane, 𝜃7 [MPa]; (iii) the Young’s modulus of the steel of the bolts, 𝜃8 [MPa]; (iv) the equivalent longitudinal stiffness of each support, 𝜃9 [N/m]; and (v) the equivalent lateral stiffness of each support, 𝜃10 [N/m]. Additionally, a search domain has been included to constrain this optimization problem and thus, to guarantee the physical meaning of the updated model. Thus, the lower, 𝛉𝒍, and the upper, 𝛉𝒖, bounds of this search domain were defined as follows, 𝛉𝒍=[1.89∙105,1.89∙105,1.89∙105,1.89∙ 105,1.89∙105,1.89∙105,0.75∙103,2.1∙105,1.38∙107,4.75∙106 ] and 𝛉𝒖=[2.31∙ 105,2.31∙105,2.31∙105,2.31∙105,2.31∙105,2.31∙105,1.5∙103,2.1∙106,1.1∙ 108,3.8∙107 ]. FEM 1st vibration mode X X Exp. Num. 2nd vibration mode 3rd vibration mode 4th vibration mode 5th vibration mode 6th vibration mode 7th vibration mode Fig. 5. FE model and comparison between the numerical, 𝝓𝑛𝑢𝑚,𝑗(𝛉), and experimental, 𝝓𝑒𝑥𝑝,𝑗, vibration modes of the laboratory footbridge (being 𝑗 the considered vibration mode). Subsequently, the FE model updating was conducted considering the two approaches and the three considered NIC algorithms. Comparison among GA, PS and HS for FE Model Updating under the SOF Approach First, the initial FE model of the footbridge was updated using the SOF approach and considering the three mentioned NIC algorithms (GA, PS and HS). The solution of each updating problem was performed via its implementation using both the FE analysis package Ansys [35] and the mathematical package Matlab [39]. Two comparison criteria were considered to assess the performance of the three NIC algorithms: (i) the computational cost; and (ii) the accuracy of the adjustment. For the three NIC algorithms, the following common design variables were adopted: (i) population size, 𝑃𝑜𝑝=100; (ii) maximum number of iterations, 𝐼𝑡𝑒𝑟=50; X X Exp. Num. X X Exp. Num. X X Exp. Num. X X Exp. Num. X X Exp. Num. X X Exp. Num. and (iii) objective function tolerance, 𝑇𝑜𝑙𝐹=10−5 (so that the algorithm stops if the average relative change of the best value of the objective function is less than or equal to 𝑇𝑜𝑙𝐹). Additionally, for each NIC algorithm, the following particular design variables were considered: (i) for the GA, as selection ratio, 𝑠𝑒𝑙𝑟=0.1, as crossover ratio, 𝑐𝑟𝑜𝑠𝑟= 0.9, and as mutation ratio, 𝑚𝑢𝑡𝑟=0.1; (ii) for the PS algorithm, as particle’s inertia, 𝑖𝑝𝑎𝑟=0.95, as self-adjustment weight, 𝑠𝑒𝑙𝑤=0.05, and as social-adjustment weight, 𝑠𝑜𝑐𝑤=0.02; and (iii) for the HS algorithm, as new population size, 𝑁𝑒𝑤_𝑃𝑜𝑝=40, as harmony memory consideration rate, 𝐻𝑀𝐶𝑅=0.9, as pitch adjustment rate, 𝑃𝐴𝑅=0.3, and as bandwidth, 𝑏𝑤=1 100(𝛉𝐮−𝛉𝐥). In order to assess the impact of the residual weights on the updated value of the single-objective function, 𝑓(𝛉), a trial-and-error criterion [5] was adopted. Table 2 shows the eleven cases established in terms of the different values considered for the residuals associated with the natural frequencies, ∑𝑤𝑗𝑓, and vibration modes., ∑𝑤𝑗𝑚, respectively. Additionally, each case was run ten times to take into account the random generation of the initial population. As result of the updating process, Table 2 shows the average values of the residuals for the three considered algorithms in terms of the different value of the weighting factors. Additionally, Fig. 6a illustrates the representation in a multi-objective functional space of the residuals associated with each case and computational algorithm. The average values of the updated parameters obtained using the three considered NIC algorithms under the SOF approach are shown in Table 3. Table 2. Average values of the residuals in terms of the weighting factors under the SOF approach for the three considered NIC algorithms (GA, PS and HS). Case ∑𝑤𝑗𝑓 ∑𝑤𝑗𝑚 GA [𝑥10−3] PS [𝑥10−3] HS [𝑥10−3] ∑𝑟𝑗𝑓(𝛉)2 ∑𝑟𝑗𝑚(𝛉)2 ∑𝑟𝑗𝑓(𝛉)2 ∑𝑟𝑗𝑚(𝛉)2 ∑𝑟𝑗𝑓(𝛉)2 ∑𝑟𝑗𝑚(𝛉)2 01 1.0 0.0 3.63 2.69 3.56 2.68 3.55 2.69 02 0.9 0.1 3.57 2.66 3.56 2.68 3.55 2.69 03 0.8 0.2 3.59 2.62 3.58 2.51 3.58 2.66 04 0.7 0.3 3.58 2.53 3.57 2.52 3.61 2.52 05 0.6 0.4 3.57 2.56 3.57 2.50 3.62 2.50 06 0.5 0.5 3.58 2.51 3.58 2.51 3.58 2.57 07 0.4 0.6 3.58 2.52 3.58 2.49 3.62 2.51 08 0.3 0.7 3.59 2.48 3.58 2.49 3.63 2.51 09 0.2 0.8 3.58 2.51 3.58 2.49 3.62 2.49 10 0.1 0.9 3.58 2.53 3.59 2.47 3.70 2.47 11 0.0 1.0 3.59 2.50 3.59 2.49 3.64 2.46 Finally, the average simulation time required to solve the updating problem for each case was: (i) 74072 s for the GA; (ii) 43857 s for the PS; and (iii) 34380 s for the HS. According to these results, two main conclusions may be drawn: (i) the accuracy of the three mentioned NIC algorithm is similar; and (ii) the HS algorithm allows reducing the simulation time required to perform the updating process. In this manner, the HS algorithm allows performing the FE model updating of this structure more efficiently than the other algorithms, since it reduces the computational cost without compromising the accuracy of the adjustment. Comparison among GA, PS and HS for the FE Model Updating under the MOF Approach Subsequently, the initial FE model of the laboratory footbridge was also updated under the MOF approach and considering the three mentioned NIC algorithms. Once more, the updating process was performed via its implementation using both Ansys [35] and Matlab [39] software packages and considering the same comparison criteria as above. For the three NIC algorithms under the MOF approach, the same design variables (common and particular) as those employed for the SOF case were adopted. Additionally, the number of elements of the Pareto front was established in 25. The updating problem was solved ten times for each NIC algorithms to take into account the random generation of the initial population. As result of the updating process, Fig. 6b illustrates the average Pareto front obtained for each considered NIC algorithm. Additionally, the average value of the updated parameters for the three considered NIC algorithms under the MOF approach is given in Table 3 (the NBI method [24] has been adopted for this purpose). 2.20E-03 2.30E-03 2.40E-03 2.50E-03 2.60E-03 2.70E-03 3.50E-03 4.00E-03 4.50E-03 5.00E-03 5.50E-03 6.00E-03 6.50E-03 7.00E-03 7.50E-03 f2 f1 SOF-GA SOF-PS SOF-HS 2.20E-03 2.30E-03 2.40E-03 2.50E-03 2.60E-03 2.70E-03 3.50E-03 4.00E-03 4.50E-03 5.00E-03 5.50E-03 6.00E-03 6.50E-03 7.00E-03 7.50E-03 f2 f1 MOF-GA MOF-PS MOF-HS a) b) Fig. 6. Comparison among the different residuals in the functional space considering the three NIC algorithms (GA, PS and HS): a) the SOF approach; b) the MOF approach; and c) both the SOF and MOF approaches. Additionally, the average simulation time required to obtain the Pareto front corresponding to each mentioned NIC algorithm has been recorded: (i) 307057 s for the GA; (ii) 240663 s for the PS; and (iii) 121676 s for the HS. Table 3. Average updated value of the physical parameters of model, 𝛉, after the updating process, considering the two approaches and the three mentioned computational algorithms. 𝛉 SOF MOF GA PS HS GA PS HS 𝜃1 2.22∙105 2.23∙105 2.23∙105 2.28∙105 2.29∙105 2.29∙105 𝜃2 2.11∙105 2.16∙105 2.09∙105 2.12∙105 2.25∙105 2.22∙105 𝜃3 2.09∙105 2.10∙105 2.09∙105 2.02∙105 2.19∙105 2.13∙105 𝜃4 2.07∙105 2.09∙105 2.05∙105 2.18∙105 2.19∙105 2.19∙105 𝜃5 2.10∙105 2.04∙105 2.12∙105 2.15∙105 2.15∙105 2.17∙105 𝜃6 2.14∙105 2.09∙105 2.12∙105 2.17∙105 2.12∙105 2.15∙105 𝜃7 1.43∙103 1.43∙103 1.42∙103 1.02∙103 1.09∙103 1.03∙103 𝜃8 4.40∙105 4.33∙105 4.29∙105 2.90∙105 2.96∙105 2.83∙105 𝜃9 4.26∙107 4.24∙107 4.24∙107 4.29∙107 4.12∙107 4.13∙107 𝜃10 1.50∙107 1.49∙107 1.49∙107 1.45∙107 1.42∙107 1.43∙107 According to these results, three main conclusions can be drawn: (i) the accuracy of both the PS and HS algorithm is greater than the one provided by the GA since the Pareto front obtained by these two algorithms allows reducing the overall value of the residuals between the numerical and experimental modal properties of the structure; (ii) the accuracy between the PS and HS algorithms is similar; and (iii) the HS algorithm allows reducing again the simulation time required to perform the updating process. In this manner, the HS algorithm is again the most efficient NIC algorithm to solve the updating problem under the MOF approach. 2.20E-03 2.30E-03 2.40E-03 2.50E-03 2.60E-03 2.70E-03 3.50E-03 4.00E-03 4.50E-03 5.00E-03 5.50E-03 6.00E-03 6.50E-03 7.00E-03 7.50E-03 f2 f1 SOF-GA SOF-PS SOF-HS MOF-GA MOF-PS MOF-HS c) Apéndice B Articulo B Naranjo-Pérez, J., Sáez, A., Jiménez-Alonso, J.F., Pachón, P., Compán, V. (2019). A hybrid UKF-MGA algorithm for finite element model updating of historical constructions. IABSE Symposium, Guimaraes 2019: Towards a Resilient Built Environment Risk and Asset Management – Report pp. 29-36. ISBN: 978-385748163-5 105 IABSE Symposium 2019 Guimarães Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 1 A Hybrid UKF-MGA Algorithm for Finite Element Model Updating of Historical Constructions Javier Naranjo-Pérez, Andrés Sáez Department of Continuum Mechanics and Structural Analysis, Universidad de Sevilla, Seville, Spain Javier F Jiménez-Alonso, Pablo Pachón, Víctor Compán Department of Building Structures and Geotechnical Engineering, Universidad de Sevilla, Seville, Spain Contacting author: [email protected] Abstract The finite element model (FE) updating is a calibration method that allows minimizing the discrepancies between the numerical and experimental modal parameters. As result, a more accurate FE model is obtained and the structural analysis can represent the real behaviour of the structure. However, it is a high computational cost process. To overcome this issue, alternative techniques have been developed. This study focuses on the use of the unscented Kalman filter (UKF), which is a local optimization algorithm based on statistical estimation of parameters taken into account the measurements. The dome of a real chapel is considered as benchmark structure. A FE model is updated applying two different algorithms: (i) the multi-objective genetic algorithm and (ii) a hybrid unscented Kalman filter-multi-objective genetic algorithm (UKF-MGA). Finally, a discussion of the results will be presented to compare the performance of both algorithms. Keywords: model updating; historical constructions; unscented Kalman filter; multi-objective genetic algorithm. 1. Introduction The structural assessment of historical constructions is an increasingly relevant issue. The most usual approach for the assessment of their structural behaviour is to design a mathematical model (e.g. finite element model). To build the Finite Element (FE) model, the material and geometrical properties are assumed and the boundary conditions implemented may not be the proper. These facts may involve large discrepancies between the numerical results obtained from these mathematical models and the experimental results obtained from field vibration tests. The model updating may be formulated as an optimization procedure where these discrepancies are minimized. Friswell et al. [1] introduced in 1995 the main techniques which may be employed for FE model updating. The procedure consists in selecting one or several physical parameters of the numerical model and changing theirs values iteratively. For each iteration, the differences between numerical and experimental results (usually modal parameters) are evaluated. In this manner, viewed from the optimization perspective, the parameter identification solved as an inverse problem can be considered as a general minimization problem. IABSE Symposium 2019 Guimarães: Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 2 The objective is to find the physical parameters that minimizes the mean square error (objective function) between actual and simulated modal parameters. Due to the high nonlinearity intrinsic to structural models, the objective functions could present many local minima. To prevent the algorithm to converge to a local minimum, the use of global optimization algorithms must be considered. Genetic algorithm (GA) has been widely used to solve a parameter identification inverse problem. The high computational cost has led to the emergence of hybrid techniques that reduce the convergence time of global optimization algorithms. In this study, a hybrid unscented Kalman filter-multi-objective genetic algorithm (UKF-MGA) is considered. The unscented Kalman filter (UKF) [6] is an algorithm used to estimate the parameters of a nonlinear system subjected to white noise. The parameters updating is based on the estimation of the Gaussian means and covariances of the quantities being estimated, minimizing the variance of the estimation error. The Kalman filter finds the local minimum around the current estimate. This study focuses on the comparison of the results obtained using both the maximum likelihood approach and a hybrid UKF-MGA algorithm. As benchmark structure, the chapel of the Würzburg Residence (Germany) is considered. This building, whose construction dates from the early eighteenth century, was declared a World Heritage Site by UNESCO in 1981. The updating procedure is based on the actual modal parameters of the structure, obtained through the accelerations recorded during an ambient vibration test. Finally, the number of iterations, the computation time and the differences between the updated and experimental modal parameters are compared for both algorithms. 2. FE Model Updating based on the maximum likelihood method 2.1 Basics of FE Model Updating The FE model updating tool aims to design a numerical model which is more adequate to the real behaviour of the structure [1]. It is an iterative process in which the values of certain pre-selected physical parameters are changed until the optimal solution is reached. The updated FE model, under the maximum likelihood method, is obtained by minimizing the objective functions defined in terms of the relative differences between experimental and numerical modal parameters. Within the framework of multi-objective optimization, the problem to be solved consists in minimizing the objectives functions defined with regard to the residual of the natural frequencies and the residual of the vibration modes. The expressions to calculate both residuals are: 𝑟𝑓,𝑗=𝑓𝑛𝑢𝑚,𝑗−𝑓𝑒𝑥𝑝,𝑗 𝑓𝑒𝑥𝑝,𝑗 , 𝑗=1,…,𝑚𝑓 (1) 𝑟𝑠,𝑗 2=(1−√𝑀𝐴𝐶𝑗)2 𝑀𝐴𝐶𝑗, 𝑗=1,…,𝑚𝑓 (2) where 𝑓𝑒𝑥𝑝,𝑗 is the experimental natural frequency 𝑗, 𝑓𝑛𝑢𝑚,𝑗 is the numerical natural frequency 𝑗 and 𝑀𝐴𝐶𝑗 is the Modal Assurance Criterion. The MAC is a parameter used to assess the correlation between the numerical and experimental natural frequency 𝑗 and is calculated by the following expression: 𝑀𝐴𝐶𝑗=(𝜙𝑛𝑢𝑚,𝑗 𝑇·𝜙𝑒𝑥𝑝,𝑗)2 (𝜙𝑛𝑢𝑚,𝑗 𝑇·𝜙𝑛𝑢𝑚,𝑗)·(𝜙𝑒𝑥𝑝,𝑗 𝑇·𝜙𝑒𝑥𝑝,𝑗) (3) being 𝜙𝑛𝑢𝑚 and 𝜙𝑒𝑥𝑝 the numerical and experimental vibration modes, respectively. Both residuals are used to define the objective functions. Therefore, the minimization problem may be defined as follows: min𝑓1(𝜃)=min1 2[∑𝑟𝑓,𝑗 2(𝜃) 𝑚𝑓 𝑗=1 ]1 2 (4) min𝑓2(𝜃)=min1 2[∑𝑟𝑠,𝑗 2(𝜃) 𝑚𝑓 𝑗=1 ]1 2 (5) IABSE Symposium 2019 Guimarães: Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 3 where 𝑓1(𝜃) and 𝑓2(𝜃) are the objective functions related to the natural frequencies and vibration modes, respectively, 𝜃 is the possible solution vector containing the updated physical parameters of the model and 𝑚𝑓 is the number of vibration modes considered. 2.2 Multi-objective Genetic Algorithm (MGA) The multi-objective global optimization algorithm (MGA) employed in this study is the natureinspired NSGA-II (Non-Dominated Sorting Genetic Algorithm) algorithm [3]. This algorithm consists of the followings steps. First, the initialization consists in creating an initial random population (solutions) and the assessment of the objective functions. The value of the objective functions are then classified to distinguish between dominated and non-dominated (Pareto front) solutions. Second, the stochastic selection function selects several classified solutions and generate a new population by using the crossover and the mutation functions. The former allows creating a new solution from two previous ones and the latter is used to randomly modify the value of one parameter of the new solution in order to search on new areas of the domain. Finally, the objective functions are evaluated for the new generation and the non-dominated solutions are obtained. The two last steps are repeated until a stop criterion is reached. As result, the Pareto’s optimal front is build where each point represents a possible solution. 3. FE Model Updating considering the hybrid UKF-MGA algorithm The Kalman filter [2] is an estimator that consists in estimating the unknown variables by considering the measurements which contain statistical white noise. The mathematical representation of a general parameter identification problem solved by using the Kalman filter is based on the state-space model: 𝜽𝑘=𝜽𝑘−1+𝒘𝑘−1 (6) 𝒛𝑘=𝒉(𝜽𝑘,𝒙𝑘)+𝒗𝑘 (7) where 𝜽 is the vector containing the model parameters and 𝒛 extracts the modelling outputs of the modelling function 𝒉( ). The estimation process noise and the modelling uncertainty are taken into account in the vectors 𝒘 and 𝒗, respectively. Both terms are assumed to be white Gaussian noise with zero-mean and covariance matrices 𝑸 and 𝑹, respectively. The matrix 𝑹 is the sum of measurements noise and modelling noise: 𝑹=𝑹𝑚𝑒𝑎𝑠+𝑹𝑚𝑜𝑑𝑒𝑙 [4]. Nevertheless, the modelling uncertainty is neglected in this study as the same model is considered for each iteration of the simulation. An extension of the Kalman filter to nonlinear systems, the so-called extended Kalman filter (EKF) was developed [5]. The EKF requires linearization (Taylor series expansion) of Equation (7). In case 𝒉 is a numerical model, there is no analytical solution of the Jacobian of the function. Instead, a numerical method such as the finite difference algorithm must be used to approximate the Jacobian. To overcome this issue, the unscented Kalman filter (UKF) was proposed by Julier et al. [6]. The UKF is an estimator for nonlinear systems where no linearization of the modelling function 𝒉 is needed. This algorithm uses the unscented transform technique to select a group of sample points (sigma points) which are propagated through the nonlinear functions to calculate the mean and the covariance of the estimated parameters [7]. It implies that the second order Gaussian approximation of the estimates is preserved. In this study, the genetic algorithm hybridized with UKF is considered due to the advantages of this algorithm over the EKF. A detailed description of the algorithm for parameter estimation problems is given below. As a Kalman filter, the UKF consists of two steps: (i), the prediction and; (ii) the correction (or update). The former consists in assessing the sigma points and the prior estimation error, modelling error and the cross covariance. The latter involves computing Kalman's Gain matrix (𝑲) and correcting the prior estimates of the parameters and the covariance, based on the measurements 𝒛𝑜𝑏𝑠. IABSE Symposium 2019 Guimarães: Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 4 The sigma points are 2𝑛+1 distributed points in the vicinity of the current estimate. They are defined from the prior estimation error covariance [8]: 𝓢=𝜂√𝑷𝑘−1 (8) (𝜒𝑘−1)0=𝜽 𝑘−1 (9) (𝜒𝑘−1)𝑖=𝜽 𝑘−1+𝓢𝑖, 𝑖=1,2,…,𝑛 (10) (𝜒𝑘−1)𝑖+𝑛=𝜽 𝑘−1−𝓢𝑖, 𝑖=1,2,…,𝑛 (11) being 𝑛 the number of parameters, 𝜂 is the scaling parameter computed as 𝜂=√𝑛+𝜆, 𝑷 is the estimation error covariance and 𝝌 are the 2𝑛+1 sigma points. The sigma points are multiplied by a weight, defined as: 𝑊0=𝜆 𝑛+𝜆 (12) 𝑊𝑖=𝑊𝑖+𝑛=1 2(𝑛+𝜆), 𝑖=1,2,…,𝑛 (13) This algorithm is integrated into the MGA. Specifically, the UKF searches around the candidate values of the parameters proposed by the genetic algorithm. Once the value of the parameters 𝜽 𝑘|𝑘 are obtained, the objective functions are assessed by minimising the mean square error between numerical and experimental modal parameters. The genetic algorithm uses these values of the objective functions to propose a new possible solution (vector of parameters). The algorithm is shown below. Initial Step: 𝜽 0=𝜽𝑝𝑟𝑖𝑜𝑟 𝑷0 𝑚=𝑷𝑝𝑟𝑖𝑜𝑟 Main loop: 𝑓𝑜𝑟 𝑘=1:𝑁𝑈𝐾𝐹 (number of iterations of the UKF) Prediction step: (𝜒𝑘−1)𝑖 (𝜒𝑘|𝑘−1)𝑖=(𝜒𝑘−1)𝑖 𝜽 𝑘|𝑘−1=∑𝑊𝑖·(𝜒𝑘|𝑘−1)𝑖 2𝑛 0 𝑷𝑘|𝑘−1 𝜃=∑(𝑊𝑖[(𝜒𝑘|𝑘−1)𝑖−𝜽 𝑘|𝑘−1]·[(𝜒𝑘|𝑘−1)𝑖−𝜽 𝑘|𝑘−1]𝑇) 2𝑛 0+𝑸 (𝓏𝑘|𝑘−1)𝑖=𝒉((𝜒𝑘|𝑘−1)𝑖) 𝒛𝑘|𝑘−1=∑𝑊𝑖·(𝓏𝑘|𝑘−1)𝑖 2𝑛 0 𝑷𝑘|𝑘−1 𝑧=∑(𝑊𝑖[(𝓏𝑘|𝑘−1)𝑖−𝒛𝑘|𝑘−1]·[(𝓏𝑘|𝑘−1)𝑖−𝒛𝑘|𝑘−1]𝑇) 2𝑛 0+𝑹 𝑷𝑘|𝑘−1 𝜃𝑧 =∑(𝑊𝑖[(𝜒𝑘|𝑘−1)𝑖−𝜽 𝑘|𝑘−1]·[(𝓏𝑘|𝑘−1)𝑖−𝒛𝑘|𝑘−1]𝑇) 2𝑛 0 Correction (update) step: IABSE Symposium 2019 Guimarães: Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 5 𝑲𝑘=𝑷𝑘|𝑘−1 𝜃𝑧 (𝑷𝑘|𝑘−1 𝑧)−1 𝜽 𝒌=𝜽 𝑘|𝑘−1+𝑲𝑘(𝒛𝑜𝑏𝑠−𝒛𝑘|𝑘−1) 𝑷𝑘 𝑚=𝑷𝑘|𝑘−1 𝑚−𝑲𝑘𝑷𝑘|𝑘−1 𝑧𝑲𝑘 𝑇 𝑒𝑛𝑑 4. Model updating of a real historical construction As benchmark structure, the domes of the chapel integrated in the Würzburg Residence (Germany) is considered. The domes are masonry surfaces of bricks grown in a single lawyer. The bricks are 30 cm thick, and at the base of the dome, they are reinforced by increasing the thickness until 45 cm. The reinforcement is also presented in the radial ribs reaching a section of 45x45 cm. In the construction process of the domes, the nerves are built alongside the rest of the sheet. a) b) Figure 1. a) Interior of the chapel and b) exterior of the dome of the chapel of the Würzburg Residence 4.1 Preliminary Numerical Model The preliminary FE model was built using the Ansys software [9]. The structural model of the dome was defined by considering 3D shell elements (SHELL63) which have 4 nodes per element. The effect of the lateral walls of the chapel on the dome was established, in a simplified manner, by means of two 1D spring elements (COMBIN14) whose stiffness represents the bending rigidity of the wall in the direction under consideration. The following mechanical properties of the material have been considered: Young’s modulus, E (GPa) = 2; Poisson’s ratio, 𝜈 = 0.2; density, 𝜌 (kg/m3) = 1700 and horizontal stiffness of the walls, 𝑘ℎ = 1000 kN/m. The latter was determined by considering that these elements were made with the same material than the rest and assimilating their behaviour to a cantilever. To adequately model the behaviour of the walls, four different zones were assumed (a wall for each lateral dome and two walls for the main dome) (see Figure 2). The numerical modal analysis led to the three natural frequencies shown in Table 1. Figure 2. FE model of the dome of the chapel 4.2 Ambient Vibration Test and Operational Modal Analysis An ambient vibration test was performed in order to experimentally identify the modal parameters of the structure. To this effect, the accelerations IABSE Symposium 2019 Guimarães: Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 6 were recorded during 10 minutes by eight uniaxial accelerometers at 51 instrumented points. The signals were processed using the time domain SSI algorithm [10], implemented in the software ARTeMIS Modal [11]. This whole procedure has been detailed in [12]. The three identified natural frequencies are shown in Table 1. Besides, the relative difference between the numerical and experimental natural frequencies and the MAC are calculated for each vibration mode. Table 1. Numerical and experimental natural frequencies Vibration Mode 𝒇𝒏𝒖𝒎 [Hz] 𝒇𝒆𝒙𝒑 [Hz] 𝚫𝒇 [%] MAC 1 (lateral) 3.362 2.764 21.635 0.809 2 (longitudinal) 3.431 3.231 6.190 0.877 3 (lateral) 3.988 4.015 -0.672 0.648 From the results in Table 1, it can be noticed that despite the use of a very detailed FE model, the natural and experimental mode shapes are not in good agreement as the errors between them are high. A FE model updating process must be conducted to guarantee that the numerical model represents more accurately the actual behaviour of the dome. 4.3 Finite Element Model Updating This process was implemented in the software programme MATLAB [13]. As global optimization algorithm, the MGA was used. This study focuses on the comparison of the two previous described techniques to find the global minimum of the objective functions: (i), the maximum likelihood method (MLM) based on MGA and; (ii), the hybrid UKF-MGA algorithm. The parameters selected for the model updating process are the same for both problems. A sensitivity study was performed to select the physical parameters that influence more the dynamic behaviour. This influence was taken into account on the basis of the modal strain energy. The seven parameters chosen were the followings: the Young’s modulus of the dome, the Young’s modulus of the ribs, the Young’s modulus of the starting wall, the longitudinal stiffness of the walls 1 and 4, the lateral stiffness of the walls 2 and 3, the lateral stiffness of the walls 1 and 4 and the longitudinal stiffness of the walls 2 and 3. The search range of the parameters were: (i) for the Young’s moduli, [1–4] GPa and; (ii) for the stiffness, [200-2200] kN/m. The objective functions were defined as expressed in Equations (4) and (5). 4.3.1 MLM-MGA method The value of the updated parameters is shown in Table 2. Table 2. Updated physical parameters (MLMMGA) Parameter Updated value 𝑬 dome [Pa] 1.164E+9 𝑬 ribs[Pa] 1.468E+9 𝑬 starting wall [Pa] 1.319E+9 𝒌𝒙 walls 1 and 4 [N/m] 7.531E+5 𝒌𝒚 walls 2 and 3 [N/m] 3.941E+5 𝒌𝒚 walls 1 and 4 [N/m] 1.729E+6 𝒌𝒙 walls 2 and 3 [N/m] 9.793E+5 The modal analysis of the updated FE model led to the three natural frequencies given in Table 3. The improvement of the relative differences between natural frequencies (except the third vibration mode) and the increase of MAC ratios exhibit the good performance of the updating procedure. Table 3. Natural frequencies of the updated FE model (MLM-MGA) Vibration Mode 𝒇𝒏𝒖𝒎,𝒖𝒑 [Hz] 𝒇𝒆𝒙𝒑 [Hz] 𝚫𝒇 [%] MAC 1 2.725 2.764 -1.411 0.861 2 3.237 3.231 -0.186 0.921 3 4.188 4.015 4.299 0.728 IABSE Symposium 2019 Guimarães: Towards a Resilient Built Environment - Risk and Asset Management March 27-29, 2019, Guimarães, Portugal 7 4.3.2 Hybrid UKF-MGA algorithm In the implementation of the UKF-MGA the following parameters were chosen [14]: initial estimation error covariance, 𝑷0 𝑚= 𝑑𝑖𝑎𝑔((𝑢𝑏−𝑙𝑏 100 )2); covariance matrix of process noise, 𝑸=0.005𝑷0; the covariance matrix of measurements noise is taken from the OMA, 𝑹= 0.001·𝐈(𝑟,𝑟), with 𝐈(𝑟,𝑟) the identity matrix having dimensions equal to the size of the measurements data, 𝑟 and 𝜆=0.0001. The measurements vector, 𝒛𝑜𝑏𝑠, was defined taking into account both the experimental natural frequencies and the mode shapes. The latter was implemented in terms of the MAC, considering the unit for the three vibration modes. In this manner, the difference 𝒛𝑜𝑏𝑠−𝒛 represent the error between observed and numerical modal parameters. The updating process gave as a result the model parameters presented in Table 4. Table 4. Updated physical parameters (UKF-MGA) Parameter Updated value 𝑬 dome [Pa] 2.848E+9 𝑬 ribs[Pa] 3.424E+9 𝑬 starting wall [Pa] 2.778E+9 𝒌𝒙 walls 1 and 4 [N/m] 6.418E+5 𝒌𝒚 walls 2 and 3 [N/m] 2.371E+5 𝒌𝒚 walls 1 and 4 [N/m] 1.862E+6 𝒌𝒙 walls 2 and 3 [N/m] 9.728E+5 The natural frequencies obtained after this optimization are shown in Table 5. This algorithm also conducts to a good adjustment of the numerical and experimental modal parameters. Table 5. Natural frequencies of the updated FE model (UKF-MGA) Vibration Mode 𝒇𝒏𝒖𝒎,𝒖𝒑 [Hz] 𝒇𝒆𝒙𝒑 [Hz] 𝚫𝒇 [%] MAC 1 2.741 2.764 -0.832 0.859 2 3.207 3.231 -0.743 0.928 3 4.423 4.015 10.161 0.695 4.4 Discussion of the results A discussion of the performance of both algorithms, GA and GA-UKF, is presented. The comparison is developed in terms of two key criteria: regarding the relative differences between numerical and experimental natural frequencies and MAC ratios and with respect to the number of iterations and computation time. First, from Table 3 and Table 5, it can be observed that the relative differences, Δ𝑓, are smaller for the MLM-GA (except the first vibration mode). In addition, the MAC ratios are larger for this algorithm (except for the second vibration mode). Second, the number of iterations and the computational cost is compared (see Table 6). The algorithms were run in a 3.60 GHz processor with a 29.3 RAM. For the MLM-GA, the number of iterations of the multi-objective genetic algorithm was established in 100 and a population of 200 solutions (vector of possible parameters) was considered. The time to complete the procedure was around 47000 seconds. However, for the UKFMGA, the number of iterations of the multiobjective genetic algorithm was set to 20, the initial population was made up of 10 vectors and 5 iterations of the unscented Kalman filter were considered. The reduced number of iterations and population of the genetic algorithm is due to the multiple (2𝑛+1, being 𝑛 the number of parameters) evaluations of the objective functions for each iteration of the UKF. The duration of the UKF-MGA is around 16000 seconds. This is 34% less than the time of the GA. Based on the two criteria, the UKF-MGA represents a useful tool to perform a model updating process, as the computational cost is considerably reduced and the results, although worst, are similar to those obtained with the MLM-MGA. Table 6. Comparison of the MLM-MGA and UKFMGA algorithms 𝑵 (MGA) Pop. 𝑵 (UKF) Time [s] MLM-MGA 100 200 - 47710 UKF-MGA 20 10 5 16068 involved in Bayesian updating may significantly be reduced by parallelization, the ease of implementation of the maximum likelihood method for such practical applications has prompted us to tackle the updating problem using these later algorithms. In order to solve the resulting optimization problem, either local or global optimization algorithms can be employed (Nocedal and Wright, 1999). As the relationship between the considered physical parameters and the modal properties of the structure is clearly nonlinear, multiple optimums are normally expected in the objective function of this problem. Provided that local optimization algorithms may have difficulties to find the global optimum under such scenario, global optimization algorithms are usually employed to solve the updating problem. To this end, computational intelligence algorithms, like genetic algorithms (GA) and artificial neural networks (ANN) have been widely used for practical engineering applications (Marwala, 2010). Both algorithms have shown their efficiency when solving nonlinear optimization problems (Koh and Perry, 2009; Levin and Lieven, 1998a). For instance, GA have been successfully employed to assess the structural behaviour of existing bridge structures (Jiménez-Alonso and Sáez, 2016) or to detect damage in civil engineering structures (Hao and Xia, 2002). ANN have as well been implemented to conduct the finite-element-model updating of a wide variety of structural systems (Lu and Tu, 2004; Hasançebi and Dumlupınar, 2013; Guo et al., 2017). Many other success stories are scattered in the scientific literature for finite element updating. The main limitation of these computational algorithms is the high simulation time required to perform the finite-element-model updating when the complexity of the structure increases. Two trends have been considered to overcome this drawback (if we keep parallelization aside in this review): (i) the implementation and validation of more recent global optimization algorithms and (ii) the hybridization between local and global optimization algorithms. As representative examples of the first trend, Levin and Lieven (1998b) performed a comparison between the performance of the GA and the simulated annealing (SA) to solve the finite-element-model updating problem of a flat plate wing; Perera et al. (2010) compared the results obtained between GA and the particle swarm optimization (PSO) for the finite-element-model updating of a one-story reinforced concrete frame. On the other hand, two outstanding examples of the second trend can be remarked: (i) the work reported by Shabbir and Omenzetter (2015), who proposed a hybrid sequential niche-PSO algorithm to perform the finite-element-model updating of a footbridge, and (ii) the study provided by Astroza et al. (2016), who implemented a hybrid unscented Kalman filter-SA algorithm to perform the finite-element-model updating of a steel frame structure. In this paper, a novel computational algorithm devised to improve the efficiency of finite-element-model updating of civil engineering structures is proposed, implemented and further validated. The design of this new algorithm employs the maximum likelihood approach and stems from the two trends mentioned above. On the one hand, as basis of this new algorithm, the Harmony Search (HS) global optimization algorithm is adopted. The harmony search algorithm is based on the improvisation process in which musicians seek harmony (Geem et al., 2001). This algorithm has been applied efficiently in multiple practical engineering applications (Yang and Koziel, 2011). However, according to the best of the authors’ knowledge, it has been rarely applied for the finite-element-model updating of civil engineering structures (Fadel et al., 2012; Kaveh et al., 2014). On the other hand, the performance of the original HS algorithm is improved via its hybridization with a local optimization algorithm: the so-called Unscented Kalman Filter (UKF) is considered herein. This algorithm belongs to the Kalman filter family (Kalman, 1960) which originally focused on estimating the means and the covariance of the constitutive parameters of a linear system through a recursive process with two steps: (i) the prediction; and (ii) the correction steps. Subsequently, the original Kalman filter was adapted to nonlinear systems (Jazwinski, 1970) through the so-called Extended Kalman Filter (EKF). The basic idea of this algorithm is to linearize locally the nonlinear function, which represents the behaviour of the system, before applying the Kalman filter. Thus, one of the main drawbacks of this algorithm is that it only allows computing a first order approximation of the nonlinear function. In order to overcome this limitation, the Unscented Kalman Filter (UKF) was subsequently proposed by Julier and Uhlmann (1997). This algorithm defines a set of sampling points which allow estimating the mean and covariance of the constitutive parameters after their propagation through the nonlinear system. The direct estimation of these statistical properties allows computing a second order approximation (or even a third order approximation in case of a Gaussian random vector is propagated) of the nonlinear function (Van Der Merwe, 2004; Nguyen and Nestorović, 2015). The resulting hybrid UKF-HS algorithm is then formulated to tackle two types of optimization problems: (i) single-objective and (ii) multi-objective approaches. The implementation details are next thoroughly presented and discussed. The performance of the new algorithm is validated and further compared with the results of other two global computational intelligence algorithms (GA and HS). To this end, the finite element model updating of a benchmark reconfigurable steel footbridge (Hudson and Reynolds, 2017) located at the laboratory of the Vibration Engineering Section of the University of Exeter (U.K.) is considered. The experimental modal properties of the structure were obtained via the signal processing of the records obtained during a forced vibration test (Maia and Silva, 1997). The paper is organized as follows. First, some basics about finite-element-model updating under the maximum likelihood method are presented. Subsequently, the three considered computational intelligence algorithms (GA, HS, UKF-HS) are described in detail, with special emphasis on the newly proposed UKF-HS algorithm. Subsequently, the performance of the three mentioned algorithms is compared when they are implemented for the finite-element-model updating of a laboratory footbridge. Finally, some concluding remarks are drawn to close the paper in the fifth section. 2. Basics of Finite-element-model Updating under the Maximum Likelihood Method As indicated above, the finite-element-model updating problem may be considered as a parameter identification problem in which the values of the most relevant physical parameters of the structure are to be estimated through the minimization of the differences between the numerical and experimental modal properties of the structure (Mottershead et al., 2011). Different estimators are normally used to solve this parameter identification problem. Among them, the maximum likelihood method has been widely used to perform the finite-element-model updating of civil engineering structures (Marwala, 2010) and it will be the approach considered herein. In this manner, the finite-element-model updating problem is transformed into an optimization problem, where the updated model follows from iterative modifications of the most relevant physical parameters, until an optimal solution is obtained. The objective function of such optimization problem may be defined in terms of the square relative differences (residuals) between the numerical and the experimental modal properties of the structure (natural frequencies and vibration modes). As design parameters of this problem, the most relevant physical parameters of the model are considered. Provided that the relationship between the residuals and the design parameters is nonlinear, global optimization algorithms are needed to perform the optimization process (Nocedal and Wright, 1999). Additionally, a search domain is usually established to guarantee that the solution obtained retains its physical meaning. Two approaches can be considered to define the objective function: (i) the singleobjective approach; and (ii) the multi-objective approach. Under the single-objective approach, the optimization problem consists in minimizing a single-objective function defined in terms of the weighted residuals between the numerical and experimental modal properties of the structure. These weights take into account the relative contribution of each residual. The value of these weights can be established either by a statistical criterion (Friswell and Mottershead, 1995) or by a correlation criterion (Teughels et al., 2002). According to the first criterion, the weights are determined in terms of the uncertainty associated with the estimates of the experimental modal properties of the structure. According to the second criterion, the weights are defined by an iterative process (trial and error) which ensures the best agreement between the numerical and experimental modal properties of the structure. The second criterion will be considered herein. Under the multi-objective approach, the optimization problem consists in minimizing a multi-objective function defined by several functional components. The overall objective function is normally defined by two functional components for practical engineering applications (Jin et al., 2014). One component may be defined in terms of the residuals associated with the natural frequencies, 𝑟𝑓,𝑗(𝛉), whilst the other component may be defined in terms of the residuals associated with the vibration modes, 𝑟𝑚,𝑗(𝛉) (𝑗 being the considered vibration mode). In this case, no weights are needed to define the objective function. However, a set of possible solutions (the socalled Pareto front) is obtained as result of the optimization process. Hence, a subsequent decision making problem must be addressed in order to select the best solution among the different elements of the Pareto front. Among the different criteria provided in literature (Infantes et al., 2019), the normal boundary intersection (NBI) method has been considered herein (Deb and Gupta, 2011). Residuals 𝑟𝑓,𝑗(𝛉) and 𝑟𝑚,𝑗(𝛉) can be defined as follows: 𝑟𝑓,𝑗(𝛉)=𝑓num,𝑗(𝛉) − 𝑓𝑒𝑥𝑝,𝑗 𝑓exp,𝑗 𝑗=1,2,3,…𝑛𝑓 (1) 𝑟𝑚,𝑗(𝛉)2=(1−√𝑀𝐴𝐶𝑗(𝛉))2 𝑀𝐴𝐶𝑗(𝛉) 𝑗=1,2,3,…𝑛𝑓 (2) where 𝛉 is a vector containing the most relevant physical parameters of the model; 𝑛𝑓is the total number of vibration modes considered to conduct the updating; 𝑓𝑛𝑢𝑚,𝑗 and 𝑓𝑒𝑥𝑝,𝑗 [Hz] are the numerical and experimental natural frequencies (𝑗 being the number of the considered vibration mode), respectively; and 𝑀𝐴𝐶𝑗 denotes the Modal Assurance Criterion (Allemang and Brown, 1982) which assesses the correlation between the numerical and experimental vibration modes. The 𝑀𝐴𝐶𝑗 ratio can be expressed as follows: 𝑀𝐴𝐶𝑗(𝛉)=(𝜙num,𝑗(𝛉)·𝜙exp,𝑗)2 (𝜙num,𝑗 𝑇(𝛉)·𝜙num,𝑗(𝛉))·(𝜙exp,𝑗 𝑇·𝜙exp,𝑗) (3) with 𝜙𝑛𝑢𝑚,𝑗 and 𝜙𝑒𝑥𝑝,𝑗 being the numerical and experimental vibration modes, respectively. On the one hand, the formulation of the finite-element-model updating problem using the single-objective approach can be formulated as follows: min𝑓(𝛉)=12[∑𝑤𝑓,𝑗·𝑟𝑓,𝑗(𝛉)2 𝑛𝑓 𝑗]1/2+12[∑𝑤𝑚,𝑗·𝑟𝑚,𝑗(𝛉)2 𝑛𝑓 𝑗]1/2 subject to{ 𝛉𝐥≤𝛉≤𝛉𝐮;∑𝑤𝑓,𝑗+𝑤𝑚,𝑗=1 𝑛𝑓 𝑗} (4) 𝑤𝑓 and 𝑤𝑚 being the weights associated with the residuals of the natural frequencies and vibrations modes; and 𝛉𝐥 and 𝛉𝐮 being the lower and upper bounds of the search domain of the considered physical parameters, respectively. On the other hand, the formulation of the finite-element-model updating problem using the multi-objective approach can be formulated as follows: min𝑓(𝛉)=min(𝑓1(𝛉) 𝑓2(𝛉)) where {𝑓1(𝛉)=12[∑𝑟𝑓,𝑗(𝛉)2 𝑛𝑓 𝑗]1/2 𝑓2(𝛉)=12[∑𝑟𝑚,𝑗(𝛉)2 𝑛𝑓 𝑗]1/2 subject to: 𝛉𝐥≤𝛉≤𝛉𝐮 (5) In summary, a typical finite-element-model updating consists of the following steps (Marwala, 2010): (i) the experimental modal properties of the structure are identified via the signal processing (either experimental or operational modal analysis may be employed to this end) of the records obtained during a vibration test (either forced vibration test or ambient vibration test); (ii) the numerical modal properties of the structure are obtained via a modal analysis based on a finite-element model; (iii) the most relevant physical parameters of the model are modified iteratively (using a global optimization algorithm); (iv) the objective function is evaluated; and (v) this iterative process is repeated until some convergence criterion is met. As result of this process, either a vector of updated parameters (single-objective approach) or a Pareto front (multi-objective approach) is obtained. In the latter case, a subsequent decision making problem must be addressed in order to select the vector of updated parameters. 3. Computational Intelligence Algorithms for the Finite-element model Updating of Civil Engineering Structures. The finite-element-model updating problem using the maximum likelihood method is usually tackled in practical engineering applications by employing computational intelligence algorithms (Marwala, 2010). These algorithms present as main advantage that the solution obtained is a global optimum of the objective function. However, they have the drawback of the high simulation time required to perform the updating process, especially when the complexity of the model increases. In order to overcome this limitation, several strategies have been adopted, as summarized in the introduction: from parallelizing the algorithms to the development of alternative global optimizers or the devising of hybrid local-global algorithms. In this paper, a novel hybrid Unscented Kalman Filter-Harmony Search (UKF-HS) algorithm is proposed and validated. This hybrid algorithm takes advantage of both the acceleration scheme provided by the unscented Kalman filter algorithm and the global search characteristics of the harmony search algorithm. The performance of this hybrid algorithm is assessed by comparing the results of the finite-element-model updating of a laboratory footbridge conducted using three different computational intelligence algorithms: (i) the classic genetic algorithms (GA); (ii) the harmony search algorithm (HS); and (iii) the novel hybrid UKF-HS algorithm. Further comparison of the performance of the new algorithm with the quite numerous optimizers available in the literature is out of the scope of this paper. However, it is relevant to point out here that previous studies (Marwala, 2010) reveal that Particle Swarm Optimizers (PSO) outperform other algorithms like Nelder-Mead (NM) or Simulated Annealing (SA). More recent studies (Jiménez-Alonso et al., 2017) show that the HS optimizer outperforms both GA and PSO when addressing the finite element model updating of civil structures. For this reason, our proposal selects HS as global optimizer and hybridizes it with the local UKF algorithm, in order to speed up the updating process. In order to make this paper as self-contained as possible, the three mentioned computational algorithms considering both single-objective and multi-objective approaches are next described. First, a general overview of GA is included. Next, the HS algorithm is summarized. Subsequently, the theoretical formulation of the UKF algorithm is presented. Finally, the newly proposed hybrid UKF-HS algorithm is described in detail. 3.1. Genetic Algorithms (GA) Genetic algorithms (GA) are nature-inspired computational algorithms based on the natural selection theory. Thus, the assessment of each chromosome (physical parameter vector) is based on the value of the objective function for each candidate. The GA, under the single-objective approach, can be summarized in the following steps (Koh and Perry, 2009): (i) an initial population of parameter vectors is randomly created; (ii) the objective function is evaluated for all the candidates; (iii) an iterative process is developed where parameter vectors are selected in order to create a new population; (iv) a new population is created using both mechanisms, the crossover, which allows obtaining a new vector from to previous ones, and the mutation, which consists in modifying the value of one component of the parameter vector to explore new areas of the search domain; (v) the objective function is evaluated for the new population at each iteration and the candidates with higher fitness are selected; (iv) the steps (iii) to (v) are repeated iteratively until a convergence criterion is met. The flowchart of the GA is shown in Figure 1a. On the other hand, the genetic algorithms, under the multi-objective approach (MGA), are addressed by the non-dominated sorting genetic algorithms (the so-called NSGA-II algorithm) proposed by Deb et al. (2002). According to this algorithm, the non-dominated solutions are classified using the non-dominated sorting method. In this manner, a new population is generated considering only the non-dominated solutions. Both mechanisms are controlled by an operator denominated crowded comparison operator. This operator is based on two attributes of each solution: (i) the rank; and (ii) the crowding distance. As result, a set of solutions, which constitutes the Pareto front, is obtained. The flowchart of the MGA algorithm is shown in Figure 1b. Figure 1. Flowchart of the genetic algorithms: a) GA and b) MGA. 3.2. Harmony Search Algorithm (HS) The HS algorithm was proposed by Geem et al. (2001). The algorithm is inspired in the musical improvisation where harmony is searched according to aesthetic requirements. As the GA, the HS finds the global minimum of the objective function by modifying iteratively the value of the considered physical parameters of the model. Although this algorithm has been implemented successfully for several practical engineering applications (Manjarres et al., 2013); however it has not been yet implemented for the finite-element-model updating of civil engineering structures according to the best of the authors’ knowledge. The HS algorithm, under the single-objective approach, can be summarized in the following steps: (i) the harmony matrix, 𝐇, is created (which contains the initial candidate solutions, parameter vector, created randomly); (ii) the objective function is evaluated for each solution; (iii) a new harmony is generated using three different mechanisms (memory consideration, pitch adjustments and randomization); (iv) the Initial Population Evaluate Selection Crossover Mutation Convergence? Solution Yes No Operations Update Population Initial Population Evaluate Update Population Evaluate Is Dominated? Update Population Yes No Convergence? Yes No Pareto Front a) b) objective function is evaluated for each new harmony, (v) the harmony matrix, 𝐇, is updated with the best generated harmonies; and (vi) the steps (iii) to (v) are repeated until a convergence criterion is met. When a new harmony is generated, each element of a new candidate vector can be defined in terms of either a previous value stored in the harmony matrix, 𝐇, or adopting a new random value. This fact is controlled by the harmony memory consideration rate, 𝐻𝑀𝐶𝑅. This ratio establishes the probability of selecting a previous component of the harmony matrix, 𝐇. Additionally, when some elements adopts the value of a previous one, it can be mutated according to the pitch adjustment rate, 𝑃𝐴𝑅. This ratio establishes the probability of mutating an element of the candidate vector. The adjustment is based on an additional parameter, the so-called bandwidth, 𝑏𝑤, which is added or subtracted to mutate the considered candidate vector. The flowchart of the HS algorithm is shown in Figure 2a. The HS algorithm, under the multi-objective approach (MHS), is an extension of the above mentioned algorithm which allows minimizing multi-objective functions. At each iteration, a set of new possible candidate solutions is created according to both the 𝐻𝑀𝐶𝑅 and 𝑃𝐴𝑅 ratios. The classification of the non-dominated solutions is performed using both the non-dominated sorting method (Deb et al., 2002) and the crowding distance. In order to restore the initial size of the harmony matrix, 𝐇, the worst solutions in terms of crowding distance are deleted. These steps are repeated iteratively until a convergence criterion is met. As result of this process, a set of possible solutions is obtained, the so-called Pareto front. The flowchart of the MHS algorithm is shown in Figure 2b. Figure 2. Flowchart of the harmony search algorithm: a) HS and b) MHS. Initialize H Convergence? Solution Yes No Operations Evaluate Selection Improvisation Randomization Update H Initialize H Evaluate Update H Evaluate Is Dominated? Update H Yes No Convergence? Yes No Pareto Front b)a) 3.3. Unscented Kalman Filter (UKF) A local minimization algorithm has been combined with the HS algorithm to reduce the simulation time needed to perform the updating process. The unscented Kalman filter (UKF) has been selected for this purpose (Julier and Uhlmann, 1997). This algorithm is a derivative free estimator, which involves that neither Jacobians nor Hessians must be computed, widely used for the state-estimation of nonlinear dynamic systems. Among the different transformations, the scaled unscented transformation proposed by Wan and Van Der Merwe (2000) has been adopted herein to define the UKF algorithm. The formulation for a parameter identification problem may be represented as follows (Wan et al., 2000): 𝛉𝑘=𝛉𝑘−1+𝐰𝑘−1 (6) 𝐳𝑘=ℎ(𝛉𝑘)+𝐯𝑘 (7) where 𝛉 is the parameter vector; ℎ is the nonlinear modelling function, 𝐳 extracts the outputs of the dynamic system; 𝐰 is the statistical noise of the identification process; and 𝐯 is the statistical noise of observation process. Both type of noise are assumed to be uncorrelated and white Gaussian noise with zero-mean and covariance matrices 𝐐 and 𝐑 respectively. It is necessary to remark that the matrix 𝐑 may be computed by means of two terms (Tarantola, 2005): (i) the measurement noise; and (ii) the modelling noise. As the same models have been employed during the simulation and estimation phases, the effects of modelling noise has not been considered herein (Astroza et al., 2016). The UKF is a straightforward extension of the unscented transformation (Julier and Uhlmann, 2004). The estimation is addressed by considering 2𝑛𝑑+1 (being 𝑛𝑑 the number of parameters) deterministic sampling points (sigma points) which can be used to evaluate the true mean and the covariance. The propagation of these sigma points through the nonlinear function ℎ leads to compute the true posterior mean and covariance of the estimated parameters up to the second order of the Taylor series expansion of the nonlinear function (third order of the Taylor series expansion for a Gaussian inputs). The determination of the sigma points is based on the square-root decomposition of the posterior covariance matrix, 𝐏. For this reason, the matrix 𝐏 must be positive semidefinite at each step. Even though the square-root decomposition can be efficiently derived using the Cholesky factorization (𝐀=√𝐏=chol(𝐏) being 𝐏= 𝐀𝐀T), the matrix 𝐏 is still updated at each iteration and numerical errors can give a nonpositive semidefinite matrix 𝐏. The square-root UKF algorithm, proposed by Van Der Merwe and Wan (2001), overcomes this issue since it avoids factorising at each step. This algorithm propagates directly the matrix 𝐀 and guarantees that the covariance matrix is positive semi-definite. The square-root UKF algorithm has two steps: (i) the prediction step; and (ii) the correction (update) step. The prediction step is carried out considering the prior model to evaluate the sigma points and predict the estimates of the estimation error covariance, 𝐀𝜃, and calculate the model outputs, 𝔃𝑘|𝑘−1. The 2𝑛𝑑+1 sigma points are computed as: (𝝌𝑘−1)0=𝛉𝑘−1|𝑘−1 (8) (𝝌𝑘−1)𝑖= 𝛉𝑘−1|𝑘−1+√(𝑛𝑑+𝜆)(𝐀𝑘−1|𝑘−1 𝜃)𝑖 𝑖=1,2,…,𝑛𝑑 (9) (𝝌𝑘−1)𝑖+𝑛= 𝛉𝑘−1|𝑘−1−√(𝑛𝑑+𝜆)(𝐀𝑘−1|𝑘−1 𝜃)𝑖 𝑖=1,2,…,𝑛𝑑 (10) where 𝛉𝑘−1|𝑘−1 are the posterior parameters estimated at the previous step and 𝜆 is a scaling parameter. The sigma points are weighted under the assumption that the sum of all the weights must be equal to the unit. The weights are defined as: 𝑊0=𝜆 𝑛𝑑+𝜆 (11) 𝑊𝑖=𝑊𝑖+𝑛=1 2(𝑛𝑑+𝜆) 𝑖=1,2,…,𝑛𝑑 (12) For parameter identification problems, the estimation error covariance can be calculated by applying an exponential weighting on past data, 𝐀𝑘 𝜃=𝛾−1/2𝐀𝑘−1 𝜃, being 𝛾 a scalar factor slightly less than the unit (Van Der Merwe and Wan, 2001). The correction (update) step consists in estimating the posterior mean and estimation error covariance in terms of the Kalman’s gain matrix, 𝐊; the measurements, 𝐳𝑜𝑏𝑠; and the model outputs. The Kalman’s gain matrix is derived from the model outputs error covariance, 𝐒𝑧, and the cross covariance between the estimation error and the model outputs error covariances, 𝐏𝜃𝑧. In order to calculate the model outputs error covariance, 𝐒𝑧, the following expression (considering that 𝑊𝑖>0 for all 𝑖≥1) can be used (Terejanu, 2011): 𝐒𝑘𝑧=∑𝑊𝑖[(𝔃𝑘|𝑘−1 i−𝒛𝑘|𝑘−1)·(𝔃𝑘|𝑘−1 i−𝒛𝑘|𝑘−1)𝑇] 2𝑛𝑑 0+𝐑= [√𝑊𝑖(𝔃𝑘|𝑘−1 i−𝒛𝑘|𝑘−1), √𝐑]·[√𝑊𝑖(𝔃𝑘|𝑘−1 i−𝒛𝑘|𝑘−1)𝑇, √𝐑𝑇]𝑇+ 𝑊0[(𝔃𝑘|𝑘−1 0−𝒛𝑘|𝑘−1)·(𝔃𝑘|𝑘−1 0−𝒛𝑘|𝑘−1)𝑇] 𝑓𝑜𝑟 𝑖=1:2𝑛𝑑 (13) The first term can be expressed by means of a 𝑞𝑟 factorization 1 as the product of an orthogonal matrix 𝐆𝑘∈ ℝ2𝑛𝑑+𝑟×𝑟 and an upper triangular matrix 𝐒𝑘𝑧∈ℝ𝑟×𝑟, being 𝑟 the number of measurements. The last term, can be taken into account performing a rank 1 update to Cholesky factorization 2 . Therefore, the matrix 𝑺𝑘𝑧 can be calculated as (Terejanu, 2011): 𝐒𝑘𝑧=𝑞𝑟([√𝑊1:2𝑛𝑑·[(𝔃𝑘|𝑘−1)1:2𝑛𝑑−𝒛𝑘|𝑘−1], √𝐑 ]) (14) 𝐒𝑘𝑧=𝑐ℎ𝑜𝑙𝑢𝑝𝑑𝑎𝑡𝑒(𝐒𝑘𝑧, (𝔃𝑘|𝑘−1)0−𝒛𝑘|𝑘−1, sgn(𝑊0)𝑊0 ) (15) Once the posterior mean of the parameter estimates, 𝛉𝒌|𝒌 and the posterior estimation error covariance, 𝐀𝑘|𝑘 𝜃, are obtained, the algorithm follows a iterative process until the maximum number of iterations is reached. The considered algorithm may be formulated as follows. 1 The 𝑞𝑟 factorization of a matrix 𝐌 allows expressing the matrix as 𝐌=𝐐𝐑. Here, 𝐌 is an 𝑚by-𝑛 matrix, 𝐑 is an 𝑚-by-𝑛 upper triangular matrix and 𝐐 is an 𝑚-by-𝑚 unitary matrix. 2 If 𝐑=𝑐ℎ𝑜𝑙(𝐀), the Cholesky factor of the rank 1 update 𝐀+𝐯𝐯T is written as 𝐒= 𝑐ℎ𝑜𝑙𝑢𝑝𝑑𝑎𝑡𝑒(𝐑,𝐯). Initial Step: 𝛉0=𝛉𝑝𝑟𝑖𝑜𝑟 (16) 𝐀0 𝜃=𝑐ℎ𝑜𝑙(𝐏0θ) (17) Main loop: 𝑓𝑜𝑟 𝑘=1:𝑁𝑈𝐾𝐹 (number of iterations of the UKF) Prediction step: Calculate the 2𝑛𝑑+1 Sigma Points: (𝝌𝑘−1)𝑖 (18) (𝝌𝑘|𝑘−1)𝑖=(𝝌𝑘−1)𝑖 (19) 𝛉𝑘|𝑘−1=∑𝑊𝑖·(𝝌𝑘|𝑘−1)𝑖 2𝑛𝑑 0 (20) 𝐀𝑘|𝑘−1 𝜃=𝛾−12 𝐀𝑘−1 𝜃 (21) (𝔃𝑘|𝑘−1)𝑖=𝒉((𝝌𝑘|𝑘−1)𝑖) (22) 𝒛𝑘|𝑘−1=∑𝑊𝑖·(𝔃𝑘|𝑘−1)𝑖 2𝑛𝑑 0 (23) Correction (update) step: 𝐒𝑘|𝑘−1 𝑧=𝑞𝑟([√𝑊1:2𝑛𝑑·[(𝔃𝑘|𝑘−1)1:2𝑛𝑑−𝒛𝑘|𝑘−1] √𝐑 ]) (24) 𝐒𝑘|𝑘−1 𝑧=𝑐ℎ𝑜𝑙𝑢𝑝𝑑𝑎𝑡𝑒(𝐒𝑘|𝑘−1 𝑧,(𝔃𝑘|𝑘−1)0−𝒛𝑘|𝑘−1,𝑠𝑔𝑛(𝑊0)) (25) 𝐏𝑘|𝑘−1 𝜃𝑧 =∑(𝑊𝑖[(𝝌𝑘|𝑘−1)𝑖−𝛉𝑘|𝑘−1]·[(𝔃𝑘|𝑘−1)𝑖−𝒛𝑘|𝑘−1]𝑇) 2𝑛𝑑 0 (26) 𝐊𝑘=(𝐏𝑘|𝑘−1 𝜃𝑧 /𝐒𝑘|𝑘−1 𝑧𝑇) /𝐒𝑘|𝑘−1 𝑧 (27) 𝛉𝒌|𝒌=𝛉𝑘|𝑘−1+𝐊𝑘(𝒛𝑜𝑏𝑠−𝒛𝑘|𝑘−1) (28) 𝐔=𝐊𝑘𝐒𝑘|𝑘−1 𝑧 (29) 𝐀𝑘|𝑘 𝜃=𝑐ℎ𝑜𝑙𝑢𝑝𝑑𝑎𝑡𝑒(𝐀𝑘|𝑘−1 𝜃,𝐔,−1) (30) 𝑒𝑛𝑑 3.4. Hybrid UKF-HS algorithm The local minimization procedure followed by the square root UKF is implemented in both the HS and the MHS algorithms to reduce the simulation time required to solve the finite-element-model updating problem. The proposed hybrid algorithm is a local-global optimization algorithm which combines the main virtues of the two component algorithms. In this sense, it should be emphasized that this proposed algorithm is derivative free; not being required the computation of the Jacobian. Besides it guarantees the positive semi-definiteness of the covariance matrix as the intrinsic properties of the square root UKF algorithm are preserved. The general layout of this new hybrid algorithm follows the general scheme of the above mentioned HS algorithm. Thus, the different steps that configure both the single-objective approach (UKF-HS) and the multi-objective approach (UKF-MHS algorithm) have been included herein. First, the new hybrid UKF-HS algorithm consists of the following steps: (i) the HS algorithm creates the initial candidate solution (parameter vector); (ii) the square-root UKF algorithm computes the mean of the