ILP Solving in CP-SAT

cp-sat.py

from ortools.sat.python import cp_model
import math
# Custom library for visualisation. Can be omitted.
from visualiser import Visualiser

NUM_CELLS = 25
CELLS_PER_SIDE  = int(math.sqrt(NUM_CELLS))
START_CELL = 5
TARGET_CELL = 23

# Weights matrix has infinity edge cost between cells to begin with
W=[]
for i in range(NUM_CELLS):
    list=[]
    for j in range(NUM_CELLS):
        list.append(1000)
    W.append(list)

# Name the cells like below (example for 5x5 grid):
#  0  1  2  3  4
#  5  6  7  8  9
# 10 11 12 13 14
# 15 16 17 18 19
# 20 21 22 23 24
# Set finite edge weight for adjacent cells
# identifies corner cases where consecutively identified cells are not adjacent (like 9&10)
for i in range(NUM_CELLS):
    if (i+1)%CELLS_PER_SIDE != 0:
        W[i][i+1] = 1 
    if i+CELLS_PER_SIDE < NUM_CELLS:
        W[i][i+CELLS_PER_SIDE] = 1
    if i%CELLS_PER_SIDE != 0:
        W[i][i-1] = 1 
    if i-CELLS_PER_SIDE >= 0:
        W[i][i-CELLS_PER_SIDE] = 1

for i in W:
    print(i)

# Source goes downwards, destination goes rightwards
# Matrix of decision variables that represent transition between cells
P=[]

# [CREATE MODEL]-------------------------------------------------------------------------------------------------------
model = cp_model.CpModel()

# [START VARIABLE]-----------------------------------------------------------------------------------------------------
# create p matrix
for i in range(NUM_CELLS):
    list=[]
    for j in range(NUM_CELLS):
        list.append(model.NewBoolVar(str(i)+','+str(j)))
    P.append(list)

# [START CONSTRAINTS]-------------------------------------------------------------------------------------------------
# ________CONSTRAINT 1________
# Flow constraints (relationship between in and out constraints)
for i in range(NUM_CELLS):
    out_edges=P[i] #The ith row. Sum gives out degree
    in_edges=[P[j][i] for j in range(NUM_CELLS)] #the jth column. Sum is in degree
    if i is START_CELL:
        # out_degree-in_degree=1 for start_cell
        model.Add(sum(out_edges) - sum(in_edges) == 1)
    elif i is TARGET_CELL:
        # out_degree-in_degree=-1 for sink_cell
        model.Add(sum(out_edges) - sum(in_edges) == -1)
    else:
        # out_degree=in_degree for all other cells
        model.Add(sum(out_edges) - sum(in_edges) == 0)

# ________CONSTRAINT 2________
# Set the limit on in degree
for i in range(NUM_CELLS):
    in_edges=[P[j][i] for j in range(NUM_CELLS)]
    if i is START_CELL:
        # in degree should be = 0 for start_cell
        model.Add(sum(in_edges) == 0)
    elif i is TARGET_CELL:
        # in degree should be = 1 for target_cell
        model.Add(sum(in_edges) == 1) #must visit cell
    else:
        # in degree should be 0 or 1 otherwise
        model.Add(sum(in_edges)<=1)

# [START OBJECTIVE]-------------------------------------------------------------------------------------------------------
# The objective function states that the cost of the entire travel should be minimum
# I think summation of P[i][j]*W[i][j] over)) for i in range(NUM_CELLS) for j in range(NUM_CELLS) all i and all j should be minimum
model.Minimize(sum(P[i][j] * W[i][j] for i in range(NUM_CELLS) for j in range(NUM_CELLS)))

# [START SOLVING]-------------------------------------------------------------------------------------------------------
solver = cp_model.CpSolver()
solver.Solve(model)
print(solver.ResponseStats())

# [SHOW SOLUTION]-------------------------------------------------------------------------------------------------------
for i in range(NUM_CELLS):
    for j in range(NUM_CELLS):
        print(int(solver.BooleanValue(P[i][j])), end=" ")
    print()

output_matrix=[]
for i in range(NUM_CELLS):
    list=[]
    for j in range(NUM_CELLS):
        list.append(solver.BooleanValue(P[i][j]))
        if solver.BooleanValue(P[i][j]):
            print(i, "->", j)
    output_matrix.append(list)

# Visualise the path. Can be omitted
visualiser = Visualiser()
visualiser.visualise(output_matrix)