How to linearise binary constraints

Question

I'm trying to find a valid previous state in the game of life using linear programming.

This is what I have so far:

board = [[0, 1, 0],
         [0, 1, 0],
         [0, 1, 0]]

width = len(board)
height = len(board[0])

rows = [str(i) for i in range(width)]
cols = [str(i) for i in range(height)]

choices = LpVariable.dicts("is_alive", (rows, cols), 0, 1, cat="Binary")
or_variables = LpVariable.dicts("OR_variable", (rows, cols), 0, 1, cat="Binary")

prob = LpProblem("Game_of_Life_reversal", LpMinimize)
prob += 0, "Arbitary Objective Function"

for r in rows:
    for c in cols:
        cell = choices[r][c]
        neighbours = getNeighbours(r,c, width, height, choices)
        if board[int(r)][int(c)] == 1:
            #If it's alive, and it was previously alive, then it must've had 2 or 3 neighbours
            #If it's alive, and it was previously dead, then it must've had 3 neighbours
            prob += lpSum(neighbours) >= 2 * cell + (1 - cell) * 3
            prob += lpSum(neigbours) <= 3
        else:
            #If it's dead, and was previously alive, then it must've had 0, 1 or 4+ neighbours
            #If it's dead, and was previously dead, then it must've had 0, 1, 2 or 4+ neighbours
            or_var = or_variables[r][c]
            prob += lpSum(neighbours) >= 4 * (1 - or_var)
            prob += lpSum(neighbours) <= cell * or_var + 9 * (1 - or_var) + or_var * (1 - cell) * 2

The last constraint prob += lpSum(neighbours) <= or_var + 9 * (1 - or_var) + or_var * (1 - cell) * 2 is throwing TypeError: Non-constant expressions cannot be multiplied.

I'm guessing this is because I'm doing or_var * (1 - cell) * 2 which violates linearity.

Is there a workaround?

Answer 1

Consider that we have a grid n x n with. Our variables are 0 <= C(i , j) <= 1 for 1 <= i <= n and 1 <= j <= n.

Let C*(i , j) = sum of the variables associated to the neighbors of the cell (i , j).

For example, if the cell (i , j) is not a border:

C*(i, j) = C(i+1, j) + C(i-1, j) + C(i, j+1) + C(i, j-1) + C(i-1, j-1) + C(i-1, j+1) + C(i+1, j-1) + C(i+1, j+1).

If it's the upper left corner:

C*(1, 1) = C(1 , 2) + C(2 , 2) + C(2 , 1)

And so on.

---

1 ) For each cell (i, j) not in the border of the grid:

1.1) If the cell is now occupied:

3 - C(i, j) <= C*(i , j) <= 3

1.2) If the cell is now not occupied:

Let 0 <= x(i , j) <= 1,

4 * x(i , j) <= C*(i , j) <= 8 * x(i , j) + 2 - C(i , j)

---

2 ) For each cell (i , j) in the border of the grid:

2.1) If the cell is a corner:

2.1.1) If the cell is now occupied:

Use the same constraint of the item (1.1).

2.1.2) If the cell is now not occupied:

0 <= C*(i , j) <= 2 - C(i , j)

---

2.2) If the cell is an edge:

Same constraints of item (1).

---

Now you got n² contraints, plus at most n² - 4 contraints for the possibly defined x(i , j), plus n² contraints of the form 0 <= C(i , j) <= 1. All the contraints are linear.