Maxime Labonne - Introduction to Constraint Programming in Python

Constraint Programming is a technique to find every solution that respects a set of predefined constraints.

It is an invaluable tool for data scientists to solve a huge variety of problems, such as scheduling, timetabling, sequencing, etc. In this article, we’ll see how to use CP in two different ways:

Satisfiability: the goal is to find one or multiple feasible solutions (i.e., solutions that respect our constraints) by narrowing down a large set of potential solutions;
Optimization: the goal is to find the best feasible solution according to an objective function, just like Linear Programming (LP).

We’ll use CP-SAT from Google OR-Tools, an excellent free and open source CP solver. Note that it is different from MPSolver, which is dedicated to Linear and Mixed Integer Programming. The difference between CP and LP is quite confusing, we’ll touch on this topic at the end of the article.

You can run the code with the following Google Colab notebook.

!python -m pip install --upgrade --user -q ortools

🪖 I. Satisfiability with the 3 scouts problem

In the previous article, we created an army to defeat our opponent. But there was one small problem: we had to guess how powerful his army was.

This time, let’s send scouts to know the exact number. Our 3 scouts observed the enemy camp, and this is what they tell us:

Scout 1: “the number of soldiers is a multiple of 13”;
Scout 2: “the number of soldiers is a multiple of 19”;
Scout 3: “the number of soldiers is a multiple of 37”;
They all agree that the number of soldiers doesn’t exceed 10,000.

Our scouts have a personal way of counting soldiers, but we can combine these three observations to make a model.

Let’s call the number of soldiers army. We can translate our problem into the following congruence system:

army \equiv 0 \mod 13 \\ army \equiv 0 \mod 19 \\ army \equiv 0 \mod 37

If you’re not familiar with this notation, this is what it means in programming terms:

army\ \% \ 13 = 0 \\ army\ \% \ 19 = 0 \\ army\ \% \ 37 = 0

Let’s implement it with OR-Tools. The first thing we need to do is to import and create the CP-SAT model and solver.

from ortools.sat.python import cp_model

# Instantiate model and solver
model = cp_model.CpModel()
solver = cp_model.CpSolver()

The modeling process is very similar to what we did in Linear Programming.

The first step to create our CP model is to declare the variables. In this example, we only have one: army, the number of soldiers.

We have to give lower and upper bounds. The lower bound is 1 since we know there’s an army, and the upper bound is 10,000 according to the scouts:

1 \leq army \leq 10\ 000

In OR-Tools, we use the NewIntVar method to create this variable.

# 1. Variable
army = model.NewIntVar(1, 10000, 'army')

The second step is to declare the constraints.

We identified three constraints in this example. Modulo is a special operator, so we need a specific function to handle it with CP-SAT: AddModuloEquality. You can find a reference guide at this address if you need other methods.

# 2. Constraints
# variable % mod = target → (target, variable, mod)
model.AddModuloEquality(0, army, 13)
model.AddModuloEquality(0, army, 19)
model.AddModuloEquality(0, army, 37)

Unlike Linear Programming, we don’t have to define an objective function here.

The reason is simple: there is nothing to optimize! We just want to find a feasible solution that satisfies our constraints, but there is no “good” or “bad” answers. This is a key feature of Constraint Programming.

Our model is complete, we can now ask OR-Tools to solve it.

# Find the variable that satisfies these constraints
status = solver.Solve(model)

# If a solution has been found, print results
if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
    print('================= Solution =================')
    print(f'Solved in {solver.WallTime():.2f} milliseconds')
    print()
    print(f'🪖 Army = {solver.Value(army)}')
    print()
    print('Check solution:')
    print(f' - Constraint 1: {solver.Value(army)} % 13 = {solver.Value(army) % 13}')
    print(f' - Constraint 2: {solver.Value(army)} % 19 = {solver.Value(army) % 19}')
    print(f' - Constraint 3: {solver.Value(army)} % 37 = {solver.Value(army) % 37}')

else:
    print('The solver could not find a solution.')

================= Solution =================
Solved in 0.01 milliseconds

🪖 Army = 9139

Check solution:
 - Constraint 1: 9139 % 13 = 0
 - Constraint 2: 9139 % 19 = 0
 - Constraint 3: 9139 % 37 = 0

We obtained our solution in less than a millisecond: there are 9,139 soldiers in the enemy army.

We limited the search space with an upper bound of 10,000, which gave us a unique solution. But is it still the case if we push this limit?

Another perk of CP is the ability to find every possible solution to a problem. This might take a long time when the search space is large because the solver has to brute force the entire space (instead of reducing it with heuristics). Let’s explore this feature by printing every possible solution with a new upper bound of 100,000.

With OR-Tools, we ask the solver to look for every possible solution thanks to the enumerate_all_solutions parameter. We then assign it a callback class that prints every solution the solver finds.

model = cp_model.CpModel()
solver = cp_model.CpSolver()

# 1. Variable
army = model.NewIntVar(1, 100000, 'army')

# 2. Constraints
model.AddModuloEquality(0, army, 13)
model.AddModuloEquality(0, army, 19)
model.AddModuloEquality(0, army, 37)


class PrintSolutions(cp_model.CpSolverSolutionCallback):
    """Callback to print every solution."""

    def __init__(self, variable):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__variable = variable

    def on_solution_callback(self):
        print(self.Value(self.__variable))

# Solve with callback
solution_printer = PrintSolutions(army)
solver.parameters.enumerate_all_solutions = True
status = solver.Solve(model, solution_printer)

We found 10 solutions! This was to be expected since we increased the upper bound tenfold: these solutions all are multiples of 9,139.

As you can see, this example has nothing to do with optimization: it’s a pure satisfiability problem. On another note, this congruence system can be solved manually with the Chinese remainder theorem. But CP is not limited to that…

🍻 II. Optimization and beer

Let’s see another problem: our army will face the enemy in a few days. In the meantime, the quartermaster has to prepare the rations that will be used during the campaign.

The space in the supply wagons is limited and some rations are more popular than others. There are three possible rations:

🥖 Bread: it takes only 1 space but soldiers don’t like it that much with a popularity of 3;
🥩 Meat: it takes 3 spaces and has a popularity of 10;
🍺 Beer: it takes 7 spaces but soldiers love it with a popularity of 26.

The supply wagons have a capacity of 19 spaces. How to select the best rations to maximize the popularity?

This is an optimization problem we’ve already seen: actually, it is a variant of the famous knapsack problem. We could reuse the code from the previous article and just change the input parameters.

This time, we’ll solve it using Constraint Programming. This paradigm is not limited to finding feasible solutions. It can also perform optimization using different algorithms to handle this overhead.

Let’s create a model of the problem. First of all, we have to declare three variables: 🥖bread, 🥩meat, and 🍺beer. It’s possible to have 0 of them, but their number cannot exceed the maximal capacity.

0 \leq bread \leq capacity \\ 0 \leq meat \leq capacity \\ 0 \leq beer \leq capacity

# Instantiate model and solver
model = cp_model.CpModel()
solver = cp_model.CpSolver()

# 1. Variables
capacity = 19
bread = model.NewIntVar(0, capacity, 'bread')
meat  = model.NewIntVar(0, capacity, 'meat')
beer  = model.NewIntVar(0, capacity, 'beer')

This time, we only have one constraint: the space occupied by the bread, the meat, and the beer cannot exceed the wagons’ capacity (19).

1 \times bread + 3 \times meat + 7 \times beer \leq 19

# 2. Constraints
model.Add(1 * bread
        + 3 * meat 
        + 7 * beer <= capacity)

We want to maximize the total popularity of the rations that are selected:

max\ 3 \times bread + 10 \times meat + 26 \times beer

# 3. Objective
model.Maximize(3  * bread
             + 10 * meat
             + 26 * beer)

The model is complete, CP-SAT can solve the problem!

# Solve problem
status = solver.Solve(model)

# If an optimal solution has been found, print results
if status == cp_model.OPTIMAL:
    print('================= Solution =================')
    print(f'Solved in {solver.WallTime():.2f} milliseconds')
    print()
    print(f'Optimal value = {3*solver.Value(bread)+10*solver.Value(meat)+26*solver.Value(beer)} popularity')
    print('Food:')
    print(f' - 🥖Bread = {solver.Value(bread)}')
    print(f' - 🥩Meat  = {solver.Value(meat)}')
    print(f' - 🍺Beer  = {solver.Value(beer)}')
else:
    print('The solver could not find an optimal solution.')

================= Solution =================
Solved in 0.01 milliseconds

Optimal value = 68 popularity
Food:
 - 🥖Bread = 2
 - 🥩Meat  = 1
 - 🍺Beer  = 2

We obtained the highest popularity (68) possible with a capacity of 19.

Is the constraint respected? Let’s quickly check it: 1×2🥖+3×1🥩+7×2🍺 = 19, which is indeed ≤ 19.

Okay, I’d like to ask another question: how many solutions to this problem are there? Once again, we can answer it with a specific callback to count them.

class CountSolutions(cp_model.CpSolverSolutionCallback):
    """Count the number of solutions."""

    def __init__(self):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__solution_count = 0

    def on_solution_callback(self):
        self.__solution_count += 1

    def solution_count(self):
        return self.__solution_count

solution_printer = CountSolutions()

# Instantiate model and solver
model = cp_model.CpModel()
solver = cp_model.CpSolver()

# 1. Variables
capacity = 19

bread = model.NewIntVar(0, capacity, 'Bread')
meat  = model.NewIntVar(0, capacity, 'Meat')
beer  = model.NewIntVar(0, capacity, 'Beer')

# 2. Constraints
model.Add(1 * bread
        + 3 * meat 
        + 7 * beer <= capacity)

# Print results
solver.parameters.enumerate_all_solutions = True
status = solver.Solve(model, solution_printer)
print(solution_printer.solution_count())

We found 121 solutions with a capacity of 19. But this number quickly increases: with a capacity of 1000, there are 8,080,104 possible solutions! And yet, CP-SAT finds the optimal solution in less than a second. How is it possible?

CP solvers do not brute force the problem with an exhaustive search, but combine heuristics and combinatorial search instead. More specifically, the three most popular techniques for constraint satisfaction problems are backtracking, constraint propagation, and local search.

CP-SAT is quite particular since it combines CP and SAT: it is part of a broader trend of merging CP, LP, SAT, and metaheuristics.

We said that the previous problem could be solved with Linear Programming, so let’s compare the code of both solutions:

As you can see, the syntax is quite similar but it’s not the same: model/solver vs. solver, NewIntVar instead of IntVar, etc. There’s a bit of translation to do, but it’s easily manageable.

These two techniques are incredibly close to each other: they both handle variables with constraints and perform optimization using math and heuristics. However, CP is limited to discrete parameters, while LP handles continuous ones. On the other hand, you can implement specialized constraints like “all different” in CP, but not in LP. Here is a summary of the main differences between these two technologies:

If you want to know more about this topic, I would recommend this article by Irvin J. Lustig and Jean-François Puget. CPLEX’s documentation also details the differences at this address, in terms of modeling and optimization.

Conclusion

Constraint Programming is another incredible technique in the mathematical optimization toolbox. It is a radically different approach compared to traditional, declarative programming. In this article,

We saw two applications of CP with satisfiability and optimization;
We implemented CP models in OR-Tools and played with the callback function;
We highlighted the differences between CP and LP.

We limited ourselves to simple problems in this introduction, but CP has amazing applications in complex scheduling and routing problems. This is a topic I’d love to address in a future article.

If you’re interested to know more about it, feel free to follow me on Twitter @maximelabonne. Thanks for your attention!

🥇 Linear Programming Course

🔎 Course overview

📝 Chapter 1: Introduction to Linear Programming

📝 Chapter 2: Integer vs. Linear Programming

📝 Chapter 3: Constraint Programming

📝 Chapter 4: Nonlinear Programming for Marketing Budget Allocation