0

I decided to learn simulated annealing as a new method to attack this problem with. It essentially asks how to fill a grid with -1, 0, or 1 so that each row and column sum is unique. As a test case, I used a 6x6 grid, for which there is definitely an optimal solution given by Neil:

1  1  1  1  1  1  6
1  1  1  1  1 -1  4
1  1  1  1 -1 -1  2
1  1  0 -1 -1 -1 -1
1  0 -1 -1 -1 -1 -3
0 -1 -1 -1 -1 -1 -5
5  3  1  0 -2 -4

My code usually doesn't reach the optimal case the majority of runs and even returns the wrong grid cost (old_cost should match count_conflict(grid)). Are my parameters incorrectly set, have I incorrectly implemented, or possibly is simulated annealing not a viable method here?

import random
from math import exp

G_SIZE = 6
grid = [[1]*G_SIZE for i in range(G_SIZE)]

def count_conflict(grid):
    cnt = [0]*(2*G_SIZE+1)
    conflicts = 0
    for row in grid:
        cnt[sum(row)] += 1
    for col in zip(*grid):
        cnt[sum(col)] += 1

    #print(cnt)
    for c in cnt:
        if c == 0: conflicts += 1
        if c > 1: conflicts += c-1
    return conflicts

def neighbor(grid):
    new_grid = grid[:]

    i = random.choice(range(G_SIZE))
    j = random.choice(range(G_SIZE))
    new_cells = [-1, 0, 1]
    new_cells.remove(new_grid[i][j])
    new_grid[i][j] = random.choice(new_cells)

    return new_grid

def acceptance_probability(old_cost, new_cost, T):
    if new_cost < old_cost: return 1.0
    return exp(-(new_cost - old_cost) / T)


# Initial guess
for i in range(1, G_SIZE):
    for j in range(0, i):
        grid[i][j] = -1

#print(grid)

old_cost = count_conflict(grid)
T = 10.0
T_min = 0.1
alpha = 0.99
while T > T_min:
    for i in range(1000):
        new_sol = neighbor(grid)
        new_cost = count_conflict(new_sol)
        ap = acceptance_probability(old_cost, new_cost, T)
        print(old_cost, new_cost, ap, T)
        if ap > random.random():
            grid = new_sol
            old_cost = new_cost

    T *= alpha

for row in grid:
    print(row)

print(count_conflict(grid))
Community
  • 1
  • 1
qwr
  • 9,525
  • 5
  • 58
  • 102
  • Somehow, I think the base problem should be posted on math question board, as it seems you can translate this problem into a group of linked equations with several variables, which transposes into matrices computation. – DainDwarf Jan 07 '16 at 07:24
  • 3
    I'm pretty sure SA doesn't provide any guarantees of finding a global minimum. The more times you roll the dice, the (possibly) more likely you are of finding it under some circumstances. Play with your cool-down profile. This is something that you have to do semi-manually, at least at first. Can you compare your routine to something known - why not run in parallel with `scipy.optimize.anneal` and compare behaviors. – uhoh Jan 07 '16 at 07:27
  • @uhoh I've fiddled with the parameters some, but the program usually stops just short of the intended goal. Since simulated annealing works with N queens problem, I hoped it would work here. – qwr Jan 07 '16 at 07:29
  • 2
    Be careful what you mean by "works". You have to look carefully at how it was run, number of starting cases, temperatures, etc. SA can *work* and *not-work* depending on the details. You could also try running 100 cases in parallel if you used a 3D NumPy array instead. Just a though. – uhoh Jan 07 '16 at 07:30
  • 2
    There is, in general, no algorithm in existence that can guarantee finding global optima. – David Heffernan Jan 07 '16 at 07:38
  • @DavidHeffernan I want maybe not global optima, but one of the best possible solutions (like that of backtracking) – qwr Jan 07 '16 at 07:41
  • Doesn't sound like you are clear on what you want at all. This is not a programming question. – David Heffernan Jan 07 '16 at 07:54

2 Answers2

1

a few things to do first, and which might quickly lead you to a working solution, without having to do anything else (eg, swap the heuristic):

  • add a line near the top outside of your main iterative loop, to calculate the cost of your t0 state (ie, your starting configuration);

  • inside the main loop, insert a single print statement just after the line that calculates the cost for the current iteration--that writes to a file, the value returned by the cost function for that iteration; just below that add a line that prints that value every 20 iterations or something like that (eg, about once each second is about as fast as we can comprehend scrolling data)

    if n % 10 == 0: print(what_cost_fn_returned_this_iteration)

  • don't call acceptance_probability; there is no natural convergence criterion in combinatorial optimization problems; the usual practice is to break out of the main loop when any of these happen:

    the max iteration count has been reached

    the current minimum value of the cost function over the past __ iterations has changed less than __%; for instance if over the last 100 iterations, the cost (by comparing a min and max using a moving window) varies less than 1%

    after reaching a minimum during iteration, the cost is now consistently increasing with iteration count

A few other observations:

  • with the diagnostics in place (see above) you will be able to determine: (i) from some initial cost, what is my solver doing? ie, is it moving in a more-or-less direct path to lower-and-lower values? Is it oscillating? Is it increasing? (If the latter, the fix is usually that you have a sign backwards)

  • a 6 x 6 matrix is very very small--that doesn't leave a lot for the
    cost function to work with

  • re-write your cost function so that a "perfect" solution returns a
    zero cost, and all others have a higher value

  • 1000 iterations is not a lot; try increasing that to 50,000

doug
  • 69,080
  • 24
  • 165
  • 199
0

new_grid = grid[:] makes a shallow copy. A deep copy or modifying the grid in place and reverting to the original solves the issue.

qwr
  • 9,525
  • 5
  • 58
  • 102