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.

Finding Almost All Paths

A question posted on Stack Overflow (and subsequently deleted) led to a blog post by Erwin Kalvelagen on how to find all paths between two nodes in a directed graph (possibly with self-loops, i.e. arcs from a node to itself) subject to two constraints: no arc can be used more than once in a path; and there is an upper limit $M$ on the number of arcs used. Note that a path might visit a *node* more than once. It just cannot repeat an arc. The original question seems to have referred to an undirected graph, but Erwin's post works with directed graphs and so will I.

Erwin explored some mixed integer linear programming (MIP) models in his post, and a subsequent post on OR Stack Exchange led to more proposals of MIP models (including one from me). I also suggested that a "brute force" approach might be faster than any of the MIP models. In what follows, I will spell out both my MIP model and the brute force approach I used. Java code for both (which requires CPLEX and the Apache Commons Collections library) are in my code repository.

In what follows $A$ is the set of arcs in the graph, $s$ is the origin node for all paths, $t$ is the destination node for all paths, and $M$ is the maximum number of arcs to include in a path.

MIP model

 

The MIP model uses the following variables. 

• $u_{a}\in\left\{ 0,1\right\} $ is 1 if and only if arc $a$ is used on the path.
• $f_{a}\in\left\{ 0,1\right\} $ is 1 if and only if arc $a$ is the first arc on the path.
• $\ell_{a}\in\left\{ 0,1\right\} $ is 1 if and only if arc $a$ is the last arc on the path.
• $y_{ab}\in\left\{ 0,1\right\} $ is 1 if and only if arc $b$ immediately follows arc $a$ on the path.
• $z_{a}\in\left[0,M\right]$ will be the number of arcs preceding arc $a$ on the path (0 if $a$ is not on the path).

Some of the variables can be eliminated (fixed at 0) at the outset.

• $f_{a}=0$ if node $s$ is not the tail of arc $a.$
• $\ell_{a}=0$ if node $t$ is not the head of arc $a.$
• $y_{ab}=0$ if the head of arc $a$ is not the tail of arc $b.$

Since we are just interested in finding feasible solutions, the objective is to minimize 0.

Use of the $z$ variables mimics the Miller-Tucker-Zemlin approach to avoiding subtours in the traveling salesman problem. A common knock on the MTZ formulation for the TSP is that it has a somewhat loose continuous relaxation. Since we have a trivial objective (all feasible solutions are optimal), that is not a concern here.

The constraints are as follows.

• There must be one first arc and one last arc.
  $$\sum_{a\in A}f_{a}=1.$$$$\sum_{a\in A}\ell_{a}=1.$$
• At most $M$ arcs can be used.$$\sum_{a\in A}u_{a}\le M.$$
• An arc is used if and only if it is either the first arc or follows another arc on the path.$$f_{a}+\sum_{b\in A}y_{ba}=u_{a}\quad\forall a\in A.$$
• The last arc must be a used arc.$$\ell_{a}\le u_{a}\quad\forall a\in A. (1)$$
• The sequence value of an unused arc is 0.$$z_{a}\le Mu_{a}\quad\forall a\in A.$$
• No arc can follow the last arc.$$\ell_{a}+\sum_{b\in A}y_{ab}\le1\quad\forall a\in A. (2)$$
• If an arc is used, either it is the last arc or another arc follows it.$$\ell_{a}+\sum_{b\in A}y_{ab}=u_{a}\quad\forall a\in A. (3)$$
• If an arc $b$ follows arc $a$, the sequence number of arc $b$ must be one higher than the sequence number of arc $a.$ $$z_{a}-z_{b}+\left(M+1\right)y_{ab}\le M\quad\forall a,b\in A,a\neq b.$$

Over on OR Stack Exchange, Rob Pratt correctly pointed out that constraint (3) implies constraints (1) and (2). I've left them in the Java code anyway. The CPLEX presolver removes them automatically.

Finding all solutions 

To find all solutions, one possible approach is to solve whatever MIP model you choose, then add a "no-good" constraint that eliminates the solution just found (and only that one) and solve again, until eventually the aggregation of "no-good" constraints makes the model infeasible. What I did instead was to use the "populate" command in CPLEX, which accumulates a pool of solutions. Besides a time limit, I used two non-default parameter settings: I cranked up the solution pool capacity (the maximum number of solutions to find) to something high enough to let it find all solutions; and I set the solution pool intensity to its highest value (4), which tells CPLEX to aggressively seek out all feasible solutions.

Brute force approach

 

The brute force approach is remarkably straightforward. We will use a queue of (incomplete) paths that I will call the "to-do list". Start by creating a length one path for each arc with tail $s$ and add them to the to-do list. We now proceed in a loop until the to-do list is empty. At each iteration, we pull a path $P$ off the to-do list, identify the arcs whose tails match the head of the last arc in $P$, remove any arcs already on $P$, and for each surviving arc $a$ create a new path $P'$ by extending $P$ with $a$. If the head of arc $a$ is $t$, $P'$ is a new $s-t$ path, which we record. Either way, if $P'$ has length less than $M$, we add it to the to-do list, and eventually try to extend it further.

 

Do they work?

 

Would I be posting this if they didn't. :-) I tested both methods on the digraph from Erwin's post, which has 10 nodes and 22 arcs (with two self-loops). The source and sink are nodes 1 and 10 respectively. With $M=3$ there are four paths (which both methods find), and with $M=4$ there are nine paths (which both methods find). In both cases, brute force took about 1 ms. on my PC (using non-optimized Java code with a single thread). CPLEX times were rather stochastic, but I think it fair to say that building the models took around 55+ ms. and solving them typically took 20 ms. or more.

When I set a maximum length ($M$) of 10, things got interesting. The brute force approach found 268 paths (in about 6 ms), while CPLEX found only 33. Neither time limit nor pool size were a factor, so I assume that this is just a limitation of the aggressive setting for pool intensity. That means that to find all possible solutions using CPLEX, the solve / add constraint / solve again approach would be necessary.