CP Approach to Integer Optimization

Constraint programming (CP) represents a different approach to optimization than that of classical optimization theory. CP is based on feasibility (finding feasible solutions to a problem), and focuses on the constraints and variables rather than the objective function. For many types of problems, CP can find optimal solutions faster than a MIP solver.

Should I use MIP or CP-SAT?

For standard integer programming problems, in which a feasible point must satisfy all constraints, the MIP solver is faster. In this case, the feasible set is convex: for any two points in the set, the line segment joining them lies entirely in the set.

On the other hand, for problems with highly non-convex feasible sets, the CP-SAT solver is often faster than the MIP solver. Such problems arise when the feasible set is defined by many constraints joined by "or," meaning that a point only needs to satisfy one of the constraints to be feasible. For an example, see Assignment with allowed groups.

Example: transforming non-integer constraints

To increase computational speed, both the CP-SAT solver and the original CP solver work over the integers. To solve a problem in which some of the constraints have non-integer terms, you must first transform those constraints by multiplying them by a sufficiently large integer. The following example and solution illustrate this.

Maximize 2x + 2y + 3z subject to the following constraints:
x + 72 y + 32 z25
3x - 5y + 7z45
5x + 2y - 6z37
x, y, z0
x, y, z integers

Since the first constraint has non-integer terms, it must be transformed into one with all integer terms. This is shown in the Constraints section below.

Solution using the CP-SAT solver

The following sections present a Python program that solves the problem above using the CP-SAT solver. We don't particularly recommend using the CP-SAT solver for a problem like this, instead of the MIP solver. The purpose of the example is just to show how to set up an integer programming problem for CP-SAT, and how to transform non-integer constraints.

Some sections of the program are quite similar to the MIP program described previously, but there are some significant differences, particularly in defining the constraints and objective. If you're planning on using both solvers, it's worthwhile to compare these two examples.

Declare the model

The following code declares the model for the problem.

from ortools.sat.python import cp_model


def main():
  model = cp_model.CpModel()

Create the variables

The following code creates the variables for the problem.

  var_upper_bound = max(50, 45, 37)
  x = model.NewIntVar(0, var_upper_bound, 'x')
  y = model.NewIntVar(0, var_upper_bound, 'y')
  z = model.NewIntVar(0, var_upper_bound, 'z')

Define the constraints

As noted above, the CP-SAT solver works over the integers. Since the first constraint,

x + 72 y + 32 z25

has non-integer terms, you must first multiply the entire constraint by a sufficiently large integer so that all terms are integers. In this case, you can multiply by 2, which results in the new constraint

2x + 7y + 3z50

This doesn't change the problem, since the solution set for the original constraint is the same as for the transformed constraint.

The following code defines the three linear constraints for the problem:

  model.Add(2*x + 7*y + 3*z <= 50)
  model.Add(3*x - 5*y + 7*z <= 45)
  model.Add(5*x + 2*y - 6*z <= 37)

Define the objective function

The following code defines the objective function for the problem and declares it to be a maximization problem:

  model.Maximize(2*x + 2*y + 3*z)

Call the solver

The following code calls the solver and displays the results.

  solver = cp_model.CpSolver()
  status = solver.Solve(model)

  if status == cp_model.OPTIMAL:
    print('Maximum of objective function: %i' % solver.ObjectiveValue())
    print()
    print('x value: ', solver.Value(x))
    print('y value: ', solver.Value(y))
    print('z value: ', solver.Value(z))

The output is shown below:

Maximum of objective function: 35

x value:  7
y value:  3
z value:  5

The entire program

The entire program is shown below.

from __future__ import print_function
from ortools.sat.python import cp_model


def main():
  model = cp_model.CpModel()
  var_upper_bound = max(50, 45, 37)
  x = model.NewIntVar(0, var_upper_bound, 'x')
  y = model.NewIntVar(0, var_upper_bound, 'y')
  z = model.NewIntVar(0, var_upper_bound, 'z')

  model.Add(2*x + 7*y + 3*z <= 50)
  model.Add(3*x - 5*y + 7*z <= 45)
  model.Add(5*x + 2*y - 6*z <= 37)

  model.Maximize(2*x + 2*y + 3*z)

  solver = cp_model.CpSolver()
  status = solver.Solve(model)

  if status == cp_model.OPTIMAL:
    print('Maximum of objective function: %i' % solver.ObjectiveValue())
    print()
    print('x value: ', solver.Value(x))
    print('y value: ', solver.Value(y))
    print('z value: ', solver.Value(z))


if __name__ == '__main__':
  main()

Solution using the original CP solver

The following sections present a Python program that solves the integer optimization problem using the original CP solver. (However, we recommend the faster CP-SAT solver.)

Declare the solver

the following code declares the solver.

from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import solver_parameters_pb2

def main():
  # Instantiate a CP solver.
  parameters = pywrapcp.Solver.DefaultSolverParameters()
  solver = pywrapcp.Solver('simple_CP', parameters)

Create the variables

The following code creates the variables.

  var_upper_bound = max(50, 45, 37)
  x = solver.IntVar(0, var_upper_bound, 'x')
  y = solver.IntVar(0, var_upper_bound, 'y')
  z = solver.IntVar(0, var_upper_bound, 'z')

The command solver.NewIntVar(0, var_upper_bound, 'x') creates an integer variable whose values are between 0 and var_upper_bound. Although the problem doesn't impose upper bounds on the individual variables, we still need to supply an upper bound to the solver. So we just choose a value that is large enough that it doesn't change the solutions. Setting the value to the maximum of the constraint upper bounds (shown below) is sufficient.

Define the constraints

The following code defines the constraints for the problem.

  solver.Add(2*x + 7*y + 3*z <= 50)
  solver.Add(3*x - 5*y + 7*z <= 45)
  solver.Add(5*x + 2*y - 6*z <= 37)

Define the objective function

The following code defines the objective for the problem and declares it to be a maximization problem:

  objective = solver.Maximize(2*x + 2*y + 3*z, 1)

Decision builders

Before calling the solver, you need to create a decision builder for the problem. The decision builder creates the search tree and determines the order in which the solver searches solutions. The following code creates the decision builder using the solver's Phase method.

  decision_builder = solver.Phase([x, y, z],
                                  solver.CHOOSE_FIRST_UNBOUND,
                                  solver.ASSIGN_MAX_VALUE)

The Phase method has three input parameters:

  • The first parameter contains the decision variables — the variables the solver uses to decide which node of the tree to visit next. In this case, the decision variables are x, y and z.
  • The second parameter specifies how the solver chooses the next variable for the search. Here the code uses CHOOSE_FIRST_UNBOUND, which means the solver chooses the first unbound variable. Other options are CHOOSE_LOWEST_MIN and CHOOSE_RANDOM.
  • The third parameter specifies how the solver chooses the next value of the current search variable to check. Here, we've chosen ASSIGN_MAX_VALUE, which means the solver chooses the largest value that hasn't already been searched. We could also have chosen ASSIGN_MIN_VALUE or ASSIGN_RANDOM_VALUE.

In this example, the decision builder doesn't play a crucial role, because the search tree is so small. But in larger problems, choosing the right parameters for the decision builder can speed up the search significantly.

Create a solution collector

The CP solver normally returns many feasible solutions to an optimization problem. A useful tool for storing these results is the solution collector. The following code shows how to create one.

  collector = solver.LastSolutionCollector()
  # Add the decision variables.
  collector.Add(x)
  collector.Add(y)
  collector.Add(z)
  # Add the objective.
  collector.AddObjective(2*x + 2*y +3*z)

In this example, the solution collector stores just the last solution returned by the solver. Since this is a maximization problem, each solution returned by the solver has a larger objective value than the previous one. As a result, the last solution returned by the solver is the optimal solution.

Call the solver

The following code calls the solver and displays the results.

  solver.Solve(decision_builder, [objective, collector])
  if collector.SolutionCount() > 0:
    best_solution = collector.SolutionCount() - 1
    print("Maximum of objective function:", collector.ObjectiveValue(best_solution))
    print()
    print('x= ', collector.Value(best_solution, x))
    print('y= ', collector.Value(best_solution, y))
    print('z= ', collector.Value(best_solution, z))
The results of the search, shown below, are the same as the results returned by the MIP solver:
$ python my_projects/int_prog.py
Optimal objective value = 23

x = 3
y = 2

If you want to see all solutions returned by the solver, rather than just the optimal one, you can use the AllSolutionCollector as follows:

  collector = solver.AllSolutionCollector()
  collector.Add(x)
  collector.Add(y)
  collector.AddObjective(obj_expr)
  solver.Solve(decision_builder, [objective, collector])
    Here's the output of the program:
Maximum of objective function: 35

x value:  7
y value:  3
z value:  5

The entire program

The entire program is shown below.

from __future__ import print_function
from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import solver_parameters_pb2

def main():
  # Instantiate a CP solver.
  parameters = pywrapcp.Solver.DefaultSolverParameters()
  solver = pywrapcp.Solver('simple_CP', parameters)
  var_upper_bound = max(50, 45, 37)
  x = solver.IntVar(0, var_upper_bound, 'x')
  y = solver.IntVar(0, var_upper_bound, 'y')
  z = solver.IntVar(0, var_upper_bound, 'z')
  solver.Add(2*x + 7*y + 3*z <= 50)
  solver.Add(3*x - 5*y + 7*z <= 45)
  solver.Add(5*x + 2*y - 6*z <= 37)
  objective = solver.Maximize(2*x + 2*y + 3*z, 1)
  decision_builder = solver.Phase([x, y, z],
                                  solver.CHOOSE_FIRST_UNBOUND,
                                  solver.ASSIGN_MAX_VALUE)
    # Create a solution collector.
  collector = solver.LastSolutionCollector()
  # Add the decision variables.
  collector.Add(x)
  collector.Add(y)
  collector.Add(z)
  # Add the objective.
  collector.AddObjective(2*x + 2*y +3*z)
  solver.Solve(decision_builder, [objective, collector])
  if collector.SolutionCount() > 0:
    best_solution = collector.SolutionCount() - 1
    print("Maximum of objective function:", collector.ObjectiveValue(best_solution))
    print()
    print('x= ', collector.Value(best_solution, x))
    print('y= ', collector.Value(best_solution, y))
    print('z= ', collector.Value(best_solution, z))

if __name__ == '__main__':
  main()

Send feedback about...