User Tools

Site Tools


Project 02: Reasoning About Constraints in a Real-World Problem

  • Released on: 10/26
  • Due on: 11/23 11:59 pm

There are three objectives for this project:

  1. Understand the basics of the Vehicle Routing Problem (VRP).
  2. Become familiar with two different constraint encodings of the VRP.
  3. Design and implement an extension to the baseline VRP model.

Obtain the Project

sudo apt-get update
sudo update_intro_to_autonomy

As this project doesn't focus on integrating together several tools (unlike Project 1). We will be returning to the Jupyter notebook for the project. You can find the project files under projects/project2.

Dependencies (Students Not Using the VM)

This project makes use of Numberjack to formulate problems, and SCIP as a solver.

Find SCIP and an academic license in chef/site-cookbooks/intro-to-autonomy/files/default/scipopt. Place the contents in /usr/local/src. (No need to decompress yet.)

Run the following script from /usr/local/src as root to compile SCIP with NumberJack flags:


set -e
# Check SCIP not compiled yet
if [ -f "$file" ]
  echo "Compiled SCIP library found. Exiting compilation of SCIP bash script successfully."
  exit 0
  echo "SCIP compiled library not found in: $file. Proceeding to compile SCIP....."
  tar xvf scipoptsuite-3.2.0.tgz
  cd scipoptsuite-3.2.0
  make scipoptlib ZIMPL=false ZLIB=false READLINE=false GAMS=false GMP=false LEGACY=true SPX_LEGACY=true -j3
  echo "SCIP compiled sucessfully. Exiting compilation of SCIP bash script successfully."

Then, set a SCIP environment variable, which will be used when Numberjack is installed:

ENV['ZIBPATH'] = '/usr/local/src/scipoptsuite-3.2.0'

Finally, install Numberjack, matplotlib, and numpy through pip. The Python files required for this project are located in the Project 2 git repository.

If You Cannot Use SCIP

If you are unable to use SCIP in the VM, please follow these instructions.

Project Overview

The second half of 16.413 focuses on constraint modelling and reasoning. So far, in the psets for this half, you've been applying these techniques to toy problems. The overall goal of this project is to introduce you to a very common (and very important) real-world problem: The Vehicle Routing Problem (VRP). There are many different ways to encode this problem. In the project, you will see two different encodings, each coming from a different research community. You will implement both encodings using Numberjack. Last, you will choose a side-constraint and add it to the basic encoding of the problem.

The VRP problem should already be vaguely familiar to you, as it is similar to the problem you were modelling in the first part of the project.


For this project, you must deliver a single tar archive of the project2 folder (see pset_submission). You must include a fully evaluated Jupyter notebook, as well as any additional files needed to run your notebook. There are no auto-graded cells in the notebook, so the final structure of the notebook is completely up to you. The notebook must answer all questions from the project and explain what is going on in your code.

It is recommended that you do not modify any .py files we've given you, in case we need to push a fix.

Part 0: Introduction

First, you must learn about the VRP problem. Read the following sections of Handbook of Constraint Programming (2006) (also posted on Stellar):

  • Chapter 23 intro
  • All of section 23.1
  • Section 23.2, except 23.2.3, 23.2.4, 23.2.5, 23.2.6
  • Section 23.3, except 23.3.3, 23.3.4, 23.3.5
  • All of section 23.6
  • And all of section 23.7.

While this paper discusses adding time windows to the problem, we will be ignoring time windows for this project.

Sections 23.3.[1-2] are very dense, especially if you're not familiar with the constraint programming community. Unfortunately, this paper mixes together the problem encoding (the constraints) and how the problem is solved (the propagation rules). We will be focusing on the encoding. Still, if you have trouble picking apart the encoding from the propagation rules, don't worry; we'll walk you through the encoding more in depth later in the project. However, it is important you understand the general gist of how they are modelling the problem as we will be using the same terminology.

Once you are done reading, look at the code we've provided you. contains most of what you need to know. It defines four classes.

  1. VRP is the class that represents the VRP problem. The important functions are capacity, demand, distance, num_vehicles, and num_customers. Additionally, we only consider problems with integer demands and capacities.
  2. Route details the route for a single vehicle in the solution.
  3. VRPSolution is a wrapper class that holds a list of Route objects. If there is no solution, routes should be None.
  4. VRPModel is the base class for models of the VRP problem. You will be extending this class in your solutions.

Important Note 1: the paper you read was written by someone in the constraint programming community. It is common in that community for arrays, variable names, etc. to be 1-indexed. However, Numberjack inherits Python's 0-indexing. In the code provided in, we use 0-indexing. Vehicles are numbered 0-(m-1), customers are numbered 0-n with 0 being treated specially as the warehouse.

Important Note 2: You will notice that VRP.cost() returns an integer instead of a real number. Normally this is not necessary when encoding VRPs, but it is useful in this case because of some Numberjack peculiarities that will be (briefly) discussed later. It will probably result in non-optimal answers, but they'll be close enough for our purposes. contains some code to fix the internal behavior of numberjack. You shouldn't need to touch it. provides some useful functionality that will be discussed in later sections.

Part 1: Mathematical Programming Formulation

The first step in this project is to encode the VRP using mathematical programming. Mathematical programming deals with integer and real-valued variables. Constraints are mathematical (in)equalities over the variables. In this case, we are using integer variables and all of the constraints are linear.

Implement the mathematical programming formulation of VRP found in Section 23.2.1 of the reading. You should finish the implementation of VRPMilpModel.create_model() and VRPMilpModel.solve(). Remember that we don't care about time windows for this project, so you can ignore those constraints.

Note that there are errata in the paper's encoding:

  1. Equation 23.7 should be ∀ S ⊆ C, k ∈ M
  2. Add the constraint Σ_{i ∈ ({0} ∪ C)} x[i, 0, k] = 0 ∀ k ∈ M (this is necessary because we aren't encoding time windows.)

Helpful notes:

  1. You may find powerset() from useful in implementing constraint 23.7.

Answer the following questions:

  1. How many constraints and variables are there in the model for m vehicles, n customers? Where m = {2, 5, 10} and n = {2, 5, 10}. How quickly does the number of constraints and variables grow with respect to m and n? (Think Big O notation)
  2. In your own words, what is the purpose of constraint 23.7? Plot an example (incorrect) solution that would be valid if Constraint 23.7 were omitted.

Part 2: Constraint Programming Formulation

Now we turn our attention to the CP formulation of the VRP.

Intro to CP

The CP community takes a much higher-level view of most problems than the mathematical programming community does. As such, they have a much wider array of constraint and variable types available to them.

For instance, in addition to integer and real-valued variables, CP solvers can reason natively about enumerated variables (think an integer variable with a name for each of its values instead of just a number), as well as arrays and matrices of variables. This opens up many intriguing possibilities. For instance, you can declare a Variable with the domain [0,9] and a 10 element VarArray with each variable having the domain [0,4] like such:

x = Variable(10)
a = VarArray(10, 5)

The real power, however comes in when you use Variables to index into arrays. The following constraint specifies that a[x] must be 2. So a solution with x=3, a=[0,0,0,2,0,0,0,0,0,0] would be valid, but a solution with x=0, a=[0,0,0,2,0,0,0,0,0,0] would not be.

a[x] == 2

Similar things work for matrices (Matrix class in Numberjack).

If we look at constraints in CP, we can see they have a wide array of constraints available to them, the most interesting of which are termed “global” constraints. Global constraints can often be represented as a series of smaller constraints and'ed together. However, CP solvers can reason about them natively and use powerful propagation techniques to speed up the search for a solution. One of the most well known global constraints is the AllDiff constraint. Assume you have a variable array with 3 elements. You can enforce that each Variable in the array takes on a different value using a series of != constraints, like so:

b = VarArray(4, 5)
b[0] != b[1]
b[0] != b[2]
b[0] != b[3]
b[1] != b[2]
b[1] != b[3]
b[2] != b[3]

Or you could use an AllDiff constraint like:

b = VarArray(4, 5)

You can see the global constraints included in Numberjack at:

For this project, you will make use of the Circuit constraint (implemented in When given a variable array, the Circuit constraint asserts that if you treat every element of the array as a pointer back into the array, the path created uses every variable in the array and makes a closed circuit. For instance a solution with a=[1, 2, 3, 4, 0] would satisfy Circuit(a), because element 0 points to element 1, which points to element 2, which points to element 3, which points to element 4, which points back to element 0. A solution with a=[1, 2, 0, 4, 3] would not satisfy Circuit(a), because there is not one circuit that covers every variable (instead, there are two circuits 0→1→2→0 and 3→4→3).

Last, a hallmark of CP encodings of problems is that they often have redundant constraints. While redundant constraints tend to be useless or even a problem for mathematical solvers (especially if the redundant constraints are over real-valued variables), redundant constraints tend to give more information to CP solvers which allows them to solve problems faster. Redundant constraints and variables in the VRP case would be having successor and predecessor variables for each visit.

A CP Encoding

As thinking about constraints the CP way is probably not familiar to most of you, we've provided more code to get you started. Additionally, CP is a fast moving field and the encoding from the paper can be vastly abbreviated using global constraints. As such, the CP model we want you to encode is below.

The successor of each vehicle's last visit should be the next vehicle's first visit.

1. s[l(k)] = f(k+1) ∀ k ∈ M - {m}
2. s[l(m)] = f(0)

Paths must be consistent. That is, the successor of a visit's predecessor should be itself. And the predecessor of a visit's successor should be itself.

3. s[p[i]] = i ∀ i ∈ V
4. p[s[i]] = i ∀ i ∈ V

The vehicles used on each path should be consistent.

5. v[i] = v[p[i]] ∀ i ∈ V - F
6. v[i] = v[s[i]] ∀ i ∈ V - L

There should be no subtours (routes that don't visit the warehouse).

7. Circuit(s)
8. Circuit(p)

Model the load of the vehicle on each route:

9. q[i] = q[p[i]] + r[i] ∀ i ∈ V - F
10. q[i] = q[s[i]] - r[s[i]] ∀ i ∈ V - L

Model the capacity of each vehicle for a given route:

11. q[i] <= Q[v[i]] ∀ i ∈ V

We are interested in minimizing the total distance travelled.

12. minimize Σ_{∀ i ∈ V - L} δ[i, s[i]]

These constraints were noted as missing by some students on Piazza. They are added at the end to maintain numbering of previous constraints.

We need to make sure that vehicles only travel the route from their start visit to their end visit.

13. v[f(k)] = k ∀ k ∈ M
14. v[l(k)] = k ∀ k ∈ M

Implementation and Questions

Implement the CP formulation of VRP described above. You should finish the implementation of VRPCpModel.create_model() and VRPCpModel.solve().


  1. If you're not careful, redundant constraints can become incompatible constraints. You might want to model without the redundant constraints and variables first, then add the redundancy in once that's working.
  2. Numberjack has the ability to interface with Mistral, a true CP solver. However, either Mistral or Numberjack's interface to Mistral doesn't support real valued variables, so we are using SCIP to solve this problem. Additionally, the distance matrix we've defined in the provided code is real-valued even though the distances returned by VRP.distance() are all integers. This is to cut down on the number of integer variables SCIP has to deal with. Finding solutions with this encoding will likely take longer with SCIP than they would a true CP solver.

Answer the following questions:

  1. How many constraints and variables are there in the model for m vehicles, n customers? Where m = {2, 5, 10} and n = {2, 5, 10}. How quickly does the number of constraints and variables grow with respect to m and n? (Think Big O notation)

Part 3: Testing

Design eight VRP instances and run them through each of your encodings to get solutions. Describe each instance, describe each solution (in enough detail to convince the grader that the solution is correct, no need to convince that it is optimal), and provide a plot of each solution.

Five of your instances should have solutions.

Three of your instances should have no solution, including two that have no solution for non-trivial reasons.

Part 4: Side Constraint Extension

Choose one of the side constraints described in section 23.1.1 of the reading and augment either your CP or MILP model to include it.

You will need to make a new subclass of VRP and VRPModel for your augmented problem.

Describe the side constraint you chose and what variables and constraints it added to your encoding.

Provide 5 instances of your augmented problem (including 2 without a solution). Describe each instance, describe each solution (in enough detail to convince the grader that the solution is correct, no need to convince that it is optimal), and provide a plot of each solution.

Part 5: Wrap-up

Answer the following questions:

  1. Which formulation do you find easier to understand? Why?
  2. Do you see any application of VRPs to your own research?
  3. This is the first time this part of the graduate student project has been given. What are your thoughts? Anything we should definitely keep or add for next year? Anything we should definitely ditch?

Parting Note

We hope you enjoyed learning about reasoning about constraints in real-world problems and can find an application of it to your own research.

If you do see an application to your own research, please be aware that while Numberjack is great for pedagogical purposes, it has some serious issues that prevent it from being fully utilized in a non-exploratory research context. If you'd like to learn more about other options for specifying constraint programs, please drop the course staff a line and we can give you some pointers.


There are two errors in the provided code for part 2:

1. Correct

q = VarArray(len(V), 0, n + 2*m - 1, "q")


q = VarArray(len(V), 0, max([self.vrp.capacity(k) for k in M]), "q")

2. Correct

v = VarArray(len(V), 0, n + 2*m - 1, "v")


v = VarArray(len(V), 0, m - 1, "v")
projects/project2.txt · Last modified: 2016/11/16 21:16 by aytonb