Robust design optimization using Kriging models: Application to the robust design optimization of truss structures

==Resumen==

El problema de diseño óptimo robusto de estructuras es una tarea computacionalmente costosa como consecuencia del acoplamiento de los procesos de cuantificación de incertidumbre y de optimización. Para hacer frente a este problema, en este artículo se propone una metodología, basada en modelos Kriging, para resolver de forma eficiente el problema de cuantificación de incertidumbre en el proceso de optimización. El modelo Kriging aproxima, de forma simultánea, la respuesta estructural en el dominio de diseño y en el dominio estocástico, permitiendo desacoplar los procesos de cuantificación de incertidumbre y de optimización. La metodología propuesta incluye un criterio de actualización de los modelos Kriging basado en la estimación del error en la predicción, que mejora la aproximación en las regiones cercanas al frente de Pareto. Se han resuelto 3 problemas para mostrar la aplicabilidad y la precisión de la metodología propuesta. Los resultados muestran que la metodología es adecuada para resolver el problema de diseño óptimo robusto con una precisión razonable y un número de evaluaciones del modelo de simulación muy inferior al que requieren los métodos convencionales.

==Abstract==

Conventional methods addressing the robust design optimization problem of structures usually require high computational requirements due to the nesting of uncertainty quantification within the optimization process. In order to address such a problem, this work proposes a methodology, based on Kriging models, to efficiently assess the uncertainty quantification in the optimization process. The Kriging model approximates the structural performance both in the design domain and in the stochastic domain, which allows to decouple the uncertainty quantification process and the optimization process. In addition, an infill criterion based on the variance of the Kriging prediction is included to update the Kriging model towards the global Pareto front. Three numerical examples show the applicability and the accuracy of the proposed methodology. The results show that the proposed method is appropriate to solve the robust design optimization problem with reasonable accuracy and a considerably lower number of function calls than required by conventional methods.

==Palabras clave==

Optimización estructural ; Diseño óptimo robusto ; Metamodelos ; Modelos Kriging

==Keywords==

Structural optimization ; Robust design optimization ; Metamodels ; Kriging models

==1. Introducción==

Tradicionalmente, las incertidumbres en las condiciones de carga, en las propiedades de los materiales, en la geometría o en las condiciones de contorno de las estructuras se han incluido en el proceso de diseño mediante hipótesis basadas en la experiencia o en criterios de ingeniería tales como la utilización de factores de seguridad. Utilizando dichas hipótesis se obtiene un modelo simplificado, basado en los valores nominales de las variables y de los parámetros de diseño. Sin embargo, las soluciones óptimas que se han alcanzado con este enfoque determinista tienen un comportamiento óptimo únicamente bajo condiciones cercanas a las fijadas en el proceso de optimización, pudiendo deteriorarse en gran medida cuando las condiciones se alejan de las de diseño.

La necesidad de incorporar las incertidumbres en el proceso de diseño ha estimulado el interés por la investigación de procedimientos capaces de proporcionar diseños más robustos y fiables. En la actualidad existen 2 formulaciones que consideran la respuesta probabilista de las estructuras en el proceso de diseño óptimo: el diseño óptimo basado en fiabilidad (''Reliability-Based Design Optimization'' , RBDO  [[#bib0005|[1]]] ) y el diseño óptimo robusto (''Robust Design Optimization'' , RDO  [[#bib0010|[2]]] , [[#bib0015|[3]]]  and [[#bib0020|[4]]] ). La primera de ellas consiste en un problema de optimización en que se incluyen los efectos de la incertidumbre por medio de probabilidades de fallo y de valores esperados. La segunda trata de determinar un diseño menos sensible a las incertidumbres de las variables y de los parámetros que intervienen en la respuesta estructural mediante la utilización de índices de robustez.      

La principal dificultad para la aplicación de los actuales métodos de diseño óptimo robusto es el elevado coste computacional, tanto para evaluar los momentos estadísticos de la respuesta estructural como para obtener las sensibilidades necesarias para el proceso de optimización. Para abordar este problema, diversos autores proponen la utilización de aproximaciones o metamodelos durante el proceso de diseño óptimo robusto [[#bib0025|[5]]] , [[#bib0030|[6]]] , [[#bib0035|[7]]] , [[#bib0040|[8]]] , [[#bib0045|[9]]]  and [[#bib0050|[10]]] . Dellino et al. [[#bib0045|[9]]]  combinaron el método Taguchi de diseño robusto con superficies de respuesta para acelerar el proceso de cuantificación de incertidumbre. Li y Kang [[#bib0025|[5]]]  utilizaron modelos Kringing para aproximar los momentos estadísticos de la respuesta estructural y un algoritmo de recocido simulado para la búsqueda de un óptimo global. Guoqi y Daugi [[#bib0050|[10]]]  presentaron una metodología basada en la utilización de ''support vector regression''  para la resolución del problema de diseño óptimo robusto de problemas de alta dimensionalidad. Martínez y Martí  [[#bib0035|[7]]]  propusieron un procedimiento adaptativo basado en modelos Kriging para la resolución del problema de diseño óptimo robusto multiobjetivo. Jin et al. [[#bib0055|[11]]]  realizaron un estudio de la aplicabilidad y la precisión de diferentes técnicas de metamodelos para la resolución del problema de optimización bajo incertidumbre, en el que concluyeron que la precisión de los resultados óptimos dependía en gran medida de la capacidad del metamodelo para capturar globalmente la respuesta estructural. Beyer y Sendhoff [[#bib0060|[12]]]  realizaron una revisión del estado de la cuestión del problema de diseño óptimo robusto, en la que pusieron de manifiesto la importancia de una mayor investigación en nuevas técnicas de metamodelos que garanticen la precisión global de la aproximación. La dificultad para obtener una aproximación global precisa ha estimulado la investigación en procedimientos adaptativos que garanticen la obtención de un diseño óptimo global. En este sentido, autores como Jurecka et al. [[#bib0040|[8]]]  adaptaron técnicas de optimización dirigidas por metamodelos al campo del diseño óptimo robusto.      

El objetivo de este artículo es proponer una metodología para la resolución del problema de diseño óptimo robusto, que utiliza modelos Kriging para resolver de forma eficiente el problema de la cuantificación de incertidumbre en el proceso de optimización. A diferencia de otros trabajos, en la metodología propuesta el modelo Kriging aproxima simultáneamente la respuesta estructural en el dominio de diseño y en el dominio estocástico. El modelo Kriging sustituye al simulador y actúa como soporte del proceso de cuantificación de incertidumbre y del proceso de optimización. Asimismo, esta metodología incluye un criterio de actualización de los modelos Kriging, basado en la estimación del error en la predicción, que mejora la aproximación en las regiones cercanas al frente de Pareto. La metodología propuesta permite: 1) evaluar la robustez del diseño utilizando un número reducido de evaluaciones en comparación con los procedimientos convencionales, tales como los métodos de simulación; 2) desacoplar los procesos de cuantificación de incertidumbre y de optimización, reduciendo el número de evaluaciones del simulador en comparación con el enfoque anidado; y 3) mejorar la precisión de la aproximación en las regiones cercanas a los óptimos de Pareto.

==2. Metamodelos==

Un metamodelo (''model of model''[[#bib0065|[13]]] ) se utiliza habitualmente para sustituir globalmente a un modelo de simulación que requiere un alto coste computacional para su evaluación. Esta aproximación puede utilizarse, entre otras aplicaciones, para facilitar el proceso de optimización, explorar el espacio de diseño o llevar a cabo análisis de fiabilidad. En la literatura especializada existe una gran variedad de metamodelos, tales como las superficies de respuesta [[#bib0070|[14]]] , los modelos Kriging  [[#bib0075|[15]]]  and [[#bib0080|[16]]] , las funciones de base radial [[#bib0085|[17]]]  o las redes neuronales  [[#bib0075|[15]]]  and [[#bib0090|[18]]] . De entre las diferentes técnicas de metamodelos, los modelos Kriging [[#bib0080|[16]]]  han alcanzado una gran popularidad en los últimos años gracias a su gran flexibilidad para aproximar respuestas con alto grado de no linealidad [[#bib0055|[11]]]  y a que proporcionan información estadística del error cometido en la predicción [[#bib0095|[19]]] .      

===2.1. Modelos Kriging===

Los modelos Kriging asumen que el modelo de simulación puede ser aproximado por una realización de un proceso estocástico ''gaussiano  ''<math display="inline">G(\boldsymbol{\mbox{x}})</math>  con una media <math display="inline">E[G(\boldsymbol{\mbox{x}})]=</math><math>\boldsymbol{\mbox{f}}{\left(\boldsymbol{\mbox{x}}\right)}^T\quad \boldsymbol{\beta }</math>  y una covarianza <math display="inline">Cov[G(\boldsymbol{\mbox{x}}),G(\boldsymbol{\mbox{x}}^{{'}})]=</math><math>{\alpha }^2R(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{x}}^{{'}},\boldsymbol{\phi })</math> , a priori desconocidas. En las expresiones anteriores, '''''β'''''  = [''β''<sub>1</sub> , …, ''β''<sub>''p''</sub> ]<sup>''T''</sup>  es un vector de parámetros desconocidos, '''''f''''' ('''x''' ) = [''f''<sub>1</sub> ('''x''' ), …, ''f''<sub>''p''</sub> ('''x''' )]<sup>T</sup>  es un conjunto conocido de funciones de <math display="inline">\boldsymbol{\mbox{x}}\in {\mathbb{R}}^n</math>  (funciones de regresión), ''α''<sup>2</sup>  es la varianza de <math display="inline">G(\boldsymbol{\mbox{x}})</math>  y <math display="inline">R(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{x}}^{{'}},\boldsymbol{\phi })</math>  es la función de correlación entre '''x'''  y '''x''' ′. El proceso estocástico <math display="inline">G(\boldsymbol{\mbox{x}})</math>  representa el conocimiento a priori de la función que aproximar, y por esta razón la selección de la función de correlación debe ser consistente con la información conocida de dicha función. El modelo de correlación más ampliamente utilizado es el modelo exponencial anisótropo generalizado:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{x}}^{{'}},\boldsymbol{\phi })=</math><math>exp\left(\sum_{i=1}^n-{\phi }_i\vert x_i-x_i^{{'}}\vert ^s\right),\quad 1\leq s\leq 2.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

Los parámetros '''''β''''' , ''α''<sup>2</sup>  y '''''ϕ'''''  son desconocidos a priori y se han de determinar a partir de las respuestas del simulador <math display="inline">Y=</math><math>\lbrace y_1,\ldots ,y_m\rbrace </math>  para un conjunto de diseños <math display="inline">X=</math><math>\lbrace \boldsymbol{\mbox{x}}_1,\ldots ,\boldsymbol{\mbox{x}}_m\rbrace </math> . Utilizando técnicas bayesianas, la distribución posterior de <math display="inline">G(\boldsymbol{\mbox{x}})</math> , condicionada a las observaciones <math display="inline">Y=</math><math>\lbrace y_1,\ldots ,y_m\rbrace </math> , es ''gaussiana''[[#bib0100|[20]]]  con media:

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\mbox{ˆ}}{y}(\boldsymbol{\mbox{x}})\equiv E[G(\boldsymbol{\mbox{x}})\vert Y]=</math><math>\boldsymbol{\mbox{f}}{\left(\boldsymbol{\mbox{x}}\right)}^T\overset{\mbox{ˆ}}{\boldsymbol{\beta }}+</math><math>\boldsymbol{\mbox{r}}{\left(x\right)}^T\boldsymbol{\mbox{R}}^{-1}(Y^T-</math><math>\boldsymbol{\mbox{F}}\overset{\mbox{ˆ}}{\boldsymbol{\beta }}),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

y covarianza:

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Cov[G(\boldsymbol{\mbox{x}}),G(\boldsymbol{\mbox{x}}^{{'}})\vert Y]=</math><math>{\alpha }^2\lbrace R(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{x}}^{{'}},\boldsymbol{\phi })-</math><math>\boldsymbol{\mbox{r}}{\left(\boldsymbol{\mbox{x}}\right)}^T\boldsymbol{\mbox{R}}^{-1}\boldsymbol{\mbox{r}}(\boldsymbol{\mbox{x}}^{{'}})+</math><math>\boldsymbol{\mbox{u}}{\left(\boldsymbol{\mbox{x}}\right)}^T{\left(\boldsymbol{\mbox{F}}^T\boldsymbol{\mbox{R}}^{-1}\boldsymbol{\mbox{F}}\right)}^{-1}\boldsymbol{\mbox{u}}(\boldsymbol{\mbox{x}}^{{'}})\rbrace ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{u}}(x)=\boldsymbol{\mbox{f}}(\boldsymbol{\mbox{x}})-</math><math>\boldsymbol{\mbox{F}}^T\boldsymbol{\mbox{R}}^{-1}\boldsymbol{\mbox{r}}(\boldsymbol{\mbox{x}}),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\mbox{ˆ}}{\boldsymbol{\beta }}={\left(\boldsymbol{\mbox{F}}^TR^{-1}\boldsymbol{\mbox{F}}\right)}^{-1}\boldsymbol{\mbox{F}}^TR^{-1}Y,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\alpha }^2=\frac{1}{m}{\left(Y-\boldsymbol{\mbox{F}}\overset{\mbox{ˆ}}{\boldsymbol{\beta }}\right)}^T\boldsymbol{\mbox{R}}^{-1}(Y-</math><math>\boldsymbol{\mbox{F}}\overset{\mbox{ˆ}}{\boldsymbol{\beta }}),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4c)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_{ij}=f_j(\boldsymbol{\mbox{x}}_i),\quad i=1,\ldots ,m,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4d)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>j=1,\ldots ,p,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4e)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_{ij}=R(\boldsymbol{\mbox{x}}_i,\boldsymbol{\mbox{x}}_j,\boldsymbol{\phi }),\quad i,j=</math><math>1,\ldots ,m,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4f)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>r_i(\boldsymbol{\mbox{x}})=R(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{x}}_i,\boldsymbol{\phi }),\quad i=</math><math>1,\ldots ,m,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4g)
|}

donde <math display="inline">\overset{\mbox{ˆ}}{y}(\boldsymbol{\mbox{x}})</math>  representa la mejor predicción del modelo de simulación en el punto '''x''' , '''R'''  es la matriz de covarianza y <math display="inline">Cov[G(\boldsymbol{\mbox{x}}),G(\boldsymbol{\mbox{x}})\vert Y]</math>  es la varianza de la predicción, la cual proporciona intervalos de confianza sobre la predicción media <math display="inline">\overset{\mbox{ˆ}}{y}(\boldsymbol{\mbox{x}})</math> .      

En [[#eq0010|(2)]]  y [[#eq0015|(3)]]<math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{\beta }}</math>  y ''α''<sup>2</sup>  se obtienen analíticamente y solo dependen de los parámetros de correlación '''''ϕ''''' . En este trabajo los parámetros de correlación se han obtenido resolviendo el problema de optimización [[#bib0100|[20]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>max\left(-\frac{m\quad ln({\alpha }^2)+ln(\vert \boldsymbol{\mbox{R}}\vert )}{2}\right).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

Las funciones de correlación y de regresión más adecuadas dependen de las características de la función original, a priori desconocida. En este trabajo se han considerado funciones de regresión del tipo constante, lineal, cuadrática y cúbica, en combinación con las funciones de correlación recogidas en la [[#tbl0005|tabla 1]] . La mejor combinación de funciones de correlación y de regresión se selecciona en base al error cuadrático medio (''root mean square error'' , RMSE) y al coeficiente de determinación (R<sup>2</sup> ). Estos indicadores se obtienen utilizando las respuestas del simulador '''y''' ('''x'''<sub>v</sub> ) para una muestra de validación '''x'''<sub>v</sub> :

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Funciones de correlación

|-

! Nombre
! R('''x''' , '''x''' ′, '''''ϕ''''' )                                                    
|-

| Exponencial
| <math display="inline">exp\left({\sum }_{i=1}^n-\right. </math><math>\left. {\phi }_i\vert x_i-x_i^{{'}}\vert \right)</math>
|-

| Exponencial-Gaussiana
| <math display="inline">exp\left({\sum }_{i=1}^n-\right. </math><math>\left. {\phi }_i\vert x_i-x_i^{{'}}\vert ^s\right),1\leq s\leq 2</math>
|-

| Gaussiana
| <math display="inline">exp\left({\sum }_{i=1}^n-\right. </math><math>\left. {\phi }_i\vert x_i-x_i^{{'}}\vert ^2\right)</math>
|-

| Lineal
| <math display="inline">m\acute{a}x\lbrace 0,1-{\phi }_i\vert x_i-</math><math>x_i^{{'}}\vert \rbrace </math>
|-

| Esférica
| <math display="inline">1-1,5{\xi }_i+0,5{\xi }_i^3,\quad {\xi }_i=</math><math>m\acute{\imath}n\lbrace 1,{\phi }_i\vert x_i-x_i^{{'}}\vert \rbrace </math>
|-

| Cúbica
| <math display="inline">1-3{\xi }_i^2+2{\xi }_i^3,\quad {\xi }_i=</math><math>m\acute{\imath}n\lbrace 1,{\phi }_i\vert x_i-x_i^{{'}}\vert \rbrace </math>
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>RMSE=\sqrt{\frac{{\sum }_{i=1}^{m_v}{\left(y_i(\boldsymbol{\mbox{x}}_v)-{\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}}_v)\right)}^2}{m_v}},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R^2=1-\frac{{\sum }_{i=1}^{m_v}{\left(y_i(\boldsymbol{\mbox{x}}_v)-{\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}}_v)\right)}^2}{{\sum }_{i=1}^{m_v}{\left(y_i(\boldsymbol{\mbox{x}}_v)-\overline{y}(\boldsymbol{\mbox{x}}_v\right)}^2},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

donde <math display="inline">{\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}}_v)</math>  y ''y''<sub>''i''</sub> ('''x'''<sub>v</sub> ) son la componente ''i-ésima  ''  del vector de respuestas obtenidas con el modelo Kriging y el simulador, respectivamente, <math display="inline">\overline{y}(\boldsymbol{\mbox{x}}_v)</math>  es la media de todos los valores obtenidos con el simulador, y <math display="inline">m_v</math>  es el tamaño de la muestra de validación. Para no tener que utilizar muestras de validación, en este trabajo se estiman estos indicadores utilizando una medida basada en errores de validación cruzada denominada PRESS ''(prediction error sum of squares)''[[#bib0105|[21]]] :

<span id='eq0070'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{PRESS}_{RMSE}=\sqrt{\frac{{\sum }_{i=1}^m{\left({\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}})-{\overset{\mbox{ˆ}}{y}}_{-i}(\boldsymbol{\mbox{x}})\right)}^2}{m}},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

donde <math display="inline">{\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}})</math>  es la predicción para el punto ''i  ''  del modelo Kriging ajustado con la muestra inicial <math display="inline">X</math> , y <math display="inline">{\overset{\mbox{ˆ}}{y}}_{-i}(\boldsymbol{\mbox{x}})</math>  se calcula eliminando el punto de diseño ''i  ''  de la muestra inicial <math display="inline">X</math> , ajustando de nuevo el modelo Kriging con los puntos restantes y obteniendo la predicción en el punto eliminado. Trabajos como los de  [[#bib0110|[22]]] , [[#bib0115|[23]]] , [[#bib0120|[24]]]  and [[#bib0125|[25]]]  han demostrado que, cuando el número de puntos de la muestra inicial es suficiente, la estimación del error cuadrático medio (RMSE) a través del PRESS es un buen indicador para la selección de metamodelos de acuerdo con la precisión de la aproximación. En este trabajo se ha utilizado el PRESS<sub>''RMSE''</sub>  para la selección de las funciones de correlación y de regresión que proporcionan el mejor ajuste. El PRESS<sub>''RMSE''</sub>  resulta un buen indicador para la selección de metamodelos; sin embargo, no resulta un indicador de utilidad para determinar el grado de precisión global del metamodelo. En este trabajo se ha utilizado una estimación del R<sup>2</sup>  mediante errores de validación cruzada:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{PRESS}_{R^2}=1-\frac{{\sum }_{i=1}^m{\left({\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}})-{\overset{\mbox{ˆ}}{y}}_{-i}(\boldsymbol{\mbox{x}})\right)}^2}{{\sum }_{i=1}^m{\left({\overset{\mbox{ˆ}}{y}}_i(\boldsymbol{\mbox{x}})-\overline{y}(\boldsymbol{\mbox{x}})\right)}^2},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

donde <math display="inline">\overline{y}(\boldsymbol{\mbox{x}})</math>  es la media de las predicciones.      

===2.2. Actualización de los modelos Kriging===

La estimación del coeficiente de determinación <math display="inline">{PRESS}_{R^2}</math>  se ha utilizado como indicador de la precisión global del modelo Kriging. Para valores de <math display="inline">{PRESS}_{R^2}\leq 0.98</math>  se actualiza el modelo Kriging con puntos situados en las zonas de peor ajuste. Los nuevos puntos se obtienen mediante la resolución del problema de optimización:

<span id='eq0080'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{\boldsymbol{\mbox{x}}}{m\acute{a}x}\in {\mathbb{R}}^n\quad MSE(\boldsymbol{\mbox{x}}),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10a)
|}

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s.\quad a:\left|\boldsymbol{\mbox{x}}-\boldsymbol{\mbox{x}}_\boldsymbol{\mbox{i}}\right|\geq \epsilon ,\quad \boldsymbol{\mbox{x}}_\boldsymbol{\mbox{i}}\in X,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10b)
|}

donde MSE('''x''' ) es el error cuadrático medio, que en este trabajo se toma como la varianza de la predicción en el punto '''x'''  definida en [[#eq0015|(3)]] . Para evitar problemas de condicionamiento de la matriz de covarianza '''R'''  se incluyen las restricciones [[#eq0085|(10b)]] , que garantizan que la distancia entre los puntos candidatos para la actualización y los existentes en la muestra actual (<math display="inline">\boldsymbol{\mbox{x}}_\boldsymbol{\mbox{i}}\in X</math> ) es superior a un valor predefinido ''ϵ'' .      

==3. Formulación del problema de diseño óptimo robusto==

El problema de diseño óptimo robusto se formula como un problema de optimización multiobjetivo donde las funciones objetivo son los 2 primeros momentos estadísticos de la respuesta estructural ''f''<sub>'''x''' ,'''z'''</sub> :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{\boldsymbol{\mbox{x}}\in {\mathbb{R}}^\boldsymbol{\mbox{n}}}{min}\lbrace {\mu }_{f(\boldsymbol{\mbox{x},\mbox{z}})}(\boldsymbol{\mbox{x}}),{\sigma }_{f(\boldsymbol{\mbox{x},\mbox{z}})}(\boldsymbol{\mbox{x}})\rbrace ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11a)
|}

<span id='eq0095'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s.\quad a:{\mu }_{g_j(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})+</math><math>{\beta }_j{\sigma }_{g_j(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})\leq 0,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11b)
|}

<span id='eq0100'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\sigma }_{h_k(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})\leq {\sigma }_k^T,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11c)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{x}}^{inf}\leq \boldsymbol{\mbox{x}}\leq \boldsymbol{\mbox{x}}^{sup},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11d)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>j=1,\cdots ,m_d,\quad {\beta }_j>0,</math>
|-
|<math>k=m_d+1,2,\cdots ,m,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

donde <math display="inline">\boldsymbol{\mbox{x}}\in {\mathbb{R}}^n</math>  es un vector que contiene los valores deterministas de las variables de diseño, o las medias y las desviaciones estándar en el caso de que las variables sean aleatorias, <math display="inline">\boldsymbol{\mbox{z}}\in {\mathbb{R}}^p</math>  es un vector de parámetros aleatorios, y '''x'''<sup>inf</sup>  y '''x'''<sup>sup</sup>  son los límites inferior y superior de las variables de diseño. Las restricciones del tipo [[#eq0095|(11b)]]  actúan como restricciones de probabilidad bajo la hipótesis de una distribución de probabilidad normal de la función ''g''<sub>''j''</sub> ('''x''' , '''z''' ). El coeficiente ''β''<sub>''j''</sub>  representa el índice de fiabilidad para la restricción ''j-ésima'' , que marca la probabilidad de que la restricción original sea satisfecha. Las restricciones del tipo  [[#eq0100|(11c)]]  limitan la desviación estándar de las restricciones de igualdad ''h''<sub>''k''</sub> ('''x''' , ''z'' ).      

En el problema de diseño óptimo robusto las 2 funciones por minimizar son habitualmente conflictivas, en el sentido de que no existe un punto de diseño que sea óptimo para ambas. En este caso, el criterio de optimalidad utilizado en la optimización mono-objetivo es sustituido por el óptimo de Pareto. Un diseño es óptimo de Pareto para el problema de diseño óptimo robusto si, y solo si, no existe ningún otro diseño factible que mejore su valor en media, sin degradar el valor de la desviación estándar o viceversa. El conjunto de óptimos de Pareto representa la solución del problema de optimización multiobjetivo y geométricamente componen el conocido como frente de Pareto.

Para la obtención del conjunto de óptimos de Pareto es práctica habitual sustituir el conjunto de funciones objetivo por una única función de escalarización. En este trabajo se ha utilizado la función de escalarización aumentada de Tchebycheff [[#bib0130|[26]]] , que para el problema de diseño óptimo robusto es:

<span id='eq0115'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{\boldsymbol{\mbox{x}}}{m\acute{\imath}n}\quad m\acute{a}x[{\omega }_1C_1(\boldsymbol{\mbox{x}})+</math><math>\Phi (\boldsymbol{\mbox{x}}),\quad {\omega }_2C_2(\boldsymbol{\mbox{x}})+</math><math>\Phi (\boldsymbol{\mbox{x}})],</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>siendo:</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C_1(\boldsymbol{\mbox{x}})=\vert {\mu }_{f(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})-</math><math>{\psi }_1^\ast \vert ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12c)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C_2(\boldsymbol{\mbox{x}})=\vert {\sigma }_{f(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})-</math><math>{\psi }_2^\ast \vert ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12d)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\Phi (\boldsymbol{\mbox{x}})=\rho \left[C_1(\boldsymbol{\mbox{x}})+\right. </math><math>\left. C_2(\boldsymbol{\mbox{x}})\right],</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12e)
|}

donde ''ρ  ''  es un número positivo suficientemente pequeño (en este trabajo 0,05), <math display="inline">{\psi }_1^\ast </math>  y <math display="inline">{\psi }_2^\ast </math>  son el valor de mínima media y de mínima desviación estándar y ''ω''<sub>1</sub>  y ''ω''<sub>2</sub>  son pesos con valores comprendidos en el intervalo [0,1] que cumplen ''ω''<sub>1</sub>  + ''ω''<sub>2</sub>  = 1. Minimizando la función de escalarización aumentada de Tchebycheff es posible obtener óptimos de Pareto incluso en la zona cóncava del frente, a diferencia de otras funciones de escalarización tales como la de sumas ponderadas  [[#bib0010|[2]]]  and [[#bib0130|[26]]] .      

==4. Metodología propuesta==

En la metodología propuesta el modelo Kriging aproxima, de forma simultánea, la respuesta estructural en el dominio de diseño y en el dominio estocástico. Esta metodología permite: 1) evaluar la robustez del diseño utilizando un número reducido de evaluaciones en comparación con los procedimientos convencionales; 2) desacoplar los procesos de cuantificación de incertidumbre y de optimización, reduciendo el número de evaluaciones del simulador en comparación con el enfoque anidado; y 3) mejorar la precisión de la aproximación en las regiones cercanas a los óptimos de Pareto.

La metodología propuesta consta de las siguientes 5 etapas:

<span id='sec0035'></span>
===4.1. Aproximación global de la respuesta en el dominio de las variables de diseño y de los parámetros aleatorios, mediante===
* La realización de un diseño de experimentos en el dominio de las variables de diseño y de los parámetros aleatorios <math display="inline">\left[\left(\boldsymbol{\mbox{x}}_1,\boldsymbol{\mbox{z}}_1),\ldots ,(\boldsymbol{\mbox{x}}_n,\boldsymbol{\mbox{z}}_n\right)\right]</math> ). En este trabajo se ha utilizado un hipercubo latino óptimo obtenido mediante la función ''lhsdesign''  de la ''toolbox''  de estadística de Matlab  [[#bib0135|[27]]]  con el criterio de optimalidad ''maximin''  (maximización de la mínima distancia entre puntos)  [[#bib0140|[28]]] .                          
* La evaluación de la respuesta estructural, con el simulador, para cada punto de diseño:<math display="inline">\left[f(\boldsymbol{\mbox{x}}_1,\boldsymbol{\mbox{z}}_1),\ldots ,f(\boldsymbol{\mbox{x}}_n,\boldsymbol{\mbox{z}}_n)\right]</math> .                          
* El ajuste de un modelo Kriging de la respuesta estructural <math display="inline">\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})</math> .                          
* La actualización de los modelos Kriging a través de la maximización de la varianza en la predicción (10).

===4.2. Cuantificación de la incertidumbre===

Los momentos estadísticos de la respuesta estructural se obtienen mediante una simulación MonteCarlo, utilizando los modelos Kriging obtenidos en la etapa [[#sec0035|4.1]] . En este trabajo el número de evaluaciones en la simulación MonteCarlo ha sido de 50.000.      

<span id='sec0045'></span>
===4.3. Aproximación global de los momentos estadísticos de la respuesta en el dominio de las variables de diseño, mediante===
* La realización de un diseño de experimentos en el dominio de las variables de diseño <math display="inline">\left[\boldsymbol{\mbox{x}}_1,\ldots ,\boldsymbol{\mbox{x}}_m\right]</math> ). En este trabajo se ha utilizado un hipercubo latino óptimo mediante el procedimiento definido en [[#sec0035|4.1]] .                          
* La evaluación de la media y de la desviación estándar de las respuestas estructurales para cada punto de diseño:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[{\mu }_{\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}}_1),\ldots ,{\mu }_{\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}}_m)\right],</math>
|-
|<math>\left[{\sigma }_{\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}}_1),\ldots ,{\sigma }_{\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}}_m)\right].</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}
* El ajuste de un modelo Kriging de los momentos estadísticos <math display="inline">{\tilde{\mu }}_{\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})</math> , <math display="inline">{\tilde{\sigma }}_{\tilde{f}(\boldsymbol{\mbox{x}},\boldsymbol{\mbox{z}})}(\boldsymbol{\mbox{x}})</math> .                          
* La actualización de los modelos Kriging mediante la maximización de la varianza en la predicción (10).

===4.4. Optimización multiobjetivo===

El problema de optimización multiobjetivo se resuelve mediante la función de escalarización aumentada de Tchebycheff [[#eq0115|(12a)]] , utilizando los modelos Kriging obtenidos en el apartado [[#sec0045|4.3]] . En este trabajo, el problema de optimización [[#eq0115|(12a)]]  se ha resuelto mediante un algoritmo genético, más concretamente con la función ''ga''  de la ''toolbox''  de optimización de Matlab  [[#bib0135|[27]]] .      

===4.5. Actualización de los modelos Kriging en el dominio estocástico===

Los modelos Kriging obtenidos en [[#sec0035|4.1]]  se actualizan en el dominio de las variables de diseño utilizando el conjunto de soluciones óptimas de Pareto '''x'''<sup>*</sup> . Para actualizar los modelos Kriging en el dominio estocástico se buscan puntos situados en zonas donde la calidad de la aproximación es pobre y que son de gran importancia para evaluar los momentos estadísticos (zonas de elevada densidad de probabilidad). Los puntos candidatos vienen dados por la resolución del problema de optimización en el dominio estocástico Ω:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{m\acute{a}x}_{\boldsymbol{\mbox{z}}\in \Omega }\left(MSE(\boldsymbol{\mbox{x}}^\ast ,\boldsymbol{\mbox{z}})p_Z(\boldsymbol{\mbox{z}})\right),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s.\quad a:\left|\left\{\boldsymbol{\mbox{x}}^\ast ,\boldsymbol{\mbox{z}}\rbrace -\right. \right. </math><math>\left. \left. \lbrace \boldsymbol{\mbox{x}}_\boldsymbol{\mbox{i}},\boldsymbol{\mbox{z}}_\boldsymbol{\mbox{i}}\right\}\right|\geq \epsilon ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>i=1,\ldots ,n,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

donde MSE('''x'''<sup>*</sup> , '''z''' ) es el error cuadrático medio, que en este trabajo se toma como la varianza de la predicción en el punto ('''x'''<sup>*</sup> , '''z''' ) definida en [[#eq0015|(3)]] , ''p''<sub>''Z''</sub> ('''z''' ) es la densidad de probabilidad de los parámetros aleatorios, y ''n''  es el número de puntos de la muestra inicial.      

Los puntos candidatos para la actualización del metamodelo ('''x'''<sup>*</sup> , '''z'''<sup>*</sup> ) se añaden a la muestra original, que se utiliza para la actualización de los modelos Kriging.      

==5. Aplicación numérica==

===5.1. Función de Branin===

El objetivo de este primer ejemplo es mostrar la aplicabilidad y la precisión de la metodología propuesta mediante la resolución del problema de diseño óptimo robusto de una función analítica, la función de Branin [[#bib0145|[29]]] . Esta función se utiliza habitualmente como test en problemas de optimización global y se ha considerado en este trabajo por presentar mínimos locales y por ser claramente multimodal en la variable ''x'' , lo que da lugar a un amplio conjunto de soluciones para el problema de diseño óptimo robusto. La función de Branin es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>f(x,z)={\left(z-\frac{5,1}{4{\pi }^2}x^2+\frac{5}{\pi }x-\right. </math><math>\left. 6\right)}^2+10\left(1-\frac{1}{8\pi }\right)cos(x)+</math><math>10.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

El espacio de diseño está limitado a −5 ≤ ''x''  ≤ 10. Se considera como fuente de incertidumbre el parámetro aleatorio ''z'' , que varía según una distribución normal con media ''μ''<sub>''z''</sub>  = 5 y desviación estándar ''σ''<sub>''z''</sub>  = 2. El problema de diseño óptimo robusto es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{x}{m\acute{\imath}n}\lbrace {\mu }_{f(x,z)}(x),{\sigma }_{f(x,z)}(x)\rbrace ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>-5\leq x\leq 10.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15b)
|}

El modelo Kriging inicial <math display="inline">\tilde{f}(x,z)</math>  se ha obtenido a partir de un hipercubo latino óptimo de 40 puntos en el dominio (''x'' , ''z'' ). La minimización de PRESS<sub>''RMSE''</sub>[[#eq0070|(8)]]  da como resultado una función de correlación ''gaussiana''  y una función de regresión constante. La muestra inicial de 40 puntos se ha incrementado con 5 puntos adicionales en las zonas de mayor error  [[#eq0080|(10a)]] . El coeficiente de determinación estimado es <math display="inline">{PRESS}_{R^2}=</math><math>0</math> ,989, mientras que el obtenido a partir de una muestra de 10.000 puntos es R<sup>2</sup>  = 0,999.      

La función real y el modelo Kriging resultante se muestran en la [[#fig0005|figura 1]] . Los 2 primeros momentos estadísticos de la respuesta se han aproximado en el dominio de diseño a partir de un hipercubo latino óptimo de 30 puntos, mejorado con 10 puntos adicionales en las zonas de más error [[#eq0080|(10a)]] .

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr1.jpg|center|375px|Función de Branin. (a) función real, (b) modelo Kriging.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Función de Branin. (a) función real, (b) modelo Kriging.

</span>
|}

El problema de optimización multiobjetivo se ha resuelto utilizando la función de escalarización aumentada de Tchebycheff, adoptando un vector de pesos de tamaño 10, cuyos valores varían uniformemente entre 0 y 1. Los parámetros utilizados para el algoritmo genético son: número de repeticiones (8), tamaño de la población (24), probabilidad de cruce (80%), elitismo manteniendo los 2 mejores diseños y número máximo de generaciones (80). El conjunto de diseños óptimos obtenidos tras 3 actualizaciones del modelo Kriging forma el frente de Pareto de la [[#fig0010|figura 2]] . Este frente permite obtener una solución de compromiso entre el diseño de mínima media (''x''<sup>*</sup>  = 9, 915) y el diseño de mínima desviación estándar (''x''<sup>*</sup>  = 0, 695). En este caso en la [[#fig0010|figura 2]]  se ha seleccionado el punto de menor distancia al origen. La [[#fig0015|figura 3]]  muestra los óptimos de Pareto sobre los modelos Kriging de la media y de la desviación estándar. Para comprobar el error cometido con la aproximación de los momentos estadísticos, se ha realizado una simulación de MonteCarlo con una muestra de 50.000 puntos, que se evaluaron con la función real y con el modelo Kriging. La [[#tbl0010|tabla 2]]  muestra los resultados para los diseños de mínima media y de mínima desviación estándar. Los errores en ambos diseños son inferiores al 0,1% en media y al 0,15% en desviación estándar.

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr2.jpg|center|372px|Función de Branin. Frente de Pareto, soluciones de mínima media, mínima ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Función de Branin. Frente de Pareto, soluciones de mínima media, mínima desviación estándar y solución de compromiso seleccionada.

</span>
|}

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr3.jpg|center|378px|Función de Branin. Óptimos de Pareto sobre los modelos Kriging de la media (a) y ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Función de Branin. Óptimos de Pareto sobre los modelos Kriging de la media (a) y de la desviación estándar (b).

</span>
|}

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2.

Errores en la aproximación de los momentos estadísticos

|-

! 
! 
! Kriging
! Real
! Error (%)
|-

| ''x''<sup>*</sup>  = 9.915                                                    
| ''μ''<sub>''f'' (''x'' ,''z'' )                                                        </sub> (''x'' )                                                    
| 9.887
| 9.895
| −0,081
|-

| 
| ''σ''<sub>''f'' (''x'' ,''z'' )                                                        </sub> (''x'' )                                                    
| 10.113
| 10.127
| −0,138
|-

| ''x''<sup>*</sup>  = 0.695                                                    
| ''μ''<sub>''f'' (''x'' ,''z'' )                                                        </sub> (''x'' )                                                    
| 21.408
| 21.410
| −0,009
|-

| 
| ''σ''<sub>''f'' (''x'' ,''z'' )                                                        </sub> (''x'' )                                                    
| 5.700
| 5.706
| −0,105
|}

===5.2. Estructura de 4 barras===

El segundo ejemplo numérico es la estructura de 4 barras [[#bib0010|[2]]]  mostrada en la [[#fig0020|figura 4]] .

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr4.jpg|center|342px|Estructura de 4 barras.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Estructura de 4 barras.

</span>
|}

En el nudo 4 de la estructura se aplica una fuerza horizontal ''P'' . Las áreas ''A''<sub>1</sub>  (barras 1 y 3) y ''A''<sub>2</sub>  (barras 2 y 4) se consideran variables de diseño deterministas. Los módulos de Young ''E''<sub>1</sub>  (barras 1 y 3) y ''E''<sub>2</sub>  (barras 2 y 4) se consideran variables aleatorias que siguen distribuciones normales con medias <math display="inline">{\mu }_{E_1}=</math><math>210</math>  y <math display="inline">{\mu }_{E_2}=</math><math>100</math>  y desviaciones estándar <math display="inline">{\sigma }_{E_1}=</math><math>21</math>  y <math display="inline">{\sigma }_{E_2}=</math><math>5</math> . La densidad del material es ''ρ''  = 1. El problema consiste en la minimización del desplazamiento horizontal (''u  '' ) del nudo 4, sujeto a una restricción de masa <math display="inline">w\leq 5</math> . La formulación del problema de diseño óptimo robusto es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{A_1,A_2}{m\acute{\imath}n}\lbrace {\mu }_u(A_1,A_2),{\sigma }_u(A_1,A_2)\rbrace ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s.\mbox{   }a:w\leq 5,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>0\leq A_{1,2}\leq 2.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16c)
|}

El modelo Kriging inicial se ha obtenido a partir de un hipercubo latino óptimo de 100 puntos, mejorado con 50 puntos adicionales en las regiones de mayor error. Las funciones de correlación y de regresión resultantes de la minimización del PRESS<sub>''RMSE''</sub>[[#eq0070|(8)]]  son ''gaussiana''  y cuadrática, respectivamente. La respuesta real y el modelo Kriging resultante se muestran en la  [[#fig0025|figura 5]] . El coeficiente de determinación estimado es <math display="inline">{PRESS}_{R^2}=</math><math>0</math> ,986, mientras que el obtenido a partir de una muestra de 10.000 puntos es R<sup>2</sup>  = 0,984.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr5.jpg|center|341px|Estructura de 4 barras. (a) desplazamiento real, (b) modelo Kriging.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Estructura de 4 barras. (a) desplazamiento real, (b) modelo Kriging.

</span>
|}

El problema de optimización multiobjetivo se ha resuelto utilizando la función de escalarización aumentada de Tchebycheff, adoptando un vector de pesos de tamaño 10, cuyos valores varían uniformemente entre 0 y 1. Los parámetros utilizados para el algoritmo genético son: número de repeticiones (8), tamaño de población (40), probabilidad de cruce (80%), elitismo manteniendo los 2 mejores diseños y número máximo de generaciones (100).

La [[#fig0030|figura 6]]  muestra el frente de Pareto obtenido a partir de una población inicial de 10 diseños, tras 10 actualizaciones del modelo Kringing inicial. Las figuras [[#fig0035|7]]  muestran los óptimos de Pareto sobre los modelos Kriging de la media y de la desviación estándar del desplazamiento. Todos los diseños obtenidos están situados en la frontera que delimita la región factible y varían desde el diseño óptimo de mínima media (''A''<sub>1</sub> =1,766, ''A''<sub>2</sub> =0,002) hasta el diseño óptimo de mínima desviación estándar (''A''<sub>1</sub> =0,3685, ''A''<sub>2</sub> =1,9789).

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr6.jpg|center|364px|Estructura de 4 barras. Frente de Pareto, soluciones de mínima media, mínima ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Estructura de 4 barras. Frente de Pareto, soluciones de mínima media, mínima desviación estándar y solución de compromiso seleccionada.

</span>
|}

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr7.jpg|center|374px|Estructura de 4 barras. Óptimos de Pareto sobre los modelos Kriging de la media ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Estructura de 4 barras. Óptimos de Pareto sobre los modelos Kriging de la media (a) y de la desviación estándar (b).

</span>
|}

La [[#fig0040|figura 8]]  muestra las funciones de densidad de probabilidad y acumulada de probabilidad, reales y aproximadas, para los diseños óptimos de mínima media y de mínima desviación estándar. Los errores cometidos en la aproximación de los momentos estadísticos se muestran en la [[#tbl0015|tabla 3]] .

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr8.jpg|center|362px|Función de densidad de probabilidad y función acumulada de probabilidad. (a) ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

Función de densidad de probabilidad y función acumulada de probabilidad. (a) Diseño de mínima media, (b) Diseño de mínima desviación estándar.

</span>
|}

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3.

Errores en la aproximación de los momentos estadísticos

|-

! 
! 
! Kriging
! Real
! Error (%)
|-

| ''A''<sub>1</sub>  = 1, 766                                                    
| ''μ''<sub>''u''</sub> (''A''<sub>1</sub> , ''A''<sub>2</sub> )                                                    
| 3.843e-3
| 3.846e-3
| 0,073
|-

| ''A''<sub>2</sub>  = 0, 002                                                    
| ''σ''<sub>''u''</sub> (''A''<sub>1</sub> , ''A''<sub>2</sub> )                                                    
| 4.007e-4
| 3.961e-4
| 1,151
|-

| ''A''<sub>1</sub>  = 0, 369                                                    
| ''μ''<sub>''u''</sub> (''A''<sub>1</sub> , ''A''<sub>2</sub> )                                                    
| 3.965e-3
| 3.967e-3
| 0,033
|-

| ''A''<sub>2</sub>  = 1, 979                                                    
| ''σ''<sub>''u''</sub> (''A''<sub>1</sub> , ''A''<sub>2</sub> )                                                    
| 1.721e-4
| 1.786e-4
| 3,629
|}

===5.3. Celosía de 25 barras con incertidumbre en cargas y material===

En este ejemplo se aplica la metodología propuesta a la resolución del problema de diseño robusto de la celosía de 25 barras mostrada en la [[#fig0045|figura 9]] . Las variables de diseño consideradas son las dimensiones de los elementos estructurales clasificadas en 7 grupos, tal y como muestra la [[#fig0045|figura 9]] . Las variables de diseño están limitadas inferiormente a un valor de 0,6452 cm<sup>2</sup> . La estructura está sometida a unas cargas verticales en los nudos del cordón inferior cuyo valor nominal es <math display="inline">\overline{F}=</math><math>44,48</math>  kN. El material de las barras tiene una densidad ''ρ'' =7.833,41 kg/m<sup>3</sup>  y un módulo de Young <math display="inline">\overline{E}</math> =206.842,8 MPa. La tensión máxima admisible a tracción y compresión es ''σ''<sub>''y''</sub> |<sub>''t'' ,''c''</sub>  = 206.843 MPa. El desplazamiento de todos los nudos en las direcciones x e y está limitado a un valor máximo de 5,08 cm. El problema está parametrizado con 7 variables de diseño (áreas) y 2 variables aleatorias (''ξ''<sub>1</sub>  y ''ξ''<sub>2</sub> ) que siguen una distribución normal. <math display="inline">{\xi }_1\sim N(\mu =</math><math>0,\sigma =0.05)</math>  está relacionada con las fuerzas aplicadas mediante la fórmula <math display="inline">F=</math><math>\overline{F}(1+{\xi }_1)</math> , mientras que <math display="inline">{\xi }_2\sim N(\mu =</math><math>0,\sigma =0.05)</math>  está relacionada con el módulo de Young mediante la fórmula <math display="inline">E=</math><math>\overline{E}(1+{\xi }_2)</math> . Siguiendo la línea de Lagaros et al. [[#bib0150|[30]]] , el problema de diseño robusto consiste en la minimización del peso de la estructura y de la varianza de la respuesta estructural representada por la desviación estándar del desplazamiento vertical del nudo 10. Se analiza: 1) la precisión de los modelos Kriging de la respuesta estructural en el dominio estocástico y de diseño; 2) la precisión de la metodología propuesta para la obtención de los momentos estadísticos, 3) la precisión de los diseños óptimos obtenidos en comparación con métodos convencionales.

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr9.jpg|center|361px|Celosía de 25 barras.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

Celosía de 25 barras.

</span>
|}

El modelo Kriging inicial se ha obtenido a partir de un hipercubo latino óptimo de 200 puntos. Las funciones de correlación y de regresión resultantes de la minimización del PRESS<sub>''RMSE''</sub>[[#eq0070|(8)]]  son exponencial-''gaussiana  ''  y cuadrática, respectivamente. El coeficiente de determinación estimado es <math display="inline">{PRESS}_{R^2}</math> =0,9973, mientras que el obtenido a partir de una muestra de 10.000 puntos es R<sup>2</sup>  = 0, 9971. La [[#fig0050|figura 10]]  muestra el modelo Kriging de la desviación estándar del desplazamiento vertical del nudo 10 en el espacio de diseño.

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr10.jpg|center|354px|Celosía de 25 barras. Modelo Kriging de la desviación estándar de uy,10.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 10.

Celosía de 25 barras. Modelo Kriging de la desviación estándar de ''u''<sub>''y'' ,10                                                </sub> .                  

</span>
|}

El problema de optimización multiobjetivo se ha resuelto utilizando la función de escalarización aumentada de Tchebycheff, adoptando un vector de pesos de tamaño 10. El frente de Pareto resultante tras 8 iteraciones se muestra en la [[#fig0055|figura 11]]  (a). El resultado obtenido se compara con el conseguido con métodos convencionales existentes en la bibliografía. En este trabajo se ha seleccionado el método ''Univariate Dimension Reduction''  (UDR)  [[#bib0155|[31]]]  and [[#bib0160|[32]]]  por ser uno de los métodos más eficientes para la obtención de los momentos estadísticos en procesos de diseño óptimo robusto [[#bib0165|[33]]] . La [[#fig0055|figura 11]]  (a) muestra el frente de Pareto resultante estimando los momentos estadísticos mediante el método UDR. En esta figura se observa un buen acuerdo entre los frentes de Pareto obtenidos mediante el método UDR y la metodología propuesta. El número de evaluaciones del modelo de simulación necesarias utilizando el método UDR fueron 3,6e5 de acuerdo con la expresión: (4 × ''N''<sub>''va''</sub>  + 1) × ''N''<sub>''opt''</sub>  × ''N''<sub>''pop''</sub>  × ''N''<sub>''gen''</sub> . ''N''<sub>''va''</sub>  es el número de variables aleatorias de problema, ''N''<sub>''opt''</sub>  el tamaño de vector de pesos, ''N''<sub>''pop''</sub>  el tamaño de la población, y ''N''<sub>''gen''</sub>  el número de generaciones. El número de evaluaciones del modelo de simulación necesarias con la metodología propuesta (280) es considerablemente inferior.

<span id='fig0055'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_336370485-1-s2.0-S0213131513000461-gr11.jpg|center|332px|(a) Frentes de Pareto obtenidos utilizando modelos Kriging y UDR (Univariate ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 11.

(a) Frentes de Pareto obtenidos utilizando modelos Kriging y UDR ''(Univariate Dimension Reduction)'' . (b) Errores cometidos en la aproximación de la desviación estándar mediante modelos Kriging para los óptimos de Pareto.                  

</span>
|}

La [[#fig0055|figura 11]]  (b) representa, para cada diseño óptimo de Pareto, el error cometido en la estimación de la desviación estándar mediante el modelo Kriging en comparación con una simulación MonteCarlo de 1e5 simulaciones. Los errores oscilan entre −3% y 5%.      

==6. Conclusiones==

Se ha propuesto una metodología para la resolución del problema de diseño óptimo robusto que utiliza modelos Kriging para resolver, de forma eficiente, el problema de cuantificación de incertidumbre en la optimización. A diferencia de otros trabajos, en la metodología propuesta el modelo Kriging aproxima, de forma simultánea, la respuesta estructural en el dominio de diseño y en el dominio estocástico. Asimismo, la metodología incluye un criterio de actualización de los modelos Kriging, basado en la estimación del error en la predicción, que permite mejorar la calidad de la aproximación en las regiones cercanas al frente de Pareto.

La metodología propuesta permite: 1) obtener los indicadores de robustez del diseño (momentos estadísticos de la respuesta estructural) utilizando un número reducido de evaluaciones en comparación con los procedimientos convencionales; 2) desacoplar los procesos de cuantificación de incertidumbre y de optimización, resolviendo el problema con un número de evaluaciones del simulador muy inferior al enfoque anidado, y 3) mejorar la precisión de la aproximación en las regiones cercanas a los óptimos de Pareto hasta alcanzar un nivel deseado.

Para demostrar la aplicabilidad y la precisión de la metodología propuesta se han resuelto un ejemplo matemático y 2 ejemplos de diseño óptimo robusto de estructuras de nudos articulados. La metodología propuesta no está limitada a problemas de diseño óptimo robusto de estructuras de nudos articulados. Su principal campo de aplicación es el diseño óptimo robusto de problemas de ingeniería que requieren un elevado coste computacional para su resolución debido a: 1) la utilizacion de software ''in house''  o comercial, que obliga en muchos casos a utilizar estrategias de cuantificacion de incertidumbre no intrusivas (MonteCarlo), y 2) la utilización de complejos modelos de simulación de alto coste computacional.      

Los resultados obtenidos muestran que los modelos Kriging son una herramienta eficaz para aproximar los momentos estadísticos de la respuesta estructural con un bajo coste computacional, y para actuar como soporte de la optimización en la búsqueda de óptimos globales. Asimismo, la reducción que se consigue en el coste computacional los hace muy adecuados para su aplicación a problemas de ingeniería que requieren modelos complejos de simulación.

Una ventaja adicional de la metodología propuesta es que las aproximaciones globales creadas durante la optimización se pueden reutilizar en nuevos procesos de optimización, en exploración de los espacios de diseño y estocástico, o en aplicaciones computacionalmente exigentes. Quedan para futuras investigaciones el estudio de nuevas técnicas de metamodelos que permitan su aplicación a problemas de alta dimensionalidad, así como su aplicación a otras formulaciones del problema de optimización bajo incertidumbre tales como el diseño óptimo basado en fiabilidad.

==Agradecimientos==

Este trabajo se ha desarrollado con el apoyo financiero del Ministerio Español de Economía y Competitividad, en el proyecto de investigación DPI2011-26394 «Sistema integrado de diseño óptimo robusto de topología de estructuras».

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] M. Valdebenito, G. Schuëller; A survey on approaches for reliability-based optimization; Struct. Multidiscip. Optim., 42 (2010), pp. 645–663 10.1007/s00158-010-0518-6</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] I. Doltsinis, Z. Kang; Robust design of structures using optimization methods; Comput. Meth. Appl. Mech. Eng., 193 (2004), pp. 2221–2237</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] J. Díaz, S. Hernández, L. Romera, A.N. Fontán; Diseño y análisis térmico bajo incertidumbe de estructuras aeronaúticas; Rev. Int. Métodos Numér. Cálc. Diseño ing., 17 (2011), pp. 95–104</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] J. Pons-Prats, G. Bugeda, F. Zárate, E. Oñate; Optimización robusta en aplicaciones aeronáuticas con la combinación de cálculo estocástico y algoritmos evolutivos; Rev. Int. Métodos Numér. Cálc. Diseño ing., 28 (2012), pp. 18–32</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] K. Lee, D. Kang; A robust optimization using the statistics based on Kriging metamodel; J. Mech. Sci. Technol., 20 (2006), pp. 1169–1182</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] J. Martínez, P. Martí, Diseño óptimo robusto de estructuras basado en modelos Kriging, en: Congress on Numerical Methods in Engineering, Coimbra, 14 to 17 June, pp. 447-448.</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] J. Martínez, P. Martí, Metamodel-based multi-objective robust design optimization of structures, en: Computer Aided Optimum Design in Engineering XII, New Forest, UK, 20 - 22 June, pp. 35-45.</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] F. Jurecka, M. Ganser, K. Bletzinger; Update scheme for sequential spatial correlation approximations in robust design optimisation; Comput. Struct., 85 (2007), pp. 606–614</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] G. Dellino, J. Kleijnen, C. Meloni; Robust optimization in simulation: Taguchi and Response Surface Methodology; Int. J. Prod. Econ., 125 (2010), pp. 52–59</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] X. Guoqi, H. Dagui; A metamodeling method based on support vector regression for robust optimization; CJME, 23 (2010)</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] R. Jin, X. Du, W. Chen; The use of metamodeling techniques for optimization under uncertainty; Struct. Multidiscip. Optim., 25 (2003), pp. 99–116</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] H. Beyer, B. Sendhoff; Robust optimization - A comprehensive survey; Comput. Meth. Appl. Mech. Eng., 196 (2007), pp. 3190–3218</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] J. Kleijnen; A comment on blannings metamodel for sensitivity analysis: The regression metamodel in simulation; Interfaces, 5 (1975), pp. 21–23</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] R. Myers, D. Montgomery; Response surface methodology;  (2nd ed.)New York, Wiley-Interscience (2002)</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] T. Simpson, J. Poplinski, P.N. Koch, J. Allen; Metamodels for computer-based engineering design: Survey and recommendations; Eng. Comput., 17 (2001), pp. 129–150 10.1007/PL00007198</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] D.G. Krige; A statistical approach to some basic mine valuation problems on the Witwatersrand, Journal of the Chemical; Metallurgical and Mining Society of South Africa, 52 (1951), pp. 119–139</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] H. Nakayama, M. Arakawa, R. Sasaki; Simulation-based optimization using computational intelligence; Optim. Eng., 3 (2002), pp. 201–214 10.1023/A:1020971504868</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] M. Papadrakakis, N.D. Lagaros, Y. Tsompanakis; Structural optimization using evolution strategies and neural networks; Comput. Meth. Appl. Mech. Eng., 156 (1998), pp. 309–333</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] D.R. Jones; A taxonomy of global optimization methods based on response surfaces; J. Global Optim., 21 (2001), pp. 345–383</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] J. Sacks, W.J. Welch, T.J. Mitchell, H.P. Wynn; Design and analysis of computer experiments; Stat. Sci., 4 (1989), pp. 409–435</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] M. Meckesheimer, A.J. Booker, R.R. Barton, T.W. Simpson; Computationally inexpensive metamodel assessment strategies; AIAA, 40 (2002), pp. 2053–2060</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] F. Viana, R.T. Haftka; Cross Validation Can Estimate How Well Prediction Variance Correlates with Error; AIAA, 47 (2009), pp. 2266–2270</li>
<li><span id='bib0115'></span>
[[#bib0115|[23]]] F. Viana, R.T. Haftka; V.S. Jr, Multiple surrogates: how cross-validation errors can help us to obtain the best predictor; Struct. Multidiscip. Optim., 39 (2009), pp. 439–457</li>
<li><span id='bib0120'></span>
[[#bib0120|[24]]] F. Viana, R. T. Haftka, V. Steffen, S. Butkewitsch, M. F. Leal, Ensemble of Surrogates: a Framework based on Minimization of the Mean Integrated Square Error, en: 49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials, Schaumburg, IL, pp. 1-27.</li>
<li><span id='bib0125'></span>
[[#bib0125|[25]]] F. Viana, V. Picheny, R.T. Haftka; Using Cross Validation to Design Conservative Surrogates; AIAA, 48 (2010), pp. 2286–2298</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] K. Miettinen; Nonlinear Multiobjective Optimization; Kluver Academic Publishers, Dordrecht (1999)</li>
<li><span id='bib0135'></span>
[[#bib0135|[27]]] Matlab, The language of technical computing. (R2011a), MathWorks Inc., 2011.</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] M. Johnson, L. Moore, D. Ylvisaker; Minimax and maximin distance designs; J. Stat. Plann. Infer., 26 (1990), pp. 131–148</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] F. Branin; Widely convergent method for finding multiple solutions of simultaneous nonlinear equations; IBM J. Res. Dev., 16 (1972), pp. 504–522</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] N.D. Lagaros, V. Plevris, M. Papadrakakis; Reliability based robust design optimization of steel structures; Int. J. Simul. Multi. Design Optim., 1 (2007), pp. 19–29</li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] S. Rahman, H. Xu; A univariate dimension-reduction method for multi-dimensional integration in stochastic mechanics; Probab. Eng. Mech., 19 (2004), pp. 393–408</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] H. Xu, S. Rahman; A generalized dimension-reduction method for multidimensional integration in stochastic mechanics; Int. J. Numer. Methods Eng., 61 (2004), pp. 1992–2019</li>
<li><span id='bib0165'></span>
[[#bib0165|[33]]] I. Lee, K.K. Choi, L. Du, D. Gorsich; Dimension reduction method for reliability-based robust design optimization; Comput. Struct., 86 (2008), pp. 1550–1562</li>
</ol>