Multiobjective integer programming

Question

I wish to use integer programming to enumerate pareto optimal solutions. I would like to implement an algorithm that uses gurobi or a similar single-objective integer programming solver to do this, but I don't know any such algorithms. Could someone please propose an algorithm for enumerating the efficient frontier?

Answer 1

In this answer, I'll address how to enumerate all pareto efficient solutions of 2-objective pure integer optimization problems of the form

min_x {g'x, h'x}
s.t.  Ax <= b
      x integer

Algorithm

We start the algorithm by optimizing for one of the objectives (we'll use g here). Since this is a standard single-objective integer optimization problem, it can be easily solved with gurobi or any other LP solver:

min_x g'x
s.t.  Ax <= b
      x integer

We initialize a set P, which will eventually contain all the pareto efficient solutions, to P = {x*}, where x* is the optimal solution of this model. To get the next point on the efficient frontier, which has the second-smallest g'x value and an improved h'x value, we can add the following constraints to our model:

1. Remove x* from the feasible set of solutions (details on these x != x* constraints are provided later in the answer).
2. Add a constraint that h'x <= h'x*

The new optimization model that we need to solve is:

min_x g'x
s.t.  Ax <= b
      x != x* for all x* in P
      h'x <= h'x* for all x* in P
      x integer

Again, this is a single-objective integer optimization model that can be solved with gurobi or another solver (once you follow the details below on how to model x != x* constraints). As you repeatedly solve this model, adding the optimal solutions to P, solutions will get progressively larger (worse) g'x values and progressively smaller (better) h'x values. Eventually, the model will become infeasible, which means no more points on the pareto frontier exist, and the algorithm terminates.

At this point, there may be some pairs of solutions x, y in P for which g'x = g'y and h'x > h'y, in which case x is dominated by y and can be removed. After filtering in this way, the set P represents the full pareto efficient frontier.

x != x* Constraints

All that remains is to model constraints of the form x != x*, where x and x* are n-dimensional vectors of integer variables. Even in one dimension this is a non-convex constraint (see here for details), so we need to add auxiliary variables to help us model the constraint.

Denote the n variables stored in the optimization model (collectively denoted x) as x_1, x_2, ..., x_n, and similarly denote the variable values in x* as x*_1, x*_2, ..., x*_n. We can add new binary variables y_1, y_2, ..., y_n to the model, where y_i is set to 1 when x_i > x*_i. Because x_i and x*_i are integer valued, this is the same as saying x_i >= x*_i + 1, and this can be implemented with the following constraints (M is a large constant):

x_i - x*_i >= 1 - M(1-y_i) for all i = 1, ..., n

Similarly, we can add new binary variables z_1, z_2, ..., z_n to the model, where z_i is set to 1 when x_i < x*_i. Because x_i and x*_i are integer valued, this is the same as saying x_i <= x*_i - 1, and this can be implemented with the following big-M constraints:

x_i - x*_i <= -1 + M(1-z_i) for all i = 1, ..., n

If at least one of the y or z variables is set, then we know that our x != x* constraint is satisfied. Therefore, we can replace x != x* with:

y_1 + y_2 + ... + y_n + z_1 + z_2 + ... + z_n >= 1

In short, each x != x* constraint can be handled by adding 2n binary variables and 2n+1 constraints to the model, where n is the number of variables in the original model.

Answer 2

PolySCIP is an academic open-source solver for multi-objective mixed integer linear programs. In case you do not want to implement your own solver or want to compare yours with another one.