constraintsor-toolspulpinteger-programmingcp-sat

How to use CP-SAT solver or other tools to find 3D arrays with constraints on rows, columns and tubes (representing school classes, terms and days)?


I will be grateful to anyone who can help me write some Python code to enumerate the 21×2×3 arrays, indexed with i, j and k, which are two-thirds filled with 0's and one-third filled with the values 'Ava', 'Bob', 'Joe', 'Mia', 'Sam', 'Tom', 'Zoe' in such a way that:

  1. fixed the index i you have exactly two empty 2-tuples and one 2-tuple with different non-zero values;

  2. fixed the index k you have exactly fourteen empty 2-tuples and seven 2-tuple with different non-zero values;

  3. fixed the indexes j and k you have a 21-tuple with fourteen zero values and exactly one occurrence of each of the non-zero values, respecting the following constraints:

    a) "Ava" can appear only in a row with index 0, 1, 4, 6, 10, 11, 13, 14, 15, 19 or 20;

    b) "Bob" can appear only in a row with index 2, 3, 5, 7, 8, 9, 12, 16, 17 or 18;

    c) "Joe" can appear only in a row with index 2, 4, 5, 7, 8, 10, 14, 15, 18 or 20;

    d) "Mia" can appear only in a row with index 0, 1, 3, 6, 9, 12, 13, 16, 17 or 19;

    e) "Sam" can appear only in a row with index 1, 2, 7, 9, 15, 17 or 20;

    f) "Tom" can appear only in a row with index 0, 3, 8, 10, 12, 16 or 19;

    g) "Zoe" can appear only in a row with index 4, 5, 6, 11, 13, 14 or 18.

As a result I would like to obtain something like this:

[ 0   0       [Tom Mia      [ 0   0
  0   0        Ava Sam        0   0
  0   0        Sam Bob        0   0
  0   0        Bob Tom        0   0
  0   0         0   0        Joe Zoe
  0   0        Joe Zoe        0   0
  0   0         0   0        Zoe Ava
 Joe Sam        0   0         0   0
  0   0         0   0        Tom Bob
  0   0         0   0        Mia Sam
 Tom Ava        0   0         0   0
 Ava Zoe        0   0         0   0
 Bob Mia        0   0         0   0
  0   0        Mia Ava        0   0
  0   0        Zoe Joe        0   0
 Sam Joe        0   0         0   0
  0   0         0   0        Bob Tom
  0   0         0   0        Sam Mia
 Zoe Bob        0   0         0   0
 Mia Tom        0   0         0   0
  0   0 ]       0   0 ]      Ava Joe]

Rows represent school classes, columns represent school terms (there are 2 of them), tubes represent class days (there are 3 of them: Monday, Wednesday and Friday). So the first horizontal slice of the above solution means that class 1A has lesson only on Wednesday, in the the first term with teacher Tom and in the second term with teacher Mia. (Teachers can only work in some classes and not in others.)

Thanks in advance!


Update n. 1

As a starting point, I tried to attack the following toy problem:
Enumerate all arrays with a given number of rows and 3 columns which are two-thirds filled with "0" and one-third filled with "1" in such a way that summing the values in each row you always get 1 and summing the values in each column you always get rows / 3.
Finally, after struggling a bit, I think I managed to get a solution with the following code, that I kindly ask you to correct or improve. (I have set rows = 6 because the number of permutations of the obvious solution is 6!/(2!*2!*2!) = 90, whereas setting rows = 21 I would have got 21!/(7!*7!*7!) = 399,072,960 solutions.)

from ortools.sat.python import cp_model

# Create model.
model = cp_model.CpModel()

# Create variables.
rows = 6
columns = 3
x = []
for i in range(rows):
  x.append([model.NewBoolVar(f'x[{i}][{j}]') for j in range(columns)])

# Add constraints.
for i in range(rows):
  model.Add(sum(x[i]) == 1)

# Uncomment the following four lines of code if you want to solve the slightly more general problem that asks to enumerate
# all boolean arrays, with a given number of rows and columns, filled in such a way that summing the values in each
# row you always get 1 and summing the values in each column you always get no more than the ceiling of (rows / columns).
# if rows % columns != 0:
#   for j in range(columns):
#     model.Add(sum(x[i][j] for i in range(rows)) <= rows // columns + 1)
# else:
  for j in range(columns):
    model.Add(sum(x[i][j] for i in range(rows)) == rows // columns)


class MyPrintedSolution():

    def __init__(self, sol, sol_number):
        self.sol = sol
        self.sol_number = sol_number

    def PrintReadableTable(self):
        print(f'Solution {self.sol_number}, printed in readable form:')
        counter = 0
        for v in self.sol:
          if counter % columns != columns-1:
            print(v, end = ' ')
          else:
            print(v)
          counter += 1
        print()

    def PrintRawSolution(self):
        print(f'Solution {self.sol_number}, printed in raw form:')
        counter = 0
        for v in self.sol:
          print(f'{v}', end = '')
          counter += 1
        print('\n')


class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):

    def __init__(self, variables, limit):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0
        self.__solution_limit = limit

    def solution_count(self):
        return self.__solution_count

    def on_solution_callback(self):
        self.__solution_count += 1
        solution = [self.Value(v) for v in self.__variables]
        myprint = MyPrintedSolution(solution, self.__solution_count)
        myprint.PrintReadableTable()
        # myprint.PrintRawSolution()
        if self.__solution_count >= self.__solution_limit:
          print(f'Stop search after {self.__solution_limit} solutions')
          self.StopSearch()


# Create solver and solve model.
solver = cp_model.CpSolver()
# solver.parameters.num_workers = 16  # Solver works better with more workers. (At least 8, 16 if enough cores.)
# solver.parameters.log_search_progress = True
solver.parameters.enumerate_all_solutions = True
# solver.parameters.max_time_in_seconds = 10.0
solution_limit = 100000
solution_printer = VarArraySolutionPrinter([x[i][j] for i in range(rows) for j in range(columns)], solution_limit)
solver.Solve(model, solution_printer)

Update n. 2

Following @Christopher Hamkins' initial roadmap and subsequent precious suggestions, I think I finally got what I wanted, using the following code (although I am of course always open to corrections or further suggestions).

from ortools.sat.python import cp_model

# Create model.
model = cp_model.CpModel()

# Create variables.
classes = 21 # indexed with "i", but one could as well have chosen "c"
terms = 2 # indexed with "j", but one could as well have chosen "t"
days = 3 # indexed with "k", but one could as well have chosen "d"
persons = 8 # indexed with "p"
persons_names = [' 0 ', 'Ava', 'Bob', 'Joe', 'Mia', 'Sam', 'Tom', 'Zoe']
classes_names = ['1A', '1B', '1C', '1D', '1E', '1F', '1G', '2A', '2B', '2C', '2D', '2E', '2F', '2G', '3A', '3B', '3C', '3D', '3E', '3F', '3G']
classes_p = [[] for _ in range(persons)]
classes_p[0] = list(range(classes))
classes_p[1] = [0, 1, 4, 6, 10, 11, 13, 14, 15, 19, 20] # list of classes in which person 1 can work
classes_p[2] = [2, 3, 5, 7, 8, 9, 12, 16, 17, 18] # list of classes in which person 2 can work
classes_p[3] = [2, 4, 5, 7, 8, 10, 14, 15, 18, 20] # list of classes in which person 3 can work
classes_p[4] = [0, 1, 3, 6, 9, 12, 13, 16, 17, 19] # list of classes in which person 4 can work
classes_p[5] = [1, 2, 7, 9, 15, 17, 20] # list of classes in which person 5 can work
classes_p[6] = [0, 3, 8, 10, 12, 16, 19] # list of classes in which person 6 can work
classes_p[7] = [4, 5, 6, 11, 13, 14, 18] # list of classes in which person 7 can work
x = {}
for i in range(classes):
  for j in range(terms):
    for k in range(days):
      for p in range(persons):
        x[i, j, k, p] = model.NewBoolVar(f'x[{i}, {j}, {k}, {p}]')

# Add constraints.
"""
For all i, j, k constrain the sum of x[i, j, k, p] over p in the range of people to be equal to 1,
so exactly nobody or one person is selected at a given slot.
"""
for i in range(classes):
  for j in range(terms):
    for k in range(days):
      model.Add(sum(x[i, j, k, p] for p in range(persons)) == 1)

"""
For all i constrain the sum of x[i, j, k, p] over all j, k, p in their respective ranges (except p = 0)
to be exactly equal to 2, so exactly two people are in a given row.
"""
for i in range(classes):
  model.Add(sum(x[i, j, k, p] for j in range(terms) for k in range(days) for p in range(1, persons)) == 2)

"""
For all i, k, and for p = 0, add the implications
x[i, 0, k, 0] == x[i, 1, k, 0]
"""
for i in range(classes):
 for k in range(days):
    model.Add(x[i, 0, k, 0] == x[i, 1, k, 0])

"""
For all i, p (except p = 0), constrain the sum of x[i, j, k, p] over all j and k
to be at most 1.
"""
for i in range(classes):
  for p in range(1, persons):
    model.Add(sum(x[i, j, k, p] for j in range(terms) for k in range(days)) <= 1)
    # for k in range(days): # Equivalent alternative to the previous line of code
    #   model.AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])

"""
For all j, k constrain the sum of x[i, j, k, p] over all i, p in their respective ranges (except p = 0)
to be exactly equal to 7, so exactly seven people are in a given column.
"""
for j in range(terms):
  for k in range(days):
    model.Add(sum(x[i, j, k, p] for i in range(classes) for p in range(1, persons)) == 7)

"""
For all j, k, p (except p = 0) constrain the sum of x[i, j, k, p] over all i
to be exactly equal to 1, so each person appears exactly once in the column.
"""
for j in range(terms):
  for k in range(days):
    for p in range(1, persons):
      model.Add(sum(x[i, j, k, p] for i in range(classes)) == 1)

"""
For all j and k, constrain x[i, j, k, p] == 0 for the row i in which each person p can't appear.
"""
for p in range(persons):
  for i in enumerate(set(range(classes)) - set(classes_p[p])):
    for j in range(terms):
      for k in range(days):
        model.Add(x[i[1], j, k, p] == 0)


class MyPrintedSolution():

    def __init__(self, sol, sol_number):
        self.sol = sol
        self.sol_number = sol_number

    def PrintReadableTable1(self):
        print(f'Solution {self.sol_number}, printed in first readable form:')
        print('    |      Mon      |      Wed      |      Fri     ')
        print(' Cl | Term1   Term2 | Term1   Term2 | Term1   Term2') 
        print('----------------------------------------------------', end='')
        q = [_ for _ in range(8)] + [_ for _ in range(24, 32)] + [_ for _ in range(8, 16)] + [_ for _ in range(32, 40)] + [_ for _ in range(16, 24)] + [_ for _ in range(40, 48)]
        r = []
        for i in range(21):
            r += [n+48*i for n in q]
        shuffled_sol = [self.sol[m] for m in tuple(r)]
        counter = 0
        for w in shuffled_sol:
          if (counter % (persons * days * terms)) == 0:
            print('\n ', classes_names[counter // (terms * days * persons)], sep='', end=' |')
          if w:
            print('  ', persons_names[counter % persons], sep='', end='   ')
          counter += 1
        print('\n')

    def PrintReadableTable2(self):
        print(f'Solution {self.sol_number}, printed in second readable form:')
        print(' Cl |      Term1             Term2       ')
        print(' Cl | Mon   Wed   Fri   Mon   Wed   Fri  ')
        print('----------------------------------------', end = '')
        counter = 0
        for v in self.sol:
          if (counter % (persons * days * terms)) == 0:
            print('\n ', classes_names[counter // (terms * days * persons)], sep = '', end = ' |')
          if v:
            print(' ', persons_names[counter % persons], sep = '', end = '  ')
          counter += 1
        print('\n')

    def PrintRawSolution(self):
        print(f'Solution {self.sol_number}, printed in raw form:')
        counter = 0
        for v in self.sol:
          print(f'{v}', end = '')
          counter += 1
        print('\n')


class VarArraySolutionPrinter(cp_model.CpSolverSolutionCallback):

    def __init__(self, variables, limit):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variables = variables
        self.__solution_count = 0
        self.__solution_limit = limit

    def solution_count(self):
        return self.__solution_count

    def on_solution_callback(self):
        self.__solution_count += 1
        solution = [self.Value(v) for v in self.__variables]
        myprint = MyPrintedSolution(solution, self.__solution_count)
        myprint.PrintReadableTable1()
        # myprint.PrintReadableTable2()
        # myprint.PrintRawSolution()
        if self.__solution_count >= self.__solution_limit:
          print(f'Stop search after {self.__solution_limit} solutions')
          self.StopSearch()


# Create solver and solve model.
solver = cp_model.CpSolver()
# solver.parameters.num_workers = 16  # Solver works better with more workers. (At least 8, 16 if enough cores.)
# solver.parameters.log_search_progress = True
solver.parameters.enumerate_all_solutions = True
# solver.parameters.max_time_in_seconds = 10.0
solution_limit = 20
solution_printer = VarArraySolutionPrinter([x[i, j, k, p] for i in range(classes) for j in range(terms) for k in range(days) for p in range(persons)], solution_limit)
status = solver.Solve(model, solution_printer)

Update n. 3

@AirSquid proposed a solution using PuLP which is to me almost as valuable as the one using CP-SAT. It provides only one solution at a time, but (it has other advantages and) one can always get around this by adding some ad hoc further constraints, for example to see a different solution with a certain person in a specific position.


Solution

  • Your "toy" problem is definitely going in the right direction.

    For your actual problem, try making a 21×2×3x8 array x indexed with i, j, k and p (for person) of BoolVar's. The last index represents the person, it will need 0 to represent "nobody" and for the rest Ava = 1, Bob = 2, etc., so its max value will be one more than the number of people. If the variable X[i,j,k,p] is true (1) it means that the given person p is present at the index i, j, k. If X[i,j,k,0] is true, it means a 0 = nobody is present at i, j, k.

    For all i, j, k, constrain the sum of x[i, j, k, p] for p in the range of people to be equal to 1, so exactly nobody or one person is selected at a given slot.

    For point 1: fixed the index i you have exactly two empty 2-tuples and one 2-tuple with different non-zero values:

    For all i constrain the sum of x[i, j, k, p] for all j, k, p in their respective ranges (except p = 0) to be exactly equal to 2, so exactly two people are in a given row.

    For all i, k, and for p = 0, add the implications

    x[i, 0, k, 0] == x[i, 1, k, 0]

    This will ensure that if one of the pair is 0, so is the other.

    For all i, k and p except p = 0, add the implications

    x[i, 0, k, p] implies x[i, 1, k, p].Not and

    x[i, 1, k, p] implies x[i, 0, k, p].Not

    (Actually one of these alone should be sufficient)

    You can directly add an implication with the AddImplication(self, a, b) method, or you can realize that "a implies b" means the same thing as "b or not a" and add the implication with the AddBoolOr method. For the first implication, with x[i, 0, k, p] as a, and x[i, 1, k, p].Not as b, therefore adding:

    AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])
    

    Note that both variables are negated with Not in the expression.

    Since the other implication assigns x[i, 1, k, p] as a, and x[i, 0, k, p].Not as b, the resulting expression is exactly the same

    AddBoolOr([x[i, 0, k, p].Not(), x[i, 1, k, p].Not])
    

    so you only need to add it once.

    This will ensure a tuple will consist of two different people.

    Alternative formulation of the last part:

    For all i and p except p = 0, constraint the sum of x[i, j, k, p] for all j and k to be exactly equal to 1.

    For point 2: fixed the index k you have exactly fourteen empty 2-tuples and seven 2-tuple with different non-zero values;

    For all j and k constrain the sum of x[i, j, k, p] for all i and p (except p=0) in their respective ranges to be exactly equal to 7, so exactly seven people are in a given column.

    For all j, k, and p (except p = 0) constrain the sum of x[i, j, k, p] over all i to be exactly equal to 1, so each person appears exactly once in the column (that is, once for each value of the indices j and k, for some value of i).

    For point 3:

    For all j and k, Constrain x[i, j, k, p] == 0 for the row i in which each person p can't appear.

    Let us know how it works.