Yet Another Math Programming Consultant: Wolf, goat and cabbage problem: MIP and network model

The Wolf, goat, and cabbage problem can be stated as [1]:

A farmer with a wolf, a goat, and a cabbage must cross a river by boat. The boat can carry only the farmer and a single item. If left unattended together, the wolf would eat the goat, or the goat would eat the cabbage. How can they cross the river without anything being eaten?

A long transcontinental flight was a good opportunity to try to attack this problem. Here are some approaches to model this. Not at all a very useful or practical model, but still interesting, I think (although I may be in a small minority here).

Inventory model

In this model, we keep track of the inventory of items (wolf, goat, cabbage). We assume the starting inventory is: all items are on the left bank of the river. The final inventory should be: all items are on the right bank. In this model, I assume the following numbering scheme:

----     31 SET trip  trips

trip1 ,    trip2 ,    trip3 ,    trip4 ,    trip5 ,    trip6 ,    trip7 ,    trip8 ,    trip9 ,    trip10


----     31 SET dir  direction of trip

L->R,    R->L


----     31 SET tripDir  trip direction combos

              L->R        R->L

trip1          YES
trip2                      YES
trip3          YES
trip4                      YES
trip5          YES
trip6                      YES
trip7          YES
trip8                      YES
trip9          YES
trip10                     YES

We don't know in advance how many trips we need. In a MIP model, however, we need to statically size the problem in advance. We hope 10 trips is more than enough. If it is too small, the model will be infeasible. If too large, the MIP model is larger and more difficult to solve than needed.

Let's first model the passengers on each trip. We have 0 or 1 of them.

$\begin{align}&\color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\in\{0,1\}\\ &{\mathit{item}}\in\{{\mathit{wolf}},{\mathit{goat}},{\mathit{cabbage}}\}\end{align}$

The capacity constraint is obvious:

$\sum_{\mathit{item}} \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}} \le 1 \>\>\forall {\mathit{trip}}$

The inventory is handled using an inventory balance equation, something we often see in production scheduling models:

$\begin{align}&\color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}} \in \{0,1\}\\ &\mathit{side}\in\{L,R\} \\ & \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}} = \color{darkred}{\mathit{inv}}_{\mathit{trip}-1,\mathit{side},\mathit{item}} + \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\cdot I_{\mathit{arrival}}(\mathit{trip},\mathit{side}) - \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\cdot I_{\mathit{departure}}(\mathit{trip},\mathit{side}) + \color{darkblue}{\mathit{initialInv}}_{\mathit{item},\mathit{trip}} \end{align}$

The $\color{darkblue}{\mathit{initialInv}}_{\mathit{item},\mathit{trip}}$ has only non-zero values for the first trip. It establishes the initial inventory (all items on the left bank). The notation $I_{\mathit{arrival}}(\mathit{trip},\mathit{side})$ is an indicator function: $I_{\mathit{arrival}}(\mathit{trip},\mathit{side})= \begin{cases} 1& \text{if this is an arrival}\\ 0&\text{otherwise}\end{cases}$ Similar for the departures. Note that these indicator functions are defined over set elements, which are completely exogenous. We can precompute this easily, forming some useful data structure. In the GAMS model, I used a 2d set for this.

To forbid leaving certain combinations without supervision (the cases are: goat/cabbage, wolf/goat), we can use some simple linear constraints: $\sum_{\mathit{item}} \color{darkblue}{\mathit{eat}}_{\mathit{case},\mathit{item}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}}\le 1 \>\>\forall \mathit{case},\mathit{trip},\mathit{side}$ This is a bit too strict: we only want this check for departures, and not for arrivals. Also, we may want to relax this constraint when we are done (when all items have arrived on the right bank). So this becomes:

$\sum_{\mathit{item}} \color{darkblue}{\mathit{eat}}_{\mathit{case},\mathit{item}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}}\le 1 + \color{darkred}{\mathit{done}}_{\mathit{trip}} \>\>\forall \mathit{case},\mathit{trip},\mathit{side}|I_{\mathit{departure}}(\mathit{trip},\mathit{side}) =1$

The $\color{darkred}{\mathit{done}}_{\mathit{trip}}\in \{0,1\}$ variable indicates when we are done. We could say $\color{darkred}{\mathit{done}}_{\mathit{trip}} = \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{wolf}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{goat}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{cabbage}}$ We don't have to add: $\color{darkred}{\mathit{done}}_{\mathit{trip}} = \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{wolf}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{goat}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{cabbage}} \cdot (1-\color{darkred}{\mathit{inv}}_{\mathit{trip},L,\mathit{wolf}}) \cdot (1-\color{darkred}{\mathit{inv}}_{\mathit{trip},L,\mathit{goat}}) \cdot (1-\color{darkred}{\mathit{inv}}_{\mathit{trip},L,\mathit{cabbage}})$ as the inventory balance takes care of the inventory at the left bank. The binary multiplication (this is nonlinear and nonconvex) can be easily linearized as: $\begin{align}& \color{darkred}y=\color{darkred}x_1 \cdot \color{darkred}x_2 \cdot \color{darkred}x_3 \iff \begin{cases} \color{darkred}y \le \color{darkred}x_1 \\ \color{darkred}y \le \color{darkred}x_2\\ \color{darkred}y \le \color{darkred}x_3 \\ \color{darkred}y \ge \color{darkred}x_1+ \color{darkred}x_2+\color{darkred}x_3-2 \end{cases}\end{align}$

We can enforce the $\color{darkred}{\mathit{done}}_{\mathit{trip}}$ to remain 1 once we are done by: $\color{darkred}{\mathit{done}}_{\mathit{trip}} \ge \color{darkred}{\mathit{done}}_{\mathit{trip}-1}$ This constraint is not really needed (the objective takes care of this), but it is conveying what we know about optimal solutions. It can also help with performance (not really an issue here). I prefer to keep it in the model.

The objective I used is: maximize the number of "done" trips: $\max \sum_{\mathit{trip}}\color{darkred}{\mathit{done}}_{\mathit{trip}}$ This will generate a tight schedule (no unnecessary trips).

The complete model is:

Inventory based model
$\begin{align}\max&\sum_{\mathit{trip}}\color{darkred}{\mathit{done}}_{\mathit{trip}} \\ & \sum_{\mathit{item}} \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}} \le 1 && \forall {\mathit{trip}} \\ &\color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}} = \color{darkred}{\mathit{inv}}_{\mathit{trip}-1,\mathit{side},\mathit{item}} + \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\cdot I_{\mathit{arrival}}(\mathit{trip},\mathit{side}) \\ & \>\>\> - \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\cdot I_{\mathit{departure}}(\mathit{trip},\mathit{side}) + \color{darkblue} {\mathit{initialInv}}_{\mathit{item},\mathit{trip}} && \forall \mathit{trip},\mathit{side},\mathit{item} \\ & \sum_{\mathit{item}} \color{darkblue}{\mathit{eat}}_{\mathit{case},\mathit{item}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}}\le 1 + \color{darkred}{\mathit{done}}_{\mathit{trip}} && \forall \mathit{case},\mathit{trip},\mathit{side}\|I_{\mathit{departure}}(\mathit{trip},\mathit{side}) =1 \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \le \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{item}} && \forall \mathit{trip},\mathit{item} \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \ge \sum_{\mathit{item}} \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{item}} - 2 && \forall {\mathit{trip}} \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \ge \color{darkred}{\mathit{done}}_{\mathit{trip}-1} && \forall \mathit{trip} \gt 1 \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \in \{0,1\} \\ & \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}} \in \{0,1\} \\ & \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}} \in \{0,1\} \end{align}$

Inventory based model

$\begin{align}\max&\sum_{\mathit{trip}}\color{darkred}{\mathit{done}}_{\mathit{trip}} \\ & \sum_{\mathit{item}} \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}} \le 1 && \forall {\mathit{trip}} \\ &\color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}} = \color{darkred}{\mathit{inv}}_{\mathit{trip}-1,\mathit{side},\mathit{item}} + \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\cdot I_{\mathit{arrival}}(\mathit{trip},\mathit{side}) \\ & \>\>\> - \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}\cdot I_{\mathit{departure}}(\mathit{trip},\mathit{side}) + \color{darkblue} {\mathit{initialInv}}_{\mathit{item},\mathit{trip}} && \forall \mathit{trip},\mathit{side},\mathit{item} \\ & \sum_{\mathit{item}} \color{darkblue}{\mathit{eat}}_{\mathit{case},\mathit{item}} \cdot \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}}\le 1 + \color{darkred}{\mathit{done}}_{\mathit{trip}} && \forall \mathit{case},\mathit{trip},\mathit{side}|I_{\mathit{departure}}(\mathit{trip},\mathit{side}) =1 \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \le \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{item}} && \forall \mathit{trip},\mathit{item} \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \ge \sum_{\mathit{item}} \color{darkred}{\mathit{inv}}_{\mathit{trip},R,\mathit{item}} - 2 && \forall {\mathit{trip}} \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \ge \color{darkred}{\mathit{done}}_{\mathit{trip}-1} && \forall \mathit{trip} \gt 1 \\ & \color{darkred}{\mathit{done}}_{\mathit{trip}} \in \{0,1\} \\ & \color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}} \in \{0,1\} \\ & \color{darkred}{\mathit{inv}}_{\mathit{trip},\mathit{side},\mathit{item}} \in \{0,1\} \end{align}$

The results of the model are:

----    116 VARIABLE pax.L  passengers: items taken on each trip

              wolf        goat     cabbage

trip1                        1
trip3                                    1
trip4                        1
trip5            1
trip7                        1


----    116 VARIABLE inv.L  inventory just after trip

            L.wolf      L.goat   L.cabbage      R.wolf      R.goat   R.cabbage

trip1            1                       1                       1
trip2            1                       1                       1
trip3            1                                               1           1
trip4            1           1                                               1
trip5                        1                       1                       1
trip6                        1                       1                       1
trip7                                                1           1           1
trip8                                                1           1           1
trip9                                                1           1           1
trip10                                               1           1           1


----    116 VARIABLE done.L  if 1 we are done

trip7  1,    trip8  1,    trip9  1,    trip10 1


----    116 VARIABLE z.L                   =        4.000  objective

We can generate a slightly more meaningful report:

----    132 PARAMETER trace  trip report

                          L.wolf      L.goat   L.cabbage      R.wolf      R.goat   R.cabbage

initial.    .                  1           1           1
trip1  .L->R.goat              1                       1                       1
trip2  .R->L.                  1                       1                       1
trip3  .L->R.cabbage           1                                               1           1
trip4  .R->L.goat              1           1                                               1
trip5  .L->R.wolf                          1                       1                       1
trip6  .R->L.                              1                       1                       1
trip7  .L->R.goat                                                  1           1           1

On the left we see the passengers taken along each trip. On the right we have the inventory just after each trip. Only after trip 7, we are done.

This is not the only solution. Basically, we can reverse the schedule of the wolf and the cabbage. A no-good cut can be added to the model to forbid a previous solution (but allow all other solutions). For an optimal binary solution $\color{darkblue}x^*_i$ , such a cut looks like $\sum_{i|\color{darkblue}x^*_i=1} x_i - \sum_{i|\color{darkblue}x^*_i=0} x_i \le \sum_{i|\color{darkblue}x^*_i=1} 1 - 1$ In our case, we can apply this to just the variable $\color{darkred}{\mathit{pax}}_{\mathit{trip},\mathit{item}}$ . The variable $\color{darkred}{\mathit{inv}}$ is completely determined once we know $\color{darkred}{\mathit{pax}}$ , so we can ignore it here.

If we add this no-good cut, we indeed see there are two different passenger schedules:

----    116 VARIABLE pax.L  passengers: items taken on each trip

              wolf        goat     cabbage

trip1                        1
trip3                                    1
trip4                        1
trip5            1
trip7                        1

----    160 VARIABLE pax.L  passengers: items taken on each trip

              wolf        goat     cabbage

trip1                        1
trip3            1
trip4                        1
trip5                                    1
trip7                        1

All in all, this problem was not super easy to model. Of course, that allows us to discuss quite a few formulation tricks for a rather simple problem.

GAMS Model

$onText

A farmer with a wolf, a goat, and a cabbage must cross a river by boat. The boat can carry
only the farmer and a single item. If left unattended together, the wolf would eat the goat,
or the goat would eat the cabbage. How can they cross the river without anything being eaten?

https://en.wikipedia.org/wiki/Wolf,_goat_and_cabbage_problem

$offText


*-------------------------------------------------------------------------
* data
*-------------------------------------------------------------------------

sets
** basic sets
   dummy 'for ordering in displays' /initial/
   trip 'trips' /trip1*trip10/
   item /wolf,goat,cabbage/
   side 'left/right bank' /L,R/
   dir 'direction of trip' /'L->R','R->L'/
   case 'forbidden item combinations' /case1,case2/

** derived sets
   tripDir(trip,dir) 'trip direction combos'
   arrival(trip,side) 'this is an arrival'
   departure(trip,side) 'this is a departure'
;

* trips 1,3,5,.. are L->R, trips 2,4,6,..  are R->L
tripDir(trip,'L->R') = mod(ord(trip),2) = 1;
tripDir(trip,'R->L') = mod(ord(trip),2) = 0;
display trip,dir,tripDir;

* arrival or departure
arrival(trip,'L') = tripdir(trip,'R->L');
arrival(trip,'R') = tripdir(trip,'L->R');
departure(trip,side) = not arrival(trip,side);
display arrival,departure;

* all items are on the left bank initially
parameter InvInitial(item,side,trip) 'start inventory before trip1' /
    (wolf, goat, cabbage).L.trip1   1
/;

* the final state is to have all items on the right bank
parameter InvTarget(item,side) 'target inventory' /
    (wolf, goat, cabbage).R   1
/;
    
* these combinations can not be left alone
table Eat(case,item) 'combinations not allowed'
        wolf  goat  cabbage
case1     1     1     
case2           1      1
;
 
*-------------------------------------------------------------------------
* Model
*-------------------------------------------------------------------------

variables
   pax(trip,item)      'passengers: items taken on each trip'
   inv(trip,side,item) 'inventory just after trip'
   done(trip)          'if 1 we are done'
   z                   'objective'
;
binary variable pax,inv,done;

equations
   countPassengers(trip)       'max number of passengers/items on a trip'
   invBal(side,item,trip)      'inventory balance after trip'
   eating(trip,side,case)      'some combinations of items are forbidden when unattended'
   isDone1(trip,item)          "done = prod(item,inv(trip,'R',item))" 
   isDone2(trip)               "done = prod(item,inv(trip,'R',item))"
   ordering(trip)              'done(trip)>=done(trip-1)' 
   obj                         'objective'
;

* capacity constraint 
countPassengers(trip).. sum(item,pax(trip,item)) =l= 1;

* inventory just after trip
invBal(side,item,trip)..
    inv(trip,side,item) =e= inv(trip-1,side,item) + pax(trip,item)$arrival(trip,side) - pax(trip,item)$departure(trip,side) + invInitial(item,side,trip); 

* forbidden combinations
eating(departure(trip,side),case)..
    sum(item, Eat(case,item)*inv(trip,side,item)) =l= 1+done(trip);

* done = prod(item,inv(trip,'R',item))
isDone1(trip,item).. done(trip) =l= inv(trip,'R',item);
isDone2(trip).. done(trip) =g= sum(item,inv(trip,'R',item))-2;

* done(trip)>=done(trip-1)
* ordering(trip-1).. makes sure we only generate the constraints that are needed
ordering(trip-1).. done(trip) =g= done(trip-1);

* maximize number of "done" trips
obj.. z =e= sum(trip, done(trip));

* optional
done.fx(trip)$(ord(trip)=card(trip)) = 1;
inv.fx(trip,side,item)$(ord(trip)=card(trip)) = invTarget(item,side);

model m /all/;
solve m maximizing z using mip;
abort$(m.modelstat <> %modelStat.optimal% and m.modelstat <> %modelStat.integerSolution%) "No solution";

*-------------------------------------------------------------------------
* raw results
*-------------------------------------------------------------------------

option pax:0,inv:0:1:2,done:0;
display pax.l,inv.l,done.l,z.l;

*-------------------------------------------------------------------------
* more meaningful trip report
*-------------------------------------------------------------------------

alias(item,item2);

parameter trace(*,*,*,side,item) 'trip report';
trace('initial','','',side,item) = invInitial(item,side,'trip1');
loop(trip,
   trace(tripDir(trip,dir),item,side,item2)$(pax.l(trip,item)>0.5) = inv.l(trip,side,item2);
   trace(tripDir(trip,dir),'',side,item2)$(sum(item,pax.l(trip,item))<0.5) = inv.l(trip,side,item2);
   break$(done.l(trip)>0.5);
);
option trace:0:3:2;
display trace;

*-------------------------------------------------------------------------
* find second solution
* for binary variable x(i):
*   ecut(cut).. sum(i$sol(i,cut),x(i)) - sum(i$(sol(i,cut)=0),x(i)) =l= sum(i$sol(i,cut),1) - 1;
*-------------------------------------------------------------------------

sets
   cut 'static set' /cut1*cut2/ 
   dcut(cut) 'dynamic subset'
   ti(trip,item) 'trip-item combos'
;
ti(trip,item) = yes;

parameter sol(trip,item,cut) 'earlier solutions';

equation ecut(cut) 'no-good cut';

ecut(dcut).. sum(ti$sol(ti,dcut),pax(ti)) - sum(ti$(sol(ti,dcut)=0),pax(ti)) =l= sum(ti$sol(ti,dcut),1) - 1;

* find second solution with same objective value
model m2/m,ecut/;
dcut('cut1') = yes;
sol(ti,'cut1') = round(pax.l(ti));
z.fx = round(z.l);
solve m2 maximizing z using mip;
abort$(m2.modelstat <> %modelStat.optimal% and m2.modelstat <> %modelStat.integerSolution%) "No second solution";
display pax.l;

* there should be no third solution: this model is infeasible 
dcut('cut2') = yes;
sol(ti,'cut2') = round(pax.l(ti));
solve m2 maximizing z using mip;
abort$(m2.modelstat = %modelStat.optimal% or m2.modelstat = %modelStat.integerSolution%) "Expected model to be infeasible";

Shortest path formulation

Here, we use a very different formulation. We'll try to formulate this as a shortest-path network model. Let's start with the network nodes. The nodes form a state that has two parts:

The location of the node (left or right). We can also interpret this as where the boat (or the farmer) is at this point in time.
The inventory. I store both the left and right inventory, but we could save a bit of space by just storing one of them.

Here is an initial list of nodes:

----     52 SET state  state: node location + inventory

              L.wolf      L.goat   L.cabbage      R.wolf      R.goat   R.cabbage

node1 .L         YES         YES         YES
node2 .L         YES         YES                                             YES
node3 .L         YES                     YES                     YES
node4 .L         YES                                             YES         YES
node5 .L                     YES         YES         YES
node6 .L                     YES                     YES                     YES
node7 .L                                 YES         YES         YES
node8 .L                                             YES         YES         YES
node9 .R         YES         YES         YES
node10.R         YES         YES                                             YES
node11.R         YES                     YES                     YES
node12.R         YES                                             YES         YES
node13.R                     YES         YES         YES
node14.R                     YES                     YES                     YES
node15.R                                 YES         YES         YES
node16.R                                             YES         YES         YES

Not all nodes are allowed. E.g. node 9 is illegal: the farmer is on the right bank, so we can't have wolf+goat or goat+cabbage on the left. When we check all nodes, we are only left with the following list:

----     83 SET rnode  reduced node set

node1 ,    node2 ,    node3 ,    node5 ,    node6 ,    node11,    node12,    node14,    node15,    node16

The set of arcs is not the cartesian product of the nodes. We can only ship 0 or 1 item. So, we have a somewhat sparse network. Here is what I work with:

----    133 SET pax  passengers along arc

                     wolf        goat     cabbage       empty

node1 .node11                     YES
node2 .node12                     YES
node2 .node14         YES
node3 .node11                                             YES
node3 .node12                                 YES
node3 .node15         YES
node5 .node14                                 YES
node5 .node15                     YES
node6 .node14                                             YES
node6 .node16                     YES
node11.node1                      YES
node11.node3                                              YES
node12.node2                      YES
node12.node3                                  YES
node14.node2          YES
node14.node5                                  YES
node14.node6                                              YES
node15.node3          YES
node15.node5                      YES
node16.node6                      YES

I use an LP model to solve the shortest path problem:

Shortest path LP Model
$\begin{align}\min& \sum_{(i,j)\in A}\color{darkred}f_{i,j}\\ & \sum_{j\|(j,i)\in A}\color{darkred}f_{j,i} + \color{darkblue}{\mathit{supply}}_i = \sum_{j\|(i,j)\in A}\color{darkred}f_{i,j} + \color{darkblue}{\mathit{demand}}_i && \forall i\\ & \color{darkred}f_{i,j} \ge 0 \\ & \color{darkblue}{\mathit{supply}}_i = \begin{cases} 1 & \text{if $i$ is the source node} \\ 0 & \text{otherwise}\end{cases}\\ & \color{darkblue}{\mathit{demand}}_i = \begin{cases} 1 & \text{if $i$ is the sink node} \\ 0 & \text{otherwise}\end{cases} \end{align}$

Shortest path LP Model

$\begin{align}\min& \sum_{(i,j)\in A}\color{darkred}f_{i,j}\\ & \sum_{j|(j,i)\in A}\color{darkred}f_{j,i} + \color{darkblue}{\mathit{supply}}_i = \sum_{j|(i,j)\in A}\color{darkred}f_{i,j} + \color{darkblue}{\mathit{demand}}_i && \forall i\\ & \color{darkred}f_{i,j} \ge 0 \\ & \color{darkblue}{\mathit{supply}}_i = \begin{cases} 1 & \text{if $i$ is the source node} \\ 0 & \text{otherwise}\end{cases}\\ & \color{darkblue}{\mathit{demand}}_i = \begin{cases} 1 & \text{if $i$ is the sink node} \\ 0 & \text{otherwise}\end{cases} \end{align}$

The raw results are:

----    160 VARIABLE f.L  flow

             node2       node3       node6      node11      node12      node14      node16

node1                                            1.000
node2                                                                    1.000
node3                                                        1.000
node6                                                                                1.000
node11                   1.000
node12       1.000
node14                               1.000

This is a bit difficult to interpret. With a bit of postprocessing, we can show:

----    188 SET trace  results from LP model

                                   L.wolf      L.goat   L.cabbage      R.wolf      R.goat   R.cabbage

initial.node1 .      .                YES         YES         YES
trip1  .node1 .node11.goat            YES                     YES                     YES
trip2  .node11.node3 .                YES                     YES                     YES
trip3  .node3 .node12.cabbage         YES                                             YES         YES
trip4  .node12.node2 .goat            YES         YES                                             YES
trip5  .node2 .node14.wolf                        YES                     YES                     YES
trip6  .node14.node6 .                            YES                     YES                     YES
trip7  .node6 .node16.goat                                                YES         YES         YES

The part on the right is the inventory immediately after a trip (i.e. traveling along the arc). As before, we need 7 trips to move all items from the left bank to the right bank.

The LP model itself is trivial and standard, but designing and populating the network was again not super easy. Compared to the first model, we moved complexity from the model equations to data manipulation.

GAMS Model

$onText

   Same problem as WolfGoatCabbage.gms, but now as shortest path problem

$offText


*----------------------------------------------------------------------
* nodes
*----------------------------------------------------------------------

* for ordering of displays we want to have this first
set dummy /'initial'/;


sets
   side 'left or right bank' /L,R/
   item /wolf,goat,cabbage/
   node /node1*node16/
   left(node) /node1*node8/
   right(node) /node9*node16/
;

alias (side,nodeloc);

* enumerate states: location of node + inventory 
set state(node,nodeloc,side,item)  'state: node location + inventory' /
* node (or boat) is on the left bank
       (node1.L).(L.wolf,L.goat,L.cabbage)     
       (node2.L).(L.wolf,L.goat,R.cabbage)    
       (node3.L).(L.wolf,R.goat,L.cabbage)    
       (node4.L).(L.wolf,R.goat,R.cabbage)    
       (node5.L).(R.wolf,L.goat,L.cabbage)    
       (node6.L).(R.wolf,L.goat,R.cabbage)    
       (node7.L).(R.wolf,R.goat,L.cabbage)    
       (node8.L).(R.wolf,R.goat,R.cabbage)
       
* node (or boat) is on the right bank
       (node9.R).(L.wolf,L.goat,L.cabbage)     
       (node10.R).(L.wolf,L.goat,R.cabbage)    
       (node11.R).(L.wolf,R.goat,L.cabbage)    
       (node12.R).(L.wolf,R.goat,R.cabbage)    
       (node13.R).(R.wolf,L.goat,L.cabbage)    
       (node14.R).(R.wolf,L.goat,R.cabbage)    
       (node15.R).(R.wolf,R.goat,L.cabbage)    
       (node16.R).(R.wolf,R.goat,R.cabbage)    
    /
;

* display state  
option state:0:2:2;
display state;

* drop node location index to form inventory parameter
parameter inv(node,side,item) 'inventory';
inv(node,side,item) = sum(state(node,nodeloc,side,item),1); 
option inv:0:1:2;
display inv;


* not all nodes are allowed
* we could have done this by hand by dropping nodes from the
* state set. Instead, let's see how this looks like in code.
set
   combo  'combinations not allowed'  /combo1,combo2/
   forbidden(combo,item) /
         combo1.(wolf,goat)
         combo2.(goat,cabbage)
   /
   nodex(node) 'forbidden node when left unattended'
 ;

* populate nodex with forbidden nodes
* check inventory of side where node is not located 
loop((combo,nodeloc,side)$(ord(side)<>ord(nodeloc)),
  nodex(node)$(sum(forbidden(combo,item)$state(node,nodeloc,side,item),1) = 2) = yes;
);
display nodex;

set rnode(node) 'reduced node set';
rnode(node) = yes;
rnode(nodex) = no;
display rnode;



*----------------------------------------------------------------------
* arcs
* we have a directed arc from:
*  1. L->R and R->L
*  2. with 0 or 1 passengers
*----------------------------------------------------------------------

alias (rnode,n1,n2);

set
  otherside(node,node) 'nodes are at opposite side'
  arc(node,node)       'edges'
  pax(node,node,*)     'passengers along arc' 
;

otherside(left(n1),right(n2)) = yes;
otherside(right(n1),left(n2)) = yes;
display otherside;

Parameter
    diff(node,node,item) 'differences in L inventory between node1->node2'
    numdiff(node,node) 'number of differences in L inventory'
    sumdiff(node,node) 'sum of differences in L inventory'
;
diff(otherside(n1,n2),item) = inv(n1,'L',item)-inv(n2,'L',item);
numdiff(otherside(n1,n2)) = sum(item,abs(diff(n1,n2,item)));
sumdiff(otherside(n1,n2)) = sum(item,diff(n1,n2,item));

sets
   diff0(node,node) 'no difference in inventory: trip with zero passengers'
   diff1(node,node) 'one difference in inventory: trip with one passenger'
;   
diff0(otherside) = numdiff(otherside) = 0;
diff1(left(n1),right(n2)) = numdiff(n1,n2)=1 and sumdiff(n1,n2)=1;
diff1(right(n1),left(n2)) = numdiff(n1,n2)=1 and sumdiff(n1,n2)=-1;
display diff0,diff1;

* zero passengers: inventory stays the same
arc(otherside)$(numdiff(otherside)=0) = yes;

* one passenger: move item from node1->node2
arc(diff1) = yes;

* pax moved
pax(arc(diff1),item) = abs(diff(diff1,item))=1;
pax(arc(diff0),'empty') = yes;
display pax;


*----------------------------------------------------------------------
* shortest path model
*----------------------------------------------------------------------

Parameters
   supply(node)  / node1 1 /
   demand(node)  / node16 1 /
;

binary variable f(node,node) 'flow';
variable z 'objective';

equations
   flowbal(node) 'flow balance'
   obj           'objective'
;

obj.. z =e= sum(arc,f(arc));

flowbal(n1)..
   sum(arc(n2,n1),f(n2,n1)) + supply(n1) =e= sum(arc(n1,n2),f(n1,n2)) + demand(n1);

model shortestpath /all/;
solve shortestpath minimizing z using rmip;
display f.l;


*----------------------------------------------------------------------
* reporting
*----------------------------------------------------------------------

set
   trip /trip1*trip10/
   trace(*,*,*,*,side,item) 'results from LP model'
;
singleton sets
  n(node) 'current node' /node1/
  next(node) 'next node'
;

alias(item,item2);

trace('initial','node1','','',side,item)$inv('node1',side,item) = yes;
loop(trip$card(n),
     next(node) = f.l(n,node)>0.5;
* with passenger     
     trace(trip,n,next,item,side,item2)$(pax(n,next,item) and inv(next,side,item2))= yes;
* without passenger
     trace(trip,n,next,'',side,item2)$(sum(pax(n,next,item),1)=0 and inv(next,side,item2)) = yes;
     n(node) = next(node);
);
option trace:0:4:2;
display trace;

References

Wolf, goat, and cabbage problem, https://en.wikipedia.org/wiki/Wolf,_goat_and_cabbage_problem

Yet Another Math Programming Consultant

Friday, March 21, 2025

Wolf, goat and cabbage problem: MIP and network model

Inventory model

Shortest path formulation

References

No comments:

Post a Comment