Recocido simulado (Simulated Annealing)

El recocido simulado (Simulated Annealing, SA) es un algoritmo aleatorizado, que aproxima el óptimo global de una función. Se llama algoritmo aleatorizado porque emplea una cierta cantidad de azar en su búsqueda y por lo tanto su salida puede variar para la misma entrada.

El problema

Se da una función $E(s)$ , que calcula la energía del estado $s$ . Se nos pide hallar el estado $s_{best}$ en el que $E(s)$ se minimiza. SA es adecuado para problemas donde los estados son discretos y $E(s)$ tiene múltiples mínimos locales. Tomaremos el ejemplo del problema del viajante (Travelling Salesman Problem, TSP).

Problema del viajante (TSP)

Se da un conjunto de nodos en el espacio bidimensional. Cada nodo se caracteriza por sus coordenadas $x$ e $y$ . La tarea es hallar un ordenamiento de los nodos que minimice la distancia a recorrer al visitar estos nodos en ese orden.

Motivación

El recocido es un proceso metalúrgico, en el que se calienta un material y se lo deja enfriar, para permitir que los átomos del interior se reordenen en un arreglo con energía interna mínima, lo que a su vez hace que el material tenga propiedades distintas. El estado es el arreglo de átomos y la energía interna es la función que se minimiza. Podemos pensar el estado original de los átomos como un mínimo local de su energía interna. Para hacer que el material reordene sus átomos, necesitamos motivarlo a atravesar una región donde su energía interna no está minimizada para alcanzar el mínimo global. Esta motivación se da calentando el material a una temperatura más alta.

El recocido simulado, literalmente, simula este proceso. Empezamos con algún estado aleatorio (material) y fijamos una temperatura alta (lo calentamos). Ahora, el algoritmo está listo para aceptar estados que tienen una energía mayor que el estado actual, porque está motivado por la alta temperatura. Esto evita que el algoritmo se quede atascado dentro de mínimos locales y se mueva hacia el mínimo global. Conforme avanza el tiempo, el algoritmo se enfría y rechaza los estados con mayor energía y se mueve hacia el mínimo más cercano que ha hallado.

La función de energía E(s)

$E(s)$ es la función que hay que minimizar (o maximizar). Mapea cada estado a un número real. En el caso del TSP, $E(s)$ devuelve la distancia de recorrer un círculo completo en el orden de nodos del estado.

Estado

El espacio de estados es el dominio de la función de energía, $E(s)$ , y un estado es cualquier elemento que pertenece al espacio de estados. En el caso del TSP, todos los caminos posibles que podemos tomar para visitar todos los nodos es el espacio de estados, y cualquiera de estos caminos se puede considerar como un estado.

Estado vecino

Es un estado del espacio de estados que está cerca del estado anterior. Esto usualmente significa que podemos obtener el estado vecino a partir del estado original usando una transformación simple. En el caso del problema del viajante, un estado vecino se obtiene eligiendo al azar 2 nodos e intercambiando sus posiciones en el estado actual.

Algoritmo

Empezamos con un estado aleatorio $s$ . En cada paso, elegimos un estado vecino $s_{next}$ del estado actual $s$ . Si $E(s_{next}) < E(s)$ , entonces actualizamos $s = s_{next}$ . En caso contrario, usamos una función de aceptación por probabilidad $P(E(s),E(s_{next}),T)$ que decide si debemos movernos a $s_{next}$ o quedarnos en $s$ . $T$ aquí es la temperatura, que inicialmente se pone a un valor alto y decae lentamente con cada paso. Cuanto más alta la temperatura, más probable es moverse a $s_{next}$ . Al mismo tiempo también llevamos un registro del mejor estado $s_{best}$ a través de todas las iteraciones. Procedemos hasta convergencia o hasta que se acaba el tiempo.

Una representación visual del recocido simulado, buscando los máximos de esta función con múltiples máximos locales.

Temperatura (T) y decaimiento (u)

La temperatura del sistema cuantifica la disposición del algoritmo a aceptar un estado con una energía mayor. El decaimiento es una constante que cuantifica la “tasa de enfriamiento” del algoritmo. Se sabe que una tasa de enfriamiento lenta ( $u$ más grande) da mejores resultados.

Función de aceptación por probabilidad (PAF)

$P(E,E_{next},T) =$

Aquí, $\mathcal{U}$ es un valor aleatorio uniforme continuo en $[0,1]$ . Esta función toma el estado actual, el siguiente estado y la temperatura, y devuelve un valor booleano, que le dice a nuestra búsqueda si debe moverse a $s$ {next} $s_{n e x t}$ o quedarse en $s$ . Nótese que para $E_{next} < E$ , esta función siempre devolverá True; en caso contrario todavía puede hacer el movimiento con probabilidad $\exp(-\frac{E_{next}-E}{T})$ , que corresponde a la medida de Gibbs .


bool P(double E,double E_next,double T,mt19937 rng){
    double prob =  exp(-(E_next-E)/T);
    if(prob > 1) return true;
    else{
        bernoulli_distribution d(prob); 
        return d(rng);
    }
}

Plantilla de código


class state {
    public:
    state() {
        // Generate the initial state
    }
    state next() {
        state s_next;
        // Modify s_next to a random neighboring state
        return s_next;
    }
    double E() {
        // Implement the energy function here
    };
};
 
 
pair<double, state> simAnneal() {
    state s = state();
    state best = s;
    double T = 10000; // Initial temperature
    double u = 0.995; // decay rate
    double E = s.E();
    double E_next;
    double E_best = E;
    mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
    while (T > 1) {
        state next = s.next();
        E_next = next.E();
        if (P(E, E_next, T, rng)) {
            s = next;
            if (E_next < E_best) {
                best = s;
                E_best = E_next;
            }
            E = E_next;
        }
        T *= u;
    }
    return {E_best, best};
}

Cómo usarlo:

Rellenar las funciones de la clase state según corresponda. Si se está intentando hallar un máximo global y no un mínimo, hay que asegurarse de que la función $E()$ devuelva el negativo de la función que se está maximizando e imprimir $-E_{best}$ al final. Fijar los parámetros de abajo según la necesidad.

Parámetros

$T$ : Temperatura inicial. Ponerla a un valor más alto si se quiere que la búsqueda corra durante más tiempo.
$u$ : Decaimiento. Decide la tasa de enfriamiento. Una tasa de enfriamiento más lenta (valor más grande de $u$ ) usualmente da mejores resultados, a costa de correr durante más tiempo. Asegurarse de que $u < 1$ .

El número de iteraciones que correrá el bucle está dado por la expresión

$N = \lceil -\log_{u}{T} \rceil$

Consejos para elegir $T$ y $u$ : si hay muchos mínimos locales y un espacio de estados amplio, poner $u = 0.999$ , para una tasa de enfriamiento lenta, que permitirá al algoritmo explorar más posibilidades. Por otro lado, si el espacio de estados es más estrecho, $u = 0.99$ debería bastar. Si no se está seguro, conviene ir a lo seguro poniendo $u = 0.998$ o más. Calcular la complejidad temporal de una sola iteración del algoritmo, y usar esto para aproximar un valor de $N$ que evite TLE, luego usar la fórmula de abajo para obtener $T$ .

$T = u^{-N}$

Implementación de ejemplo para TSP


 
class state {
    public:
    vector<pair<int, int>> points;
	std::mt19937 mt{ static_cast<std::mt19937::result_type>(
		std::chrono::steady_clock::now().time_since_epoch().count()
		) };
    state() {
        points = {%raw%} {{0,0},{2,2},{0,2},{2,0},{0,1},{1,2},{2,1},{1,0}} {%endraw%};
    }
    state next() {
        state s_next;
        s_next.points = points;
        uniform_int_distribution<> choose(0, points.size()-1);
        int a = choose(mt);
        int b = choose(mt);
        s_next.points[a].swap(s_next.points[b]);
        return s_next;
    }
 
    double euclidean(pair<int, int> a, pair<int, int> b) {
        return hypot(a.first - b.first, a.second - b.second);
    }
    
    double E() {
        double dist = 0;
        int n = points.size();
        for (int i = 0;i < n; i++)
            dist += euclidean(points[i], points[(i+1)%n]);
        return dist;
    };
};
 
int main() {
    pair<double, state> res;
    res = simAnneal();
    double E_best = res.first;
    state best = res.second;
    cout << "Length of shortest path found : " << E_best << "\n";
    cout << "Order of points in shortest path : \n";
    for(auto x: best.points) {
        cout << x.first << " " << x.second << "\n";
    }
}

Modificaciones adicionales al algoritmo:

Añadir una condición de salida basada en tiempo al bucle while para evitar TLE
El decaimiento implementado arriba es un decaimiento exponencial. Siempre se puede reemplazar esto por una función de decaimiento según las necesidades.
La función de aceptación por probabilidad dada arriba prefiere aceptar estados que tienen menor energía por el factor $E_{next} - E$ en el numerador del exponente. Se puede simplemente quitar este factor, para hacer la PAF independiente de la diferencia de energías.
El efecto de la diferencia de energías, $E_{next} - E$ , sobre la PAF se puede aumentar/disminuir aumentando/disminuyendo la base del exponente como se muestra abajo:


bool P(double E, double E_next, double T, mt19937 rng) {
    double e = 2; // set e to any real number greater than 1
    double prob =  pow(e,-(E_next-E)/T);
    if (prob > 1)
        return true;
    else {
        bernoulli_distribution d(prob); 
        return d(rng);
    }
}