A Bicriterion IP Model for a Game Strategy

A question on Mathematics Stack Exchange asks how to solve an optimization problem related to a game apparently called "speedrunner". The problem boils down to selecting integer values for variables $n_a,$ $n_c$ and $n_{b,i}$ ($i=1,\dots,n_c$) so as to minimize elapsed time $$t_{\textrm{total}} =t_{a}n_{a}+\sum_{i=1}^{n_{c}}t_{b}n_{b,i}+t_{c}n_{c},$$

where  $t_a$, $t_b$ and $t_c$ are parameters with specified values. The sole constraint specified is that winnings, defined by

$$x_{\textrm{total}} = x_{0}n_{a}\prod_{i=1}^{n_{c}}n_{b,i}$$(where $x_0$ is another parameter), must be at least $\$1$ billion. We can infer from the requirement of positive earnings that $n_a \ge 1,$ $n_c \ge 1,$ and $n_{b,i} \ge 1$ for $i=1,\dots,n_c.$ For modeling purposes, we will allow the index $i$ to exceed $n_c$ and require that $n_{b,i}=0$ for $i \gt n_c.$ An integer linear programming model would be a no-brainer were it not for two things: the product form of the expression for winnings is nonlinear; and it involves the product of a variable number of variables. 

 

The approach I chose was to express the original variables in terms of binary variables. To do that, we first need upper bounds on the original variables. We get those by guessing an upper bound $T$ on $t_{\textrm{total}}.$ (If we guess too low, we will either get a solution that exceeds that upper bound or the model will become infeasible, either of which will tip us off to try a larger guess.) Once we have that upper bound, we get an upper bound for each variable by setting the other two variables as small as possible. That leads to the following upper bounds:$$N_{a} =\left\lfloor \frac{T-t_{b}-t_{c}}{t_{a}}\right\rfloor$$

$$N_{b} =\left\lfloor \frac{T-t_{a}-t_{c}}{t_{b}}\right\rfloor$$

and 

$$N_{c} =\left\lfloor \frac{T-t_{a}}{t_{c}+t_{b}}\right\rfloor .$$ 

The denominator in the third equation arises from the fact that adding 1 to $n_c$ not only directly costs $t_c$ units of time but also adds another $n_{b,i}$ variable with minimum value 1, costing $t_b$ units of time.

 

Armed with those upper bounds, we can introduce binary variables $y_j$ $(j=1,\dots,N_a),$ $w_j$ $(j=1,\dots, N_c)$ and $z_{i,j}$ $(i=1,\dots,N_c$ and $j=0,\dots,N_b)$ along with constraints making them type 1 special ordered sets (i.e., each set contains exactly one variable with value 1) and expanding the original variables in terms of them. Specifically:

• $n_{a}=\sum_{j=1}^{N_{a}}j\cdot y_{j}$ with constraint $\sum_{j=1}^{N_{a}}y_{j}=1;$ 

• $n_{c}=\sum_{j=1}^{N_{c}}j\cdot w_{j}$ with constraint $\sum_{j=1}^{N_{c}}w_{j}=1;$ and

• $n_{b,i}=\sum_{j=0}^{N_{b}}j\cdot z_{i,j}$ for $i=1,\dots,N_{c}$

  with the following constraints for all $i\in\left\{ 1,\dots,N_{c}\right\}:$

  • $\sum_{j=0}^{N_{b}}z_{i,j}=1$ (the value of $n_{b,i}$ is uniquely defined); and
    
• $z_{i,0}+\sum_{k=i}^{N_{c}}w_{k}=1$ ($n_{b,i}=0$ if and only if $n_{c}<i$).

Armed with all that, the original objective function can be expressed as minimizing$$t_{a}n_{a}+t_{b}\sum_{i=1}^{N_{c}}n_{b,i}+t_{c}n_{c},$$where the only tweak is that we now sum over a constant number $N_c$ of terms rather than a variable number $n_c,$ relying on the fact that $n_{b,i}=0$ for $i>n_c.$

 

That brings us to the nonlinear formula for winnings. The original constraint is $ x_{\textrm{total}} \ge 10^9,$ which we can rewrite as $\log(x_{\textrm{total}}) \ge \log(10^9)$ using whatever flavor logarithm you like. Expanding the logs of $n_a$ and $n_{b,i}$ in terms of the binary variables, we get$$\log(x_{\textrm{total}}) = \log(x_{0})+\sum_{j=1}^{N_{a}}\log(j)\cdot y_{j} +\sum_{i=1}^{N_{c}}\sum_{j=1}^{N_{b}}\log(j)\cdot z_{i,j}.$$So the constraint that log winnings equal or exceed the log of $10^9$ is linear in the binary variables.

 

The solution in the Math Stack Exchange post $(n_a = 7,$ $n_c = 3$ and $n_{b,1}=n_{b,2}=n_{b,3}=18)$ turns out to be optimal given the parameter values in the question, and the IP model here will correctly achieve the minimum time (approximately 52.3 minutes) ... but it will not necessarily reproduce the best winnings within that time (approximately $\$1.02$ billion). That is because there are multiple feasible solutions that hit the correct time value, with $n_a = 7$ and with $n_c = 3$ but with the 54 cumulative $n_b$ units allocated differently (for instance, $\{19, 19, 16\}$ rather than $\{18, 18, 18\}.$ To get the best possible solution, we can use a bicriterion model with lexicographically ordered objectives, first minimizing time and then maximizing winnings (by minimizing the negative of the log of winnings).

I tested that model using the Java API to CPLEX 22.1. My code is available from my university repository.

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

Wolf, Goat and Cabbage

On Operations Research Stack Exchange, someone asked about a possible connection between the "wolf, goat and cabbage" logic puzzle and Monge's optimal transport problem. In the logic puzzle, a farmer has to get a wolf, a goat and a cabbage across a river using a boat that can only accommodate one of them (plus the farmer) at a time. If you leave the wolf and goat alone together at any point, the wolf eats the goat. If you leave the goat and cabbage alone together at any point, the goat eats the cabbage. Fortunately, the wolf has no appetite for cabbage and the cabbage does not seem to want to eat anything, else the problem would be infeasible.

Neither Monge's transport problem nor the more commonly taught (in my opinion) Hitchcock transportation problem directly apply, although you can (almost) treat the puzzle as a multiperiod commodity flow with an "inventory" of each item (wolf, goat, cabbage) on each side of the river. The "almost" part is that you need some of the variables to be integer, for two reasons. One is that the LP relaxation of the logic constraints (e.g., inventory of wolf + inventory of goat $\le 1$ on this side of the river at this time) will result in fractional values (we'll leave half a wolf and half a goat here and carry the other halves across the river) (which would greatly diminish the values of both wolf and goat). The other is that the objective is to minimize the number of trips made. It would be tempting to just assign a cost of 1 to each movement of an object in either direction, but the catch is that you will occasionally cross with nothing in the boat (besides the farmer). Those "deadheading" trips count toward the objective, but it's tricky to assign a cost to a zero flow.

To fill in some idle time, I coded up a MIP model. Mind you, I am not advocating MIP as a way to solve problems like this; I just wanted to confirm my thinking (in particular, that an LP commodity flow model would have fractional solutions). Assume that the farmer arrives at the left bank at time 0 with all three items, and that each trip (in either direction) counts as one time unit, with the first trip occurring at time 1. We need to set an upper bound $T$ on the number of trips. Since the state of the system is described by the location of four things (counting the farmer), with each have two possible locations (left bank, right bank), $T=2^4 =16$ works. The set of items will be denoted $I=\lbrace W, G, C\rbrace.$ My formulation uses the following variables.

• $L_{i,t}\in [0,1]$ and $R_{i,t}\in [0,1]$ are the inventories of item $i\in I$ at time $t\in \lbrace 0,\dots,T$ on the left and right banks respectively, after any trip occurring in that period.
• $x_{i,t}\in \lbrace 0,1 \rbrace$ and $y_{i,t}\in \lbrace 0,1 \rbrace$ are the amount of item $i$ crossing the river at time $t$ from left to right or right to left respectively.
• $z_t\in \lbrace 0,1 \rbrace$ is 1 if transport is ongoing and 0 if we are done (the farmer and all three items are on the right bank).

It would be perfectly fine (but unnecessary) to make the inventory variables integer-valued, and we could also make the inventory variables integer-valued and drop the integrality restrictions on the cartage variables ($x$ and $y$).

Some of the variables can be fixed at the outset.

• We start with all inventory on the left bank, so $L_{i,0}=1$ and $R_{i,0}=0$ for all $i\in I.$ 
• There is no trip at time 0, so $z_0=0$ and $x_{i,0}=y_{i,0}$ for all $i\in I$.
• Trips alternate left-to-right (odd time periods) and right-to-left (even time periods), so $x_{i,t}=0$ for all $i\in I$ and for all even $t$ and $y_{i,t}=0$ for all $i\in I$ and for all odd $t$.
  

The objective function is to minimize the number of trips required. $$\min \sum_{t=1}^T z_t.$$

The constraints are rather straightforward.

• The inventory on either bank in any period is the inventory on that bank from the preceding period, plus any arriving inventory and minus any departing inventory. So for $t\ge 1$ $$L_{i,t} = L_{i, t-1} - x_{i,t} + y_{i,t}\quad \forall i\in I$$ and $$R_{i,t} = R_{i,t-1} + x_{i,t} - y_{i,t}\quad \forall i\in I.$$
• In an odd numbered period (where the farmer ends up on the right bank), neither wolf and goat nor goat and cabbage can be on the left bank. So for odd $t$ $$L_{W,t} + L_{G,t} \le 1$$ and $$L_{G,t} + L_{C,t}\le 1.$$
• In an even numbered period (where the farmer ends up on the left bank), neither wolf and goat nor goat and cabbage can be on the right bank unless the problem is completed ($z_t = 0$), in which case the farmer remains on the right bank. So for even $t$ $$R_{W,t}+R_{G_t} + z_t \le 2$$ and $$R_{G,t} + R_{C,t} + z_t \le 2.$$
• Transport continues until the left bank is empty. $$3z_t \ge \sum_{i\in I} L_{i,t - 1}\quad \forall t\ge 1.$$

It does indeed produce a correct solution, using seven trips (see the Wikipedia page for the solution) ... and with integrality conditions dropped it does indeed produce a nonsense solution with fractions of items being transported.

Java code for this model (using CPLEX) is in my Git repository.