RECOCIDO SIMULADO (sudoku solver en python)

     El recocido simulado es una técnica de búsqueda aleatoria que sirve para resolver problemas de optimización global, dicho método imita el proceso de recocido de materiales, cuando un metal se enfría y posteriormente alcanza un estado cristalino de menor energía con un tamaño de cristal más grande para reducir los defectos en las estructuras metálicas. El proceso de recocido implica un control cuidadoso de la temperatura y tasa de enfriamiento.




Los pioneros en la aplicación del recocido simulado en problemas de optimización fueron Kirkpatrick, Gelatt y Vecchi en el año de 1983. Desde ese entonces ha sido estudiado extensamente. A diferencia de los métodos del descenso de gradiente y otros métodos de búsqueda deterministas  cuya principal desventaja es que tienden a estancarse en mínimos locales, la principal ventaja del recocido simulado es la habilidad de evitar el estancamiento en mínimos locales, de hecho, se ha probado que el recocido simulado  convergera al optimo global si existe suficiente aleatoriedad combinada con un enfriamiento lento.

Metafóricamente hablando, el proceso de recocido simulado es equivalente a arrojar pelotas sobre un acantilado, conforme estas bajan, rebotan y pierden energía a tal punto de llegar a una posición sin energía (mínimo local). Si a las pelotas se les permite rebotar suficientes veces con pérdidas de energía bajas, eventualmente estas caerán a posiciones globalmente bajas (mínimo global).


Modelo metropoli

     El fundamento del recocido simulado se basa en el trabajo de Metrópolis (1953) en el campo de la termodinámica estadística.

Básicamente, Metrópolis modeló el proceso de enfriamiento simulando los cambios energéticos en un sistema de partículas conforme decrece la temperatura, hasta que converge a un estado estable (congelado). Las leyes de la termodinámica dicen que a una temperatura t la probabilidad de un incremento energético de magnitud δE se puede aproximar por


siendo k una constante física denominada Boltzmann

En el modelo de Metrópolis, se genera una perturbación aleatoria en el sistema y se calculan los cambios de energía resultantes: si hay una caída energética, el cambio se acepta automáticamente; por el contrario, si se produce un incremento energético, el cambio será aceptado con una probabilidad indicada por la anterior expresión (Ec. 1).

Así, el recocido simulado puede decirse que es una versión iterada del modelo de Metropolis considerando k=1, donde el criterio que se toma en cuenta para considerar una solución de mayor costo o no es:


Donde:

𝛿: Es la diferencia de costes entre una solución vecina y la solución actual (Se ha tomado 𝛿𝐸 únicamente como 𝛿).
𝑇: Es el parámetro de temperatura.


Algoritmo Básico

   -Partiendo de una solución inicial, el algoritmo se ejecuta en varias iteraciones.
   
   -En cada iteración se genera un vecino aleatorio.
   
   -Movimientos que mejoran la función objetivo se aceptan siempre.

   -Si la solución vecina (candidato) es peor, entonces se selecciona una probabilidad dada (Ec. 2) que depende de la temperatura actual (𝑇) y de cuanto se degrade la función objetivo (𝛿).

   -Conforme avanza el algoritmo, la probabilidad decrece.

   -Se utiliza un parámetro de control llamado temperatura (𝑇), para poder determinar la aceptación de soluciones peores.

   -En cada nivel de temperatura, se explora un cierto numero de soluciones.

   -Cada vez que se ha explorado cierto numero de soluciones por nivel  de temperatura, se decrece la temperatura actual y se itera de nuevo.

   




Implementación del recocido simulado

Para implementar el recocido simulado se han definido 4 puntos los cuales van a depender del tipo de problema que se este resolviendo y la forma en el que este se quiera solucionar, a continuación se muestran dichos puntos:

• Representación
• Solución inicial
• Mecanismo de solución entre transiciones
• Secuencia de enfriamiento

     1. Representación

Depende del tipo de problema que se este resolviendo. A continuación se presentan algunos casos:

     -Vector ordenado de números enteros. Ejemplo: Problema del Viajante de Comercio (1 8 2 7 3 5 6 4)
     -Vector binario. Ejemplo: Problema de la Mochila (0 1 1 0 1 1 1 0)
     -Vector de Números Reales. Ejemplo: Problemas de Optimización con Parámetros Continuos (1.64 3.56 2.53)

Los casos anteriores son ejemplos de problemas clásicos de optimización. En los siguientes links se da mas información acerca de los problemas anteriormente enunciados.

No tienes permitido ver los links. Registrarse o Entrar a mi cuenta
No tienes permitido ver los links. Registrarse o Entrar a mi cuenta

     2. Solución inicial

     -Aleatoria: Se puede proponer una solución creada de forma aleatoria. Ej: vector permutado para el caso del problema del viajante de comercio.
     -Solución previa: Se puede obtener una solución de recocido simulado que se solución inicial de otro recocido simulado.
     -Solución obtenida con otras heurísticas: Se puede proponer como solución inicial una solución obtenida por otra heurística, como por ejemplo algoritmos genéticos o algoritmos voraces.

     3. Mecanismo de solución entre transiciones

     -Generación de una nueva solución: Se puede dividir en dos puntos:
      

           
• Definición del conjunto de vecinos: Es decir, obtener un conjunto soluciones candidatas a partir de la actual. Ej. Si se tiene un vector solución del problema del viajante de comercio So=[1,5,2,6,4,3], un conjunto de dos posibles soluciones candidatas podría ser Sc={[1,2,5,6,4,3],[1,5,3,6,4,2]}, se puede observar que las soluciones candidatas han sido obtenidas simplemente aplicando un intercambio entre dos elementos de la solución inicial del problema.

           
• Selección de un elemento de dicho conjunto: Se selecciona un elemento del conjunto de soluciones candidatas.

     -Cálculo de la diferencia de costos entre la solución actual y la vecina: Se obtiene la diferencia de costos entre la solución candidata seleccionada y la actual. Cabe destacar que el costo depende del problema que se este resolviendo, en el problema del viajante, el costo correspondería a la distancia entre ciudades, entre mas distancia mayor coste.

     -Aplicación del criterio de aceptación: Si la solución candidata es mejor (menor coste), se considera como la solución actual para la siguiente iteración, en el caso contrario, si la solución candidata es peor, se aplica el criterio dado por la ecuación 2, es decir, si generando un numero aleatorio uniforme entre 0 y 1 resulta ser menor que la probabilidad obtenida por la ecuación, se considera como solución actual a la candidata.

     4. Secuencia de enfriamiento

     -Valor inicial de temperatura: El valor de temperatura puede ser arbitrario, es decir, se pueden seleccionar temperaturas iniciales altas o bajas sin embargo una estrategia de selección de temperatura inicial puede ser el coste de la solución actual multiplicado por un parámetro entre 0 y 1 (usualmente se selecciona 0.4).

     -Mecanismo de enfriamiento: Existen varios mecanismos de enfriamiento:

• Enfriamiento basado en sucesivas temperaturas descendentes fijadas por el usuario.
• Lineal (i denota iteración) 𝑇=𝑇−𝛽, 𝑇_𝑖=𝑇_0−𝑖×𝛽
• Geométrica o exponencial 𝛼𝜖[0.8,0.99]   𝑇=𝛼𝑇
• Logaritmica o de Boltzman: Lenta pero existe prueba de convergencia global 𝑇(𝑡)=𝑇_0/ln⁡(1+𝑡)  t=1,2,---n
• Cauchy: 𝑇_𝑖=𝑇_0/(1+𝑖)
• Cauchy modificado utilizando balanceo o número fijo de iteraciones: 𝑇_(𝑖+1)=𝑇_𝑖/(1+𝛽𝑇_𝑖 ), donde 𝛽=𝑇_0−𝑇_𝐹/(𝐿−1) 𝑇_0 𝑇_𝐹

     -Criterio de parada: En teoría, el algoritmo debería finalizar cuando T=0. En la práctica, se para:

• Cuando T alcanza o está por debajo de un valor final 𝑇_𝑓 (0.01) , fijado previamente, o después de un número fijo de iteraciones.
• En ocasiones suele ser difícil dar valor de 𝑇_𝑓, se suele usar un número fijo de iteraciones, una buena opción es parar cuando no se haya aceptado ningún vecino de los L(T) generados en la iteración actual (num_éxitos=0). En ese caso, es muy probable que el algoritmo se haya estancado y no vaya a mejorar la solución obtenida.
  

A continuación se presenta la aplicación del recocido simulado para la solución de Sudokus. Nota: el código no ha sido optimizado, simplemente se ha realizado con el fin de aplicar dicha técnica:

import numpy as np
import statistics
class Sudoku_Solver:
    
    def __init__(self,So,To,alpha,L,Tf):
        self.So=So
        self.c_index=np.array(np.where(So==0))
        self.S=np.copy(self.So)
        self.To=To
        self.alpha=alpha
        self.L=L
        self.Tf=Tf
        self.T=To
    
    #Imprimir en consola el Sudoku
    def PrintSudoku(self,sudoku):
        print("\n")
        for i in range(len(sudoku)):
            line = ""
            if i == 3 or i == 6:
                print("-----------------------\n-----------------------")
            for j in range(len(sudoku[i])):
                if j == 3 or j == 6:
                    line += "|| "
                line += str(sudoku[i,j])+" "
            print(line)
            
    # Rellenar las casillas vacias del Sudoku        
    def fill_So(self):
        for k in range(3):
            for l in range(3):
                Aux=np.zeros((3,3))
                V=np.random.permutation(9)+1
                for i in range(0,3):
                    for j in range(0,3):
                        Aux[i][j]=self.S[i+k*3][j+l*3] 
                V=np.setdiff1d(V,Aux.reshape(9,))
                for i in range(len(V)):
                    e_index=np.where(Aux==0)
                    val=np.random.randint(0,len(e_index[1]))
                    Aux[e_index[0][val]][e_index[1][val]]=V[i]
                for i in range(0,3):
                    for j in range(3):
                        self.S[i+k*3][j+l*3]=Aux[i][j] 
        return self.S
    
    #Calcular el numero total de errores, en este caso la cantidad de digitos que se repiten por columna y fila
    def TotalError(self,sudoku):
        NErrors = 0 
        for i in range (0,9):
            NErrors += self.PartialError(i ,i ,sudoku)
        return(NErrors)

    #Calcular unicamente el error por fila y columna tomando como pivote a la casilla la cual se encuentre en dicha fila y columna
    def PartialError(self,row, column, sudoku):
        nErrors = (9 - len(np.unique(sudoku[:,column]))) + (9 - len(np.unique(sudoku[row,:])))
        return(nErrors)
    
    #Aplicar swap entre los digitos de dos casillas distintas 
    def Swap(self,S):
        Li=np.random.randint(0,3)
        Lj=np.random.randint(0,3)
        self.idx_iA=np.random.randint(Li*3,(Li+1)*3)
        self.idx_jA=np.random.randint(Lj*3,(Lj+1)*3)
        self.idx_iB=np.random.randint(Li*3,(Li+1)*3)
        self.idx_jB=np.random.randint(Lj*3,(Lj+1)*3)
        while (self.So[self.idx_iA][self.idx_jA]!=0) or (self.So[self.idx_iB][self.idx_jB]!=0):
            self.idx_iA=np.random.randint(Li*3,(Li+1)*3)
            self.idx_jA=np.random.randint(Lj*3,(Lj+1)*3)
            self.idx_iB=np.random.randint(Li*3,(Li+1)*3)
            self.idx_jB=np.random.randint(Lj*3,(Lj+1)*3)
        temp=S[self.idx_iA][self.idx_jA]
        S[self.idx_iA][self.idx_jA]=S[self.idx_iB][self.idx_jB]
        S[self.idx_iB][self.idx_jB]=temp
        return S
    
    #Aplicación del recocido simulado
    def sim_annealing(self):
         S=self.fill_So()
         best=np.copy(S)
         puntuacion=self.TotalError(best)
         while (self.TotalError(best)):
             Ts = []
             auxSudoku = best
             for i in range(1,10):
                auxSudoku = self.Swap(best)
                Ts.append(self.TotalError(auxSudoku))
             t=(statistics.pstdev(Ts))
             regulizer=0
             self.T=1000
             while self.T>self.Tf:
                 puntuacion_anterior=puntuacion
                 for i in range(self.L):
                     Scand=self.Swap(np.copy(S))
                     EScand=self.PartialError(self.idx_iA,self.idx_jA,Scand)+self.PartialError(self.idx_iB,self.idx_jB,Scand)
                     ES=self.PartialError(self.idx_iA,self.idx_jA,S)+self.PartialError(self.idx_iB,self.idx_jB,S)
                     delta=EScand-ES
                     if np.random.rand()<np.exp(-delta/t) or (delta<0):
                         S=Scand
                     puntuacion+=delta
                 t*=self.alpha
                 self.T=self.alpha*self.T
                 if puntuacion >= puntuacion_anterior: regulizer += 1
                 else: regulizer = 0
                 if (regulizer > 50): t += 2
                 if self.TotalError(S)<self.TotalError(best):
                     best=np.copy(S)
                     print('Best: '+str(self.TotalError(best)))
                 if self.TotalError(best)==0:
                     print(self.PrintSudoku(best))
                     break
                 
if __name__=='__main__':
    So=np.array([[0,0,0,0,9,0,1,0,5],
                 [0,0,8,0,5,4,3,2,7],
                 [0,0,5,0,1,0,0,8,6],
                 [0,8,3,7,0,6,0,0,9],
                 [5,0,6,0,3,1,0,0,0],
                 [1,2,0,0,8,0,4,0,0],
                 [6,3,1,0,0,9,0,5,2],
                 [2,7,0,5,0,0,0,0,0],
                 [8,0,0,0,7,0,0,0,0]])
    
    To=1000
    alpha=0.9
    L=50
    Tf=0.01
    sudoku_solver=Sudoku_Solver(So,To,alpha,L,Tf)
    sudoku_solver.sim_annealing()                    

Al correr el código se obtiene:



Que es la solución al sudoku introducido.


REFERENCIAS

[1] Yang, Xin-She, (2010). "Engineering optimization : an introduction with metaheuristic applications". Cambridge, United Kingdom. WILEY (347 pgs.)

Otras referencias: