Un planteamiento probabilista de los criterios optimizantes mediante el principio de la máxima entropía

Resumen

Este trabajo describe brevemente el debate científico que se produjo en su momento en el ámbito de la optimización de estructuras entre los partidarios de los métodos matemáticos de optimización de problemas no lineales y los de los criterios optimizantes. En ese sentido recuerda el planteamiento que investigadores como A. B. Templeman llevaron a cabo y que conceptualmente consistía en establecer una propuesta probabilista de los criterios optimizantes, cuyo mérito radicaba en ser completamente racional y no prejuzgar las condiciones que iban a ser activas en el diseño óptimo. Tras revisar el estado de la cuestión, se ha desarrollado una nueva formulación que simplifica y mejora las anteriores y transforma el problema de optimización en otro consistente en resolver sucesivamente sistemas de ecuaciones de la misma dimensión que el número de variables del problema. El algoritmo requiere la obtención de las sensibilidades de primer orden, con lo que también puede ser considerado como un nuevo método de gradiente. Tras la descripción teórica, el método se ha aplicado para resolver 2 ejemplos de optimización de estructuras de celosía muy comunes en la literatura de esta disciplina científica. Los resultados numéricos que se muestran permiten concluir que esta formulación obtiene la solución óptima con gran precisión.

Abstract

This work reviews A. B. Templeman's research in the maximum entropy principle applied to structural optimization, which consists in defining a probabilistic approach for the optimality criteria. This new reasoning does not assume aprioristically which conditions could be active in the optimum point, giving a completely rational point of view to the optimization problem. After reviewing the state of the art, a new formulation, which simplifies and improves previous methodologies, is developed. This new method turns the optimization problem into solving successively systems of equations with the same dimension as the number of design variables. Since this algorithm requires obtaining first order sensitivities, it can be considered as a new gradient method. After the theoretical description, the algorithm is tested with two examples of truss structures. Results show a good performance of the algorithm, based on comparisons with other alternatives found in the literature.

1

Introducción

Está comúnmente aceptado que el inicio del planteamiento moderno de la optimización estructural lo marca la publicación de L. Schmit [1] en 1960 de un trabajo sobre lo que entonces se denominaba síntesis estructural. Sin embargo, los comienzos de esta disciplina —y este periodo de tiempo puede abarcar más de 2 décadas— estuvieron marcados por la disputa entre 2 líneas de pensamiento respecto al concepto de diseño óptimo de estructuras.

Una de ellas se basaba en el citado trabajo de Schmit, que formulaba el problema como la optimización de una función F(X) condicionada por un conjunto de condiciones g_j(X) (j = 1,…,m), dependientes todas ellas del vector de variables de diseño X, y escrito como sigue:

sujeto a:

Las condiciones se correspondían con respuestas estructurales como tensiones, movimientos, frecuencias de vibración, etc., que se obtenían por medio de modelos estructurales, y el problema definido por (1) se resolvía mediante métodos matemáticos de optimización de problemas no lineales [2–5].

La otra línea de pensamiento se sustentaba en la capacidad de intuición y el sentido común ingenieril para identificar las causas que hacen que un diseño estructural sea óptimo, y en sus inicios esgrimía que un proyectista, gracias a su experiencia en el problema en estudio, podía conocer de antemano los criterios que permitirían obtener el deseado diseño óptimo, es decir, conocer las condiciones que serían activas en la solución del problema (1). Por este motivo recibían el nombre de criterios optimizantes (optimality criteria)[6], y entre ellos los más conocidos son los denominados como modos simultáneos de colapso (simultaneous failure modes)[7] y criterio de agotamiento (fully stressed design)[8], aunque también se formularon otros criterios para tratar los movimientos de las estructuras [9] o el fenómeno de pandeo [10]. Las limitaciones intrínsecas a estos métodos se intentaron solucionar adoptando planteamientos con más soporte matemático, a fin de hacerlos más versátiles [11].

En el debate científico entre ambas líneas de investigación había quedado demostrado teóricamente que el apriorismo de los optimality criteria propuestos hasta entonces les restaba validez general, pero el hecho de que en determinados tipos de problemas acertasen con el diseño óptimo suponía un acicate para quienes persistían en utilizarlos. De hecho, algunos investigadores han sido fervientes partidarios de esa línea de trabajo durante todo su periodo de actividad académica [12].

Sin embargo, en ocasiones el debate intentó superarse con planteamientos unificadores. A este respecto, cabe resaltar que los métodos duales, propuestos por Fleury [13], ponen de manifiesto que un planteamiento general de los criterios optimizantes no era otra cosa que una linealización del problema de optimización no lineal descrito en (1), y ello explicaba su éxito en algunos tipos de problemas.

Un planteamiento más ambicioso fue abordado por Templeman, y en él se intentaba eliminar todo tipo de sesgo inicial en la asunción de criterios optimizantes y, en consecuencia, permitir una formulación más racional de esta metodología. Estaba basado en el principio de la máxima entropía [14] y se desarrollará en el siguiente apartado.

2

Identificación probabilística de las condiciones activas mediante el principio de máxima entropía

2.1

Procedimiento de doble fase de optimización

El planteamiento elaborado por Templeman para el algoritmo de optimización consistía en aprovechar la condición de Kuhn-Tucker [15], que expresa que, en el punto óptimo, el vector opuesto al gradiente de la función objetivo puede descomponerse entre los vectores gradientes de las condiciones activas, sin que ninguna de las componentes sea negativa. Es decir:

siendo m_a el conjunto de condiciones activas en el punto óptimo. La expresión anterior puede generalizarse en la forma:

siendo:

donde m_p es el conjunto de condiciones pasivas, que es complementario de m_a. Como consecuencia de lo anterior, en el diseño óptimo para cualquier condición g_j(X) y componente λ se cumplirá:

dado que, o bien la condición es activa, y por tanto g_j(X) = 0, o bien la condición es pasiva, y por tanto λ = 0. A partir de ello, Templeman [16] sustituía la formulación del problema de optimización (1) por el siguiente:

sujeto a:

En esta formulación las incógnitas están constituidas por el vector de variables de diseño X y el de componentes λ. Para resolver la formulación, se establecía un proceso en 2 fases. En la primera de ellas se fijaban los valores de λ y se optimizaban las variables de diseño, mientras que en la segunda se actualizaba el vector λ manteniendo los valores de X obtenidos en la etapa anterior. El proceso anterior se repite el número de veces necesario hasta conseguir la convergencia. A continuación se describen estas 2 fases en detalle.

• Fase 1. Optimización de las variables de diseño X

El problema por resolver es:

sujeto a:

Para elegir el valor inicial, cada λ_j se interpreta como la probabilidad de que una condición sea activa. Ello es coherente con el hecho de que cuando una condición g_j (X) sea pasiva su valor asociado λ_j será nulo. Dado que al comienzo del proceso no existe ningún indicio basado en información racional que permita asignar valores diferentes a las componentes del vector λ, se adoptan valores iguales para todas ellas, es decir:

Una vez adoptado este criterio, cualquier algoritmo de optimización permitirá obtener los valores óptimos del resto de variables de diseño X1*, que harán activas algunas condiciones y mantendrán otras pasivas. Esto es inconsistente con los valores de λ_j definidos inicialmente, y por ello la siguiente fase consiste en actualizar el vector λ.

• Fase 2. Actualización de los componentes del vector λ

Para actualizar el vector λ de forma neutral, es decir, sin favorecer a los valores que corresponden a unas condiciones del problema frente a otros, se utilizaba el principio de la máxima entropía, que indica que los nuevos valores de λ_j serán los que maximicen la siguiente función:

sujeta a las condiciones:

La condición (8b) se introduce para normalizar los valores de λ_j, y la expresión (8c) describe el hecho de que al no corresponder las componentes de λ y las condiciones g_j(X) a la misma iteración, el sumatorio correspondiente no se anula, sino que produce un residuo ¿_k que irá decreciendo en las sucesivas iteraciones. Por ello, para ¿_k deberá adoptarse en cada iteración una expresión que sea decreciente. Una muy habitual es ¿_k = e^−k.

La solución de (8) proporciona los nuevos valores de λ y, por tanto, con ellos puede reiniciarse la fase 1 y, tras ella, la fase 2, prosiguiendo sucesivamente de esa manera hasta concluir el proceso cuando exista la adecuada convergencia entre los valores de las variables de diseño X y las componentes del vector λ.

2.2

Procedimiento con una sola fase de optimización

El problema definido por (8) es una optimización condicionada con restricciones de igualdad que se puede resolver mediante el clásico método de multiplicadores de Lagrange, definiendo la función lagrangiana siguiente:

y estableciendo las condiciones de punto estacionario como:

La solución del sistema de ecuaciones resultante permite expresar los valores de λ_j en la forma:

y el multiplicador μ puede obtenerse de la expresión:

Al tener de forma explícita la expresión de los elementos del vector λ, el doble bucle de optimización, uno para las variables de diseño y otro para λ, puede reducirse a uno solo. Es decir, el problema definido por (2) puede expresarse como:

sujeta a:

El valor de μ debería haberse obtenido de la expresión (12), pero dado que su solución no es simple, el procedimiento seguido hasta ahora era resolver el problema (13) repetidamente con diferentes valores de μ hasta alcanzar la convergencia (16).

3

Nuevo planteamiento para la resolución del problema de optimización

El problema definido por (13) se puede resolver por cualquier método de optimización numérica, aunque su implementación computacional ha mostrado que el planteamiento producía inestabilidades numéricas de cierta importancia.

Templeman [17] puso de manifiesto que el planteamiento evolucionaba muy rápidamente en las primeras iteraciones pero se ralentizaba posteriormente y podían aparecer problemas numéricos que aparentemente se asociaban al crecimiento del parámetro de control μ. Asimismo, Erlander [18,19] usó esta metodología en problemas de transporte y encontró las mismas dificultades.

No obstante, si se repara en que se trata de una optimización condicionada de igualdad, para resolverlo se puede plantear la creación de una función lagrangiana:

y establecer las condiciones de punto estacionario como:

tras llevar a cabo los cálculos correspondientes se llega a:

si se linealizan las condiciones del problema:

resulta:

Las expresiones anteriores se pueden escribir de una forma mucho más condensada si se define la terminología siguiente:

Como consecuencia, el problema queda definido de la siguiente manera:

La expresión (20a) representa un sistema de ecuaciones lineales de la misma dimensión que el número de las variables de diseño. Dado un diseño inicial X₀, la solución del problema descrito por la expresión (20) permite encontrar el nuevo punto X como:

Resolviendo el sistema de ecuaciones de (20a) de forma independiente para ▿F y para el otro sumando, se puede escribir:

entrando en (20b):

tras alguna operación algebraica:

sustituyendo ese valor de γ en (22) y a través de (21), se obtiene el nuevo diseño X.

En resumen, las etapas del proceso son las siguientes:

• 1)

  Se define un valor inicial de μ.

• 2)

  Se define el valor inicial X₀.

• 3)

  Se evalúan las condiciones g_j, las sensibilidades ∂gj∂xi y las probabilidades λj=eμgj∑k=1meμgk.

• 4)

  Se forman los vectores y matrices λ, Λ, g, G, ∇ G, ∇ F.

• 5)

  Se resuelve el sistema de ecuaciones para obtener ΔX y con ello el nuevo diseño X.

• 6)

  Si todos los componentes de ΔX no son lo suficientemente pequeños, se vuelve al paso 3.

• 7)

  Cuando ya lo son, se compara el vector X con el obtenido con el valor anterior de μ y, si no se ha conseguido la convergencia, se vuelve al paso 3.

• 8)

  Si se ha alcanzado la convergencia, se ha obtenido el diseño óptimo.

El procedimiento descrito permite ser interpretado desde distintos puntos de vista:

• a)

  Desde la perspectiva de λ, es un planteamiento probabilista de los criterios de optimización.

• b)

  Desde el punto de vista de ΔX, es un algoritmo de gradiente, ya que ese vector representa un avance a lo largo de una dirección y se obtiene utilizando sensibilidades de primer orden.

• c)

  Desde el punto de vista de ΔG, el planteamiento consiste en una secuencia de problemas lineales.

• d)

  Para los problemas lineales, dado que las condiciones son expresiones de primer orden, el método no requiere iteraciones en γ. Se trata, por tanto, de una variante de los métodos de punto interior, inaugurados por Karmarkar [20] y, en consecuencia, una alternativa al método simplex.

• e)

  Desde una perspectiva general, es interesante observar que el planteamiento adoptado reduce la resolución del problema de optimización a la tarea de calcular repetidamente sistemas de ecuaciones lineales cuya dimensión es la misma que la del número de variables de diseño que tenga el problema, siendo indiferente a este respecto el número de condiciones existentes.

En consecuencia, el planteamiento presenta varias novedades, como unificar los métodos de criterios optimizantes (optimality criteria) y de programación matemática —que históricamente eran 2 líneas separadas de investigación— y ser un nuevo método de resolución de problemas lineales mediante las denominadas técnicas de punto interior, ya que la aproximación lineal que se lleva a cabo no supone ninguna simplificación en ese tipo de problemas. Por otra parte, para los problemas no lineales supone un nuevo método de gradiente mucho más simple que el método de las direcciones eficientes, que para la obtención de la dirección de búsqueda requiere la formulación y la resolución de un problema de optimización lineal. Por el contrario, en el método formulado por los autores solo es necesario resolver un sistema de ecuaciones lineales.

El método planteado abre una serie de líneas de investigación, como por ejemplo, entre otras: el estudio del rango de magnitudes del parámetro μ, que representa una ponderación del tipo de las usadas en los métodos de función penalty, puesto que los valores de ese parámetro pueden resultar dependientes del problema.

Se puede definir una segunda línea en relación con el sistema de ecuaciones que debe resolverse en cada iteración; debe recordarse que su rango es el mismo que el del número de variables de diseño, pero requiere conocer las sensibilidades de todas las condiciones del problema. A este respecto sería interesante utilizar esquemas de constraint screening a fin de considerar únicamente en cada etapa las condiciones que superen un umbral que permita pensar que puedan ser activas. Esta técnica, que es muy habitual en otros métodos, puede ser implementada muy sencillamente, ya que los valores de los elementos del vector λ proporcionan dichas probabilidades.

4

Ejemplos de aplicación

El algoritmo anterior se ha aplicado a la optimización de 2 modelos estructurales de celosías utilizados en la literatura técnica.

4.1

Estructura plana en celosía

Se trata de una estructura plana en celosía (fig. 1) con cordón superior curvo con 25 barras cuyas áreas transversales son las variables de diseño [21]. La geometría de la misma aparece en la figura, y las acciones exteriores están constituidas por una carga concentrada P = 300 kN que puede actuar en cualquiera de los nudos del cordón inferior, por lo que existen 5 casos diferentes de carga.

Figura 1.

Estructura plana en celosía.

El módulo de elasticidad del material es E = 210 GPa. La función objetivo del problema es el volumen del material que se desea minimizar, y las condiciones son solamente de tipo tensional. Con todo ello, la formulación del problema queda de la siguiente forma:

sujeto a:

Se ha considerado la simetría del problema a la hora de resolverlo y, por ello, se ha asignado una variable de diseño por cada 2 barras simétricas.

En las figuras 2–4 y en la tabla 1 se presentan los resultados obtenidos por el algoritmo basado en el principio de la máxima entropía propuesto por los autores. En la figura 2 aparece la evolución de la función objetivo. El hecho de que la función objetivo aumente debe interpretarse teniendo en cuenta que se parte de un diseño inicial inválido, es decir, un diseño que no cumple las condiciones de limitación de tensión. Por lo tanto, al proceder con la optimización es necesaria una aportación de material para obtener un diseño válido en la convergencia al mínimo local de la función objetivo.

Figura 2.

Evolución de la función objetivo.

Tabla 1.

Área y tensión en las barras en el diseño óptimo mediante el algoritmo propuesto

Variable	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11	x12	x13
Barras	1,11	2,12	3,9	4,10	5,7	6,8	13,17	14,16	15	18,24	19,25	20,23	21,22
Área (cm^2)	43,34	16,56	36,33	5,99	33,18	3,86	11,81	12,15	13,14	5,81	16,60	10,85	6,81
σ (MPa)	−104	130	−104	130	−104	−104	130	130	130	−104	121	130	130

La figura 3 informa de la tensión máxima que se obtiene en cada iteración en cada una de las barras para el conjunto de casos de carga, y se observa que las barras a compresión convergen al valor de –104 MPa y las barras a tracción al de 130 MPa, excepto las barras 19 y 25.

Figura 3.

Evolución de la tensión en las barras.

En la figura 4 se presenta la evolución de las variables de diseño, que son las áreas de las secciones transversales de cada barra, y en la tabla 1 aparecen sus valores finales.

Figura 4.

Evolución del área de las barras.

Para proceder a la validación del algoritmo propuesto, este ejemplo se ha resuelto también mediante el algoritmo SQP implementado en el software de optimización Altair Optistruct [22] y mediante el método modificado de las direcciones eficientes, implementado en el código Genesis [23].

Los resultados obtenidos por los algoritmos anteriores difieren, ya que el algoritmo SQP converge al fully stressed design, es decir, el diseño en el que todas las barras trabajan a la máxima tensión. Sin embargo, el método modificado de las direcciones eficientes alcanza el diseño óptimo, donde las barras 19 y 25 no necesitan trabajar a su máxima capacidad y, pese a ello, se obtiene un diseño con un volumen de material inferior al fully stressed design. El algoritmo propuesto por los autores converge a este mismo resultado, que a su vez coincide con el resultado clásico de la literatura [21].

En la tabla 2 se presenta una comparación de los valores de las variables de diseño en el resultado final alcanzado por cada algoritmo. Los valores correspondientes de la función objetivo son 0,238 m³ para el método modificado de las direcciones eficientes y para el algoritmo de la máxima entropía, y 0,239 m³ para el método de secuencia de problemas cuadráticos.

Tabla 2.

Comparación del área de las barras en el diseño óptimo mediante el algoritmo propuesto, el método modificado de las direcciones eficientes (MMFD) y el de secuencia de problemas cuadráticos (SQP)

Variable	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11	x12	x13
Barras	1,11	2,12	3,9	4,10	5,7	6,8	13,17	14,16	15	18,24	19,25	20,23	21,22
Propuesto	43,34	16,56	36,33	5,99	33,18	3,86	11,81	12,15	13,14	5,81	16,60	10,85	6,81
MMFD	43,36	16,58	36,28	5,88	33,18	3,83	11,81	12,17	13,13	5,61	16,91	10,85	6,81
SQP	43,37	15,87	37,11	7,35	32,99	5,98	13,90	12,01	13,45	8,77	12,57	10,53	6,59

4.2. Estructura espacial de celosía

Se trata de una estructura espacial en celosía de 25 barras cuyas áreas transversales son las variables de diseño [24], que toman como valor inicial 1 cm² y como límite inferior 0,1 cm². La geometría de la misma aparece en la figura 5, con las coordenadas de los nudos indicadas en la tabla 3. Las acciones exteriores están constituidas por 2 cargas concentradas P = 128 kN y otras 2 cargas de 126 kN dispuestas como se muestra en la tabla 4 y en la figura 5.

Figura 5.

Estructura espacial en celosía.

Tabla 3.

Coordenadas de los nudos de la estructura

Nudo	X (m)	Y (m)	Z (m)
1	−0,95	0,00	5,08
2	0,95	0,00	5,08
3	−0,95	0,95	2,54
4	0,95	0,95	2,54
5	0,95	−0,95	2,54
6	−0,95	−0,95	2,54
7	−2,54	2,54	0,00
8	2,54	2,54	0,00
9	2,54	−2,54	0,00
10	−2,54	−2,54	0,00

Tabla 4.

Cargas actuantes sobre la estructura

Nudo	Cargas (kN)
	X	Y	Z	Resultante
1	−88,9	88,9	−22,6	128
2	88,9	−88,9	−22,6	128
3	−88,9	88,9	0	126
5	88,9	−88,9	0	126

El módulo de elasticidad del material es E = 70,6 GPa. La función objetivo del problema es el volumen del material que se desea minimizar, y las condiciones son de tipo tensional. Con todo ello, la formulación del problema queda de la siguiente forma:

sujeto a:

En las figuras 6–8 aparece la evolución de la función objetivo, de la tensión en las barras y del valor de su área a lo largo de las iteraciones, y en la tabla 5 aparecen los valores finales de aquellas. Al igual que en el ejemplo anterior, la función objetivo crece dado que el diseño inicial es inválido y se necesita la aportación de material para satisfacer las condiciones de diseño en tensión. Los resultados obtenidos se han verificado mediante el algoritmo SQP implementado en el software de optimización Altair Optistruct [22] y mediante el método modificado de las direcciones eficientes implementado en el código Genesis [23]. El diseño óptimo proporcionado por estos códigos coincide con el obtenido utilizando la formulación presentada por los autores, arrojando únicamente diferencias en cifras decimales no significativas.

Figura 6.

Evolución de la función objetivo.

Figura 7.

Evolución de la tensión en las barras.

Figura 8.

Evolución del área de las barras.

Tabla 5.

Área y tensión máxima en las barras en el diseño óptimo

Variable	x1	x2	x3	x4	x5	x6	x7	x8	x9	x10	x11	x12	x13
Barra	1	2,5	3,4	6,9	7,8	10,11	12,13	14,17	15,16	18,21	19,20	22,25	23,24
Área (cm^2)	0,58	0,44	5,04	0,10	4,67	0,58	0,10	2,27	0,59	0,10	0,10	0,10	2,33
σ (MPa)	276,0	−276,0	276,0	81,7	−276,0	276,0	126,8	276,0	−276,0	31,1	−95,8	51,9	−276,0

5

Conclusiones

Se pueden extraer las siguientes conclusiones de esta investigación:

• 1.

  Se ha revisado un planteamiento existente que proponía de forma probabilista las técnicas de diseño optimizantes (optimality criteria), haciendo uso del principio de la máxima entropía.

• 2.

  Se ha generado un nuevo algoritmo que actualiza y mejora las formulaciones anteriores y requiere un único bucle de optimización para la resolución del problema.

• 3.

  El algoritmo puede entenderse como una versión de las formulaciones basadas en secuencia de problemas lineales.

• 4.

  Para los problemas de optimización lineal el algoritmo requiere muy pocas iteraciones y es una alternativa al método simplex.

• 5.

  La obtención del nuevo diseño en cada iteración requiere el uso de derivadas de primer orden, por lo que la técnica presentada puede entenderse como una variante de los métodos de gradiente.

• 6.

  El nuevo algoritmo transforma el problema de optimización en otro consistente en resolver sucesivos sistemas de ecuaciones, cuya dimensión es la del número de variables del problema y es independiente del número de condiciones.

• 7.

  La formulación descrita se ha aplicado a la resolución de problemas no lineales de optimización estructural y los resultados obtenidos reproducen con gran precisión los que se pueden encontrar en las referencias existentes.