Guía de optimización en ingeniería

Este artículo contiene una guía básica para plantear y resolver un problema de optimización en ingeniería. Para fijar ideas nos centraremos en un ejemplo sencillo e iremos explicando el proceso sobre ese ejemplo.

1 Ejemplo: definición del problema

Consideramos un trapecio construido de un material homogéneo del que podemos elegir la inclinación del lado derecho (ver la figura). Buscamos la inclinación de ese lado derecho que proporcione menor área, pero manteniendo la componente horizontal del centro de gravedad dentro de la base para darle estabilidad.

Buscamos el trapecio con menor área cuyo centro de gravedad caiga en la base, cuando movemos el lado derecho (en rojo) de manera que cambie el punto a donde se apoya

Es importante identificar en este punto dos cosas:

1. Coste que queremos optimizar
2. Restricciones que tenemos en el problema

En este caso sencillo, el coste consiste en el área total del trapecio mientras que la restricción es que el centro de gravedad caiga en la base.

2 Planteamiento matemático

Vamos a plantear este problema matemáticamente. Los pasos a seguir son:

1. Identificar los grados de libertad y asociarle uno o varios parámetros para describir todas las posibles geometrías. Llamaremos a los parámetros [math] (a_1,a_2,...,a_n)[/math]
2. Escribir en función de los parámetros anteriores: la función coste [math] f_{coste}(a_1,a_2,...,a_n)[/math], los intervalos [math] [L_i,M_i] [/math] en los que se pueden mover los parámetros [math] a_i, \; i=1,...,n,[/math] y las restricciones en términos de los parámetros [math] g_{rest}(a_1,a_2,...,a_n)[/math].

Una vez hecho esto el planteamiento matemático es: Buscamos los valores [math] (a_1,a_2,...,a_n)[/math] que minimicen el coste

[math] f_{coste} (a_1,a_2,...,a_n) [/math]

cuando los parámetros toman valores en los intervalos

[math] L_i \leq a_i \leq M_i, \quad i=1,2,...,n, [/math]

y además verifican las restricciones

[math] \begin{array}{l} g_{rest1} (a_1,a_2,...,a_n)\leq 0\\ g_{rest2} (a_1,a_2,...,a_n)\leq 0\\ ...\\ g_{restk} (a_1,a_2,...,a_n)\leq 0\\ \end{array} [/math]

Vamos a seguir estos dos pasos en nuestro ejemplo:

2.1 Identificar los grados de libertad y asociarles parámetros

En el ejemplo que tratamos sólo hay un grado de libertad que es la pendiente del lado derecho del trapecio. En general existen varios parámetros que podemos elegir para representar este grado de libertad: la distancia a, la propia pendiente del lado en rojo, el ángulo que forma ese lado con la horizontal, etc. Conviene elegir el que nos haga más fáciles los cálculos.

En este caso elegimos la longitud de la base a. En situaciones más complicadas habrá más parámetros para describir toda la variabilidad del problema. Por ejemplo, podemos incluir la altura h o la longitud del lado superior b. Sin embargo, para fijar ideas, en lo sucesivo estos dos parámetros los tomaremos constantes e iguales a 1.

[math] h=b=1. [/math]

2.2 Escribir la función coste y las restricciones en términos de los parámetros

La función coste es el área del trapecio, que viene dada por

[math] f_{coste}(a)=h\frac{b+a}{2}=\frac{1+a}{2}. [/math]

El parámetro [math]a[/math] supondremos varía en el intervalo [math][0,1][/math] puesto que otros valores bien no tienen sentido (por ejemplo [math]a\lt0[/math]) o bien podemos descartarlos por ser claramente poco eficientes (por ejemplo [math]a\gt1[/math]).

En cuanto a las restricciones, en este caso sólo tenemos una: que la coordenada x del centro de gravedad caiga en la base. Si llamamos [math]S[/math] a la superficie que ocupa el trapecio, la coordenada [math]x_G[/math] del centro de gravedad viene dada por

[math] x_{G}=\frac{\iint_S x dS}{\iint_S dS}=\frac{\int_0^a\int_0^1 x \; dydx+\int_a^1\int_{\frac{x-a}{1-a}}^1x\;dydx}{\frac{1+a}{2}} =\frac{a^2+a+1}{3(a+1)}. [/math]

Por lo tanto, la restricción se escribe [math] x_{G}=\frac{a^2+a+1}{3(a+1)} \leq a [/math]

El planteamiento matemático del problema es entonces: Buscamos [math]a[/math] que minimice el coste

[math] f_{coste}(a)=\frac{1+a}{2} [/math]

en el intervalo

[math] 0\leq a \leq 1, [/math]

con la restricción

[math] g_{rest}(a)=\frac{a^2+a+1}{3(a+1)}-a \leq 0 [/math]

3 Solución en Matlab: Fuerza bruta

En esta sección vamos a ver cómo hacer que Matlab evalúe la función coste y las restricciones para cada posible parámetro. De esta forma, evaluando en un número suficientemente elevado de parámetros y eligiendo el de coste mínimo obtenemos un sencillo método de optimización.

Empezaremos incluyendo en Matlab la función coste. Existen varias formas de definir una función en Matlab. La más sencilla tiene la siguiente sintaxis:

% Función coste
f_coste = @(x) (1+x)/2;

Observamos que aparece en primer lugar el nombre de la función, f_coste. Las variables de las que depende la función se colocan después de la @ y entre paréntesis. Por último se escribe la expresión de la función en términos de las variables.

Una vez definida esta función en Matlab podemos calcular el coste de cualquier parámetro simplemente escribiendo f_coste(0.1), f_coste(0.4), etc.

A continuación añadimos funciones para cada una de las restricciones. Como sólo tenemos una, tendríamos:

% Función coste
f_coste = @(x) (1+x)/2;
% Restricciones
g_rest = @(x) (x^2+x+1)/(3*(x+1))-x;

Una vez definidas las funciones vamos a implementar nuestro primer método de optimización. La idea es elegir un buen número de parámetros que aproximen todos los posibles valores del intervalo en el que buscamos el mínimo y evaluemos tanto la función objetivo como las restricciones. Si al probar un candidato obtenemos un valor menor que el que teníamos, nos quedamos con ese candidato. El código en Matlab sería:

%% Código para optimizar el ejemplo
% función coste
f_coste=@(x) 1+x/2;

% Restricción
g_rest=@(x) (x^2+x+1)/(3*(x+1))-x;

%% Optimización por fuerza bruta
% Intervalo de parámetros:
h=0.01;                 % distancia entre dos valores del parámetro
L=0;                    % extremo izquierdo
M=1;                    % extremo derecho
a=L:h:M;                % vector con los parámetros
% Inicialización
Optimo=10000;           % inicializamos con un valor muy alto
indice_optimo=0;        % inicializamos el indice donde está el óptimo con 0 
% Probamos todos los valores a(i) con i>=1
for i=1:length(a)
    restriccion=g_rest(a(i));     % evaluamos la restricción
    if restriccion<0              % solo en este caso a(i) es un parámetro admisible y sigo adelante
        coste = f_coste(a(i));    % Calculo el coste
        if coste<Optimo          % en este caso mejoramos!
            Optimo = coste;       % actualizo el coste con el nuevo parámetro
            indice_optimo=i;      % actualizo el índice en el que está el óptimo
        end
    end        
end
sprintf('El óptimo se alcanza en a = %f',a(indice_optimo))

Conviene tener en cuenta que el programa anterior no calcula el óptimo del problema de optimización. Tan solo calcula el mejor de los candidatos probados. Si tomamos una división más fina o más gruesa del conjunto de parámetros, el valor óptimo obtenido cambia. Es también habitual tomar una serie de valores aleatorios dentro del conjunto de parámetros. En este caso también el óptimo que da el programa cambia.

Si hemos elegido suficientes puntos para describir el conjunto de parámetros se suele tomar el óptimo como una aproximación del problema de optimización original.