The Weight Problem of Bachet de Meziriac

In a comment in this post there was a question about the Weight Problem of Bachet de Meziriac.

From this description: A merchant had a forty pound weight that broke into four pieces as a result of a fall. When the pieces were subsequently weighed, it was found that the weight of each piece was a whole number of pounds and that the four pieces could be used to weigh every integer weight between 1 and 40 pounds. What were the weights of the pieces? This problem comes from the French mathematician Claude Gaspard Bachet de Méziriac (1581-1638) who solved it in his famous book Problèmes plaisants et délectables qui se font par les nombres, published in 1624. We assume a balance scale where we can place weights in the  left or right pan like: Photo by Nikodem Nijaki (Own work) [CC BY-SA 3.0], via Wikimedia Commons

MINLP formulation

A first simple MINLP formulation can look like:

The variable \(\delta_{i,k} \in \{-1,0,+1\}\) indicates whether we place the weight \(i\) on the left or the right scale during trial \(k\). Of course the equation Check is non-linear, so we need to solve this with a MINLP solver. The non-convexity we introduced requires a global MINLP solver such as Baron or Couenne. This formulation is probably similar to how one would formulate things for use with a Constraint Programming solver.

CP (Constraint Programming) formulation

Indeed a CP formulation in Minizinc can be found on Paul Rubins blog. CP solvers typically have no problem with many types of non-linearities as long everything stays integer-valued (CP solvers have limited support for continuous variables).

MIP formulation

It is not too difficult to make a MIP model from this. The linearization of the multiplication of a binary variable times a continuous variable is well known: 

\[ \boxed{\begin{align}&y= \delta \cdot x\\ &\delta\in\{0,1\}, x\ge 0 \end{align}}\]

\[\>\>\Longleftrightarrow\>\>\]

\[\boxed{\begin{align} & y\le \delta M\\ & y\le x\\ & y\ge x-M(1-\delta)\\ & \delta \in\{0,1\},x,y \ge 0 \end{align}}\]

Here \(M\) is a tight upper-bound on \(x\).

Instead of \(\delta_{i,k} \in \{-1,0,+1\}\) we now use \(\delta_{i,k,lr} \in \{0,1\}\) where \(lr=\{left,right\}\). We don’t want to place the same weight on the left and right scale, so we add the constraint \(\delta_{i,k,’left’}+\delta_{i,k,’right’}\le 1\). Note that this constraint is not really needed: the final result would not change if we left this constraint out.  Here are the changes to the model:

We also added the Order constraint: \(x_{i+1}\le x_i\). This gives a nicer display of the results but also reduces symmetry. This can make large models (much) easier to solve.

Results

The final result is:

----     40 VARIABLE x.L  the four weights w1 27,    w2  9,    w3  3,    w4  1

The MINLP solvers Baron and Couenne do actually pretty good on the MINLP version of this model.