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

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.

ReplyDelete# 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 payoff[s] := LoadValue.val;

let PACKED[s] := {i in ITEMS, k in KNAPSACKS : Packed[i,k] > 0.5};

}

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

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[1] 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).

ReplyDeleteTallys: 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".

DeleteI 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 [1] entry of that vector does not exist. Can you post your working version (just the portion you changed)?