Algoritmos genéticos

En este artículo presentamos una idea de cómo funcionan los algoritmos genéticos para minimizar problemas con restricciones y mostramos la implementación práctica en un caso sencillo. Supondremos que la función coste y las restricciones dependen sólo de dos variables para simplificar, pero todo se puede trasladar fácilmente a casos más generales.

1 Planteamiento del problema de optimización

Consideramos el problema siguiente: Buscamos minimizar la función coste [math] f_{coste} (a_1,a_2) [/math]

cuando los parámetros toman valores en los intervalos

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

y además verifican las restricciones (supondremos también sólo 2 para simplificar)

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

2 Idea general de un algoritmo genético

Este tipo de algoritmos está basado en un proceso iterativo en el que se pasa de una muestra aleatoria de puntos a otra simulando un ciclo vital:

1. Selección del mejor o mejores.
2. Combinación de los mejores con el resto.
3. Añadir un factor aleatorio.

El proceso es como sigue:

• Elegimos una muestra aleatoria inicial de nuestros parámetros, [math] \{a_i \}_{i=1}^N [/math]
• Elegimos el número de iteraciones J: veces que vamos a modificar la muestra
• Aplicamos el ciclo vital: para j = 1, ..., J; [math] \{a_i \}_{i=1}^N \to \{b_i \}_{i=1}^N[/math]

A. Selección:

1. Evaluamos la función coste en los puntos de la muestra [math] \{a_i \}_{i=1}^N [/math].
2. Determinamos los M mejores, supongamos [math] \{m_i \}_{i=1}^M [/math].
3. Eliminamos de nuestra muestra los P peores.
4. El resto los llamamos restantes [math] \{r_i \}_{i=1}^R [/math]

B. Combinación: Combinamos los restantes con los M mejores y los mejores entre si;

[math]c_{ij}=(m_i+r_j)/2; [/math]

[math]d_{ij}=(m_i+m_j)/2; [/math]

C. Efecto aleatorio: Añadimos nuevos puntos aleatorios [math] \{l_i \}_{i=1}^L [/math]

D. Construimos la nueva muestra:

[math]\{b_i \}_{i=1}^N = \{m_i \}_{i=1}^M \cup \{c_{ij} \}_{i,j=1}^{M,R} \cup \{d_{ij} \}_{i,j=1}^{M,M} \cup \{l_i \}_{i=1}^L [/math]

• Fin del ciclo vital: Renombramos [math]\{b_i \}_{i=1}^N \to \{a_i \}_{i=1}^N[/math] y repetimos el ciclo vital
• Evaluamos la función en los puntos de la última muestra y nos quedamos con el mejor

Observaciones:

1. El algoritmo presentado no tiene por qué dar una aproximación del mínimo. No existe nigún resultado matemático general que garantice este hecho si no se imponen condiciones a la función objetivo y las restricciones. Sin embargo, a menudo si que proporciona una buena aproximación y por ello es un algoritmo bastante popular.
2. Este método es muy costoso pero significativamente menos que la fuerza bruta.
3. En el algoritmo anterior hay muchos parámetros que pueden modificarse dando lugar a variantes.

Para ver un ejemplo de cómo aplicar con Matlab un algoritmo genético ver el siguiente enlace: Ejemplo de algoritmo genético

3 Programa en MATLAB

A continuación escribimos nuestro algoritmo de optimización de presas de gravedad usando el algoritmo genético. Observamos que primero de todo debemos crear una función con las restricciones, que debe tener como salida las restricciones de igualdad y desigualdad en un formato específico. Para ello, creamos en un fichero aparte la función simple_constraint.m con las siguientes líneas:

function [c, ceq] = simple_constraint(n,g_rest1,g_rest2,g_rest3);
c=[];
ceq=[g_rest1(n);g_rest2(n);g_rest3(n)];

Una vez creada esta función la usamos en nuestro programa como sigue:

% Este programa calcula el talud óptimo de aguas abajo de una presa de gravedad 
% manteniendo la estabilidad y las condiciones de tensiones en cimentación 
% en la hipótesis de tensión plana. Supondremos que el perfil es un triángulo del cual 
% sólo podemos elegir el talud de aguas abajo que llamaremos n. El talud de 
% aguas arriba supondremos que es vertical. Buscaremos el perfil con área mínima.

% Parámetros específicos del problema
NMN       = 170;  % Nivel máximo normal de embalse (msnm)
NAP       = 174;  % Cota del vértice resistente (se considera en este caso igual al nivel de avenida de proyecto) (msnm)
Cota_cim  = 108;  % Cota de la base de cimentación, que en este caso es horizontal (msnm)
Ter_abajo = 118;  % Cota del terreno natural aguas abajo (se considera que el nivel del río aguas abajo está a esa cota) (msnm)
mu_h      = 2.4;  % peso específico del hormigón en t/m3
mu_a      = 1.0;  % peso específico del agua en t/m3
Phi       = pi/4; % ángulo de rozamiento en el contacto presa-cimiento (radianes)
F_Phi     = 1.2;  % coeficiente de seguridad al rozamiento exigido por la normativa en la hipótesis de drenes ineficaces
S_Amin    = 5;  % compresión mínima en el pie de aguas arriba (t/m2)
S_Bmax    = 84;  % compresión máxima en el pie de aguas abajo (t/m2) 

% Variables de optimización: n = pendiente aguas abajo

% variables y funciones útiles para el cálculo (las acciones se evalúan en toneladas-fuerza por metro lineal de sección)
h_NMN = NMN-Cota_cim;                       % altura de agua a embalse lleno sobre la base de cimentación (m)
H_a   = NAP-Cota_cim;                       % altura de la presa (m)
h_CE  = Ter_abajo-Cota_cim;                 % subpresión en el pie de aguas abajo (mca)
base  = @(n) H_a*n;                         % ancho de la base del triángulo (m)
P     = @(n) 0.5*base(n)*H_a*mu_h;          % peso de la sección de hormigón (t/m)
E     = 0.5*h_NMN^2*mu_a;                   % empuje hidrostático con embalse lleno (no depende de n)
S     = @(n) 0.5*(h_NMN+h_CE)*base(n)*mu_a; % empuje debido a la subpresión en cimentación (t/m)
N     = @(n) P(n)-S(n);                     % Normal, obtenida por el equilibrio de fuerzas verticales
% Necesitaremos también las tensiones en los pies de aguas arriba y aguas abajo: S_A y S_B
% Para ello usaremos lo siguiente: 
%    1. Supondremos que la distribución de las tensiones verticales en la base
%    varía linealmente desde el pie de aguas arriba (donde toma el valor S_A) hacia el de
%    aguas abajo (donde toma el valor S_B).
%    2. El equilibrio de momentos en el pie de aguas abajo
% Estas dos condiciones se escriben:
% Crec. lineal de ten. nor.    N = 1/2*(S_A+S_B)*base(n)
% 0 = Momento del peso         0 = - P(n)*base(n)*2/3
% + momento empuje agua            + E*h_NMN/3
% + momento subpresion S_2         + h_CE*base(n)*(base(n)/2)*mu_a
% + momento subpresion S_1         + 0.5*base(n)*(h_NMN-h_CE)*2/3*base(n)  
% + momento de la normal N_2       + S_A*10*base(n)*(base(n)*0.5) 
% + momento de la normal N_1       + 0.5*base(n)*(S_B-S_A)*(base(n)/3)
% Tenemos 2 ecuaciones con dos incógnitas. Despejando sacamos el valor de
% S_A y S_B
S_A = @(n) (P(n)*base(n)*2/3 - E*h_NMN/3 ...
    - h_CE*base(n)*(base(n)/2)*mu_a ...
    - 0.5*base(n)*(h_NMN-h_CE)*2/3*base(n)...
    - 0.5*base(n)*(2*N(n)/base(n))*base(n)/3)/(base(n)^2*(1/2-1/3));
S_B = @(n) 2*N(n)/base(n)-S_A(n);
% Una vez definidas todas las funciones que necesitamos, calculamos la
% función coste y las restricciones

% Función coste: área del triángulo
f_coste = @(n) 0.5*base(n)*H_a;

% Calculamos ahora las restricciones:
% %% 1. Estabilidad frente al deslizamiento: E-T<0
% donde T es la fuerza necesaria para vencer el rozamiento:
% T = N*tan(Phi)/F_Phi (despreciando la cohesión y añadiendo el coef. seg. al rozamiento)
g_rest1 = @(n) E-N(n)*tan(Phi)/F_Phi; 
% la estabilidad requiere g_rest1(n) < 0

% %% 2. Compresión en el pie de aguas abajo menor que el máximo admitido
% por el cimiento. 
g_rest2 = @(n) S_B(n)-S_Bmax;
% la estabilidad requiere g_rest2(n) < 0

% %% 3. El pie de aguas arriba debe mantenerse comprimido, con un valor superior a S_Amin.
g_rest3 = @(n) S_Amin-S_A(n);
% la estabilidad requiere g_rest3(n) < 0

% %% Planteamiento matemático del problema de optimización:
% Buscamos n en el intervalo [0.1,80] que minimice la función coste:
% f_coste(n)
% con las restricciones:
% g_rest1(n)<0, g_rest2(n)<0, g_rest3(n)<0.

%% Optimización usando el algoritmo genético de Matlab
% Intervalo de parámetros:
LB  = 0.1;         % extremo izquierdo de los intervalos
UB  = 80;          % extremo derecho de los intervalos

% Escribimos la función de restricciones
Constraint = @(n) simple_constraint(n,g_rest1,g_rest2,g_rest3);

[n,fval] = ga(f_coste,1,[],[],[],[],LB,UB,Constraint);

sprintf('El óptimo con algoritmo genético se alcanza en n = %f',n)

%% Optimización usando la función fmincon de Matlab (programación no lineal) 

n1 = fmincon(f_coste,1,[],[],[],[],LB,UB,Constraint);

sprintf('El óptimo con la función fmincon se alcanza en n = %f',n1)