Simulated annealing-n queen problem

Simulated annealing  

Here is just a quick understanding , Used to solve common problems .

Take a brief look at simulated annealing ：

Simulated annealing is a simulated physical annealing method , adopt N Sub iteration （ Anneal ）, A on an approximation function The most value （ Maximum or minimum ）.

Simulated annealing feasibility ：

Let's imagine a solid in there , If the temperature of the solid is very high at this time , Use our high school knowledge , At this time, the internal energy of the object is large , Because the kinetic energy of the particles moving irregularly inside is large , Molecular heat moves faster . At this point as a very disorder The state of , As it cools down , As the temperature drops to normal room temperature , The particles inside it slowly ‘ Orderly ’ , At this time, it is Stable .

Pay attention to the content in the text ：

1. The temperature is high -- Fast movement , Low temperature -- Slow motion .
2. The temperature drops slowly .
3. When the temperature returns to normal room temperature, it tends to be orderly ( Near optimal solution ).

therefore , The general direction of simulated annealing is ：

Simulated annealing is a kind of Loop algorithm .

1. Let's set an initial temperature T（ The temperature should be higher , such as 2000）
2. Every cycle Annealing once .（ How to operate specifically I'll explain later ）
3. Then it goes down T Temperature of , Let's get through T And a “ Cooling coefficient ”ΔT（ A close approach 1 Decimals of , such as 0.990.99) Multiply , achieve Lower slowly The effect of temperature , Until it's close to 0（ We use it eps To represent a close 0 Number of numbers ( such as 0.00001), as long as T<eps You can exit the loop ）

The simple code is as follows ：

T = 2000  #  Represents the starting temperature 
dT = 0.99  #  Is the coefficient delta T
eps = 1e-14  #  amount to 0.0000000000000001
while T > eps:
    '''
     Here is the operation of each annealing 
    '''
    T = T * dT  #  The temperature drops a little at a time , T * 0.99

Take the following function as an example ：

Such as   

The answers we ask for are nothing more than two ： The independent variables xx And the corresponding function Maximum f(x)

1. Find some at random x0( Define random in domain ), And the corresponding f(x0).

 

2.  Start annealing

As mentioned above x0 amount to A particle , So we will have a Disordered movement , That is, move randomly to the left or right .

But remember one key point ： The mobile Range and The current temperature T of .

temperature T The bigger it is , The more you move . temperature T The smaller it is , The less you move . This is simulating the disordered motion of particles .

 

3. Accept better

In the picture above, we moved to x2 It's about , obvious f(x2) be better than f(x0).

So we will answer the question to update : x0=x2,f(x0)=f(x2), It's also a greedy Thought .

 

4. With certainty probability Accept Worse state

This is also the best part of annealing .

Why should we accept a worse state ？ Because there may be a... Next to a poor state Higher The mountains of .

If we only pay attention to x0 The first random point area ( Image left ), With A drop in temperature 、 Decrease of left-right jump range , We will end up stuck here , Only the locally optimal answer is obtained .

 

And if we find the low point of the mountain on the right , With A certain probability accepts it （ The probability is related to the temperature and the criticality of the current value ）, Before the jump amplitude decreases , Find the best possible .

So how many probability do we accept it ？ We use a formula to express （ We just need to remember this formula , This is derived by scientists ）：

Explain the meaning of each part separately ： 

 ： Because we use functions To represent a probability value , So we just need to focus on x by negative The part of ： 

The value range of the negative part is in (0,1) In the open section ,x The smaller it is , The closer the 0, The bigger the closer 1.

Because in 0 To 1 Between , So this value is equivalent to the probability . such as =0.97, So the probability we accept is 97%,

The value range of the positive part will be greater than 1, In other words, the probability will exceed 100%, So I will definitely choose （ In fact, it is the last way to find a better situation ）

k： It's actually a physical constant , We're in code It won't use .

T： It's simple , It's the current temperature . So actually the denominator is T,k treat as 1 Use .

Δf： Let's focus on what is Δf.

Actually, from the previous function Can be found in ,Δf Must be a negative number .

Δf Is the difference between the function value of the current solution and the function value of the objective solution ,Δf=−|f(x0)−f(x1)|, And it must be a negative number .

For example, now let's find the of a function Maximum , So if f(x0)<f(x1) 了 , That means the result It's getting better ( Specific analysis of specific problems ), We sure Choose it （ See the first 3 spot ： Accept better ）.

If f(x0)>f(x1), That means the result It's getting worse , We need probability to choose it , therefore Δf=−(f(x0)−f(x1)).

 

We can use    And (0,1) Compare random numbers generated between , If it is greater than, it will be accepted; otherwise, it will not be accepted

 

ps: It is not difficult to see that the higher the temperature , The greater the probability of acceptance .Δf The smaller it is , Explain that the worse the answer (f(x1) And f(x2) There's a lot of difference , But also accept the probability ), Then the probability of acceptance is smaller .

 

So the conclusion is ：

• If the result of the random function is better , We must choose it ( namely x0=x1,f(x0)=f(x1))
• If the result of the random function is worse , We use The probability of accepting it

 

Example 1： By simulated annealing Value

We use functions f(x) in x Beat to find the result . Might as well set ：

import math
import random


# x Represents the sum of the squares of the numbers we randomly generate n The closeness of 
def func(x, n):
    return abs(x * x - n)


#  Simulated annealing 
def SA(T, dT, eps, n):
    #  Radical sign 10 about 3.1622  Here, the scope is artificially limited   Can also be directly random   Anyway, it will bounce back and forth under the influence of temperature and is random 
    x = random.uniform(1.6, 5.4)
    #  Work out x Sum of squares n The gap between f(x0)
    f = func(x, n)
    while T > eps:
        #  The number under the root sign must be nonnegative   So we can directly random nonnegative numbers 
        new_x = random.randint(-13, 20) * T #  The temperature effect bounces back and forth and is random   The lower range is arbitrarily limited here 
        new_f = func(new_x, n)
        if new_f < f:
            x = new_x  #  Replace   If there are multiple variables, replace multiple 
            f = new_f
        elif math.exp((f - new_f) / T) > random.random():  #  Probability acceptance  ex stay (0,1) Between   actually random.random() stay [0,1)
            x = new_x
            f = new_f
        #  cool 
        T *= dT
    #  Output results   Approximate value   If I am  3.1201885471764994
    print(x)


if __name__ == '__main__':
    n = 10  # n Represents the value to be approximated by our last function 
    T = 20000  #  initial temperature , The initial temperature mainly depends on the range of jump required at the beginning of the topic .
    dT = 0.993  #  Rate of change , It needs to be a little slower here , Write 0.995  perhaps  0.997 Fine , But the closer you get 1 The slower the speed 
    eps = 1e-14  # 10 Of -14 The power is already very small , It's no problem to write this probability 

    SA(T, dT, eps, n)

Example 2： The sum of point distances is the smallest

On the given plane N（N<=100） A little bit , You need to find a point like this , Make this point to N The sum of distances of points shall be as small as possible . Output the minimum distance and （ Round to the nearest whole number ）.

Input Input

first line N, Next N Two integers per line , Express N A little bit

4
0 0
0 10000
10000 10000
10000 0

Output Output

One positive integer per line , Minimum distance and .

28284

Their function func() Namely ： Find out one Random point A=(x0,y0) To NN The sum of distances between two points DD, That's the point A And all N Add the distances of points , namely ：

n queens problem

In short, it's a n*n On the chessboard of , There is only one chess piece in each line , And the column where the chess piece is located 、 Left slash 、 There are no other pieces in the right slash .

import copy
import random
import time
import math
import numpy as np


def init():
    cache = {}
    m = np.zeros((8, 8), dtype=int)
    for i in range(0, 8):
        temp = random.randrange(0, 8)
        m[temp, i] = 1
        cache["queen" + str(i)] = [temp, i]
    return m, cache


#  Calculate the number of no collisions in the current state 
def compute_weight(coord_cache):
    weight = 0
    for i in range(0, 8):
        x, y = coord_cache["queen" + str(i)]
        for j in range(i + 1, 8):
            x2, y2 = coord_cache["queen" + str(j)]
            if x2 - x == j - i or x2 - x == i - j:
                weight += 1
            if x2 == x:
                weight += 1
    return 28 - weight


#  Generate a new solution randomly 
def random_adjust(coord_cache):
    temp = copy.deepcopy(coord_cache)
    row = random.randrange(0, 8)
    column = random.randrange(0, 8)
    temp["queen" + str(column)] = [row, column]  #  Adjust the Queen's position 
    return temp


def draw(coord_cache):
    m = np.zeros((8, 8), dtype=int)
    for i in range(8):
        row, column = coord_cache["queen" + str(i)]
        row, column = int(row), int(column)
        m[row][column] = 1
    return m


def cool(T, eps, dt, L):
    m, coord_cache = init()
    print(" The initialization status of the eight queens is ：\n", m)
    while T > eps:  #  Temperature cycling 
        for i in range(L):  #  Each temperature cycle L Time 
            weight = compute_weight(coord_cache)
            print(" The collision free degree of the current state is ：", weight)
            if weight == 28:  #  The non collision degree is 28, It shows that the optimal solution is found 
                return True
            new_coord_cache = random_adjust(coord_cache)  #  A new solution is obtained by random adjustment 
            new_weight = compute_weight(new_coord_cache)  #  Calculate the collision degree of the new solution 
            print(" The new solution produced by random adjustment is ：\n", draw(new_coord_cache))
            print(" The collision free degree of the new solution generated by random adjustment is ：", new_weight)
            if new_weight >= weight:  #  If the collision degree of the new solution is smaller, this solution is accepted 
                coord_cache = new_coord_cache
                print(" This is a better solution , Direct reception ：\n", draw(coord_cache))
            else:
                if random.random() < math.exp((new_weight - weight) / T):  #  Otherwise, the probability of simulated annealing is accepted as a new solution 
                    coord_cache = new_coord_cache
                    print(" The current receiving probability is ：", math.exp((new_weight - weight) / T))
                    print(" This is a worse solution , But it was accepted ：\n", draw(coord_cache))

        T = T * dt


def Cool(T, eps, dt, L, num):
    t1 = time.time()
    success = 0
    fail = 0
    for i in range(num):
        if cool(T, eps, dt, L):
            success += 1
            print(" The first {0} Examples to find the optimal solution ".format(i))
        else:
            fail += 1
            print(" The first {0} Examples failed ".format(i))
    t2 = time.time()
    print("{} The successful solution of the three examples is ：{}".format(num, success))
    print("{} The percentage of cases successfully resolved is ：{}".format(num, success / num))
    print("{} The example of failure among the examples is ：{}".format(num, fail))
    print("{} The percentage of examples that failed is ：{}".format(num, fail / num))
    print("{} The time required to run the algorithm in the three examples is ：{} second ".format(num, t2 - t1))


Cool(5, 0.001, 0.98, 150, 1000)

n Queen other solutions ：

An operation ( Only one feasible solution is output here )：

def draw(n, ini):
    for i in ini:
        print('+---' * n + '+', end='')
        print('')
        print('| ' + ' | '.join(['#' if j == '0' else 'Q' for j in list(i[2:].rjust(n, '0'))]) + ' | ')
    print('+---' * n + '+', end='')


def DFS(row, colomn, left, right):
    global ini
    global result
    global flage
    bits = ((1 << n) - 1) & ~(colomn | left | right)  #  Available columns for the current row 
    if bits == 0 and len(ini) != 0:
        ini.pop()
    while bits:
        pos = bits & -bits  #  Get the position of the row 
        bits = bits & (bits - 1)  # bits ^= pos
        if row == n - 1:
            result += 1
            if flage:
                ini.append(bin(pos))
                draw(n,ini)
                flage-=1
        else:
            if flage:
                ini.append(bin(pos))
            DFS(row + 1, colomn | pos, (left | pos) >> 1, (right | pos) << 1)
            if bits == 0 and len(ini) != 0:
                ini.pop()

n = 6
result = 0
flage = 1  #  To preserve a single solution 
ini = []
DFS(0, 0, 0, 0)
print('\n','{} In this case '.format(result))

Another mathematical idea ：

#  Judgment point [row, col] Whether the queen can be placed 
def check(matrix, row, col):
    for i in range(row):
        if abs(matrix[i] - col) == 0 or abs(matrix[i] - col) == abs(row - i):
            return False
    return True


def queen(matrix, row):
    global result
    n = len(matrix)
    if row == n:
        for col in matrix:
            print(' # ' * col + ' Q ' + ' # ' * (n - (col + 1)))
        print('')
        result += 1
    for col in range(n):
        if check(matrix, row, col):
            matrix[row] = col
            queen(matrix, row + 1)


n = 10
result = 0
matrix = [0] * n
queen(matrix, 0)
print(' share {} Seed result '.format(result))