## Saturday, February 26, 2011

### 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];
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
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 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;
}
}
}


1. Self-replies are usually the sign of a disturbed mind, but in this case I'll take my chances. Bob Fourer pointed out a more elegant way of coding the exclusion constraints on the AMPL user group. I trust he won't mind my reposting it here. The model is unchanged; I'll repeat the extra code, suitably modified, below.

# we now create some structure to record multiple solutions
param want := 5; # number of solutions we want
param nfound default 0; # number of solutions found
param payoff {1..nfound}; # objective values of recorded solutions
set PACKED {1..nfound} in {ITEMS, KNAPSACKS}; # records which items were packed (and where) in each solution

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

# solve until either the problem becomes infeasible or the right number of solutions have been found
for {s in 1..want} {
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 PACKED[s] := {i in ITEMS, k in KNAPSACKS : Packed[i,k] > 0.5};
}

This version is more readable. (Thanks, Bob.)

2. Firstly: I want to thank Paul and Bob for posting this. This is exactly what I was looking for today. Secondly: the code didn't work for me as posted. I had to move the "let nfound := s" command after the "let PACKED[s] := ..." command. Otherwise, AMPL would give me an error message saying that PACKED was undefined. I believe that as soon as the value of nfound is updated, AMPL tries to update the Exclude constraint, which doesn't work because PACKED[s] hasn't been updated yet (just a guess; Bob may have a better explanation for this).

1. Tallys: Thanks for reporting the bug. I thought I ran the updated code before posting it, but obviously not. I just noticed a typo: "Packed[i,k] > 0.5" should be "Pack[i,k]>0.5".

I can confirm the error message you report, but no permutation of those last three lines works for me. If either of the last two lines occurs before nfound is incremented from its default value of 0 to 1, I (properly) get an error that the  entry of that vector does not exist. Can you post your working version (just the portion you changed)?

Due to intermittent spamming, comments are being moderated. If this is your first time commenting on the blog, please read the Ground Rules for Comments. In particular, if you want to ask an operations research-related question not relevant to this post, consider asking it on Operations Research Stack Exchange.