# Introduction to Constraint Programming in Python

The Programming Paradigm to Find One Solution Among 8,080,104 Candidates

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**:

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.

```
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.

```
# 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.

```
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.')
```

We obtained our solution in less than a millisecond: there are **9,139 soldiers** in the enemy army. Huzzah, we can now fire the scouts!

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.

```
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$$

```
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$$

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

The model is complete, CP-SAT can **solve the 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.')
```

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

📝 Chapter 1: Introduction to Linear Programming