« A modulo propagator for ECLiPSe/ic | Main | Two new tools for MiniZinc, and a paper »

Optimizing shopping baskets: The development of a MiniZinc model

Yesterday (Thursday) was a holiday in Sweden and I had a busy day. Apart from writing a modulo propagator for ECLiPSe/ic, I also spend a lot of time with the Stack Overflow question How to use constraint programming for optimizing shopping baskets?. It asked for some tips (in a Java constraint programming system) for solving a problem of optimizing a shopping basket where different shops has different prices for items, as well as a delivery cost for orders below a certain limit.

You can read what I wrote and the developments of the models in my answer (I am hakank).

Update 2010-05-16: Added a better version. See below, Version 6.

First version

Even though the question was about a Java constraint system (I mentioned JaCoP, though), I rather quickly throw together a MiniZinc model instead as some kind of a prototype for a Java model. I'm a very lazy person and MiniZinc is excellent for prototyping.

This first version is shopping_basket.mzn (which actually include the next requirement: to handle shops which don't have all items). It solved the very simple example described in the problem section without any problem.

This model is quite straightforward: x is a matrix of 0/1 decision variables which shop to buy which item (the "knapsack" part). The more tricky part is to calculate the delivery costs if the total is below a certain limit (delivery_costs_limits). Here is the constraint part of the model:
constraint
   % we must buy all items (from some shop)
   forall(i in 1..num_items) (
       sum([x[i,j] | j in 1..num_shops]) = 1
   )
   /\ 
   total = sum(i in 1..num_items) (
       % costs of the products
       sum(j in 1..num_shops) (
          x[i,j]*costs[i,j]
       )
   ) + sum(j in 1..num_shops) (
       %  and delivery costs if total for a shop < limit
       let {
         var int: this_cost = sum([costs[i,j]*x[i,j] | i in 1..num_items])
       } in
       delivery_costs[j]*bool2int(this_cost > 0 /\ this_cost < delivery_costs_limits[j])
   );

The N/A problem: Not all shops sells all items

The next question was how to handle the "N/A" cases, i.e. for shops that don't have some item. It was solved by setting the price to a large number (999999 or 99999). However this is not very nice to the constraint solver. I really wanted to set these N/A to 0 (zero) but didn't get this correct (see below for a wrong turn on this).

As mentioned above the model shopping_basket.mzn includes the N/A problem.

Another representation

Almost always there are different approached how to model a problem. In the next version, shopping_basket2.mzn, I represented x - instead of a 0/1 matrix - as an array of length 1..num_items with the domains 1..num_shops to represent which shop to buy an item. The idea was that the constraint solvers now had much less work to do.

The constraint section is now simpler in some part, but calculating the delivery costs required a little more code:
constraint
 total = sum(i in 1..num_items) (
     costs[i,x[i]]
 ) + 
 sum(j in 1..num_shops) (
   let {
   var int: this_cost = 
     sum([costs[i,j]*bool2int(x[i]==j) | i in 1..num_items])
   } in
   delivery_costs[j]*bool2int(
      this_cost > 0 /\ 
      this_cost < delivery_costs_limits[j])
 );

Real data

The simple example data was easily solved for these two models, but the real challenge began when handling real data: 29 items and 37 shops. The two simplistic models was, however, not strong enough to handle this, so I started to think harder about the problem.

A wrong turn

The first attempt was continuing with the matrix approach and trying to be clever to represent N/A by 0 (zero), see above. This was done in shopping_basket3.mzn, but this model is plain wrong! The reason for this rather embarrissing bug was that the test case I used was not correct (completely mea culpa of course).

Final(?) version, part I: limiting the domains

OK, back to the drawing board. After some deep thoughts (i.e. an afternoon nap) I realized that the first models was too much of a proof-of-concept, too much prototype. In these I had broken one of the cardinality rules of constraint programming: not thinking about the domains of all decision variables.

In the first models the total cost (total) was defined as
var int: total :: is_output;
and the temporary variables this_cost was also defined as var int. This means that there is no limits of these variables (not even that they should be positive) which gives little hints for the constraint solvers.

The remedy is shown in the final(?) version shopping_basket5.mzn (shopping_basket4.mzn is basically the same thing for the matrix approach, and was not published).

Here the maximum total for any shop is first calculated, which can be quite large given the 99999 for N/A, but it is still limited and positive.
% total price
int: max_total :: is_output = 
      sum(i in 1..num_items) (
         max([costs[i,j] | j in 1..num_shops] )
      );
total is then defined with this limit:
var 0..max_total: total :: is_output; 
The temporary variable this_cost is handled accordingly.

Final(?) version, part II: labeling

This model was now in better shape, but still too slow. For all these versions I tried a lot of different labelings (search heuristics) and with many solvers (including the two MIP solvers, MiniZinc/mip and ECLiPSe/eplex, but they didn't accept the models).

I first saw a real improvement with the combination of Gecode/fz and the largest,indomain_max labeling. It solved the problem in 25 seconds (and with 69964 failures) :
total = 42013
x = [13, 20, 17, 18, 18, 13, 17, 17, 20, 13, 17, 17, 13, 13, 
     18, 17, 13, 13, 17, 8, 13, 36, 10, 17, 13, 13, 17, 20, 13] 
Just a minute or so after publishing this result (which seemed to please the person asking the question), I tested a different labeling, largest,indomain_split, and it solved the problem in 12 seconds (51013 failures). It gave a slightly different solution of where to buy item 12 (shown as bold):
total = 42013
x = [13, 20, 17, 18, 18, 13, 17, 17, 20, 13, 17, 13, 13, 13, 
     18, 17, 13, 13, 17, 8, 13, 36, 10, 17, 13, 13, 17, 20, 13]. 
These two solutions are the only two solutions for the (optimal) total of 42013. The test for all solutions takes about 1 second.

Some comments

This was a fun problem, quite simple but instructive what will happen when ignoring declaration of proper domain variables in the first versions of the models.

Even though the "real problem" was of the magnitude of 29 items and 37 shops (29 * 37 =1073), I am somewhat surprised that the problem was so hard to solve. This reminds me of the coins_grid problem which MIP solvers solve very fast, but the constraint solver have problems with it.

Update 2010-05-16: Version 6: Now faster

Well, version 5 was not the final version. Maybe this version 6 (shopping_basket6.mzn) is?

By two simple changes the problem is now solved in about 2 seconds and with 4460 failures (compared with 12 seconds and 51013). Here are the changes:
* The first change was the representation of N/A:s from 99999 to 0 (zero) . As noted above the former representation is not a good since it make all calculations, and the domains, unnecessary large.
* And finally, the following constraint states that all prices to consider must be larger than 0:
   forall(i in 1..num_items) (
      costs[i,x[i]] > 0
   ) 
Sigh. Of course! No defense of my silliness there is.

This also lessens the surprise considerable in Some comment above.