From IP to CP

Someone asked on OR Stack Exchange how to convert an integer programming model into a constraint programming model. I think you can reasonably say that it involves a "paradigm shift", for a couple of reasons.

The first paradigm shift has to do with how you frame the problem, mainly in terms of the decision variables. Math programmers are trained to turn discrete decisions with a logical flavor into binary variables. Discrete quantities, such as how many bins of a certain type to use or how many workers to assign to a task, are expressed as general integer variables, but most other things end up turning into a slew of binary variables. The problem being solved in the ORSE question illustrates this nicely.

The problem is as follows. You have $N$ participants in a tournament involving some kind of racing. Importantly, $N$ is guaranteed to be an even number. There is one track with two lanes, and races are spread over $N-1$ days. Every participant races head to head with every other participant exactly once, and nobody races twice in the same day. For whatever reason, the left lane is preferable to the right lane, and so there is a "fairness" constraint that nobody is assigned the left lane on more than $M$ consecutive days. For some reason, the author also imposed a second fairness constraint that nobody be assigned to the right lane on more than $M$ consecutive days. Dimensions for the author's problem were $N=20$ and $M=2.$

The model has to assign participant pairs (races) to days and also make lane assignments. To decide against whom I must race on a given day, someone building an IP model will use binary variables to select my opponent. Similarly, they will use binary variables to select my lane assignment each day. So the author of the question had in his IP model a variable array opp[Competitors][Competitors][Tracks][Days] taking value 1 "if competitor 'c1' races with 'c2' on track 't' on day 'd'".

CP models are more flexible in their use of variables, and in particular general integer variables. So to decide my opponent on a given day, I can just an integer variable array indexed by day where the value is the index number of my opponent on the given day. Similarly, I could (and would) use an integer variable indexed by day to indicate my lane assignment that day, although in this case that variable does turn out to be binary, since there are only two lanes.

The second paradigm shift has to do with constraints, and it ties to what solver you are using. IP models have a very limited constraint "vocabulary". They all understand linear equalities and inequalities, and some understand some combination of SOS1, SOS2, second order cone and implication constraints. That's pretty much it. CP solvers have a richer "vocabulary" of constraints, but with the caveat that not many of those constraints are universal. I would wager that every CP solver has the "all different" constraint, and they must have the usual arithmetic comparisons ($=,\neq,\lt,\le,\gt,\ge$). Beyond that, it pays to check in advance.

I wrote a CP model (in Java) using IBM's CP Optimizer (CPO) to solve the scheduling problem. Details of the model can be sussed out from the Java code, but I will mention a few pertinent details here.

• I did use an integer variable array to determine, for each combination of participant and day, the participant's opponent that day, as well as an integer array giving the lane assignment (0 or 1) for each combination of participant and day.
• To make sure that, on any day, the opponent of X's opponent is X I used CPO's  inverse constraint. The constraint inverse(f, g) says that f[g[x]] = x and g[f[x]] = x for any x in the domain of the inner function.
• To ensure that nobody raced the same opponent twice, I used allDiff, which is CPO's version of the all different constraint.
• We have to do something to force opponents in a race to be in different lanes. Let $x_{i,d}$ and $y_{i,d}$ denote respectively the opponent and lane assignment for participant $i$ on day $d.$ In mathematical terms, the constraint we want is $y_{x_{i,d},d} \neq y_{i,d}.$ Indexing a variable with another variable is impossible in an IP model. In CPO, I used the element constraint to do just that.

I added an objective function, namely to minimize the difference between the most and fewest times any participant gets assigned the preferred left lane. I also added one constraint to mitigate symmetry. Since any solution remains a solution (with the same objective value) under any permutation of the participant indices, I froze the first day's schedule as $1\  v.\  N$, $2\  v.\  N-1$, $3\  v.\  N-2$ etc.

On my decent but not screamingly fast PC, CPO found a feasible solution almost instantly and a solution with objective value 1 in under a second. In that solution, every participant gets the left lane either nine or ten times out of the 19 racing days. It's not hard to prove that 1 is the optimal value (you cannot have everybody get exactly the same number of left lane assignments), but don't tell CPO that -- it was still chugging along trying when it hit my five minute time limit.

My Java code is available from my repository under a Creative Commons 4.0 open source license.

The "Set" Card Game

Someone asked a question (actually a pair of questions) about the "Set" card game on OR Stack Exchange. This was the first I had heard of the game. For purposes of this post, the game play is irrelevant. The key is that it uses a deck of 81 cards, each of which has four "features" with three possible values for each feature. Every combination of feature values appears on exactly one card. Three cards form a "set" if, for each of the four features, either all three cards have the same value or all three cards have different values. As noted in the Wikipedia page, there are 1,080 possible "sets". (I'm putting "set" in quotes to avoid confusion with, you know, sets.) The OR SE question asked if an integer program (IP) could be used to generate all 1,080 sets.

The answer is yes, with the caveat that while you can do it, an IP model is not the most efficient way to do it. I won't repeat my IP model here, since (a) it's a bit of a space hog and (b) you can read it in my answer on ORSE. The model is pretty small (97 binary variables) and finds one "set" quickly. There are two ways to get it to find all possible "sets". Both involve adding a "no good" constraint for each solution found, i.e., a constraint that eliminates that solution but no others. The more tedious approach solves the progressively growing model 1,081 times (the last time producing a declaration of infeasibility, meaning all "sets" have been found). The less tedious approach requires a solver that supports a callback function letting you reject candidate solutions by adding a constraint cutting them off (the same "no good" constraint used in the first approach).

Since there is no objective and the constraints are mostly logical in nature, constraint programming (CP) strikes me as a more intuitive approach than IP. The specifics of a CP model will depend on what types of constraints your CP solver supports, but I think the model I tried is likely to work with a lot of solvers. Since I coded the models in Java, I'll use 0-based indexing here. My variables (all integer-valued) are: $x_s\in \lbrace 0,\dots, 80\rbrace,$ the index of the card in slot $s\in \lbrace 0,1,2\rbrace$ of the "set"; and $y_{f,s}\in \lbrace 0,1,2\rbrace,$ the value of feature $f$ in the card in slot $s.$ The constraints are as follows:

• allDifferent($x$) (no two slots can contain the same card);
• $x_0 \lt x_1 \lt x_2$ (to avoid counting the same "set" more than once, we require that the cards be listed in ascending index order);
• $y_{f,s} = a_{x_s,f},$ where $a_{c,f}$  is the value of feature $f$ in card $c$ (the value of a feature in a slot comes from the card in the slot); and
• allDifferent$(\lbrace y_{f,0}, y_{f,1}, y_{f,2}\rbrace) \vee (y_{f,0}=y_{f,1} \wedge y_{f,1}=y_{f,2})\,\forall f$ (for every feature, either all three slots have different values or all three slots have the same value).

This exploits two constraints that CP solvers typically have but IP solvers do not. One is the "all different" constraint, which forces a set of variables to take values that do not repeat. The other is the ability to use a decision variable as an index to either a constant vector or another variable vector. That comes into play in the third constraint, where the variable $x_s$ (the index of the card in slot $s$) is used as a subscript of $a$ to get the feature value of that card.

Even CP is probably overkill for this problem, though. We can compute the 1,080 possible solutions by brute force. I used a recursive function for this.

Here is how each method stacks up in terms of computation time. As I said, I coded it in Java, using CPLEX 22.1.1 for the IP models and CP Optimizer 22.1.1 for the CP model.

• brute force took between 0.1 and 0.2 seconds;
• the CP model also took between 0.1 and 0.55 seconds;
• solving the IP model once using a callback took between 1.4 and 3 seconds; and
• solving the IP model 1,081 times (without a callback) took around 56 seconds, give or take.

Solution times for all methods are a bit random (they change when I rerun the program), so treat them as approximate. In the interest of full disclosure, I throttled CPLEX to a single thread for the callback approach. Using multiple threads would have required both synchronizing the code (a pain) and also checking for repeated solutions (since hypothetically two parallel threads could independently locate the same "set"). The moral, to the extent there is one, is that sometimes brute force beats more sophisticated approaches.

You can find my Java code here if interested.

A Matrix Puzzle

A question on Mathematics Stack Exchange, "Placing number blocks so that the resulting matrix is symmetric", revolves around a square matrix whose entries are partitioned into blocks (submatrices). The submatrices need to be rearranged so that the resulting matrix is symmetric around its main diagonal. The author of the question asked if it were possible to do this using linear programming. I think the answer is no, but we can certainly do it with either an integer linear program or a constraint program. The author posts two examples ($4\times 4$ and $5\times 5$). I cobbled together Java code for a MIP model (using CPLEX) and a CP model (using CP Optimizer) and had no trouble solving both problems.

There are a couple of generalizations worth noting. First, while the example matrices contain integers, there is nothing preventing the two models from being used with other content types (real numbers, complex numbers, text). Some variable definitions would need modification, but conceptually both methods would work. Second, the author of the post used only one-dimensional blocks (row vectors or column vectors), but my code allows for arbitrary rectangular blocks. The code assumes that the "parent matrix" is square, but that would be easy enough to relax, so the same models (with minor adjustments) would work with arbitrary rectangular matrices.

I think that the puzzle makes a good vehicle for comparing MIP and CP applications to logic problems. In optimization problems, I suspect that MIP models often do a better job of computing objective bounds than do CP models. That is a non-issue here, since the problem is just to find a feasible solution. For logic problems and similar things (particularly scheduling), I think CP models tend to be more expressive, meaning certain types of constraints or relationships can be expressed more naturally with CP than with MIP (where the relationships turn into rather arcane and complicated adventures with binary variables). That applies here, where the CP model exploits the ability to use integer variables as subscripts of other variables.

As described in the PDF file, though, that subscripting ability has its limits. CP Optimizer will let you index a one-dimensional vector of variables using an integer variable, but not a two-dimensional array. In other words, x[y] is fine but x[y, z] is not (where x, y and z are all variables). The workaround I used is to flatten a 2-D matrix into a 1-D matrix. So if $N\times N$ matrix $x$ is flattened into $N^2\times 1$ vector $\hat{x},$  $x_{y,z}$ becomes $\hat{x}_{N(y - 1) + z}.$

The models are a bit too verbose to squeeze into a blog post, so I wrote them up in a separate PDF file. My Java code (which requires recent versions of CPLEX and CP Optimizer) can be had from my code repository under a Creative Commons license.

Scheduling a Round-Robin Tournament

Once in a while I come across a scheduling problem, and of course my first reaction is to think "integer program". A while back, for example, I used an IP model to work out a schedule for a rotating duplicate bridge game (in which both teams and hosts rotated) at the request of a buddy. Most recently, I saw a question on OR Stack Exchange ("Doubles Round Robin Sorting Algorithm") related to tournament scheduling. 

The problem involves six players playing pickleball (so two players on each team, with two sitting out, for each game). Given six players, there are 15 possible teams (pairings). Twelve games are scheduled each week for four weeks, with two sets of six games per week. The author of the question correctly computed that there are 45 possible game configurations. Over the course of the four weeks, he wanted every possible game played at least once (which leaves the final three game slots to be filled arbitrarily). The following conditions must be met:

• within each set of six games, every player plays four games and sits out two;
• no player sits out two consecutive games within the same set;
• no two players are partners more than once per set.

The problem poses no criterion for selecting among feasible schedules; we just want to find any one schedule that works.

My IP formulation can be found in my answer to the OR SE question. I will just mention that it uses binary variables $x_{g,s}=1$ if game $g$ goes in slot $s$ and 0 if not. Games are enumerated up front and indexed from 1 to 15. Slots refer to positions in the schedule and are numbered from 1 to 48, with six slots per set and 12 slots per week. Since no criterion is specified, we just use the default (minimize 0).

I also tried a constraint programming model, using general integer variables $s_1,\dots,s_{48},$ where $s_j \in {1,\dots,15}$ is the index of the game played in slot $j$ of the schedule. My CP formulation uses the "all different" constraint (fairly ubiquitous among CP solvers) to ensure that the first 45 slots each contain a different game. It uses the CP Optimizer "distribute" constraint to ensure that no game is played more than twice and the CP Optimizer "count" constraint to ensure that each player plays exactly four times per set (and thus sits twice) and that no team plays more than once in a set. I'm not sure how common those types of constraints are among CP solvers because I don't have much experience with CP solvers.

As it turns out, both CPLEX (for the IP model) and CP Optimizer (for the CP model) produced feasible schedules in about a second or so on my desktop PC. My intuition was that a CP model would be better suited to this type of problem, because all variables are integer and a CP solver would not be doing much if any matrix arithmetic. As it turns out, it would take a much larger test case to tell whether that intuition has any merit.

If anyone wants to play with the models, my Java code is available from my university GitLab repository. Running it would require CPLEX and CP Optimizer (and not necessarily the most recent versions of either).

"Block Party" Puzzle

A question posed on the OR Discord channel by a doctoral student led me to discover the existence of the Jane Street puzzle page. The student was asking about building a MILP model for the June 2022 puzzle, called "Block Party 4". The puzzle involves inserting numbers into a grid, with some cells already filled in. It bears a superficial resemblance to sudoku, but with a few key differences. Where a sudoku is divided into nine square regions of nine cells each, the block party puzzle grid is divided into connected regions of varying sizes and shapes. Within a region of $k$ cells, the numbers 1 through $k$ must be filled in. Finally, rather than requiring that rows and columns contain no repeated numbers, the rules require that, for each possible value $K$, if $K$ is inserted into a cell then the nearest instance of $K$ must be at distance exactly $K$ in the $L_1$ norm. So to use a 1, there must be a 1 in an adjacent cell. To use a 2, there must be a 2 in a cell two moves away but no 2 in any adjacent cell.

Since this is a problem with logic constraints and integer decisions, my instinct was to think that constraint programming would be faster than integer programming. To test this, I coded both an IP model and a CP model in Java, using CPLEX and CP Optimizer as the respective solvers. I assumed that the grid would be square, since both the June puzzle (10 x 10) and a smaller example provided (5 x 5) were. Both models can easily be adjusted for non-square grids.

Assume an $N\times N$ grid partitioned into regions, and let $M$ be the size of the largest region (and thus the largest value that can be used in the puzzle). Number the cells 1 through $N^2$ in any order. (I used a left-to-right raster scan.) For the IP model, I use binary variables $x_{i,j}$ $(i=1,\dots,N^2$, $j=1,\dots,M)$ to indicate whether value $j$ is inserted into cell $i$. For cells with known values, I fix $x_{i,j}$ to either 0 or 1 as appropriate while building the model. Also, if cell $i$ lies in a region of size $K$, then I can fix $x_{i,j}=0$ for $j>K.$

Since we just want a feasible solution, I let the IP objective function default to minimizing zero. The most obvious constraint is $$\sum_{j=1}^M x_{i,j} = 1 \quad \forall i,$$which forces a single value to be selected for each cell. Similarly, if $B$ is a block with size $K,$ then $$\sum_{i\in B}x_{i,j}=1 \quad \forall j=1,\dots,K$$forces every value from 1 to $K$ to be used exactly once in the block. Finally, for each block $B,$ each cell $i\in B$ and each legal value $j\in \lbrace 1, \dots, \vert B\vert\rbrace$ for that cell, we add these constraints: $$x_{i,j} \le \sum_{k\in N_j(i)} x_{k,j} $$ and $$x_{i,j} + x_{k,j} \le 1\quad \forall k\in N^-_j(i),$$ where $N_j(i)$ is the set of all cells at distance exactly $j$ from cell $i$ and $N^-_j(i)$ is the set of all cells at distance less than $j$ from cell $i$ (excluding cell $i$ itself). These enforce the rule that, for value $j$ to be used in cell $i,$, it must also be used in at least one cell at distance $j$ from $i$ and in no closer cell.

The CP model is a bit more straightforward to articulate. Again, there is no objective function, since we are just solving for a feasible solution. For each cell $i$, there is a single integer variable $x_i$ with domain $1,\dots,\vert B \vert$ where $B$ is the block containing cell $i.$ If we know that cell $i$ is fixed to value $k,$ we just declare $x_i$ to have domain $\lbrace k \rbrace.$ For each block, we use an "all different" constraint to enforce the requirement that the cells in the block take distinct values. For each cell $i$ and legal value $j$ for it, the implication constraint $$(x_i = j) \implies \bigvee_{k\in N_j(i)} (x_k = j)$$ where $\bigvee$ denotes disjunction ("or"), forces at least one cell at distance $j$ from $i$ to take value $j$ if cell $i$ does, while the constraints $$(x_i = j) \implies (x_k \neq j) \quad \forall k\in N^-_j(i)$$ prohibit any closer cell from using that value. (These constraints could be condensed into a conjunction on the right hand side. For reasons I have since forgotten, I did not bother to do so.)

Both models solved the 10x10 puzzle easily. My expectation was that the CP model would be faster, for several reasons. First, it has 100 general integer variables, whereas the IP model started out with 1,100 binary variables (which the presolver whittled down to 119 binary variables, compared to 90 variables for the CP model after presolving). Second, the "all different" CP constraint seems to be a more efficient way than a steaming pile of inequality constraints to enforce the rule that no two cells in the same block take the same value. Third, CP Optimizer would be doing integer arithmetic while CPLEX would be doing double precision arithmetic, and on a per-operation basis integer arithmetic should be faster. Lastly, my experience in the past has been that the one edge IP models tend to have over CP models is tighter bounds, but that has no effect in a feasibility problem (when you are not optimizing anything).

As it turns out, I was in for a surprise. Actually, make that two surprises. First, the IP model after presolving had 211 constraints, whereas the CP model after presolving had 7,399 constraints. Note that, in the implication constraints, the left side and each equality on the right side count as a constraint. I'm not sure how comparable constraints are between the two models, but I was not expecting the CP model to have so many more. Second, while both model solved in negligible time, the IP model was faster. CPLEX solved the IP model at the root node (no branching) in about 0.01 seconds on my fairly average desktop PC. CP Optimizer needed 2,678 branches and about 0.12 seconds to solve the CP model, of which 0.05 seconds was spent in the "engine" (i.e., solving) and the rest was spent in "extraction" (turning the model into something suitable for the engine).

My Java code (which requires both CPLEX and CP Optimizer but nothing else) can be found in my GitLab repository.

Wolf, Goat and Cabbage Part II

This continues my previous post about the "Wolf, Goat and Cabbage" logic puzzle. Using a MIP model to solve it strikes me as somewhat analogous to hammering in a nail with the butt end of a screwdriver handle: it works, but maybe it's not the ideal tool. To me, a constraint solver would make more sense, so I coded up a model using IBM's CP Optimizer (henceforth "CPO", part of the CPLEX Optimization Suite). Unsurprisingly, it worked (after some debugging) and produced the same solution.

MIP solvers largely share the same "vocabulary" for building models (real / integer / binary variables, linear and maybe quadratic expressions, equality and inequality constraints). From my limited exposure to constraint solvers, there is significant variation in both what things the solver does and does not understand. Some of that is variable types. A very basic solver may understand logical (true/false) and integer variables, and maybe real-valued variables (although I'm not sure handling real variables is anywhere near universal). CPO (and presumably some other solvers) understand "interval variables", which as the name suggests represent intervals of discrete values (usually time intervals, though possibly location or some other aspect). Similarly, different solvers will understand different types of global constraints. I suspect that every CP solver worth mentioning understands "all different" (no two variables in a bunch of integer variables can take the same value), but some solvers will implement "no overlap" constraints (the time interval during which I am eating and the time interval during which I am showering cannot overlap) or precedence constraints (this job has to end before this other job can start). Those kinds of constraints make certain scheduling problems easier to solve with CP than with a MIP model.

Anyway, I'm not entirely new to CPO, though far from proficient, and I tripped over a few "features" while coding the puzzle. I wanted to use boolean (true/false) variables for certain things, such as whether an item had made it to the far side of the river (true) or was stuck on the near side (false). CPO lets you declare a boolean variable but then treats it as an integer variable, meaning that you need to think in terms of 0 and 1 rather than false and true. So you can't say "if $x$ then ..."; instead, you need to say "if $x = 1$ then ..." (and trust me, the syntax is clunkier than what I'm typing here). When you go to get the value of your boolean variable $x$ after solving the model, CPO returns a double precision value. CPLEX users will be used to this, because in a MIP model even integer variables are treated as double-precision during matrix computations. CP solvers, though, like to do integer arithmetic (as far as I know), so it's a bit unclear why my boolean variable has to be converted from double precision to integer (or boolean). Even odder is that, at least in the Java API, there is a method that returns the value of an integer variable as an integer if you pass the name of the variable as the argument, but if you pass the actual variable you are going to get a double. (Did a federal government panel design this?)

In any event, logic of the CPO model is moderately straightforward, with constraints like "you can't carry something in the boat if it isn't on the bank the boat departs from" and "if the wolf and goat end up in the same place at the end of a period, the farmer better end up there too". Some things are bit less clunky with CPO than with CPLEX. For instance, to figure out what (if anything) is in the boat at a given time, the MIP model requires binary variables subscripted by time and item index (1 if that item is in the boat on that trip 0 if not). The CPO model just needs an integer variable for each time period whose value is either the index of the thing in the boat or a dummy value if the boat is empty. Furthermore, the nature of the variable automatically takes care of a capacity constraint. Since there is only one variable for what's in the boat, at most one thing (whatever that variable indicates) can be in the boat.

Some (most?) constraint solvers, including CPO, provide a way to use a variable as an index to another variable. In my code, the integer variable indicating what's in the boat at time $t$ is used to look up the location variable (near or far bank) for that item at time $t$ from a vector of location variables for all items at time $t$.

Anyway, the code in my repository has been updated to include the CPO model, and it's heavily commented in case you wanted to compare it to the MIP model.

Tournament Scheduling: CPLEX v. CP Optimizer (Part 3)

This is the third and final chapter in my opus about scheduling "pods" (players) in a tournament. Reading the first post will be essential to understanding what is going on here. Reading the second post (in which I formulate an integer programming model for the problem) is definitely optional.

I will stick with the notation introduced in the prior post. As before, there are $N$ pods (players) being organized into teams of two to play in $G$ games (two teams per game, one in white, one in black). Each game has two slots corresponding to the two teams in the game.

• $t$ indexes a team;
• $p$ and $q$ index pods;
• $g$ indexes a game;
• $s$ indexes a slot (and $s'$ is the other slot of the same game);
• $c$ is a jersey color (0 for white, 1 for black).

For a slot $s$, I will use $G_s$ and $C_s$ to denote respectively the game in which the slot exists and the jersey color (0 or 1) of the team in the slot. For a team $t$, $P_{t,0}$ and $P_{t,1}$ will be the two pods on the team. (The order they are listed is immaterial.)

The model variables are almost (but not exactly) what they were in the IP model.

• $x_{t}\in \lbrace 0, \dots, 2G-1\rbrace$ is the index of the slot in which team $t$ plays (where slots 0 and $G$ are the white and black teams in the first game and slots $G-1$ and $2G-1$ are the white and black teams in the final game).
  
• $y_{pqg}\in \lbrace 0,1\rbrace$ is 1 if pods $p$ and $q\neq p$ are playing in game $g$ on opposing teams.
  
• $z_{pg}\in \lbrace 0,1\rbrace$ is the color jersey worn by pod $p$ at the conclusion of game $g$ (regardless of whether or not $p$ played in game $g$).
• $w_{pg}\in \lbrace 0,1\rbrace$ is 1 if pod $p$ changes jersey color going into game $g.$
• $v_{pg} \in \lbrace 0,1\rbrace$ is 1 if pod $p$ is playing in game $g.$

The changes from the IP model are that $x$ is general integer rather than binary (with dimension 1 rather than 2), the third index of $y$ is the game rather than the slot, and $v$ is a new variable. As in the IP model, the objective is to minimize $\sum_p \sum_g w_{pg}.$

The constraints are as follows. Note that the requirement that every team play exactly once is captured by the assignment of a single slot value to each team (via $x$).

• Two different teams cannot occupy the same slot: $$\textrm{allDiff}(x).$$
• Occupying a slot implies playing in the game: $$x_t = s \implies (v_{P_{t,1},G_s} = 1 \wedge v_{P_{t,2},G_s} = 1) \quad \forall t,s.$$
• Occupying a slot determines the color of the jersey worn in the game: $$x_t = s \implies  (z_{P_{t,1},G_s} = C_s \wedge z_{P_{t,2},G_s} = C_s) \quad \forall t,s.$$
• Jersey color for a pod remains constant from game to game unless a change is recorded: $$ z_{p,g} \neq z_{p,g-1}  \iff w_{p,g} = 1 \quad \forall p, \forall g > 0.$$
• Playing in consecutive games precludes changing jerseys: $$(v_{p,g - 1} = 1 \wedge v_{p,g} = 1) \implies w_{p,g} = 0 \quad \forall p, \forall g > 0.$$
• Two pods are opponents if and only if they are playing in the same game with different colors: $$y_{pqg} = 1 \iff (v_{pg} = 1 \wedge v_{qg} = 1 \wedge z_{pg} \neq z_{qg} \quad \forall p, \forall q\neq p, \forall g.$$
• Every pair of pods faces each other exactly twice: $$\sum_g y_{pqg} = 2 \quad \forall p, \forall q\neq p.$$

Here "allDiff()" is the CPLEX notation for the "all different" constraint, a global constraint common to most constraint solvers. Given a collection of variables with the same range, the all different constraint says that no two variables in the collection can take the same value.

To put it mildly, I would not be shocked if there is a more efficient (easier to solver) way to write the problem as a CP model.

 

Tournament Scheduling: CPLEX v. CP Optimizer (Part 1)

 A recent question on Mathematics Stack Exchange asked about scheduling $N$ "pods" (players) in a tournament under the following conditions.

• There will be $N(N-1)/4$ games played sequentially (no games in parallel).
• Games pit teams of two pods each against each other. One team wears white and one team wears black.
  
• Each pod partners with every other pod once and plays against every other pod twice.
• No pod can play in consecutive games wearing different jersey colors.

The objective is find a schedule that minimizes the total number of times pods have to change their jersey color.

The question specifies $N=9$. Note that, since there are $N(N-1)/4$ games, a necessary condition for the problem to be feasible is that either $N$ or $N-1$ must be divisible by $4$. That condition is not, however, sufficient. An obvious counterexample is when $N=4$ and there are three games to be scheduled. The first game uses all four teams, as does the second game. Since no pod can be forced to wear different jerseys in consecutive games, all four teams would go into the second game wearing the same colors as in the first game ... which means being in the same teams, violating the rule about pods being teamed together exactly once.

I coded both an integer programming model (using CPLEX) and a constraint programming model (using CP Optimizer) in Java and tried to solve the tournament problem with $N=9$. Since everything is inherently integer and the constraints are more logical than arithmetic (the only real algebra is summing up the number of jersey changes), I had two initial expectations: that CPLEX would provide tighter lower bounds than CPO (because MIP models tend to do better with bounds than CP models); and that CPO would find feasible (and possibly optimal) schedules faster than CPLEX would, since the problem is inherently logic-based (and CP models often do better than MIP models on scheduling problems). I was half right. CPLEX did provide tighter lower bounds than CPO, at least within the time limits I was willing to use, although I don't think either came remotely close to a "tight" lower bound. CPO, however, struggled massively to find feasible schedules while CPLEX did not have too much trouble doing so.

Before going any further, I need to issue three caveats. First, the longest I ran the IP model was 30 minutes and the longest I ran the CP model was an hour. Second, while I am pretty familiar with formulating IP models, I am much less familiar running CP models, so I may have missed opportunities to make the CP model more efficient. Third, while CPLEX gives the user a fair amount of indirect control over the search process (by setting branching priorities, MIP emphasis, how frequently to try certain kinds of cuts and so on), CP Optimizer may offer the user even more flexibility in how to tailor searches I am not yet familiar enough to do anything beyond trying to indicate which variables should be the first candidates for branching (and I'm not sure I got that right).

I'll end this installment with a synopsis of some runs.

• In a 30 minute run using the new setting 5 for the MIP emphasis parameter (stressing use of heuristics to improve the incumbent, at the cost of not paying too much attention to the lower bound), CPLEX found a feasible schedule with 14 jersey changes and a lower bound of 2.97 (a 79% gap). It evaluated somewhere around 18,500 nodes.
  
• In a 30 minute run, CP Optimizer branched over 957 million times but never found a feasible schedule and ended with a best bound of 0.
• Finally, I tried running CPLEX for three minutes (in which it found an incumbent with 19 jersey changes) and then used that incumbent as a starting solution for a one hour run of CP Optimizer. CPO accepted the incumbent solution, did a bit over 1.55 billion branch operations, and failed to improve on it. The best bound was again 0. 

If you want to look at my code (which will make slightly more sense after my next couple of posts, where I will describe the models), it is available from my university GitLab repository.

A CP Model for Toasting Bread

A question on Mathematics Stack Exchange deals with a problem (apparently from the book "Thinking Mathematically") about toasting bread on a grill. The grill can hold two slices at a time and can only toast one side of each at a time. You have three slices to toast, and the issue is to figure out how to do it in the minimum possible amount of time (what operations management people refer to as the makespan).

The questioner had a solution that I was able to prove is optimal, using a constraint programming (CP) model that I coded using the Java API to IBM's CPOptimizer (part of the CPLEX Studio product). I won't swear my model is elegant or efficient, since I'm pretty new to CPO, but I think it is correct. If anyone getting started with CPO and the Java API wants to see the source code, it is available in my repository. I'll describe a few key aspects below.

I assumed in my model that there is a single cook. The fundamental components of the model are CPO "interval variables" (instances of IloIntervalVar) for each task (inserting a slice, toasting one side, removing a slice, flipping a slice) along with a dummy task for being done and a placeholder task I called "reversing". Interval variables represent time spans during which tasks are done.

In the problem, there are two ways to get from toasting the front of a slice to toasting the back: you can leave it on the grill and flip it; or you can remove it and (later) replace it with the other side up. Since I didn't know a priori which slices will be handled which way, I created interval variables for removing each slice after the first side, replacing each slice with the second side up, and flipping each slice. Those variables are declared optional, meaning each interval may or may not show up in the solution. For each slice, there is an interval variable for the "reversing" task that is not optional. Each slice has to be reversed, one way or the other. The tasks for replacing a slice (after removing it) and for flipping the slice are listed as alternatives for the reversing task for that slice, which means exactly one of flipping or replacing must be in the solution. Separate constraints ensure that a slice is reinserted if and only if it was removed after the first side toasted. Those constraints use the IloCP.presenceOf function to test whether the remove and reinsert intervals are present, and set the results equal (so both present or neither present).

The sequencing of operations (insert, toast first side, reverse, toast second side, remove) is enforced through a bunch of constraints that use IloCP.endBeforeStart (which says the first argument has to end before the second argument starts). The dummy "done" task is sequenced to start only after each slice has been removed for the final time. I'm pretty sure the objective value (the time everything is done) could be handled other ways, but I tend to think of completion as a zero-length task.

The cook can only do one thing at a time. This is handled using the IloCP.noOverlap function. It is passed a list of every interval that requires the cook's attention (everything but the actual toasting tasks), and prevents any of those intervals from overlapping with any other.

Finally, I need to prevent more than two slices from occupying the grill at any one time. The noOverlap function is no help here. Instead, I use an instance of IloCumFunctionExpr, which represents a function over time that changes as intervals begin and end. In this case, the function measures occupancy.  This is handled by treating the usage as a combination of step functions (IloCP.stepAtStart and IloCP.stepAtEnd). Usage steps up by one at the start of the task of inserting a slice and steps down by one at the end of the task of removing a slice. (Toasting and flipping have no effect on occupancy.) The Javadoc for the relevant functions a bit, um, sparse, essentially saying nothing about the step height argument. Thus I discovered the hard way (through error messages) that adding a step with height -1 when a slice was removed was not acceptable. Instead, I had to subtract a step of height +1.

Although it was not really necessary on a problem this small, I removed some symmetry resulting from the arbitrary ordering of the slices by setting precedence constraints saying that slice 1 is started before slice 2, which in turn is started before slice 3.

It is possible to model the problem as an integer program, and in fact I initially leaned that direction. The IP model, however, would be bulkier and less expressive (which would make it more prone to logic errors), and quite possibly would be slower to solve. CPOptimizer is designed specifically with scheduling problems in mind, so it is the better tool for this particular job.

Of Typewriters and Permutations (IV)

I'm continuing the recent theme of solving a quadratic assignment problem that lays out the 26 letters of the English alphabet on a one-dimensional "keyboard" for an 18th century typewriter. I thought this would be the last post, but something new turned up, so there will likely be one more.

In the blog post that started all this, Hardmath123 found a solution (via a heuristic) with cost 5,499,341. Using frequency data from Nate Brixius, I get a slightly higher cost (5,510,008), which is the value I will use for it (because, hey, it's my blog). Hardmath123 suspected but could not prove that this layout is optimal. I still can't verify optimality (but maybe in the next post I will ... or maybe not).

In my previous epistles on this subject, I tried out three MIP models and a quadratic (integer) program. In five minute runs using a beta copy of the next version of CPLEX, the best I was able to do was a solution with objective value 5,686,878.

Out of curiosity, I tried a constraint programming (CP) model (using CP Optimizer, the constraint solver in the CPLEX Studio suite). Constraint programming is best suited to models with integer variables (which we have here), although it can handle floating-point variables to some extent. It is well suited to logic constraints, and in particular is popular in scheduling, but it's not my go-to tool in part because it tends to have very weak bounding.

Defining the constraint programming model for the keyboard problem (denoted "CP" in my code) is fairly straightforward, but rather different from any of the MIP/QP models. We use general integer variables rather than binary variables to determine the position assignments. So $p_i\in \{0,\dots,25\}$ is the position of symbol $i$ in the layout, $s_j \in \{0,\dots,25\}$ is the symbol at position $j$, and $d_{ik}\in \{0,\dots,25\}$ is the distance between symbols $i$ and $k$. (Recall that we are using zero-based indexing.) As before, the objective is$$\min \sum_i \sum_j f_{ij}d_{ij},$$where $f_{ij}$ is the frequency with which symbol $j$ follows symbol $i$.

Where it gets interesting is in specifying the constraints. Although CP is designed for logic constraints ("or", "and", "if-then"), it's real power comes from what they call "global constraints", specialized constraints that tie multiple variables together. Arguably the most common global constraint is "all-different", meaning that in a group of variables no two of them can take the same value. It is the key to defining permutations, and in our model we add all-different constraints for both $p$ (no two symbols have the same position) and $s$ (no two slots have the same symbol). One or the other of those is technically redundant, but redundant constraints in a CP model can help the solver, so it's probably worth including both.

Another global constraint provided by CP Optimizer is "inverse". CP models allow variables to be used to index other variables, something we cannot do directly in math programming modes. The inverse constraint, applied to two vector variables $u$ and $v$, says that $u_{v_i}=i$ and $v_{u_i}=i$. In our model, we can apply the inverse constraint to $p$ and $s$. In effect, it says that the symbol in the position occupied by symbol $i$ is $i$, and the position of the symbol in position $j$ is $j$.

Finally, we can define the distance variables using the position variables:$$d_{ij} = |p_i - p_j|\quad \forall i,j.$$ How does the CP model perform? As I expected, the bounding is pretty bad. In a five minute run, the lower bound only makes it to 447,795 (91.87% gap). The good news is that the CP model finds Hardmath123's solution with objective value 5,510,008 ... and finds it in about 22 seconds on my PC! This is using the default search strategy; I did not assign branching priorities.

The take-away here, from my perspective, is that in some problems a CP model can be a great way to generate an optimal or near-optimal incumbent, which you can use as a starting solution in an optimization model if you need to prove optimality.

A CP Model for the Weight Problem of Bachet de Meziriac

Fellow OR blogger Erwin Kalvelagen just posted a MINLP model for a puzzle apparently posed by French mathematician Claude Gaspard Bachet de Méziriac. You can read the puzzle statement on Erwin's blog; a brief synopsis is that you are looking to find four integer weights that sum to 40, such that any object with integer weight from 1 to 40 can be weighed using them. I assume (and Erwin apparently agrees) that means weight on a balance scale, where weights can be placed either in the same pan as the object being weighed or in the opposite pan.

Just for amusement, here's a MiniZinc constraint programming model for the problem. It yields (unsurprisingly) the same unique solution that Erwin's MINLP model produces.

% Bachet de Meziriac problem

include "globals.mzn";

% sizes of the pieces
array[1..4] of var 1..40: x;

% pieces used in subtotals
array[1..40, 1..4] of var -1..1: y;

% order constraints
constraint increasing(x);

% overall sum
constraint sum(x) = 40;

% partial sums
constraint forall (i in 1..40) (sum([ x[j]*y[i,j] | j in 1..4 ]) = i);

% search
solve::int_search(x, largest, indomain_max, complete) satisfy;

Printing a CP Optimizer Model in Java

Update: I've converted the code to a static method (so that you do not need to extend the IloCP class), and added to a utility package I wrote. Please see this post for details on the software and how to get it.

The IBM ILOG CP Optimizer constraint programming solver has a Java API that currently seems to provide no easy way to print out a model. Here's a hack that prints the key ingredients (objectives and constraints). It does not explicitly list variables (they're embedded in the constraints and objectives), but that would be an easy fix. It also does not print out parameter settings. The code as written overrides the .toString() method for a class containing an instance of IloCP named 'solver'.

  @Override
  public String toString() {
    if (solver == null) {
      return "No model yet!";
    } else {
      StringBuilder sb = new StringBuilder();
      Iterator it = solver.iterator();
      while (it.hasNext()) {                           // iterate over the model
        IloAddable object = (IloAddable) it.next();    // get the next object
        if (object instanceof IloObjective) {
          // For an objective, .toString() returns the body of the objective.
          sb.append(object.getName())
            .append(": ")
            .append(object)
            .append("\n");
        } else if (object instanceof IloConstraint) {
          /*
            For some constraints, .toString() returns the body; for others
            (notably inequalities), .toString() returns the name if there is
            one, the body if the name is null.
          */
          IloConstraint ct = (IloConstraint) object;
          String name = ct.getName();             // get the name
          ct.setName(null);                       // temporarily remove the name
          sb.append(name)
            .append(": ")
            .append(ct)                           // this now gets the body
            .append("\n");
          ct.setName(name);                       // restore the name
        } else {
          sb.append(object.getClass().toString()).append("\n");
        }
      }
      return sb.toString();
    }
  }

Indexing an Array with a Variable

Sometimes, in the context of a mathematical program, someone wants to use a variable $x$ to index a vector $y$, as in \[z = y[x].\] As a starting point, we should assume that $x$ is an integer variable whose domain equals, or at least is a subset of, the index set of $y$. If you try to set $z=y[2.71828]$, or $z=y[5]$ when $y$ is indexed by ${1,\dots,4}$, you should expect Bad Things To Happen.

With that stipulation, $z=y[x]$ poses no problem in a constraint program. It cannot, however, be expressed directly in a mathematical program. If the domain of $x$ is not too large, it can be implemented somewhat obliquely in a mathematical program using binary variables.

Let's assume that $y$ is indexed over $1,\dots,N$ (with $N$ known at the outset), and that $y$ is a parameter. (We'll treat the case where $y$ is a variable in a minute.) Create $N$ binary variables $x_1,\dots,x_N$ and add the constraint \[x_1+\dots+x_N=1,\] which makes $\{x_1,\dots,x_N\}$ a type 1 special ordered set (SOS1). Then define $z$ via \[z=\sum_{i=1}^N y[i]*x_i.\]

You can do exactly this when $y$ is a vector of variables, but it adds a nonconvex quadratic constraint to the model, which likely prevents you from finding a guaranteed optimum, and greatly restricts what algorithms/software you can use. If we know a priori lower and upper bounds for $y$, say \[L_i \le y[i] \le U_i \forall i\in {1,\dots,N}\] with $L_i$ and $U_i$ parameters, we can add the $x_i$ as above and use the following set of constraints to define $z$: \[\begin{eqnarray*}y[i] -z \le (U_i - \hat{L})(1-x_i)  \forall i\in {1,\dots,N}\\ y[i] - z \ge (L_i - \hat{U})(1-x_i)  \forall i\in {1,\dots,N}\end{eqnarray*}\]where $\hat{L}=\min_{i=1,\dots,N} L_i$ and $\hat{U}=\max_{i=1,\dots,N} U_i$.

Numberjack

A response over at OR-Exchange led me to discover Numberjack, an open-source (LPGL) constraint programming platform written in Python.  I haven't had time to look at it carefully, and won't for a while (among other things, I need to learn more Python first), but I had to give them a shout-out just for the slogan on their home page (superimposed on a photograph of a forest): "Cuts your exponential search tree into logs."

Time Limits in Choco Solver

I'm introducing myself to constraint programming as part of a research project. Since I program mainly in Java (I swore off C++ after too much time swearing at it), I went looking for an open-source CP library written in Java.  My second attempt was Choco, which seems to work quite well for my purposes. Like many non-commercial software projects, though, the documentation lags a bit behind the code. For some reason, there's no Javadoc available (at least short of compiling the source code myself, which to date I've been too lazy to do).

What prompts me to write this reminder to myself is the adventure I had with time limits. The Choco documentation (a PDF file dated 01/20/2010, so it's pretty current as of this writing) mentions two methods, setTimeLimit and setCpuTimeLimit. NetBeans refuses to believe that the setCpuTimeLimit method (or anything starting with 'setCpu') exists. That's fine; setTimeLimit should work for me.

Unfortunately, my code started generating an unlabelled exception, which I could not catch even if I surrounded the solve() method (which lists no thrown exceptions) with a try ... catch block. Running the solve() method inside a SwingBackgroundTask probably didn't enhance my debugging efforts, either. I eventually deduced the exception was a consequence of hitting the time limit, but had a heck of time figuring out what was throwing it.

The solve() method returns a Boolean (not boolean) result. Eventually I paid attention to that and figured it out: if the solver hits its time limit, it returns null rather than true or false. I was using the return value in a conditional, as in

if (solver.solve()) { // do something with solution }

and that was causing a null pointer exception if the return value was null.  Live and learn.