Finding Multiple Solutions in a Binary MIP

I'm going to post some AMPL code here that illustrates how to find more than one solution to a mixed integer program in which all integer variables are binary.  The problem I'll use is a multiple knapsack where the supply of each item is one.  The AMPL code should be pretty self-explanatory, with the possible exception of the exclusion constraints.  The notion for those is as follows.  Suppose that my binary variables are $x\in \{0,1\}^d$.  I solve the original model and obtain a solution with $x=\tilde{x}$.  I now add the following constraint:$$\sum_{i : \tilde{x}_i = 0} x_i + \sum_{i : \tilde{x}_i =1}(1-x_i)\ge 1.$$Note that this constraint forces at least one of the binary variables to change its value. (Exclusion constraints for general integer variables are harder to come by.  The only way I know off-hand involves binary expansions of the general integer variables.)

The AMPL code solves the model in a loop, recording each solution and adding an exclusion constraint.  The small example I created happens by chance to have multiple optimal solutions, enough so that the five solutions generated by the code are all optimal (but distinct).  More generally, you could get a sequence of increasingly suboptimal solutions.

Three things are perhaps noteworthy in the code.  The first is the check for infeasibility:  if you've already enumerated all possible solutions (and excluded them), the model will become infeasible.  The second is that I gave default values for the exclusion constraint coefficients (and right hand sides).  AMPL otherwise solved the base model (with no exclusion constraints) but then nagged me when I tried to record the result, because the first exclusion constraint used undefined coefficients.  Apparently recording the result (or even trying to display it) triggered a reevaluation of the model.  I think I could have ducked that by waiting to update nfound, but default values were an easy fix.  Third, note that I did not check whether a binary variable equaled 0 or 1; that might bump into rounding issues.  Saying that anything above (below) 0.5 is a 1 (0) is a very safe test.  (If your rounding error is that bad, you need to clean the beads on your abacus.)

One last comment is that this is not the most efficient way to find multiple solutions, in that it requires the (expanded) model to be solved ab initio each time.  If the solver supports callbacks, a more efficient approach is to use callbacks to record incumbents, reject them, and add exclusion constraints on the fly.  CPLEX has a "solution pool" feature that lets you track multiple solutions (and if CPLEX has it, my guess is that other commercial solvers do or soon will).  Using that is probably even  more efficient than using callbacks.  If the model is easy enough to solve, though, or if your solver lacks those features, the method demonstrated below may be a good choice.

Here is the code:

# Example of how to generate and record multiple solutions to a binary MIP.
#
# The sample problem is a 0-1 multiple knapsack.
#
set ITEMS;  # items available for the knapsacks
set KNAPSACKS;  # available knapsacks
param value {ITEMS} > 0;  # item values
param size {ITEMS} > 0;  # item sizes
param capacity {KNAPSACKS} > 0;  # knapsack capacities
var Pack {ITEMS, KNAPSACKS} binary;
  # Pack[i,k] = 1 iff item i is packed in knapsack k
maximize LoadValue: sum {i in ITEMS, k in KNAPSACKS} value[i]*Pack[i,k];
  # maximize value of load
s.t. Capacity {k in KNAPSACKS}:
  sum {i in ITEMS} size[i]*Pack[i,k]  <= capacity[k];  # capacity limits
s.t. Supply {i in ITEMS}:
  sum {k in KNAPSACKS} Pack[i,k]  <= 1;  # only one of each item available

option solver cplexamp;  # to use CPLEX as the solver

data;  # ordinarily I'd use a separate file for the data
set ITEMS := 1 2 3 4 5 6 7 8 9;
set KNAPSACKS := 1 2;
param capacity := 1 6 2 4;
param:  value  size  :=
   1       6     3
   2       5     2
   3       4     4
   4       3     1
   5       2     1
   6       4     5
   7       7     6
   8       2     3
   9       1     1  ;

model;

# we now create some structure to record multiple solutions
param want := 5;  # number of solutions we want
set SOLUTIONS := 1..want;  # index set for solutions
param packed {ITEMS, KNAPSACKS, SOLUTIONS};
  # packed[i,k,s] will be 1 if solution s packed item i in knapsack k
param payoff {SOLUTIONS};  # objective values of recorded solutions
param coef {ITEMS, KNAPSACKS, SOLUTIONS} default 0;
  # coefficients of solution exclusion constraints
param rhs {SOLUTIONS} default 0;  # right hand sides of exclusion constraints
param nfound default 0;  # number of solutions found

# we add constraints to the previous model to exclude
# solutions we've already seen
s.t. Exclude {s in 1..nfound}:
  sum {i in ITEMS, k in KNAPSACKS} coef[i,k,s]*Pack[i,k] >= rhs[s];

# solve until either the problem becomes infeasible or
# the right number of solutions has been found
for {s in SOLUTIONS} {
  solve;  # solve the model
  if solve_result = 'infeasible' then break;
    # if the model has become infeasible, give up
  let nfound := s;  # bump the counter for solutions found
  let payoff[s] := LoadValue.val;
  let rhs[s] := 1;
  # create an exclusion constraint that forces at least
  # one binary variable to change value
  for {i in ITEMS, k in KNAPSACKS} {
    if Pack[i,k] > 0.5 then {
      # treat this as a 1
      let packed[i,k,s] := 1;
      let coef[i,k,s] := -1;
      let rhs[s] := rhs[s] - 1;
    }
    else {
      # treat this as a 0
      let packed[i,k,s] := 0;
      let coef[i,k,s] := 1;
    }
  }
}