You do not have permission to edit this page, for the following reason:
You can view and copy the source of this page.
x
1
==Optimización de topología robusta de estructuras continuas usando el método de Monte Carlo y modelos Kriging==2
3
4
'''A. Cordero''', '''P. Martí''', '''M. Victoria'''5
6
Departamento de Estructuras y Construcción<br />7
Universidad Politécnica de Cartagena, Cartagena, Murcia, España8
9
10
11
<span style="font-size: 85%;">'''Received''': 2 May 2016<br />12
'''Accepted''': 13 December 2016<br />13
'''Published''': July 2017</span>14
15
16
'''Resumen'''17
18
El objetivo de este trabajo es presentar una nueva metodología eficiente y precisa llamada Monte Carlo y Kriging (MCK) para la optimización de topología robusta. El objetivo es minimizar el valor esperado de la ''compliance'' considerando la existencia de incertidumbre con cargas concentradas. La incertidumbre en la carga puede presentarse en la magnitud, en la dirección y/o en la posición. La evaluación de la función objetivo se realiza utilizando el método de simulación de Monte Carlo en combinación con un modelo Kriging. Para estimar el valor esperado de la ''compliance'', se transforma el problema probabilístico en otro determinístico sujeto a múltiples estados de carga mediante el Método de Monte Carlo pero empleando un reducido número de evaluaciones del modelo de simulación. Para ello es necesario construir un modelo Kriging del modelo de simulación a partir de una pequeña muestra obtenida con un hipercubo latino del espacio de diseño y predecir la ''compliance'' en cada uno de los puntos utilizados por la simulación de Monte Carlo. Dos ejemplos demuestran la precisión y eficiencia del algoritmo. Para verificar el algoritmo propuesto, los problemas también se resuelven mediante el método de Monte Carlo estándar.19
20
'''Palabras clave''': Optimización de topología robusta, incertidumbre de la carga, método de Monte Carlo, modelos Kriging21
22
23
==Robust Topology Optimization of continuum structures using Monte Carlo method and Kriging models==24
25
'''Abstract'''26
27
The aim of this paper is to introduce an efficient and accurate new approach called Monte Carlo and Kriging (MCK) to robust topology optimization. The objective is to minimize the expected value of ''compliance'' under concentrated loading uncertainty. The loading uncertainty may occur in magnitude, direction and/or position. The Monte Carlo simulation method and Kriging model are used to evaluate the objective function. To evaluate the expected value of ''compliance'' the probabilistic problem is transformed into a multiple loading deterministic one using of Monte Carlo method but with a reduced evaluations number of simulation model. A small sample obtained with a Latin Hypercube is used to build a Kriging model of the simulation model. This is utilized to estimate the ''compliance'' in those points used by Monte Carlo simulation method. Two problems are solved to demonstrate the efficiency and accuracy of the approach. The examples are solved again using a standard Monte Carlo simulation to check the proposed approach.28
29
'''Keywords''': Robust topology optimization, loading uncertainty, Monte Carlo method, Kriging models30
31
32
==1. Introducción==33
34
La optimización estructural ha experimentado un importante crecimiento en las últimas décadas debido a aplicaciones numéricas que permiten automatizar el proceso de diseño. La optimización estructural se puede dividir en tres categorías: tamaño, forma y topología. La optimización de topología es la que mejores resultados aporta, ya que proporciona una mayor libertad a la hora de obtener nuevos diseños conceptuales.35
36
La optimización de topología tiene como objetivo encontrar la mejor distribución de material a partir de la forma y la localización de las cavidades para un dominio de diseño dado. Desde el trabajo original de Bendsoe y Kikuchi [1], la optimización de topología ha experimentado un rápido crecimiento, siendo los métodos más utilizados: el Método de homogeneización [1], el Método SIMP (''Solid Isotropic Material with Penalization'') [2], el método ESO (''Evolutionary Structural Optimization'') [3], y el método SA (''Simulated Annealing'') [4]. En los últimos años han surgido nuevas metodologías como el método LSM (''Level Set Method'') [5] y el método ITD (''Isolines Topology Design'') [6]. Una revisión actual sobre los avances realizados en optimización de topología se puede encontrar en Deathon y Grandhi [7]. En el campo industrial, la optimización de topología ha demostrado ser una herramienta potente, aplicándose a un extenso rango de disciplinas físicas tales como la transferencia de calor [8], la acústica [9], la mecánica de fluidos [10], el diseño de materiales [11], etc. También, la optimización de topología se ha convertido en una parte fundamental de muchos softwares como ABAQUS, Genesis, OptiStruct, Permas, ANSYS, etc.37
38
En la mayoría de los casos, el problema de optimización se estudia desde un punto de vista determinista, lo que se denomina como optimización de topología determinista (''Deterministic Topology Optimization'', DTO), donde los diseños son obtenidos sin considerar de forma explícita la influencia de las diversas fuentes de incertidumbre (o variabilidad) presentes en la realidad, tales como variaciones en las cargas, en las propiedades de material, en la geometría, en las condiciones de contorno, etc. Así pues, una formulación determinista puede afectar al rendimiento y operatividad de la estructura, ya que los diseños óptimos obtenidos solo tienen un comportamiento óptimo bajo condiciones operativas cercanas a los valores utilizados en el proceso de optimización. Tradicionalmente, el efecto de las incertidumbres se ha tenido en cuenta mediante la incorporación de factores de seguridad dentro de la formulación; sin embargo, esta técnica heurística puede conducir a estructuras conservadoras o poco eficientes.39
40
Con objeto de obtener estructuras menos sensibles, más racionales y con un mayor rendimiento ante condiciones reales, en los últimos años ha tomado un gran interés el desarrollo de formulaciones capaces de incluir explícitamente aquellas fuentes de incertidumbre que puedan provocar cambios significativos en el comportamiento y en el rendimiento de la estructura. De forma general, se pueden distinguir dos formulaciones, la optimización de topología basada en fiabilidad y la optimización de topología robusta.41
42
En la optimización de topología basada en fiabilidad (''Reability-Based Topology Optimization'', RBTO) el problema de optimización se formula como la minimización de una función objetivo determinista sujeta a restricciones de tipo incierto, donde las restricciones son expresadas a través de probabilidades o modos de fallo. Existe una extensa literatura sobre RBTO, pero aquí solo se recogen algunos de los modos de fallo que han sido considerados como estados críticos en deformación [12, 13], en rigidez [14, 15], y en frecuencias naturales [13, 14]. También se han empleado técnicas diferentes para estimar de modo preciso y eficaz los modos de fallo, tales como los métodos FORM y SORM (''First and Second Order Reliability Methods'') [16], el método de simulación de Monte Carlo [14], el método PMA (''Performance Measure Approach'') [13], el método RSM (''Response Surface Method'') [15] y el método SRSM (''Stochastic Response Surface Method'') [17].43
44
La optimización de topología robusta (''Robust Topology Optimization'', RTO), se formula como la minimización de una función objetivo probabilística sujeta a una serie de restricciones deterministas. La función objetivo probabilística más popular es la suma ponderada de la media y la desviación estándar de la ''compliance'' [18-25], aunque también se utiliza el valor esperado de la ''compliance'' [26, 27]. Otra alternativa consiste en tratar el problema como una optimización multiobjetivo [28], donde se optimiza de forma simultanea los objetivos en conflicto (el valor esperado y la desviación estándar) mediante un frente de Pareto (conjunto de óptimos de Pareto, entendiendo como óptimo una solución tal que no existe otra que mejore el valor del objetivo sin deteriorar el resto).45
46
Aunque las dos formulaciones incluyen la aleatoriedad o la incertidumbre dentro del proceso de optimización, difieren básicamente en que RBTO se centra en satisfacer criterios de fiabilidad de eventos extremos mientras que RTO se centra en reducir la variación de la respuesta estructural ante las fluctuaciones diarias. De modo que RBDO se asocia con el daño inducido por fallo y RTO con el gasto derivado de una pobre calidad en el rendimiento estructural.47
48
En la actualidad, la mayoría de los trabajos publicados en RTO, tratan las incertidumbres en la carga [18, 20, 23-27], aunque algunos de estos consideran también las incertidumbres en las propiedades del material [29, 30] y en la geometría [19, 22, 30].49
50
En este trabajo se utiliza la formulación RTO centrada en las incertidumbres en la carga, ya que pequeñas variaciones en sus parámetros pueden provocar grandes cambios en la respuesta estructural. Conti et al. [31] proponen una metodología donde se combina LSM con técnicas de programación estocástica para la optimización de estructuras con incertidumbre en cargas concentradas. Chen et al. [18] proponen una metodología donde integran LSM con una formulación robusta, no solo para considerar la incertidumbre en cargas concentradas, sino también para cargas distribuidas. Las cargas concentradas se modelan como variables aleatorias estadísticamente independientes, mientras que las cargas distribuidas se modelan como campos aleatorios. La alta dimensionalidad del campo aleatorio se transforma en un conjunto finito de variables aleatorias estadísticamente independientes mediante un método espectral como el desarrollo de Karhunen-Loeve. Finalmente, los momentos estadísticos de la ''compliance'' se evalúan mediante el método UDR (''Univariate Dimension Reduction'') en los puntos de cuadratura de Gauss. Dunning y Kim [26] presentan un método analítico para minimizar el valor esperado de la ''compliance'', con incertidumbre en la magnitud y en la dirección de la carga, usando complejas derivadas numéricas. En el caso más desfavorable, el problema RTO se transforma en otro de 3''n''+1 estados de carga, siendo ''n ''el número de cargas con incertidumbre. Sin embargo, la metodología se limita a variables aleatorias estadísticamente independientes definidas mediante distribuciones normales. Dunning y Kim [20] amplían la metodología anterior incorporando la desviación estándar de la ''compliance'' en la función objetivo. Zhao y Wang [27] proponen una metodología para evaluar el valor esperado de la ''compliance'', que permite utilizar cualquier tipo de función de probabilidad continua para describir la incertidumbre. Para obtener los momentos estadísticos utilizan el método de Monte Carlo y la descomposición de matrices. Zhao et al. [25] presentan un método para la optimización de topología considerando incertidumbre en la carga, utilizando el método de colocación estocástico para obtener los momentos estadísticos. Los puntos de colocación y sus respectivos pesos se determinan mediante la cuadratura del producto tensorial completo y la cuadratura dispersa de Smolyak, respectivamente.51
52
Todas las metodologías descritas anteriormente transforman el problema RTO en otro sujeto a múltiples estados de carga deterministas en donde los momentos estadísticos se calculan evaluando repetidamente la respuesta para cada uno de los estados de carga. A menudo, el coste computacional del análisis es elevado, por lo que la evaluación de la función objetivo de forma precisa y eficiente es crucial en este tipo de problemas.53
54
Actualmente se usan metamodelos (modelos aproximados), tales como el método RSM [32], los modelos RBF (''Radial Basis Function'') [33], el método ANNs (''Artificial Neural Networks)'' [34] o los modelos Kriging [35], para reducir el elevado coste computacional asociado con el cálculo de los momentos estadísticos y de los modelos de simulación. Debido a la capacidad para interpolar exactamente los valores de la respuesta en los puntos de muestreo y la flexibilidad para capturar cualquier tipo de comportamiento no lineal, los modelos Kriging son especialmente adecuados para aproximar simulaciones con un elevado coste computacional. Los modelos Kriging se han utilizado con éxito tanto en optimización basada en fiabilidad [36, 37] como en optimización robusta [38, 39], aunque todavía no se han utilizado para RTO de estructuras continuas. La aplicación de modelos Kriging para el diseño óptimo robusto de estructuras articuladas puede verse en Martínez y Martí [40].55
56
En este trabajo se propone una nueva metodología RTO, eficiente y precisa, que permite minimizar el valor esperado de la ''compliance'' considerando la existencia de incertidumbre en cargas concentradas. La incertidumbre en la carga puede presentarse en la magnitud, la dirección, la posición o una combinación de ellas, pudiendo ser descrita por cualquier tipo de función de densidad de probabilidad (''probability density function'', pdf). La evaluación de la función objetivo se realiza utilizando el método de simulación de Monte Carlo en combinación con un modelo Kriging. El objetivo es poder estimar el valor esperado de la ''compliance'' de forma precisa, transformando el problema RTO en otro sujeto a múltiples estados de carga mediante el método de simulación de Monte Carlo, pero empleando un número reducido de evaluaciones del modelo de simulación. Para ello, se debe construir un modelo Kriging del modelo de simulación a partir de una pequeña muestra obtenida usando un hipercubo latino del espacio aleatorio y obtener la respuesta estructural en los puntos de muestreo utilizados por la simulación de Monte Carlo. La metodología propuesta en este trabajo permite reducir el coste computacional, sobre todo en los casos donde existe un elevado número de cargas con incertidumbre. La metodología se ha aplicado a problemas bidimensionales, utilizando el método SIMP, aunque se pueden utilizar otros métodos de optimización de topología.57
58
El trabajo se ha organizado en cinco apartados: en el apartado 1 se realiza una revisión del estado del arte sobre el problema RTO. En el apartado 2 se describen los fundamentos teóricos usados para el desarrollo de la metodología propuesta: formulación del problema RTO usando el método SIMP, incertidumbres consideradas, método de simulación de Monte Carlo y modelos Kriging. En el apartado 3 se presenta la metodología y el algoritmo propuesto. En el apartado 4 se muestran los resultados obtenidos para dos ejemplos, y en el apartado 5 se presentan las conclusiones.59
60
==2. Fundamentos teóricos==61
62
=== 2.1. Formulación RTO con el método SIMP ===63
64
El objetivo principal de este trabajo es considerar la existencia de incertidumbre en la magnitud, la dirección y/o la posición de cargas concentradas en el problema de optimización de topología utilizando el método SIMP.65
66
El método SIMP fue presentado por Bendsøe [41] y posteriormente desarrollado por Zhou y Rozvany [42]. Actualmente es una de las técnicas más extendidas dentro de los métodos de optimización de topología basados en densidad. La idea básica del método consiste en determinar la mejor distribución de material. Para ello, se discretiza el dominio inicial mediante un conjunto de elementos finitos y se utiliza la densidad de cada elemento <math display="inline">({\rho }_{e})</math> como variable de diseño para indicar su estado, si es vacío <math display="inline">{\rho }_{e}</math>= 0 y si es material <math display="inline">{\rho }_{e}</math>=1. Para evitar densidades intermedias y alcanzar soluciones del tipo 0/1, el método SIMP utiliza una ley de interpolación con penalización:67
68
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 69
|-70
| 71
{| style="text-align: center; margin:auto;width: 100%;"72
|-73
| <math>{E}_{e}={E}_{min}+\left( E-{E}_{min}\right) {\rho }_{e}^{p}</math>74
|}75
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(1)76
|}77
78
donde <math display="inline">E</math> es el módulo de elasticidad para el material sólido, <math display="inline">{E}_{min}</math> es el módulo de elasticidad para el material vacío y el exponente <math display="inline">p</math> es un factor de penalización, siendo habituales valores <math display="inline">p\geq 1</math>. Para evitar inestabilidades numéricas, normalmente se considera que <math display="inline">{E}_{min}</math>= <math display="inline">\, E\times ({10}^{-3}\, a\, {10}^{-6})</math>.79
80
La relación entre la fuerza aplicada a una estructura y el desplazamiento que induce en la dirección de esta es el término conocido como la ''compliance'' <math display="inline">C</math>.81
82
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 83
|-84
| 85
{| style="text-align: center; margin:auto;width: 100%;"86
|-87
| <math>C={\boldsymbol{u}}^{\boldsymbol{T}}\boldsymbol{f}</math>88
|}89
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(2)90
|}91
92
donde <math display="inline">\boldsymbol{u}</math> es el vector de desplazamientos y <math display="inline">\boldsymbol{f}</math> el vector de cargas en los nodos. La aplicación del teorema de Clayperon proporciona una relación entre la ''compliance'' <math display="inline">\left( C\right)</math> y la energía elástica de deformación <math display="inline">U</math>93
94
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 95
|-96
| 97
{| style="text-align: center; margin:auto;width: 100%;"98
|-99
| <math>C=2U</math>100
|}101
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(3)102
|}103
104
Por lo que, el valor de la ''compliance'' se puede utilizar para determinar el nivel de energía de deformación contenido en la estructura. La ''compliance'' proporciona un valor matemático útil y preciso del grado de optimalidad de una estructura, con respecto a la rigidez del conjunto. De forma que minimizar la energía elástica de deformación total de la estructura, minimiza la ''compliance'' y de este modo se maximiza la rigidez estructural.105
106
En este trabajo, para poder obtener una estructura óptima robusta, se ha utilizado el valor esperado de la ''compliance'' [26, 27], aunque también se puede utilizar una suma ponderada del valor esperado y de la desviación estándar de la ''compliance''. De este modo, el problema RTO se puede formular como minimizar el valor esperado de la ''compliance'' <math display="inline">({\mu }_{C})</math> (Ec. 4.1) sujeto a una ecuación de equilibrio (Ec. 4.2), una restricción de desigualdad (Ec. 4.3) y una de doble desigualdad (Ec. 4.4):107
108
{| style="width: 100%;border-collapse: collapse;" 109
|-110
| Minimizar:111
| <math>{\mu }_{C}=E\left[ C\left( \boldsymbol{\rho },\boldsymbol{z}\right) \right] =\int_{}^{}C\left( \boldsymbol{\rho },\boldsymbol{z}\right) P\left( \boldsymbol{z}\right) d\boldsymbol{z}</math>112
| style="text-align: center;vertical-align: top;"|(4.1)113
| rowspan='4' style="text-align: right;"|(4)114
|-115
| sujeto a:116
| <math>\boldsymbol{K}\left( \boldsymbol{\rho }\mathit{\boldsymbol{,}}\boldsymbol{z}\right) \boldsymbol{\, u}\left( \boldsymbol{\rho }\mathit{\boldsymbol{,}}\boldsymbol{z}\right) =</math><math>\boldsymbol{f}(\boldsymbol{z})</math>117
| style="text-align: center;vertical-align: top;"|(4.2)118
|-119
| 120
| <math display="inline">V-{V}_{max}</math> ≤ 0121
| style="text-align: center;vertical-align: top;"|(4.3)122
|-123
| 124
| <math>\boldsymbol{0}\leq \boldsymbol{\rho }\leq \boldsymbol{1}</math>125
| style="text-align: center;vertical-align: top;"|(4.4)126
|}127
128
donde:129
130
{| style="width: 100%;border-collapse: collapse;" 131
|-132
| style="vertical-align: top;"| <math>\boldsymbol{\rho }</math>133
| style="vertical-align: top;"|es el vector de las variables de diseño (densidades de los elementos),134
|-135
| style="vertical-align: top;"| <math>\boldsymbol{z}</math>136
| style="vertical-align: top;"|es el vector de las variables aleatorias,137
|-138
| style="vertical-align: top;"| <math>C\left( \boldsymbol{\rho },\boldsymbol{z}\right)</math> 139
| style="vertical-align: top;"|es la ''compliance'' de la estructura,140
|-141
| style="vertical-align: top;"| <math>P\left( \boldsymbol{z}\right)</math> 142
| style="vertical-align: top;"|es la función de densidad de probabilidad que caracteriza la variable aleatoria <math display="inline">\boldsymbol{z}</math>,143
|-144
| style="vertical-align: top;"| <math>\boldsymbol{K}\left( \boldsymbol{\rho }\mathit{\boldsymbol{,}}\boldsymbol{z}\right)</math> 145
| style="vertical-align: top;"|es la matriz de rigidez global,146
|-147
| style="vertical-align: top;"| <math>\boldsymbol{u}\left( \boldsymbol{\rho }\mathit{\boldsymbol{,}}\boldsymbol{z}\right)</math> 148
| style="vertical-align: top;"|es el campo de desplazamientos para un vector de fuerzas <math display="inline">\boldsymbol{f}\left( \boldsymbol{z}\right) ,</math>149
|-150
| style="vertical-align: top;"| <math>V</math>151
| style="vertical-align: top;"|es el volumen, 152
|-153
| style="vertical-align: top;"| <math>{V}_{max}</math>154
| style="vertical-align: top;"|es el volumen máximo admisible.155
|}156
157
El problema de optimización (Ec. 4) se puede resolver mediante OC (''Optimality Criteria''), SLP (''Sequencial Linear Progamming''), MMA (''Method of Moving Asymptotes''), etc. [43]. Por simplicidad, el problema se ha resuelto utilizando el Método OC, donde la densidad de cada elemento <math display="inline">{\rho }_{e}</math> se actualiza mediante un esquema heurístico [44]:158
159
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 160
|-161
| 162
{| style="text-align: center; margin:auto;width: 100%;"163
|-164
| <math>{\rho }_{e}^{nuevo}=\left\{ \begin{matrix}\mathrm{max}\,\left( 0,{\rho }_{e}-m\right) \\\mathrm{min}\,\left( 1,{\rho }_{e}+m\right) \\{\rho }_{e}{B}_{e}^{\eta }\end{matrix}\right.</math> <math>\begin{matrix}si\quad {\rho }_{e}{B}_{e}^{\eta }\, \leq \mathrm{max}\,\left( 0,{\rho }_{e}-m\right) \, \, \\si\quad {\rho }_{e}{B}_{e}^{\eta }\, \geq \, \mathrm{min}\,\left( 1,{\rho }_{e}+m\right) \\otro\, caso\end{matrix}</math>165
|}166
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(5)167
|}168
169
donde <math display="inline">m</math> es un límite de movimiento, <math display="inline">\eta =</math><math>0,5</math> es un coeficiente de amortiguamiento y <math display="inline">{B}_{e}</math> es la condición de optimalidad:170
171
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 172
|-173
| 174
{| style="text-align: center; margin:auto;width: 100%;"175
|-176
| <math>{B}_{e}=\frac{-\frac{\partial {\mu }_{C}}{\partial {\rho }_{e}}}{\lambda \frac{\partial V}{\partial {\rho }_{e}}}</math>177
|}178
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(6)179
|}180
181
siendo <math display="inline">\lambda</math> el multiplicador de Lagrange que garantiza que se satisface la restricción de volumen, cuyo valor se determina usando un algoritmo de bisección.182
183
La sensibilidad de la función objetivo <math display="inline">({\mu }_{C})</math> y del volumen ( <math display="inline">V</math>) con respecto a las variables de diseño se han determinado como:184
185
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 186
|-187
| 188
{| style="text-align: center; margin:auto;width: 100%;"189
|-190
| <math>\frac{\partial {\mu }_{C}}{\partial {\rho }_{e}}=\frac{\partial E\left[ C\left( \boldsymbol{\rho },\boldsymbol{z}\right) \right] }{\partial {\rho }_{e}}=</math><math>E\left[ \frac{\partial C\left( \boldsymbol{\rho },\boldsymbol{z}\right) }{\partial {\rho }_{e}}\right]</math> 191
|}192
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(7)193
|-194
| 195
{| style="text-align: center; margin:auto;width: 100%;"196
|-197
| <math>\frac{\partial V}{\partial {\rho }_{e}}={v}_{e}</math>198
|}199
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(8)200
|}201
202
De este modo, la sensibilidad del valor esperado de la función objetivo es igual a la media de la sensibilidad de la ''compliance'', y la sensibilidad del volumen es igual al volumen del elemento <math display="inline">\left( {v}_{e}\right)</math> .203
204
Se conoce que los diseños obtenidos usando el método SIMP, pueden presentar soluciones dependientes de la malla (una malla fina puede producir diseños nuevos, en lugar de optimizar los obtenidos con una malla más gruesa), así como una alternancia de elementos sólidos y vacíos ordenados en forma de tablero de ajedrez [45] si no se utiliza algún esquema de regularización.205
206
En general, los métodos para la restricción de la densidad en problemas de optimización de topología se dividen en tres grupos [46-52]: (1) métodos para la independencia de la malla, que se subdividen en filtros de sensibilidades y de densidades; (2) métodos de restricción como el control de perímetro, el control por gradiente local, y la penalización regularizada; y (3) otros métodos como la parametrización de ondas, el campo fase, y las curvas de nivel. Estos esquemas tienen ventajas e inconvenientes y probablemente todavía el esquema ideal no ha sido desarrollado.207
208
Los diseños obtenidos introduciendo un filtro de densidad pueden contener densidades intermedias, lo que dificulta la interpretación del diseño. Por ello, las densidades intermedias, también, se pueden tratar mediante técnicas de proyección a 0/1 [48-51]. Para evitar soluciones dependientes de la malla y patrones del tipo tablero de ajedrez, se ha utilizado en este trabajo el filtro propuesto por Sigmund [49]. Este filtro modifica las sensibilidades de los elementos a partir del promedio ponderado de las sensibilidades de los elementos situados en el entorno cercano, Ec. 9.209
210
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 211
|-212
| 213
{| style="text-align: center; margin:auto;width: 100%;"214
|-215
| <math>\frac{\hat{{\partial \mu }_{C}}}{\partial {\rho }_{e}}=\frac{1}{{\rho }_{e}\sum _{i=1}^{N}\hat{{H}_{i}}}\sum _{i=1}^{N}\hat{{H}_{i}}{\rho }_{i}\frac{\partial {\mu }_{C}}{\partial {\rho }_{i}}</math>216
|}217
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(9)218
|}219
220
donde el operador convolución (factor de ponderación) <math display="inline">\hat{{H}_{i}}</math> se define como:221
222
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 223
|-224
| 225
{| style="text-align: center; margin:auto;width: 100%;"226
|-227
| <math>\hat{{H}_{i}}={r}_{min}-d\left( e,i\right)</math>228
229
<math>\left\{ i\in N\right\} ,k=1,\cdots ,N</math>230
|}231
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(10)232
|}233
234
donde el operador <math display="inline">d\left( e,i\right)</math> se define como la distancia entre el centro del elemento <math display="inline">k</math> y el centro del elemento <math display="inline">i</math>. El operador convolución <math display="inline">\hat{{H}_{i}}</math> es cero fuera del área del filtro. El operador convolución para el elemento <math display="inline">i</math> disminuye linealmente con la distancia al elemento <math display="inline">e</math>, lo que significa que en lugar de utilizar las sensibilidades reales <math display="inline">\frac{\partial {\mu }_{C}}{\partial {\rho }_{i}}</math> se utilizan las sensibilidades filtradas <math display="inline">\frac{\hat{{\partial \mu }_{C}}}{\partial {\rho }_{e}}</math>. Merece la pena resaltar que: (1) cuando el <math display="inline">{r}_{min}</math> (radio de filtro con centro en el elemento <math display="inline">e</math>) tiende a cero la sensibilidad filtrada converge a la real; (2) que todas las sensibilidades son iguales (caso de una distribución uniforme de material) cuando el valor de <math display="inline">{r}_{min}</math> tiende a infinito; y (3) que este filtro es puramente heurístico pero produce resultados muy similares a los obtenidos usando control por gradiente local, su implementación es sencilla y solo conlleva un pequeño gasto computacional extra.235
236
Apuntar que en este trabajo para la consecución de la optimización de topología se ha utilizado una versión modificada de la implementación del método SIMP “99 líneas en Matlab” realizada por Sigmund [44].237
238
===2.2. Incertidumbres===239
240
En este trabajo se consideran las incertidumbres de cargas concentradas. Las cargas se definen por medio de cuatro variables: magnitud (''F''), dirección (''θ''), posición respecto al eje ''x'' (''S''<sub>x</sub>) y posición con respecto al eje ''y'' (''S''<sub>y</sub>), cada una de las cuales pueden presentar o no incertidumbre. Cualquiera de las variables anteriores que presente incertidumbre se caracteriza mediante una pdf.241
242
La totalidad de los trabajos sobre RTO de estructuras continuas encontrados en la bibliografía [18, 20, 23-27] utilizan métodos que requieren adoptar la hipótesis de que las variables aleatorias consideradas son independientes. Aunque esta hipótesis no se cumple en una gran cantidad de casos, los problemas resueltos hasta el momento se han realizado con variables aleatorias independientes. Esta decisión ha permitido, entre otras cosas, la comparación de sus resultados con los existentes en la bibliografía. En trabajos posteriores se analizará el efecto que la dependencia entre las variables aleatorias puede introducir en los parámetros del modelo Kriging.243
244
===2.3. Método de simulación de Monte Carlo===245
246
El método de simulación Monte Carlo es un método numérico no determinístico, usado para aproximar expresiones matemáticas que resultan costosas si se evalúan con exactitud. El método se llamó así en referencia al Casino de Monte Carlo, por ser la ruleta un generador simple de números aleatorios.247
248
El uso del método de Monte Carlo, como herramienta de investigación, proviene del trabajo realizado en el desarrollo de la bomba atómica durante la Segunda Guerra Mundial en el Laboratorio Nacional de Los Álamos en EE. UU. Su origen se establece en 1949, cuando por primera vez aparece en el título de un artículo de Metropolis y Ulam [53], en el cual se denomina como un método “estadístico numérico para aproximar expresiones matemáticas complejas y costosas de evaluar con exactitud”. Hoy en día, el método de Monte Carlo se utiliza en el desarrollo de reactores nucleares, cromodinámica cuántica, radioterapia, comportamiento de radiación en la atmósfera terrestre, flujos de tráfico, evolución estelar, cálculos y predicciones económicas, búsqueda de petróleo, etc.249
250
La idea del método de Monte Carlo consiste en crear un modelo matemático del sistema que se quiere analizar, identificando aquellas variables que tienen un carácter aleatorio. Una vez determinadas estas variables, se lleva a cabo un diseño de experimentos consistente en generar una muestra para dichas variables aleatorias, y analizar el modelo de simulación para los valores generados.251
252
En este trabajo se utiliza el método de Monte Carlo para evaluar el valor esperado de la ''compliance'' (Ec. 4.1). Para ello, el primer paso consiste en identificar las variables aleatorias. Como se ha mencionado, la incertidumbre puede darse en la magnitud, la dirección, la posición respecto al eje ''x'', y la posición respecto al eje ''y'' de la carga. Estas variables se consideran estadísticamente independientes. A partir de ellas se genera aleatoriamente una muestra de diseño del espacio multidimensional <math display="inline">\boldsymbol{Z}=</math><math>{\left( {\boldsymbol{z}}_{\mathit{\boldsymbol{1}}}{\boldsymbol{,\, z}}_{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{,\ldots ,\, }}{{\boldsymbol{\, z}}_{\mathit{\boldsymbol{i}}}\mathit{\boldsymbol{,\ldots ,\, }}\boldsymbol{z}}_{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{MC}}}}\right) }^{T}</math>, uniformemente distribuida según la pdf con <math display="inline">{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{i}}}\, \in \, {R}^{s}</math>, donde la dimensión ''s ''es igual al número de variables aleatorias consideradas. Así pues, si se considera que hay una única carga con magnitud y dirección aleatoria, la muestra será <math display="inline">{\boldsymbol{z}}_{\mathit{\boldsymbol{i}}}=</math><math>\left( {F}_{i},\, {\theta }_{i}\right)</math> ''' '''para ''i ''= 1,…, ''N<sub>MC</sub>'', siendo ''N<sub>MC</sub>'' el tamaño de la muestra. A partir de la evaluación de la muestra se estima el valor esperado de la ''compliance'', transformando la integral definida en la Ec. 4.1 en la expresión:253
254
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 255
|-256
| 257
{| style="text-align: center; margin:auto;width: 100%;"258
|-259
| <math>{\mu }_{C}\cong \frac{1}{{N}_{MC}}\sum _{i=1}^{{N}_{MC}}{C}_{i}\left( \boldsymbol{\rho },{\boldsymbol{z}}_{\mathit{\boldsymbol{i}}}\right)</math> 260
|}261
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(11)262
|}263
264
Esta ecuación permite transformar el problema robusto en otro sujeto a múltiples estados de carga deterministas, siendo necesario evaluar repetidamente el modelo de simulación de la estructura para cada una de las ''N<sub>MC</sub>'' muestras.265
266
Según Rubinstein y Kroese [54] el método de simulación de Monte Carlo es la metodología más precisa para la estimación de momentos estadísticos de alto orden dimensional, lo que permite: (1) resolver problemas cuya solución es difícil o imposible de resolver, (2) estudiar la interacción entre las diferentes variables aleatorias del problema, y (3) evaluar los momentos estadísticos sin necesidad de modificar el modelo de simulación. Sin embargo, su mayor inconveniente es el elevado coste computacional, el cual es proporcional al tamaño de la muestra. La precisión de la estimación depende del tamaño de la muestra; en concreto, según el teorema del límite central, para <math display="inline">{N}_{MC}</math> grandes, el promedio de la muestra tiene distribución aproximadamente normal con media igual a la media de la variable original y desviación estándar la de la variable dividida por la raíz cuadrada de <math display="inline">{N}_{MC}</math>. Así, la desviación estándar de la variable se puede estimar por medio de la desviación estándar de la muestra (Ec. 12):267
268
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 269
|-270
| 271
{| style="text-align: center; margin:auto;width: 100%;"272
|-273
| <math>\sigma =\sqrt{\frac{1}{{N}_{MC}}\left( \sum _{i=1}^{n}{{C}_{i}\left( \boldsymbol{\rho },{\boldsymbol{z}}_{\mathit{\boldsymbol{i}}}\right) }^{2}-{{\mu }_{C}}^{2}\right) }</math>274
|}275
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(12)276
|}277
278
Actualmente se han desarrollado alternativas para mejorar la eficiencia numérica, tales como el muestreo por importancia, el muestreo estratificado, etc. [55].279
280
En este trabajo se ha utilizado el método de simulación de Monte Carlo, ya que permite resolver el problema como una suma ponderada de múltiples estados de carga deterministas, donde la precisión del método solo depende del tamaño de la muestra y no de la dimensionalidad del problema, lo que resulta muy útil en problemas con una alta dimensionalidad.281
282
===2.4. Metamodelos===283
284
Los metamodelos (modelos de modelos) se utilizan para sustituir a un modelo de simulación cuya evaluación tiene un alto coste computacional. Estos metamodelos pueden utilizarse, entre otras aplicaciones, para acelerar el proceso de optimización y para realizar visualizaciones y exploraciones del espacio de diseño óptimo. En la literatura especializada existe una gran variedad de metamodelos, tales como las superficies de respuesta [32], los modelos Kriging [56], las funciones de base radial [57], o las redes neuronales [58]. De entre las diferentes técnicas de metamodelos, los modelos Kriging han alcanzado una gran popularidad en los últimos años, debido a: 1) su gran flexibilidad para aproximar respuestas con alto grado de no linealidad [59]; 2) ser especialmente adecuados para aproximar simulaciones con un elevado coste computacional; 3) que en su formulación se imponen las condiciones para que el modelo sea un estimador lineal insesgado de los valores interpolados, y 4) proporcionan información estadística del error cometido en la predicción [60, 61].285
286
====2.4.1. Modelos Kriging====287
288
Los modelos Kriging surgieron en los años 50 para estimar las reservas en las minas de oro de Sudáfrica por el ingeniero minero Daniel Krige y el estadista Herbert Shichel. Posteriormente, Matheron adoptó el trabajo original realizado por Daniel Krige y formuló matemáticamente la mayor parte de los conceptos teóricos [61].289
290
Los modelos Kriging están basados en la interpolación estocástica, utilizando un modelo de regresión que considera la correlación existente entre éste y las observaciones para construir un modelo aproximado <math display="inline">\hat{y}(\mathit{\boldsymbol{x}})</math>, computacionalmente menos costoso que el modelo de simulación <math display="inline">y\left( \boldsymbol{x}\right) \,</math> [35].291
292
A continuación se presentan los aspectos fundamentales de los modelos Kriging para una sola respuesta ''y''('''x'''). La generalización para múltiples respuestas ( <math display="inline">{y}_{l}\left( \boldsymbol{x}\right) ,\, l=</math><math>1,\ldots q)</math>, y otros detalles de la formulación, puede verse en [56].293
294
Dados un conjunto de ''m'' puntos de diseño <math display="inline">S={\left[ {s}_{1}\ldots {s}_{m}\right] }^{T}\, con\, \, {s}_{i}\in {R}^{n}</math> y sus respuestas <math display="inline">Y=</math><math>{\left[ {y}_{1}\ldots {y}_{m}\right] }^{T}</math>, el modelo Kriging <math display="inline">\hat{y}\left( \boldsymbol{x}\right)</math> expresa la respuesta ''y''('''x''') para un vector de entrada n-dimensional <math display="inline">\boldsymbol{x}\in {R}^{n}\,</math> como una realización de un modelo de regresión <math display="inline">F</math> y una función aleatoria (proceso estocástico):295
296
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 297
|-298
| 299
{| style="text-align: center; margin:auto;width: 100%;"300
|-301
| <math>\hat{y}\left( \boldsymbol{x}\right) =F\left( \boldsymbol{\beta ,x}\right) +</math><math>z\left( \boldsymbol{x}\right)</math> 302
|}303
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(13)304
|}305
306
Como modelo de regresión se utiliza una combinación lineal de ''p'' funciones elegidas previamente <math display="inline">{f}_{j}</math>:307
308
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 309
|-310
| 311
{| style="text-align: center; margin:auto;width: 100%;"312
|-313
| <math>F\left( \boldsymbol{\beta ,x}\right) ={\beta }_{1}{f}_{1}\left( \boldsymbol{x}\right) \boldsymbol{+\ldots +\, }{\beta }_{p}{f}_{p}\left( \boldsymbol{x}\right) =</math><math>\boldsymbol{f}{\left( \boldsymbol{x}\right) }^{T}\boldsymbol{\beta }</math>314
|}315
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(14)316
|}317
318
siendo <math display="inline">{\boldsymbol{f}\left( \boldsymbol{x}\right) }^{T}=</math><math>\left[ {f}_{1}\left( \boldsymbol{x}\right) \mathit{\boldsymbol{\ldots }}{f}_{p}\left( \boldsymbol{x}\right) \right]</math> ''' '''las funciones base y <math display="inline">\, \boldsymbol{\beta }\mathit{\boldsymbol{\, }}</math> el vector de parámetros de regresión.319
320
El proceso aleatorio ''z'' se supone que tiene media cero y covarianza:321
322
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 323
|-324
| 325
{| style="text-align: center; margin:auto;width: 100%;"326
|-327
| <math>E\left( z\left( \boldsymbol{w}\right) ,z\left( \boldsymbol{x}\right) \right) ={\sigma }^{2}R\left( \theta ,\boldsymbol{w},\boldsymbol{x}\right)</math> 328
|}329
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(15)330
|}331
332
entre z('''w''') y z('''x'''), donde <math display="inline">{\sigma }^{2}</math> es la varianza del proceso para la respuesta y <math display="inline">R\left( \theta ,\boldsymbol{w},\boldsymbol{x}\right)</math> el modelo de correlación con parámetros <math display="inline">\theta</math> .333
334
El valor verdadero de la función para un punto '''x''' puede ponerse en la forma:335
336
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 337
|-338
| 339
{| style="text-align: center; margin:auto;width: 100%;"340
|-341
| <math>y\left( \boldsymbol{x}\right) =F\left( \boldsymbol{\beta ,x}\right) +</math><math>\alpha \left( \boldsymbol{\beta ,x}\right)</math> 342
|}343
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(16)344
|}345
346
donde <math display="inline">\alpha \left( \boldsymbol{\beta ,x}\right)</math> ''' '''es''' '''el error de la aproximación. La hipótesis es que con una elección adecuada de <math display="inline">\boldsymbol{\beta }</math>''' '''el error es del tipo “ruido blanco”''' '''en la región de interés de la aproximación.347
348
Para el conjunto ''S'' de puntos de diseño se tiene el vector '''F''' de respuestas:349
350
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 351
|-352
| 353
{| style="text-align: center; margin:auto;width: 100%;"354
|-355
| <math>\boldsymbol{F}={\left[ f({s}_{1})\ldots f({s}_{m})\right] }^{T}</math>356
|}357
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(17)358
|}359
360
con ''f''(x) definido en la Ec. 14. Se define '''R''' como la matriz de correlaciones del proceso estocástico entre los ''z'' en los puntos de diseño:361
362
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 363
|-364
| 365
{| style="text-align: center; margin:auto;width: 100%;"366
|-367
|<math>{R}_{ij}=R\left( \theta ,{s}_{i},{s}_{j}\right) , \quad i,j = 1,..., m</math>368
|}369
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(18)370
|}371
372
Para un punto ''x'' no introducido, sea:373
374
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 375
|-376
| 377
{| style="text-align: center; margin:auto;width: 100%;"378
|-379
| <math>\boldsymbol{r}\left( x\right) ={\left[ R\left( \theta ,{s}_{1},x\right) \ldots \, R\left( \theta ,{s}_{m},x\right) \right] }^{T}</math>380
|}381
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(19)382
|}383
384
el vector de correlaciones entre las z en los puntos de diseño y el punto ''x''.385
386
Se considera el estimador lineal:387
388
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 389
|-390
| 391
{| style="text-align: center; margin:auto;width: 100%;"392
|-393
| <math>\hat{y}\left( \boldsymbol{x}\right) =\, {\boldsymbol{c}}^{T}\boldsymbol{Y}</math>394
|}395
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(20)396
|}397
398
con''' c = c(x)''' <math display="inline">\in {R}^{m}</math>. El error es:399
400
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 401
|-402
| 403
{| style="text-align: center; margin:auto;width: 100%;"404
|-405
| <math>\hat{y}\left( \boldsymbol{x}\right) -y\left( \boldsymbol{x}\right) =\, {\boldsymbol{c}}^{T}\boldsymbol{Y-}y\left( \boldsymbol{x}\right) =</math><math>{\boldsymbol{c}}^{T}\left( \boldsymbol{F\beta +Z}\right) -\left( {\mathit{\boldsymbol{f}}\left( \boldsymbol{x}\right) }^{T}\boldsymbol{\beta +z}\right) =</math><math>{\boldsymbol{c}}^{T}\boldsymbol{Z-z+}{\left( {\boldsymbol{F}}^{T}\boldsymbol{c}\mathit{\boldsymbol{-f(}}\boldsymbol{x}\mathit{\boldsymbol{)}}\right) }^{T}\boldsymbol{\beta }</math>406
|}407
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(21)408
|}409
410
donde <math display="inline">\boldsymbol{Z}={\left[ {z}_{1}\ldots \, {z}_{m}\right] }^{T}</math> son los errores en los puntos de diseño. Para mantener el estimador insesgado se exige la condición:411
412
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 413
|-414
| 415
{| style="text-align: center; margin:auto;width: 100%;"416
|-417
| <math>{\boldsymbol{F}}^{\mathit{\boldsymbol{T}}}\boldsymbol{c}\left( \boldsymbol{x}\right) \mathit{\boldsymbol{=f(x)}}</math>418
|}419
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(22)420
|}421
422
Con esta condición el error cuadrático medio del estimador de la Ec. 20 es:423
424
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 425
|-426
| 427
{| style="text-align: center; margin:auto;width: 100%;"428
|-429
| <math>\varphi \left( \mathit{\boldsymbol{x}}\right) =E\left[ {\left( \hat{y}\left( \boldsymbol{x}\right) -y\left( \boldsymbol{x}\right) \right) }^{2}\right] =</math><math>E\left[ {\left( {\boldsymbol{c}}^{T}\boldsymbol{Z-z}\right) }^{2}\right] =E\left[ {\boldsymbol{z}}^{2}+\right. </math><math>\left. {\boldsymbol{c}}^{T}\mathit{\boldsymbol{Z}}{\boldsymbol{Z}}^{T}\boldsymbol{c}-\right. </math><math>\left. 2{\boldsymbol{c}}^{T}\boldsymbol{Zz}\right] ={\sigma }^{2}\left( 1+{\boldsymbol{c}}^{T}\boldsymbol{Rc-2}{\boldsymbol{c}}^{T}\boldsymbol{r}\right)</math> 430
|}431
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(23)432
|}433
434
La función Lagrangiana para el problema de minimizar el error <math display="inline">\varphi</math> con respecto a '''c''', sujeto a las restricciones de la Ec. 21 es:435
436
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 437
|-438
| 439
{| style="text-align: center; margin:auto;width: 100%;"440
|-441
| <math>L\left( \boldsymbol{c}\mathit{\boldsymbol{,}}\boldsymbol{\lambda }\right) =</math><math>{\sigma }^{2}\left( 1+{\boldsymbol{c}}^{T}\boldsymbol{Rc-2}{\boldsymbol{c}}^{T}\boldsymbol{r}\right) -</math><math>{\boldsymbol{\lambda }}^{T}\left( {\boldsymbol{F}}^{T}\boldsymbol{c}\mathit{\boldsymbol{-f}}\right)</math> 442
|}443
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(24)444
|}445
446
El gradiente de la Ec. 24 con respecto a '''c''' es:447
448
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 449
|-450
| 451
{| style="text-align: center; margin:auto;width: 100%;"452
|-453
| <math>{L}^{'}\left( \boldsymbol{c}\mathit{\boldsymbol{,}}\boldsymbol{\lambda }\right) \mathit{\boldsymbol{=2}}{\sigma }^{2}\left( \boldsymbol{Rc-r}\right) -</math><math>\boldsymbol{F\lambda }</math>454
|}455
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(25)456
|}457
458
y a partir de las condiciones de primer orden necesarias para la optimalidad, se obtiene el sistema de ecuaciones:459
460
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 461
|-462
| 463
{| style="text-align: center; margin:auto;width: 100%;"464
|-465
| <math>\left[ \begin{matrix}\boldsymbol{R}&\boldsymbol{F}\\{\boldsymbol{F}}^{T}&\mathit{\boldsymbol{0}}\end{matrix}\right] \left[ \begin{matrix}\boldsymbol{c}\\\tilde{\boldsymbol{\lambda }}\end{matrix}\right] =</math><math>\left[ \begin{matrix}\boldsymbol{r}\\\boldsymbol{f}\end{matrix}\right]</math> 466
|}467
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(26)468
|}469
470
donde <math display="inline">\tilde{\boldsymbol{\lambda }}=-\frac{\boldsymbol{\lambda }}{2{\sigma }^{2}}</math>.471
472
La solución del sistema de ecuaciones de la Ec. 26 es:473
474
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 475
|-476
| 477
{| style="text-align: center; margin:auto;width: 100%;"478
|-479
| <math>\begin{matrix}\tilde{\boldsymbol{\lambda }}={\left( {\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{F}\right) }^{-1}\left( {\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{r-f}\right) \\\boldsymbol{c}={\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\left( \boldsymbol{r-F}\tilde{\boldsymbol{\lambda }}\right) \end{matrix}</math>480
|}481
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(27)482
|}483
484
Sustituyendo la Ec. 27 en la Ec. 20, y teniendo en cuenta que la matriz '''R''' es simétrica, y por lo tanto también lo es su inversa '''R'''<sup>-1</sup>, se obtiene:485
486
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 487
|-488
| 489
{| style="text-align: center; margin:auto;width: 100%;"490
|-491
| <math>\begin{matrix}\tilde{y}(\boldsymbol{x})={\left( \boldsymbol{r-F}\tilde{\boldsymbol{\lambda }}\right) }^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{Y}\\={\boldsymbol{r}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{Y-}{\left( {\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{r-f}\right) }^{T}{\left( {\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{F}\right) }^{-1}{\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{Y}\end{matrix}</math>492
|}493
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(28)494
|}495
496
La solución de mínimos cuadrados para el problema de regresión <math display="inline">\boldsymbol{F\beta \cong Y}</math> es:497
498
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 499
|-500
| 501
{| style="text-align: center; margin:auto;width: 100%;"502
|-503
| <math>{\boldsymbol{\beta }}^{\mathit{\boldsymbol{\ast }}}={\left( {\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{F}\right) }^{-1}{\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{Y}</math>504
|}505
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(29)506
|}507
508
Al sustituir la Ec. 29 en la Ec. 28 se obtiene la expresión final del estimador:509
510
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 511
|-512
| 513
{| style="text-align: center; margin:auto;width: 100%;"514
|-515
| <math>\begin{matrix}\tilde{y}(\boldsymbol{x})={\boldsymbol{r}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{Y-}{\left( {\boldsymbol{F}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\boldsymbol{r-f}\right) }^{T}{\boldsymbol{\beta }}^{\mathit{\boldsymbol{\ast }}}\\={\boldsymbol{f}}^{T}{\boldsymbol{\beta }}^{\mathit{\boldsymbol{\ast }}}\mathit{\boldsymbol{+}}{\boldsymbol{r}}^{T}{\boldsymbol{R}}^{\mathit{\boldsymbol{-1}}}\left( \boldsymbol{Y}\mathit{\boldsymbol{-}}\boldsymbol{F}{\boldsymbol{\beta }}^{\mathit{\boldsymbol{\ast }}}\right) \end{matrix}</math>516
|}517
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(30)518
|}519
520
Si se sustituye uno de los puntos de diseño (''s<sub>i</sub>'') en la Ec. 30 se obtiene que <math display="inline">\tilde{y}\left( {s}_{i}\right) =</math><math>y\left( {s}_{i}\right)</math> , lo que demuestra que los modelos Kriging pasan por todos los puntos de diseño.521
522
'''Modelos de regresión'''. Los modelos de regresión considerados son polinomios de orden 0, 1 y 2 (tabla 1).523
524
525
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;"font-size: 75%;">526
<span style="text-align: center; font-size: 75%;">'''Tabla 1.''' Modelos de regresión con polinomios de orden 0, 1 y 2.</span></div>527
528
{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 529
|-530
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">Orden polinomio</span>531
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">Número de funciones</span>532
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">Funciones base</span>533
|-534
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">Constante (orden 0)</span>535
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">''p'' = 1</span>536
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">''f''<sub>1</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = 1</span>537
|-538
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">Lineal (orden 1)</span>539
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">''p'' = ''n''+1</span>540
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">''f''<sub>1</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = 1, ''f''<sub>2</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = ''x''<sub>1</sub></span><span style="text-align: center; font-size: 85%;">,…, ''f<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;"><sub>+1</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = ''x<sub>n</sub>''</span>541
|-542
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">Cuadrática (orden 2)</span>543
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">''p'' = 0,5(''n''+1)(''n''+2)</span>544
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;vertical-align: top;"|<span style="text-align: center; font-size: 85%;">''f''<sub>1</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = 1, ''f''<sub>2</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = x<sub>1</sub></span><span style="text-align: center; font-size: 85%;">,…, ''f<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;"><sub>+1</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = ''x<sub>n</sub>''</span>545
546
<span style="text-align: center; font-size: 85%;">''f<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;"><sub>+2</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = ''x''<sub>1</sub></span><span style="text-align: center; font-size: 85%;"><sup>2</sup></span><span style="text-align: center; font-size: 85%;">,…, ''f''<sub>2</sub></span><span style="text-align: center; font-size: 85%;">''<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;"><sub>+1</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = ''x''<sub>1</sub></span><span style="text-align: center; font-size: 85%;">''x''<sub>2</sub></span>547
548
<span style="text-align: center; font-size: 85%;">''f''<sub>2</sub></span><span style="text-align: center; font-size: 85%;">''<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;"><sub>+2</sub></span><span style="text-align: center; font-size: 85%;">('''x''') = ''x''<sub>2</sub></span><span style="text-align: center; font-size: 85%;"><sup>2</sup></span><span style="text-align: center; font-size: 85%;">,…, ''f''<sub>3</sub></span><span style="text-align: center; font-size: 85%;">''<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;">('''x''') = ''x''<sub>2</sub></span><span style="text-align: center; font-size: 85%;">''x<sub>n</sub>''</span>549
550
<span style="text-align: center; font-size: 85%;">''f<sub>p</sub>''</span><span style="text-align: center; font-size: 85%;">('''x''') = ''x<sub>n</sub>''</span><span style="text-align: center; font-size: 85%;"><sup>2</sup></span>551
|}<br />552
553
554
'''Modelo de correlación'''. El modelo de correlación utilizado es el exponencial generalizado:555
556
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 557
|-558
| 559
{| style="text-align: center; margin:auto;width: 100%;"560
|-561
| <math>R\left( \theta ,\boldsymbol{w},\boldsymbol{x}\right) =\mathrm{exp}\,\left( -\sum _{l=1}^{q}{\theta }_{l}{\left| {\boldsymbol{x}}_{il}-{\boldsymbol{x}}_{jl}\right| }^{{p}_{l}}\right) ,\quad 0\leq {p}_{l}\leq 2</math>562
|}563
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(31)564
|}565
566
Este modelo tiene la propiedad que si <math display="inline">{\boldsymbol{x}}_{i}=</math><math>{\boldsymbol{x}}_{\boldsymbol{j}}</math> la correlación es 1, y si | <math display="inline">{\boldsymbol{x}}_{i}-</math><math>{\boldsymbol{x}}_{j}</math>| tiende a ∞, la correlación tiende a cero. El parámetro ''θ<sub>l</sub>'' determina la rapidez con la que la correlación “cae” a medida que se mueve en la dirección coordenada ''l'', y el exponente <math display="inline">{p}_{l}\, determina\, la\, suavidad\, de\, la\, funci\acute{o}n\, en\, la\, direcci\acute{o}n\, coordenada\, l.</math>567
568
Las funciones de correlación y de regresión más adecuadas dependen de las características de la función a aproximar. En este trabajo se han considerado los modelos de regresión de orden cero, uno y dos de la Tabla 1, en combinación con el modelo de correlación exponencial generalizado. En general, la elección de las funciones de correlación y regresión se puede realizar “a priori”, en base al conocimiento que se tenga de las funciones a aproximar, o en base al error cuadrático medio (''Root Mean Square Error'', RMSE) y al coeficiente de determinación (R<sup>2</sup>). Estos dos indicadores se obtienen utilizando las respuestas del modelo de simulación (y('''x''')) para una muestra de validación ('''x'''<sub>v</sub>) de tamaño ''m<sub>v</sub>'':569
570
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 571
|-572
| 573
{| style="text-align: center; margin:auto;width: 100%;"574
|-575
| <math>RMSE=\sqrt{\frac{\sum _{i=1}^{{m}_{v}}{\left( {y}_{i}\left( {\boldsymbol{x}}_{v}\right) -{\hat{y}}_{i}\left( {\boldsymbol{x}}_{v}\right) \right) }^{2}}{{m}_{v}}}</math>576
|}577
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(32)578
|-579
| 580
{| style="text-align: center; margin:auto;width: 100%;"581
|-582
| <math>{R}^{2}=1-\frac{\sum _{i=1}^{{m}_{v}}{\left( {y}_{i}\left( {\boldsymbol{x}}_{v}\right) -{\hat{y}}_{i}\left( {\boldsymbol{x}}_{v}\right) \right) }^{2}}{\sum _{i=1}^{{m}_{v}}{\left( {y}_{i}\left( {\boldsymbol{x}}_{v}\right) -{\overline{y}}_{i}\left( {\boldsymbol{x}}_{v}\right) \right) }^{2}}</math>583
|}584
| style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"|(33)585
|}586
587
donde <math display="inline">{\overline{y}}_{i}\left( {\boldsymbol{x}}_{v}\right)</math> es la media de los valores obtenidos, con el modelo de simulación, para la muestra de validación.588
589
En este trabajo, el modelo Kriging se ha utilizado para estimar el valor esperado de la ''compliance'' (Ec. 11) a partir de un diseño de experimentos con un tamaño de población mucho menor que el de la muestra original (''N<sub>K </sub>''<sub><< </sub>''N<sub>MC</sub>''). Para ello, primero se genera una población mediante un hipercubo latino <math display="inline">\boldsymbol{S}=</math><math>{\left( {\boldsymbol{s}}_{\mathit{\boldsymbol{1\, }}}{\boldsymbol{s}}_{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{\ldots \, }}{\boldsymbol{s}}_{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{k}}}}\right) }^{T}</math>, donde la dimensión ''s ''es igual al número de variables aleatorias consideradas''.'' A continuación, se analiza el modelo de simulación, obteniendo la ''compliance'' para cada uno de los puntos <math display="inline">\boldsymbol{C}=</math><math>{\left( {C}_{1}\, {C}_{2}\ldots \, {C}_{{N}_{K}}\right) }^{T}</math>. A partir del espacio de diseño y la respuesta estructural, se ajustan los parámetros del modelo Kriging <math display="inline">\hat{y}(\boldsymbol{S,C})</math> de modo que aproxime de la mejor forma posible el modelo de simulación. Una vez se ha generado el modelo Kriging, se predice la ''compliance'' en los puntos de muestreo del método de Monte Carlo<br/> <math display="inline">\hat{\boldsymbol{C}}=</math><math>{\left( {\hat{C}}_{1}\, {\hat{C}}_{2}\ldots \, {\hat{C}}_{{N}_{MC}}\right) }^{T}</math> y a partir de estos, se evalúa el valor esperado de la ''compliance'' <math display="inline">{\, \hat{\mu }}_{C}</math>:590
591
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 592
|-593
| 594
{| style="text-align: center; margin:auto;width: 100%;"595
|-596
| <math>{\hat{\mu }}_{C}\cong \frac{1}{n}\sum _{i=1}^{n}{\hat{C}}_{i}\left( \boldsymbol{\rho },{\boldsymbol{z}}_{\mathit{\boldsymbol{i}}}\right)</math> 597
|}598
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(34)599
|}600
601
Para generar el modelo Kriging se ha utilizado la toolbox DACE [62], muy utilizada como así se recoge en la bibliografía [63-65].602
603
==3. Algoritmo de optimización==604
605
En este apartado se describe el algoritmo para la optimización de topología robusta de estructuras continuas bidimensionales considerando la existencia de incertidumbre en cargas concentradas. La incertidumbre de la carga se puede presentar en la magnitud, la dirección y/o la posición de la misma. Estas variables aleatorias se consideran estadísticamente independientes. La función objetivo a minimizar es el valor esperado de la ''compliance''. Para su evaluación de forma precisa, pero con un reducido número de evaluaciones, se combina el método de Monte Carlo con el modelo Kriging. A continuación se describen los pasos en los que se ha dividido el algoritmo:606
607
1. Definir las características del dominio inicial de la estructura (número de elementos en la dimensión ''x'', número de elementos en la dimensión ''y'', propiedades del material, regiones optimizables y no optimizables, etc.).608
609
2. Definir los parámetros de la optimización (fracción de volumen, factor de penalización, radio de filtro, tamaño de la muestra (''N<sub>MC</sub>''), tamaño del diseño de experimentos (''N<sub>K</sub>''), etc.).610
611
3. Definir las características de las cargas con incertidumbre (número de cargas con incertidumbre, identificar variables de la carga con aleatoriedad, pdf para cada una de las variables).612
613
4. Generar una muestra uniformemente distribuida de modo aleatorio según las pdfs <math display="inline">\boldsymbol{Z}=</math><math>{\left( {\boldsymbol{z}}_{\mathit{\boldsymbol{1\, }}}\mathit{\boldsymbol{\, }}{\boldsymbol{z}}_{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{\ldots \, }}{\boldsymbol{z}}_{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{MC}}}}\right) }^{T}</math> de tamaño ''N<sub>MC</sub>''.''' '''614
615
5. Generar el modelo Kriging utilizando la toolbox DACE.616
617
5.1. Diseño de experimentos mediante hipercubo latino con un tamaño de población de ''N<sub>k</sub>'' puntos, <math display="inline">\boldsymbol{S}=</math><math>{\left( {\boldsymbol{s}}_{\mathit{\boldsymbol{1\, }}}\mathit{\boldsymbol{\, }}{\boldsymbol{s}}_{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{\ldots \, }}{\boldsymbol{s}}_{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{k}}}}\right) }^{T}</math>.618
619
5.2. Obtener la ''compliance'', usando el modelo de simulación, <math display="inline">\boldsymbol{C}=</math><math>{\left( {C}_{1}\, {C}_{2}\ldots \, {C}_{{N}_{k}}\right) }^{T}</math> para la muestra generada <math display="inline">{\left( {\boldsymbol{s}}_{\mathit{\boldsymbol{1\, }}}\mathit{\boldsymbol{\, }}{\boldsymbol{s}}_{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{\ldots \, }}{\boldsymbol{s}}_{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{K}}}}\right) }^{T}</math>.620
621
5.3. Ajustar los parámetros del modelo Kriging a partir de la muestra generada y su ''compliance'' <math display="inline">\hat{y}(\mathit{\boldsymbol{S}}\boldsymbol{,C})</math>.622
623
6. Estimar la respuesta estructural <math display="inline">\hat{\boldsymbol{C}}=</math><math>{\left( {\hat{C}}_{1}\, {\hat{C}}_{2}\ldots \, {\hat{C}}_{{N}_{MK}}\right) }^{T}</math> para cada uno de los puntos de la muestra Monte Carlo <math display="inline">\boldsymbol{Z}=</math><math>{\left( {\boldsymbol{z}}_{\mathit{\boldsymbol{1\, }}}\mathit{\boldsymbol{\, }}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{\ldots \, }}{\boldsymbol{z}}_{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{MC}}}}\right) }^{T}</math>''' '''a partir del modelo Kriging.624
625
7. Evaluar el valor esperado de la ''compliance'' <math display="inline">\left( {\hat{\mu }}_{C}\right)</math> y la sensibilidad de la ''compliance'' usando las Ec. 11 y Ec. 7, respectivamente.626
627
8. Actualizar la distribución de las densidades mediante la Ec. 5.628
629
9. Si no se cumple el criterio de parada volver al punto 5. En caso contrario parar. Como criterio de parada se ha utilizado que la diferencia de ''compliance'' entre la iteración actual y la anterior sea inferior a un nivel mínimo de cambio.630
631
En la Figura 1 se muestra un diagrama de flujo del algoritmo de optimización.632
633
{| style="width: 100%;border-collapse: collapse;" 634
|-635
| style="text-align: center;vertical-align: top;"|<span id='_Toc433618559'></span>636
[[Image:Draft Samper 142890191-image1.png|center|288px]]637
638
|- style="text-align: center;vertical-align: top; font-size: 75%;"639
| |'''Figura 1.''' Diagrama de flujo. Algoritmo para la RTO considerando incertidumbre en la carga con el método de Monte Carlo con Kriging.640
|}641
642
==4. Ejemplos==643
644
En este apartado se pretende mostrar la validez del algoritmo propuesto, al que se denomina MCK (Monte Carlo con Kriging), mediante la resolución de dos ejemplos, una viga en voladizo y una T-invertida. Para los dos ejemplos el módulo de Young es de 1 para una región sólida y de 10<sup>-6 </sup>para una región vacía, el coeficiente de Poisson ''ν'' = 0,3. El tipo de elemento finito utilizado es cuadrado lineal de tamaño unidad. Para la optimización con SIMP se utiliza un factor de penalización ''p'' = 3 y un radio de filtro ''r'' = 2. El sistema de unidades es coherente.645
646
Con objeto de validar el algoritmo MCK, las soluciones obtenidas se han comparado con las obtenidas utilizando el método de simulación de Monte Carlo, que denominaremos MC. La precisión del método se determina a partir del error relativo <math display="inline">{E}_{r}</math>:647
648
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 649
|-650
| 651
{| style="text-align: center; margin:auto;width: 100%;"652
|-653
| <math>{E}_{r}=\frac{\left( {{\mu }_{C}}_{\left( MC\right) }-{{\mu }_{C}}_{\left( MCK\right) }\right) }{{{\mu }_{C}}_{\left( MC\right) }}\, \, 100</math>654
|}655
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(35)656
|}657
658
donde <math display="inline">{{\mu }_{C}}_{(MC)}</math> y <math display="inline">{{\mu }_{C}}_{(MCK)}</math> es el valor esperado de la ''compliance'' obtenida con MC y MCK, respectivamente. Por otro lado, la eficiencia del método se determina mediante ahorro en tiempo computacional <math display="inline">{E}_{t}</math>:659
660
{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 661
|-662
| 663
{| style="text-align: center; margin:auto;width: 100%;"664
|-665
| <math>{E}_{t}=\frac{\left( {t}_{\left( MC\right) }-{t}_{\left( MCK\right) }\right) }{{t}_{\left( MC\right) }}\, \, 100</math>666
|}667
| style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(36)668
|}669
670
donde <math display="inline">{t}_{(MC)\, }</math> y <math display="inline">{t}_{(MCK)}</math> son los tiempos medios por iteración con MC y MCK, respectivamente.671
672
En los dos ejemplos se han utilizado funciones base de orden cero para el modelo de regresión. Las funciones de base de orden cero son las más utilizadas en la bibliografía y constituyen el denominado ''Ordinary Kriging''. Para el modelo de correlación se ha adoptado el modelo exponencial generalizado, con ''p<sub>l</sub>'' = 1 y adoptando el mismo parámetro <math display="inline">{\theta }_{l}\,</math> para todas las direcciones coordenadas, lo que lo convierte en un modelo exponencial isótropo.673
674
===4.1. Viga en voladizo===675
676
El primer ejemplo para mostrar el efecto de considerar la existencia de incertidumbre en la magnitud, dirección y/o posición de la carga, es una viga en voladizo (Figura 2), la cual es optimizada mediante 6 configuraciones de carga (Tabla 2). La posición de la carga se define mediante dos variables ''S''<sub>x</sub> y ''S''<sub>y</sub> con respecto a la longitud y al ancho del dominio inicial. La carga se define mediante el módulo ''F'' y la dirección ''θ''. En la optimización se utiliza una restricción de volumen de <math display="inline">({V}_{f}/{V}_{0})=</math> 0,3, siendo <math display="inline">{V}_{f}</math> el volumen final y <math display="inline">{V}_{0}</math> el volumen inicial.677
678
{| style="width: 100%;border-collapse: collapse;" 679
|-680
| style="text-align: center;vertical-align: top;"|681
[[Image:Draft_Samper_142890191-image2.png|center|296px]]682
683
|- style="text-align: center;vertical-align: top; font-size: 75%;"684
| |'''Figura 2.''' Viga en voladizo. Dominio de diseño y carga aplicada.685
|}<br />686
687
Para las variables con incertidumbre que definen la carga (''F, θ'', ''S''<sub>x</sub>, ''S''<sub>y</sub>), se han considerado pdf de variables aleatorias normales N(ð; ð) de media ð y desviación estándar ð, y con truncamiento por la izquierda y por la derecha para las variables ''S''<sub>x</sub> y ''S''<sub>y</sub>. Las variables sin incertidumbre se han definido mediante su valor nominal.688
689
La configuración de carga 1 hace referencia a la optimización determinista; las configuraciones de carga 2, 3 y 4 hacen referencia a RTO donde solo se tiene en cuenta la incertidumbre en la magnitud, la dirección y la posición de la carga, respectivamente; la configuración de carga 5 hace referencia a una RTO donde se considera la existencia simultánea de incertidumbre en la magnitud y dirección de la carga; y la configuración de carga 6 hace referencia a una RTO donde se considera la existencia simultánea de incertidumbre en la magnitud, la dirección y la posición de la carga. Los parámetros utilizados para describir las pdfs que caracterizan la incertidumbre en la magnitud y la dirección de la carga son los utilizados por Dunning et al. [26]. Lógó [66] describe la incertidumbre de la posición de la carga mediante una pdf normal, pero no detalla los parámetros para su caracterización. En este trabajo la incertidumbre en la posición de la carga se define mediante una función normal, cuyos parámetros han sido elegidos de manera que la carga pueda situarse a lo largo de todo el borde lateral, pero con una mayor probabilidad de situarse sobre el centro.690
691
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size: 75%;">692
'''Tabla 2.''' Viga en voladizo. Parámetros de las configuraciones analizadas.</div>693
694
{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 695
|-696
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">Configuración</span>697
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">''F''</span>698
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">''θ ''(rad)</span>699
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">''S''<sub>x</sub></span><span style="text-align: center; font-size: 75%;"><br/>0 ≤ </span><span style="text-align: center; font-size: 75%;">''S''<sub>x</sub></span><span style="text-align: center; font-size: 75%;"> ≤ 1</span>700
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">''S''<sub>y</sub></span><span style="text-align: center; font-size: 75%;"><br/>0 ≤ </span><span style="text-align: center; font-size: 75%;">''S''<sub>y</sub></span><span style="text-align: center; font-size: 75%;"> ≤ 1</span>701
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">''N''<sub>MC</sub></span>702
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">''N''<sub>k</sub></span>703
|-704
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>705
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>706
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0 rad</span>707
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>708
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0,5</span>709
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">-</span>710
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">-</span>711
|-712
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">2</span>713
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N(1,0; 0,5)</span>714
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0 rad</span>715
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>716
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0,5</span>717
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10000</span>718
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10</span>719
|-720
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">3</span>721
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>722
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,25)</span>723
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>724
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0,5</span>725
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10000</span>726
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10</span>727
|-728
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">4</span>729
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>730
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0 rad</span>731
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>732
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N<sup>T</sup></span><span style="text-align: center; font-size: 75%;">(0,5; 0,17)</span>733
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10000</span>734
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10</span>735
|-736
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">5</span>737
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N(1,0; 0,5)</span>738
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,25)</span>739
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>740
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">0,5</span>741
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10000</span>742
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10</span>743
|-744
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">6</span>745
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N(1,0; 0,5)</span>746
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,25)</span>747
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1</span>748
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">N<sup>T</sup></span><span style="text-align: center; font-size: 75%;">(0,5; 0,17)</span>749
| style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10000</span>750
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">10</span>751
|}<br />752
753
''Configuración 1: Optimización determinista''754
755
La configuración 1 corresponde al caso determinista, donde todas las variables se han caracterizado por su valor nominal. La solución es una barra recta horizontal (Figura 3), obteniendo un valor para la ''compliance'' de 16,67. Sin embargo, si se analiza la topología determinista considerando las incertidumbres de los casos 2, 3, 4, 5 y 6, el valor esperado de la ''compliance'' aumenta a 15,45; 275,73; 450,77; 926,13 y 1763,91, respectivamente. Estos valores demuestran que la topología óptima determinista se puede comportar de forma poco fiable ante la presencia de cargas inciertas. Los valores esperados de la ''compliance'' para los casos 2, 3 y 4 revelan como la incertidumbre en la posición tiene una mayor influencia en el comportamiento estructural que la incertidumbre en la dirección, y como éstas tienen una mayor influencia que la incertidumbre en la magnitud.756
757
{| style="width: 100%;border-collapse: collapse;" 758
|- style="text-align: center;vertical-align: top; font-size: 75%;"759
| |'''760
[[Image:Draft_Samper_142890191-image3.png|center|240px]]761
'''762
|- style="text-align: center;vertical-align: top; font-size: 75%;"763
| |'''Figura 3.''' Viga en voladizo. Diseño óptimo determinista.764
|}<br />765
766
''Configuración 2: Incertidumbre en la magnitud de la carga''767
768
Para la configuración 2 se realiza la RTO considerando la existencia de incertidumbre en la magnitud de la carga. En la Figura 4, se muestra que la solución coincide con la obtenida para el diseño determinista ya que esta disposición es la que proporciona una mayor rigidez para cargas horizontales. El valor esperado de la ''compliance'' usando MC y MCK es de 14,19 y 13,78 respectivamente, lo que supone un error relativo del 2,88%. El tiempo computacional medio por iteración con MC y MCK es de 11389,4 y 6,81 segundos, respectivamente, lo que supone un ahorro del 99,94%.769
770
{| style="width: 100%;border-collapse: collapse;" 771
|-772
| style="text-align: center;vertical-align: top;"|'''773
[[Image:Draft_Samper_142890191-image4.png|center|246px]]774
'''775
| style="text-align: center;vertical-align: top;"|'''776
[[Image:Draft_Samper_142890191-image5.png|center|240px]]777
'''778
|-779
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>780
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>781
|- style="text-align: center;vertical-align: top; font-size: 75%;"782
| colspan='2' |'''Figura 4.''' Viga en voladizo. Diseño óptimo robusto considerando incertidumbre en la magnitud de la carga. (a) MC, (b) MCK.783
|}<br />784
785
''Configuración 3: Incertidumbre en la dirección de la carga''786
787
Para la configuración 3 se realiza la RTO considerando la existencia de incertidumbre en la dirección de la carga. Los diseños obtenidos (Figura 5) pasan de ser una barra horizontal, como en el caso de la optimización determinista, a ser dos barras inclinadas. El ángulo entre las barras depende del valor de la desviación estándar de la dirección de la carga, ya que cuanto mayor sea, mayor es la componente vertical de la fuerza. El valor esperado de la ''compliance'' con MC y MCK es de 27,68 y 29,96, respectivamente, lo que supone un error relativo del 2,60%. El tiempo computacional medio por iteración usando MC y MCK es de 20638,5 y 12,99 segundos, respectivamente, lo que supone un ahorro del 99,93%.788
789
{| style="width: 100%;border-collapse: collapse;" 790
|-791
| style="text-align: center;vertical-align: top;"|792
[[Image:Draft_Samper_142890191-image6.png|center|246px]]793
794
| style="text-align: center;vertical-align: top;"|795
[[Image:Draft_Samper_142890191-image7.png|center|246px]]796
797
|-798
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>799
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>800
|- style="text-align: center;vertical-align: top; font-size: 75%;"801
| colspan='2' |'''Figura 5.''' Viga en voladizo. Diseño óptimo robusto considerando incertidumbre en la dirección de la carga. (a) MC, (b) MCK.802
|}<br />803
804
''Configuración 4: Incertidumbre en la posición de la carga''805
806
Para la configuración 4 se realiza la RTO considerando la existencia de incertidumbre en la posición de la carga. Los diseños obtenidos (Figura 6) presentan una amplia zona con material en el extremo derecho del dominio, debido a las diferentes posiciones en las que puede situarse la carga. El valor esperado de la ''compliance'' es de 25,58 y 24,93 usando MC y MCK, respectivamente, lo que supone un error relativo del 2,54%. El tiempo computacional medio por iteración usando MC y MCK es de 10763,3 y 6,64 segundos, respectivamente, lo que supone un ahorro del 99,93%.807
808
{| style="width: 100%;border-collapse: collapse;" 809
|-810
| style="text-align: center;vertical-align: top;"|'''811
[[Image:Draft_Samper_142890191-image8.png|center|246px]]812
'''813
| style="text-align: center;vertical-align: top;"|'''814
[[Image:Draft_Samper_142890191-image9.png|center|246px]]815
'''816
|-817
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>818
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>819
|- style="text-align: center;vertical-align: top; font-size: 75%;"820
| colspan='2' |'''Figura 6.''' Viga en voladizo. Diseño óptimo robusto considerando incertidumbre en la posición de la carga. (a) MC, (b) MCK.821
|}<br />822
823
''Configuración 5: Incertidumbre en la magnitud y en la dirección de la carga''824
825
Para la configuración 5 se realiza la RTO considerando la existencia simultánea de incertidumbre en la magnitud y la dirección de la carga. Los diseños obtenidos (Figura 7) están constituidos por dos barras inclinadas, debido a que la influencia de la incertidumbre en la dirección de la carga es mucho más importante que la incertidumbre en la magnitud. El valor esperado de la ''compliance'' con MC y MCK es de 25,04 y 24,82, respectivamente, lo que supone un error del 0,87%. El tiempo computacional medio por iteración usando MC y MCK es de 16646,8 y 7,12 segundos, respectivamente, suponiendo un ahorro del 99,95%.826
827
{| style="width: 100%;border-collapse: collapse;" 828
|-829
| style="text-align: center;vertical-align: top;"|830
[[Image:Draft_Samper_142890191-image10.png|center|246px]]831
832
| style="text-align: center;vertical-align: top;"|833
[[Image:Draft_Samper_142890191-image11.png|center|246px]]834
835
|-836
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>837
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>838
|- style="text-align: center;vertical-align: top; font-size: 75%;"839
| colspan='2' |'''Figura 7.''' Viga en voladizo. Diseño óptimo robusto considerando incertidumbre en la magnitud y dirección de la carga. (a) MC, (b) MCK.840
|}<br />841
842
''Configuración 6: Incertidumbre en la magnitud, en la dirección y en la posición de la carga''843
844
Para la configuración 6 se realiza la RTO considerando de forma simultánea la existencia de incertidumbre en la magnitud, en la dirección y en la posición de la carga. En la Figura 8 se muestran las soluciones obtenidas. Las topologías presentan una disposición que viene determinada por la incertidumbre en la posición de la carga al ser ésta la que mayor influencia tiene sobre el problema. El valor esperado de la ''compliance'' con MC y MCK es de 52,13 y 49,96, respectivamente, lo que supone un error del 4,16%. El tiempo computacional medio por iteración usando MC y MCK es de 16162,5 y 7,02 segundos, respectivamente, lo que supone un ahorro del 99,95%.845
846
{| style="width: 100%;border-collapse: collapse;" 847
|-848
| style="text-align: center;vertical-align: top;"|849
[[Image:Draft_Samper_142890191-image12.png|center|246px]]850
851
| style="text-align: center;vertical-align: top;"|852
[[Image:Draft_Samper_142890191-image13.png|center|246px]]853
854
|-855
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>856
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>857
|- style="text-align: center;vertical-align: top; font-size: 75%;"858
| colspan='2' |'''Figura 8.''' Viga en Voladizo. Diseño óptimo robusto considerando incertidumbre en la magnitud, en la dirección y en la posición de la carga. (a) MC, (b) MCK.859
|}<br />860
861
A través de las configuraciones estudiadas, se puede concluir que: (1) las soluciones óptimas obtenidas usando RTO se comportan de forma más eficiente que la topología determinista ante la presencia de incertidumbres en la carga, (2) el algoritmo MCK proporciona soluciones igual de precisas que usando MC, y (3) el ahorro en tiempo computacional es como mínimo del 99,9% siendo el error en todos los casos estudiados inferior al 5%.862
863
===4.2. T-invertida===864
865
El segundo ejemplo es una estructura en forma de T invertida, las dimensiones del dominio de diseño se muestran en la Figura 9. La estructura está fija en su extremo superior, y sujeta a dos cargas concentradas ''F''<sub>1</sub> y ''F''<sub>2</sub>, de magnitud y dirección inciertas, aplicadas en los puntos medios de los extremos laterales. La estructura se ha optimizado bajo condiciones de carga determinista y considerando incertidumbres en la magnitud y en la dirección de la carga. Tanto para la optimización determinista, como para la optimización robusta, se ha utilizado una fracción de volumen de (''V''<sub>f </sub>/''V''<sub>o</sub>) = 0,5. Con este ejemplo se pretende demostrar el efecto de variar los parámetros de las pdfs que caracterizan las cargas. Para ello, se han estudiado dos configuraciones, donde las variables aleatorias se describen mediante las distribuciones normales recogidas en la Tabla 3. En la configuración 1 se han utilizado los parámetros usados por Dunning et al. [26], y en la configuración 2 se ha aumentado el valor de la desviación estándar de las direcciones de las cargas.866
867
{| style="width: 100%;border-collapse: collapse;" 868
|-869
| style="text-align: center;vertical-align: top;"|870
[[Image:Draft_Samper_142890191-image14.png|center|400px]]<br />871
872
|- style="text-align: center;vertical-align: top; font-size: 75%;"873
| |'''Figura 9. '''Estructura T-invertida. Dominio de diseño y cargas aplicadas.874
|}<br />875
876
<div class="center" style="width: auto; margin-left: auto; margin-right: auto;font-size: 75%;">877
'''Tabla 3.''' T-invertida. Parámetros de las configuraciones estudiadas.</div>878
879
{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 880
|-881
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Configuración</span>882
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">''F''<sub>1</sub></span>883
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">''θ''<sub>1</sub></span>884
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">''F''<sub>2</sub></span>885
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">''θ''<sub>2</sub></span>886
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">''N''<sub>MC</sub></span>887
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">''N''<sub>k</sub></span>888
|-889
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1</span>890
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(5,0; 0,5)</span>891
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,25)</span>892
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(5,0; 0,5)</span>893
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,25)</span>894
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10000</span>895
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10</span>896
|-897
| style="border-top: 1pt solid black;border-bottom: 1pt solid black;border-right: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2</span>898
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(5,0; 0,5)</span>899
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,5)</span>900
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(5,0; 0,5)</span>901
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">N(0,0; 0,5)</span>902
| style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10000</span>903
| style="border-top: 1pt solid black;border-left: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10</span>904
|}<br />905
906
La optimización determinista produce un diseño con un valor de ''compliance'' igual a 2786 (Figura 10). Si la topología óptima determinista obtenida se analiza considerando la existencia de incertidumbres en la magnitud y en la dirección de ambas cargas, el valor esperado de la ''compliance'' aumenta hasta 6792 y 7451, para la configuración 1 y 2, respectivamente. En el caso de realizar una RTO usando MCK se obtienen unos valores esperados para la ''compliance'' de 3328 y 3677, para la configuración 1 y 2, respectivamente. Al igual que en el primer ejemplo, se demuestra que los diseños obtenidos usando una formulación RTO se comportan mejor ante la existencia de incertidumbre que los diseños obtenidos mediante una formulación determinista.907
908
{| style="width: 100%;border-collapse: collapse;" 909
|-910
| style="text-align: center;vertical-align: top;"|911
[[Image:Draft_Samper_142890191-image15.png|center|240px]]<br />912
913
|- style="text-align: center;vertical-align: top; font-size: 75%;"914
| |'''Figura 10.''' Estructura T-invertida. Diseño óptimo determinista.915
|}<br />916
917
En las figuras 11 y 12 se muestran los diseños obtenidos usando MC y MCK, para la configuración 1 y 2, respectivamente. Estos diseños son muy similares a los obtenidos en la solución determinista salvo por las barras superiores centrales, las cuales pasan de estar inclinadas a verticales y reforzadas con crucetas. Esta disposición se debe a la componente horizontal de la fuerza, originada al considerar incertidumbre en la dirección de la misma. Además, se puede observar como al aumentar la desviación estándar en la dirección de la carga, el refuerzo de las crucetas también aumenta para proporcionar mayor rigidez lateral (Figuras 11b y 12b).918
919
Para determinar la precisión y eficiencia del algoritmo MCK, las soluciones se han comparado con las obtenidas usando MC. El valor esperado de la ''compliance'', usando MC para las configuraciones 1 y 2, es de 3341 y de 3742, lo que supone un error relativo de 0,38 y 1,71%, respectivamente. Por otro lado, el tiempo medio por iteración para la configuración 1 usando MCK es de 48 segundos, mientras que usando MC es de 99387, lo que supone un ahorro de tiempo del 99,95%. Para la configuración 2, el tiempo medio por iteración usando MCK es de 53 segundos, mientras que usando MC es de 114411, lo que supone un ahorro de tiempo del 99,95%.920
921
{| style="width: 100%;border-collapse: collapse;" 922
|-923
| style="text-align: center;vertical-align: top;"|924
[[Image:Draft_Samper_1428901913-image16.png|center|240px]]925
926
| style="text-align: center;vertical-align: top;"|927
[[Image:Draft_Samper_142890191-image17.png|center|240px]]928
929
|-930
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>931
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>932
|- style="text-align: center;vertical-align: top; font-size: 75%;"933
| colspan='2' |'''Figura 11.''' Estructura T-invertida. Diseño óptimo robusto considerando incertidumbre en la magnitud y dirección de la carga para la configuración 1 ( <math display="inline">{\sigma }_{\theta }=</math><math>0,25).\,</math> (a) MC, (b) MCK.934
|-935
| style="text-align: center;vertical-align: top;"|936
[[Image:Draft_Samper_142890191-image18.png|center|240px]]937
938
| style="text-align: center;vertical-align: top;"|939
[[Image:Draft_Samper_142890191-image19.png|center|240px]]940
941
|-942
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(a)</span>943
| style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">(b)</span>944
|- style="text-align: center;vertical-align: top; font-size: 75%;"945
| colspan='2' |'''Figura 12.''' Estructura T-invertida. Diseño óptimo robusto considerando incertidumbre en la magnitud y dirección de la carga para la configuración 2 ( <math display="inline">{\sigma }_{\theta }=</math><math>0,5)</math>. (a) MC, (b) MCK.946
|}947
948
==5. Conclusiones==949
950
En este trabajo se presenta un algoritmo, denominado Monte Carlo con Kriging (MCK), para la optimización de topología robusta considerando la existencia de incertidumbres en las cargas. Para ello, el problema de RTO se ha formulado de forma probabilista, tomando como función objetivo el valor esperado de la ''compliance''. En cada iteración, el valor esperado de la ''compliance'' se evalúa por medio del método de simulación de Monte Carlo en combinación con un modelo Kriging. El algoritmo propuesto puede generalizarse a otros tipos de optimización de topología determinista y a problemas tridimensionales.951
952
La precisión y eficiencia del método se ha demostrado mediante dos ejemplos. En estos ejemplos, las cargas se definen por medio de cuatro variables, la magnitud (''F''), la dirección (''θ''), la posición respecto al eje ''x'' (''S''<sub>x</sub>) y la posición con respecto al eje ''y'' (''S''<sub>y</sub>), que pueden ser aleatorias. Las variables no aleatorias se caracterizan por su valor nominal, y las variables aleatorias se caracterizan mediante funciones de densidad de probabilidad normales, aunque el algoritmo permite utilizar cualquier otro tipo de pdf. Los ejemplos estudiados demuestran que:953
954
1. El algoritmo MCK propuesto proporciona diseños casi tan precisos como los obtenidos usando el método de simulación de Monte Carlo estándar, pero de forma mucho más eficiente. Esto es debido a que el método MCK propuesto solo necesita realizar ''N''<sub>K </sub>análisis del modelo de simulación, en vez de los ''N''<sub>MC </sub>necesarios con el método de Monte Carlo (''N''<sub>K</sub> = 10 << ''N''<sub>MC</sub> = 10000). Además, el número de evaluaciones con el método propuesto es independiente del número de cargas inciertas, siendo muy útil cuando el problema presenta una alta dimensionalidad.955
956
2. El valor esperado de la ''compliance'' para los diseños obtenidos mediante una formulación robusta es menor que el valor para el diseño obtenido mediante una formulación determinista ante la existencia de incertidumbre en la carga. Esto demuestra también la importancia de incorporar la existencia de incertidumbre en el proceso de optimización de topología.957
958
Los ejemplos también muestran el efecto de considerar diferentes fuentes y niveles de incertidumbre:959
960
3. En el ejemplo 1 se demuestra que la incertidumbre en la posición de la carga tiene una mayor influencia que la incertidumbre en la dirección, y que estas a su vez tienen mayor influencia que la magnitud.961
962
4. En el ejemplo 2, se muestra como el diseño robusto se adapta al nivel de incertidumbre.963
964
==Agradecimientos==965
966
Los autores agradecen el apoyo financiero del Ministerio de Economía y Competitividad a través del proyecto de investigación “Sistema integrado de diseño óptimo robusto de topología de estructuras (DPI2011-27939-C02-01)”, el cual ha permitido la realización de esta investigación.967
968
==Referencias==969
970
[1] Bendsøe M.P., Kikuchi N. Generating optimal topologies in structural design using a homogenisation method. Comp. Meth. Appl. Mech. Engng., 71:197-224, 1988.971
972
[2] Rozvany G.I.N., Zhou M., BirkerT. Generalized shape optimization without homogenization. Struct. Optim., 4:250-254, 1992.973
974
[3] Xie Y.M., Steven G.P. A simple evolutionary procedure for structural optimization. Comput. Struct., 49(5):885-896, 1993.975
976
[4] Kirkpatrick S., Gellat C.D., Vecchi M.P. Optimization by simulated annealing. Science, 220(4598):671-680, 1983.977
978
[5] Sethian J.A. Level set method and fast marching methods. 2nd Ed. New york, Cambridge University Press, 1999.979
980
[6] Victoria M., Martí P., Querin O.M. Topology design of two-dimensional continuum structures using isolines. Comput. Struct., 87:101-09, 2009.981
982
[7] Deaton J.D., Grandhi R.V. A survey of structural and multidisciplinary continuum topology optimization: post 2000. Struct. Multidisc. Optim., 49:1-38, 2014.983
984
[8] Anflor C.T.M., Marczak R.J. Topological sensitivity analysis for two-dimensional heat transfer problems using the boundary element method. Optimization of Structures and Components, 43:11-33, 2013.985
986
[9] Dühring M.B., Jensen J.S., Sigmund O. Acoustic design by topology optimization. J. Sound Vib., 317(3-5):557-575, 2008.987
988
[10] Abdelwahed M., Hassine M., Masmoudi M. Optimal shape design for fluid flow using topological perturbation technique. J. Math. Anal. Appl., 356(2):548-563, 2009.989
990
[11] Radman A., Huang X., Xie Y.M. Topological optimization for the design of microstructures of isotropic cellular materials. Eng. Optim., 45(11):1331-1348, 2013.991
992
[12] Jung H.S., Cho S. Reability-based topology optimization of geometrically nonlinear structures with loading and material uncertainties. Finite Elem. Anal. Des., 41(3):311-331, 2004.993
994
[13] Kim C., Wang S., Bae K-r, Moon H., Choi K.K. Reliability-based topology optimization with uncertainties. J. Mech. Sci. Tech., 20(4):494-504, 2006.995
996
[14] Mogami K., Nishiwaki S., Izui K., Yoshimura M., Kogiso N. Reliability-based structural optimization of frame structures for multiple failure criteria using topology optimization techniques. Struct. Multidisc. Optim., 32(4):299-311, 2006.997
998
[15] Yoo K.S., Eom Y.S., Park J.Y., Im M.G., Han S.Y. Reliability-based topology optimization using successive standard response surface method. Finite Elem. Anal. Des., 47(7):843-849, 2011.999
1000
[16] Kharmanda G., Olhoff N., Mohamed A., Lemaire M. Reliability-based topology optimization. Struct. Multidis. Optim., 26(5):295-307, 2004.1001
1002
[17] Zhao Q., Chen X., Ma Z.D., Lin Y. Reliability-based topology optimization using stochastic response surface method with sparse grid design. Math. Probl. Engng., 2015:13, 2015.1003
1004
[18] Chen S., Chen W., Lee S. Level set based robust shape and topology optimization under random field uncertainties. Struct. Multidisc. Optim., 41:507524, 2010.1005
1006
[19] Chen S., Chen W. A new level-set based approach to shape and topology optimization under geometric uncertainty. Struct. Multidisc. Optim., 44:1-18, 2011.1007
1008
[20] Dunning P.D., Kim H.A. Robust topology optimization: minimization of expected and variance of compliance. AIAA Journal, 51(11):2656-64, 2013.1009
1010
[21] Schevenels M., Lazarov B.S., Sigmund O. Robust topology optimization accounting for spatially varying manufacturing errors. Comput. Methods Appl. Mech. Engrg., 200:3613-3627, 2011.1011
1012
[22] Jansen M., Lombaert G., Diehl M., Lazarov B.S., Sigmund O., Schevenels M. Robust topology optimization accounting for misplacement of material. Struct. Multidisc. Optim., 47:317-33, 2013.1013
1014
[23] Zhao J., Wang C. Robust topology optimization under loading uncertainty based on linear elastic theory and orthogonal diagonalization of symmetric matrices. Comput. Methods Appl. Mech. Engrg., 273:204-18, 2014.1015
1016
[24] Zhao J., Wang C. Robust structural topology optimization under random field loading uncertianty. Struct. Multidisc. Optim., 50:517-22, 2014.1017
1018
[25] Zhao Q., Chen X., Ma Z.D., Lin Y. Robust topology optimization based on stochastic collocation methods under loading uncertainties. Math. Probl. Engng., 1:1-14, 2015.1019
1020
[26] Dunning P.D., Kim H.A., Mullineux G. Introducing loading uncertainty in topology optimization. AIAA Journal, 49:760-768, 2011.1021
1022
[27] Zhao J., Wang C. Robust topology optimization of structures under loading uncertainty. AIAA Journal, 52(2):398-407, 2014.1023
1024
[28] García-Lopez N.P., Sanchez-Silva M., Medaglia A.L., Chateauneuf A. An improved robust topology optimization approach using multiobjective evolutionary algorithms. Comput. Struct., 125:1-10, 2013.1025
1026
[29] Tootkaboni M., Asadpoure A., Guest J. Topology optimization of continuum structures under uncertainty: a polynomial chaos approach. Comp. Methods Appl. Mech. Engrg., 201-204:263-275, 2012.1027
1028
[30] Lazarov B.S., Schevenels M., Sigmund O. Topology optimization considering material and geometric uncertainties iusing stochastic collocations methods. Struct. Multidisc. Optim., 46:597-612, 2012.1029
1030
[31] Conti S., Held H., Pach M., Rumpf M., Schultz R. Shape optimization under uncertainty: a sctochastic programming perspective. SIAM J. Optim., 19(4):1610-32, 2009.1031
1032
[32] Myers H., Montgomery D.C. Response surface methodology: process and product optimization using designed experiments. 2nd edition, John Wiley & Sons, New York, USA, 2002.1033
1034
[33] Powell M. The theory of radial basis function approximation in 1990. In Advances in Numerical Analysis, Vol.2: Wavelets, Subdivision Algorithms and Radial Basis Functions, W. Light (Ed.). Oxford University Press, Oxford, UK, 1992.1035
1036
[34] Papadrakakis M., Lagaros N., Tsompanakis Y. Structural optimization using evolution strategies and neural networks. Comp. Meth. in Appl. Mech. Engrg., 156(1-4):309-333, 1998.1037
1038
[35] Martin J.D., Simpson T.W. Use of Kriging models to approximate deterministic computer models. AIAA Journal, 43(4):853-863, 2005.1039
1040
[36] Dubourg V., Sudret B., Bourine J.-M. Reliability-based design optimization using Kriging surrogates and subset simulation. Struct. Multidisc. Optim., 44:673-690, 2011.1041
1042
[37] Chen Z., Peng S., Li X., Qiu H., Xiong H., Gao L., Li P. An important boundary sampling method for reliability-based design optimization using Kriging model. Struct. Multidisc. Optim., 52:55-70, 2015.1043
1044
[38] Lee K.H., Kang D.H. A robust optimization using the statistics base on kriging metamodel. J. Mech. Sci. Techn., 20(8):1169-1182, 2006.1045
1046
[39] Jurecka F. Robust design optimization based on metamodeling techniques. Technische Universität München, München, 2007.1047
1048
[40] Martínez J., Martí P. Diseño óptimo robusto utilizando modelos Kriging: aplicación al diseño óptimo robusto de estructuras articuladas. Rev. Int. Mét. Num. Cálc. Diseño Ing., 30(2):97-105, 2014.1049
1050
[41] Bendsøe M.P. Optimal shape design as a material distribution problem. Struct. Optim., 1:193-202, 1989.1051
1052
[42] Zhou M., Rozvany G. The COC algorithm. Part II: topological, geometry and generalized shape optimization. Comput. Methods Appl. Mech. Engrg., 89:197-224, 1991.1053
1054
[43] Svanberg K. The method of moving asymptotes - a new method for structural optimization. Int. J. Numer. Meth. Engrg., 24:359-373, 1987.1055
1056
[44] Sigmund O. A 99 line topology optimization code written in Matlab. Struct. Multidisc. Optim., 21:120-127, 2001.1057
1058
[45] Sigmund O., Peterson J. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Struct. Multidisc. Optim., 16:68-75, 1998.1059
1060
[46] BrunsT.E., Tortorelli D.A. Topology optimization of non-linear elastic structures and compliant mechanisms. Comput. Methods. Appl. Mech. Engrg., 190(26-27):3443-3459, 2001.1061
1062
[47] Bourdin B. Filters in topology optimization. Int. J. Numer. Meth. Engrg., 50:2143-2158, 2011.1063
1064
[48] Guest J.K., Prévost J.H., Belytschko T. Achieving minimum length scale in topology optimization using nodal design variables and projection functions. Int. J. Numer. Meth. Engrg., 61(2):238-254, 2004.1065
1066
[49] Sigmund O. Morphology-based black and white filters for topology optimization. Struct. Multidisc. Optim., 33(4-5):401-424, 2007.1067
1068
[50] Xu S., Cai Y., Cheng G. Volume preserving nonlinear density filter based on heaviside functions. Struct. Multidisc. Optim., 41:495-505, 2010.1069
1070
[51] Kawamoto A., Matsumori T., Yamakasi S., Nomura T., Kondoh T., Nishiwaki S. Heaviside projection based topology optimization by a pde-filtered scalar function. Struct. Multidisc. Optim., 44(1):19-24, 2008.1071
1072
[52] Lazarov B.S., Sigmung O. Filers in topology optimization based on Helmholtz type differential equations. Int. J. Numer. Methods Engrg., 86(6):765-781, 2011.1073
1074
[53] Metropolis N., Ulam S. The Monte Carlo method. J. Amer. Statist. Assoc., 44(247):335-341, 1949.1075
1076
[54] Rubinstein R.Y., Kroese D.P. Simulation and the Monte Carlo method. 2nd Edition, John Wiley & Sons, 2008.1077
1078
[55] Schuceller G.I. Developments in stochastic structural mechanics. Arch. Appl. Mech., 75:755-773, 2006.1079
1080
[56] Krige D.G. A statistical approach to some basic mine valuation problems on the Witwatersrand. Journal of the Chemical, Metallurgical and Mining Society of South Africa, 52:119-139, 1951.1081
1082
[57] Nakayama H., Arakawa M., Sasaki R. Simulation based optimization using computational intelligence. Optimiz. Eng., 3:201-214, 2002.1083
1084
[58] Papadrakakis M., Lagaros N.D., Tsompanakis Y. Structural optimization using evolution strategies and neural networks. Comput. Methods. Appl. Mech. Engrg., 156:309-333, 1998.1085
1086
[59] Jin R., Du X., Chen W. The use of metamodeling techniques for optimization under uncertainty. Struct. Multidisc. Optim., 25:99-116, 2003.1087
1088
[60] Jones D.R. A taxonomy of global optimization methods based on response surfaces. J. Global Optimization, 21:345-383, 2001.1089
1090
[61] Matheron G. Principles of geostatistics. Economic Geology, 58(8):1246-1266, 1963.1091
1092
[62] Lophaven N., Nielsen H.B., Sondergaard J. Aspects of the Matlab toolbox DACE. Technical Report, Informatics and Mathematical Modelling, Technical University of Denmark, DTU, Richard Petersens Plads, Lyngby, 2002.1093
1094
[63] Echard B., Gayton N., Lemaire M. AK-MCS: An active learning realiability method combining Kriging and Monte Carlo simulation. Struct. Saf., 33(2):145-154, 2011.1095
1096
[64] Kaymaz I. Application of Kriging method to structural reliability problems. Struct. Saf., 27(2):133-151, 2005.1097
1098
[65] Kleijnen J.P.C. Kriging metamodeling in simulation: a review. Eur. J. Oper. Res., 192(3):707-716, 2009.1099
1100
[66] Lógó J. SIMP type topology optimization procedure considering uncertain load position. Civil Engineering, 56(2):213-219, 2012.1101
Return to Cordero et al 2017a.
Back to Top
Document information
Published on 03/01/18
Accepted on 13/12/16
Submitted on 02/05/16
Volume 34, Issue 1, 2018
DOI: 10.23967/j.rimni.2017.5.005
Licence: CC BY-NC-SA license
Share this document
claim authorship
Are you one of the authors of this document?