The $n$-Queens Problem: Constraint Satisfaction and Optimization

Introduction

First posed anonymously in 1848 in a German chess magazine, the $n$-queens problem has been a cornerstone of computational theory, later attributed to Max Bezzel. It remains a staple example in the study of Artificial Intelligence, particularly in the domain of Constraint Programming.

The problem involves placing $n$ queens on an $n \times n$ chessboard such that no two queens threaten each other—meaning no two queens share the same row, column, or diagonal. Below is a visual representation of one solution to the 8-queens problem, illustrating the complexity and intrigue of the puzzle:

A Dynamic Programming Approach

Dynamic programming offers a potent strategy for tackling the $n$-queens problem through a depth-first search algorithm that utilizes recursion and backtracking. This method systematically explores feasible placements row by row, backtracking when a placement leads to a dead end. Below is a Python implementation of this approach, which is capable of discovering all 92 distinct solutions for the 8-queens problem near-instantaneously:

 1from tabulate import tabulate
 2
 3def solveNQueens(n_queens: int) -> list[list[str]]:
 4    assert n_queens >= 4
 5
 6    board = [[" "] * n_queens for _ in range(n_queens)]
 7    result = []
 8
 9    # Infeasible locations
10    infeasible_cols = set()
11    infeasible_diag_diff = set()
12    infeasible_diag_sum = set()
13
14    def dfs(row):
15        # A solution is valid once n queens are placed
16        if row == n_queens:
17            result.append(["".join(r) for r in board])
18            return
19
20        for column in range(n_queens):
21            # Row-column-difference to identify "\"-diagonals that
22            # are running from the top left to the bottom right.
23            # E.g.: diff =  0 => the diagonal is (0,0), (1,1), ...
24            #       diff = -1 => the diagonal is (0,1), (1,2), ...
25            #       diff =  1 => the diagonal is (1,0), (2,1), ...
26            #       ...
27            diag_diff = row - column
28            # Row-column-sum to identify "/"-diagonals that
29            # are running from the top right to the bottom left.
30            # E.g.:  sum =  7 => the diagonal is (0,7), (1,6), ...
31            #        sum =  6 => the diagonal is (0,6), (1,5), ...
32            #        ...
33            #        sum =  0 => the diagonal is (0,0)
34            diag_sum = row + column
35
36            # If the current column is valid
37            if not (
38                column in infeasible_cols
39                or diag_diff in infeasible_diag_diff
40                or diag_sum in infeasible_diag_sum
41            ):
42                # Update the infeasible locations
43                infeasible_cols.add(column)
44                infeasible_diag_diff.add(diag_diff)
45                infeasible_diag_sum.add(diag_sum)
46
47                # Place the queen
48                board[row][column] = "X"
49                dfs(row + 1)
50
51                # Backtrack
52                infeasible_cols.remove(column)
53                infeasible_diag_diff.remove(diag_diff)
54                infeasible_diag_sum.remove(diag_sum)
55                board[row][column] = " "
56
57    dfs(row=0)
58    return result
59
60
61if __name__ == "__main__":
62    result = solveNQueens(8)
63    print(tabulate(result[0], tablefmt="grid"))

While the classic 8-queens problem can be solved fairly quickly, as we scale up the problem size, the problem quickly becomes intractable due to what’s known as a combinatorial explosion—essentially, the number of possible configurations skyrockets, making a straightforward approach impractical. In fact, even when ignoring diagonal attacks, there are still $n \times (n-1) \times \ldots \times 1 = n!$ possible arrangements ($n$ in the first row, $n-1$ in the second row, …, 1 in the last row).

To illustrate, here’s how the number of solutions and the computation time explodes with increasing numbers of queens in our Python implementation:

This difficulty persists even when employing specialized constraint programming solvers, like Google’s CP-SAT. However, if the objective shifts from finding every single feasible solution to merely identifying one valid placement, the situation changes drastically. With this more focused aim, we can typically determine a feasible solution in less than 0.001 seconds for moderately sized problems (such as the ones in the table).

The $n$-Queens Problem as an Optimization Problem

Expanding on the $n$-queens problem, consider each square on the board having a cost $c_{ij}$, and the goal is to minimize the total cost of placing the $n$ queens. This variation transforms a relatively simple Constraint Satisfaction Problem into a more complex Optimization Problem.

A trivial (and highly ineffective) approach of solving this optimziation problem consists in enumerating all feasible solutions (e.g., generated by the Python code above), and comparing their cost. Alternatively, the problem can be formulated as an Integer Program, which can be efficiently solved by using a Branch-and-Bound algorithm. For the cost matrix given above, the cost-minimal placement of 8 queens is the following.

An Integer Programming Formulation

To formalize the Integer Programming model, let us use the following binary variables:

The objective function and constraints are formulated as follows:

The objective function minimizes the cost associated with each placement. The first constraint accounts for the number of queens and the next two sets of constraints restrict one queen to be placed in each column and row, respectively. The next four sets of constraints account for all diagonals, located:

1. in the first row and running to the bottom right.
2. in the first column (except the first row) and running to the bottom right.
3. in the first row and running to the bottom left.
4. in the last column (except the first row) and running to the bottom left.

Finally, the domain of the decision variables is restricted to be binary. Clearly, if, e.g., $c_{ij} = 1 , \forall i, j \in 1, 2, \ldots, n$, then the problem corresponds to the standard $n$-queens problem.

Implementing the Model with Gurobi Optimizer

The Integer Programming model can be implemented (e.g., in Python) and solved swiftly by any standard solver, e.g., Gurobi Optimizer.

 1import gurobipy as gp
 2from gurobipy import GRB
 3import numpy as np
 4from tabulate import tabulate
 5
 6n = 8
 7C = np.random.randint(1, 100, (n, n))
 8
 9# Build the model
10m = gp.Model()
11x = m.addVars(n, n, vtype=GRB.BINARY, name="x", obj=C)
12
13m.addConstr(gp.quicksum(x[i, j] for i in range(n) for j in range(n)) == n, name="Queens")
14m.addConstrs((gp.quicksum(x[i, j] for j in range(n)) <= 1 for i in range(n)), name="Row")
15m.addConstrs((gp.quicksum(x[i, j] for i in range(n)) <= 1 for j in range(n)), name="Column")
16
17# First row \-diagonal
18for j in range(n):
19    expr = gp.LinExpr(x[0, j])
20    for k in range(1, n - j):
21        expr.add(x[0 + k, j + k])
22    m.addConstr(expr <= 1, f"\-Diagonal[0][{j}]")
23
24# First column \-diagonal (except first row)
25for i in range(1, n):
26    expr = gp.LinExpr(x[i, 0])
27    for k in range(1, n - i):
28        expr.add(x[i + k, 0 + k])
29    m.addConstr(expr <= 1, f"\-Diagonal[{i}][0]")
30
31# First row /-diagonal
32for j in range(n):
33    expr = gp.LinExpr(x[0, j])
34    for k in range(1, j + 1):
35        expr.add(x[0 + k, j - k])
36    m.addConstr(expr <= 1, f"/-Diagonal[0][{j}]")
37
38# Last column /-diagonal (except first row)
39for i in range(1, n):
40    expr = gp.LinExpr(x[i, n - 1])
41    for k in range(1, n - i):
42        expr.add(x[i + k, n - 1 - k])
43    m.addConstr(expr <= 1, f"/-Diagonal[{i}][{n-1}]")
44
45m.write("model.lp")
46m.optimize()
47
48# Print the solution
49result = np.zeros((n, n))
50for i in range(n):
51    for j in range(n):
52        if x[i, j].X > 0.5:
53            result[i][j] = C[i][j]
54
55print(tabulate(result, tablefmt="grid"))

To demonstrate with an example where $n=4$, the model creates a file named model.lp to be optimized. It’s straightforward to verify that this model accurately represents our unique take on the $n$-Queens problem.

\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
  x[0,0] + 8 x[0,1] + 2 x[0,2] + 9 x[0,3] + 5 x[1,0] + x[1,1] + 9 x[1,2]
   + 2 x[1,3] + 6 x[2,0] + x[2,1] + 7 x[2,2] + 4 x[2,3] + 3 x[3,0]
   + 9 x[3,1] + 8 x[3,2] + 2 x[3,3]
Subject To
 Queens: x[0,0] + x[0,1] + x[0,2] + x[0,3] + x[1,0] + x[1,1] + x[1,2] + x[1,3]
   + x[2,0] + x[2,1] + x[2,2] + x[2,3] + x[3,0] + x[3,1] + x[3,2] + x[3,3]
   = 4
 Row[0]: x[0,0] + x[0,1] + x[0,2] + x[0,3] <= 1
 Row[1]: x[1,0] + x[1,1] + x[1,2] + x[1,3] <= 1
 Row[2]: x[2,0] + x[2,1] + x[2,2] + x[2,3] <= 1
 Row[3]: x[3,0] + x[3,1] + x[3,2] + x[3,3] <= 1
 Column[0]: x[0,0] + x[1,0] + x[2,0] + x[3,0] <= 1
 Column[1]: x[0,1] + x[1,1] + x[2,1] + x[3,1] <= 1
 Column[2]: x[0,2] + x[1,2] + x[2,2] + x[3,2] <= 1
 Column[3]: x[0,3] + x[1,3] + x[2,3] + x[3,3] <= 1
 \-Diagonal[0][0]: x[0,0] + x[1,1] + x[2,2] + x[3,3] <= 1
 \-Diagonal[0][1]: x[0,1] + x[1,2] + x[2,3] <= 1
 \-Diagonal[0][2]: x[0,2] + x[1,3] <= 1
 \-Diagonal[0][3]: x[0,3] <= 1
 \-Diagonal[1][0]: x[1,0] + x[2,1] + x[3,2] <= 1
 \-Diagonal[2][0]: x[2,0] + x[3,1] <= 1
 \-Diagonal[3][0]: x[3,0] <= 1
 /-Diagonal[0][0]: x[0,0] <= 1
 /-Diagonal[0][1]: x[0,1] + x[1,0] <= 1
 /-Diagonal[0][2]: x[0,2] + x[1,1] + x[2,0] <= 1
 /-Diagonal[0][3]: x[0,3] + x[1,2] + x[2,1] + x[3,0] <= 1
 /-Diagonal[1][3]: x[1,3] + x[2,2] + x[3,1] <= 1
 /-Diagonal[2][3]: x[2,3] + x[3,2] <= 1
 /-Diagonal[3][3]: x[3,3] <= 1
Bounds
Binaries
 x[0,0] x[0,1] x[0,2] x[0,3] x[1,0] x[1,1] x[1,2] x[1,3] x[2,0] x[2,1]
 x[2,2] x[2,3] x[3,0] x[3,1] x[3,2] x[3,3]
End

Conclusion

While the classic $n$-queens puzzle is a staple in the study of combinatorial problems, the variations on this theme are virtually limitless. One intriguing twist on this classic puzzle involves not just assigning a cost to each square on the chessboard, but also allowing the rows of the board to be permuted. This variation elevates the challenge from a standard Integer Programming problem to the more complex realm of Quadratic Programming! ♟️