Written by: Paul Rubin

Primary Source: OR in an OB World, 12/16/2018.

Okay, this is the last post on the subject. I promise! If you’re coming into this movie on the last reel, you may need to skim the last few posts to see what it’s about. I’m trying not to repeat myself too much.

To summarize where we are at: Hardmath123 posted a solution (generated by a heuristic) that might or might not be optimal; a constraint programming model gets that solution very quickly and fails to improve on it to the extent of my experience (I let it run four hours); and the MIP1 model was the least worst of the various optimization models I tried, but it had rather loose bounds.

Rob Pratt, in a comment to the second post, suggested a new constraint that does improve the bound performance of MIP1. Recall that the variables in MIP1 are \(x_{i,j}\) (binary, does symbol \(i\) go in slot \(j\)), \(p_i\) (the position where symbol \(i\) gets parked), and \(d_{ij}\) (the distance between symbols \(i\) and \(j\) in the layout). Rob’s constraint is

\(\displaystyle \sum_{i} \sum_{j>i} d_{ij} = \frac{n^3 – n}{6}\)

where \(n\) is the size of the symbol set (26 in our case). I’ll skip the derivation here; it’s based on the observation that the sum of the distances between all pairs of positions is constant, regardless of who fills those positions.

That got me thinking, and I came up with a set of constraints that are similarly motivated. Let \(\delta_j\) be the sum of the distances from position \(j\) to all other positions. Since we are using zero-based indexing, there are \(j\) positions to the left of position \(j\) with distances \(1,2,\dots,j\) (looking right to left) and \(n -1-j\) positions to the right of \(j\), with distances \(1,2,\dots, n-1-j\) (looking left to right). Using the well known formula for summing a sequence of consecutive integers,

\(\displaystyle \delta_j = \frac{j(j+1)}{2}+\frac{(n-1-j)(n-j)}{2}.\)

It follows that if symbol \(i\) lands in position \(j\), \(\sum_j d_{ij} = \delta_j\). Considering all possible landing spots for symbol \(i\) leads to the following constraints:

\(\displaystyle \sum_k d_{ik} = \sum_j \delta_j x_{ij}\quad \forall i.\)

A little experimentation showed that the combination of those constraints plus Rob’s constraint improved MIP1 more than either did alone. That’s the good news. The bad news is that I still haven’t been able to prove optimality for Hardmath123’s solution. I ran MIP1 with the extra constraints (using branching priorities, which help), using the Hardmath123 solution as the start and changing the CPLEX emphasis setting to 3. That last bit tells CPLEX to emphasize pushing the bound, which makes sense when you think you have the optimal solution and you’re trying to get to proof of optimality. After 338 minutes (5.6 hours), the incumbent had not improved (so it’s looking good as a possible optimum), but the gap was still around 14.5%, and the search tree was north of 1.25 GB and growing. To get a sense of what I might call the “futility index”, it took eleven minutes to shrink the gap from 14.58% to 14.48%.

Here are a few final takeaways from all this.

- There are often multiple ways to write a MIP model for a discrete optimization problem. They’ll often perform differently, and (at least in my experience) it’s hard to anticipate which will do best.
- Constraint programming has potential for generating good incumbent solutions quickly for some discrete optimization problems, particularly those whose structure fits with the available selection of global constraints in the solver. Permutation problems will meet this criterion. CP may not be particularly helpful for proving optimality, though, since the bounding is weak.
- Seemingly redundant constraints can help … but there’s no guarantee. (Rob suggested an additional constraint in a comment to part III. I tried it, and it slowed things down.)
- Some days you get the bear, and some days the bear gets you.

#### Paul Rubin

#### Latest posts by Paul Rubin (see all)

- Using Java Collections with CPLEX - July 22, 2019
- Evaluating Expressions in CPLEX - June 27, 2019
- R v. Python - June 17, 2019