Yet Another Math Programming Consultant: Towers of Hanoi: inventory and network formulation

The towers of Hanoi problem [1] is a famous puzzle demonstrating recursion. The task is to move a stack of disks from one peg to another. We can only move one disk at a time, and we need to obey the rule that larger disks can never be on top of a smaller disk. The moves for the 3 disk problem are:

Recursive algorithm

This problem can be programmed very elegantly using recursion:

# returns number of moves
def hanoi(n,A,B,C):
   if n==0: return 0
   n1 = hanoi(n-1,A,C,B)
   print(f"move disk from peg {A} to {B}")
   n2 = hanoi(n-1,C,B,A)
   return n1+1+n2
   
N = 3
print(f"N={N}")
cnt = hanoi(N,'A','B','C')
print(f"{cnt} moves")

The three steps to move $n$ disks from peg $A \rightarrow B$ are:

Move $n-1$ disks from $A \rightarrow C$
Move single remaining disk from $A \rightarrow B$
Move $n-1$ disks from $C \rightarrow B$

The same scheme can be used to move

$n-1$ disks.

The results are:

N=3
move disk from peg A to B
move disk from peg A to C
move disk from peg B to C
move disk from peg A to B
move disk from peg C to A
move disk from peg C to B
move disk from peg A to B
7 moves

The number of moves is $2^n-1$ . Indeed, if we run this with 4 disks, we see:

N=4
move disk from peg A to C
move disk from peg A to B
move disk from peg C to B
move disk from peg A to C
move disk from peg B to A
move disk from peg B to C
move disk from peg A to C
move disk from peg A to B
move disk from peg C to B
move disk from peg C to A
move disk from peg B to A
move disk from peg C to B
move disk from peg A to C
move disk from peg A to B
move disk from peg C to B
15 moves

This solution can be visualized as:

Inventory model

Just as in [2] where we modeled the wolf, goat, cabbage puzzle, we can model this in two ways: a formulation based on keeping track of inventory, and a shortest-path network formulation. Of course, no one in their right mind would solve this as an optimization problem. Obviously, that does not stop us from trying anyway. In an optimization model, we don't prescribe the moves to make, but rather formulate the conditions that should hold. A very different way to look at the problem.

The inventory is here modeled as the state: the location of the disks just after a move. I.e. we have $\color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}}} \in \{0,1\}$ where $t$ is the time index (one move per time period). We ignore here how the disks are organized (ordered from large to small), we just keep track of which peg a disk is assigned to. (In the next paragraphs, we will worry about the constraint that a move can only involve the smallest disks). We could formulate a constraint that says that every disk must be assigned to exactly one peg, a typical assignment constraint: $\sum_{\mathit{peg}} \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}}} = 1 \>\> \forall t,{\mathit{disk}}$ However, it turns out that this is not really necessary. The inventory balance equation takes care of that (given a valid initial inventory). Adding this constraint can help with performance.

The move variable is $\color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\in \{0,1\}$ We can only move one disk at a time, so: $\sum_{{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}=1\>\>\forall t$ Let's also forbid moves from and to the same peg (that is not a real move): $\color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}},{\mathit{peg}}} = 0$

Like in [2], we allow no moves when we are done: $\sum_{{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}= 1-\color{darkred}{\mathit{done}}_t$ where $\color{darkred}{\mathit{done}}_t \in \{0,1\}$ .

We need to make sure that the variable $\color{darkred}{\mathit{done}}_t$ is 1 if and only if all disks are placed on peg B: $\sum_{\mathit{disk}} \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{pegB}}} \ge n\cdot \color{darkred}{\mathit{done}}_t$ We rely on the objective to make sure that $\color{darkred}{\mathit{done}}_t=1$ when possible.

The inventory balance equation is a difference equation: $\color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}_1}} = \color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}_1}} - \sum_{\mathit{peg}_2} \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} + \sum_{\mathit{peg}_2} \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_2},{\mathit{peg}_1}}\>\>\forall t,\mathit{disk},\mathit{peg}_1$

The only thing left is worrying about the size of the disks. We can only move the smallest disk from a peg, and this disk should be smaller than the other disks on the target peg. In our model, we assume that we have $\color{darkblue}n$ disks with sizes $1,\dots,\color{darkblue}n$ . Just to make things more concrete. We could have used any other scheme. We need to say that: the size of the disk moved is (1) the smallest of the "from" peg, and (2) smaller than any disk on the target peg. For this, we look at the state (inventory) at the end of the previous step. Let's define a continuous variable $\color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} = \begin{cases} \min_{\mathit{disk}|\color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}}=1} \color{darkblue}{\mathit{size}}_{\mathit{disk}} & \text{if peg has one or more disks} \\ \color{darkblue}n & \text{if peg has zero disks}\end{cases}$ Instead of an equality it is easier to use a bound: $\color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} \le \begin{cases} \min_{\mathit{disk}|\color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}}=1} \color{darkblue}{\mathit{size}}_{\mathit{disk}} \ & \text{if peg has one or more disks} \\ \color{darkblue}n & \text{if peg has zero disks}\end{cases}$ This can be implemented as: $\color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} \le \color{darkblue}{\mathit{size}}_{\mathit{disk}} \cdot \color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}} + \color{darkblue}n \cdot (1-\color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}})$ With this, we can make sure we move an appropriate small disk: $\begin{align} & \sum_{\mathit{disk},\mathit{peg_2}}\color{darkblue}{\mathit{size}}_{\mathit{disk}}\cdot \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} \le \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_1}}} && \forall t,{\mathit{peg}}_1 \\ & \sum_{\mathit{disk},\mathit{peg_1}}\color{darkblue}{\mathit{size}}_{\mathit{disk}}\cdot \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} \le \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_2}}} && \forall t,{\mathit{peg}}_2 \end{align}$

The last inequality preferably should have a right-hand side that is one less. However, we cannot just subtract one from $\color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_2}}}$ as $\color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_2}}}-1$ may forbid moving disks of size $\color{darkblue}n$ . It is noted that the constraint is correct, as we can never place two disks of the same size on a peg. So, I am a bit lazy here, and just reuse the same $\color{darkred}{\mathit{smallestPrev}}$ for both the source and destination peg.

The objective is simply to maximize the number of "done" iterations: $\max \sum_t \color{darkred}{\mathit{done}}_t$ This is only needed if we overestimate the number iterations we need (the set $t$ is larger than needed). If we use a size of $2^n-1$ for the set $t$ , we will be done exactly at the last iteration.

The complete model looks like:

Inventory Model
$\begin{align}\max&\sum_t \color{darkred}{\mathit{done}}_t \\ & \sum_{{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}=1-\color{darkred}{\mathit{done}}_t&& \forall t\\& \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}},{\mathit{peg}}} = 0 && \forall t,{\mathit{disk}},{\mathit{peg}} \\ & \sum_{\mathit{disk}} \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{pegB}}} \ge n\cdot \color{darkred}{\mathit{done}}_t && \forall t\\ & \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}_1}} = \color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}_1}} - \sum_{\mathit{peg}_2} \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} + \sum_{\mathit{peg}_2} \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_2},{\mathit{peg}_1}} && \forall t,\mathit{disk},\mathit{peg}_1 \\ & \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} \le \color{darkblue}{\mathit{size}}_{\mathit{disk}} \cdot \color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}} + \color{darkblue}n \cdot (1-\color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}}) && \forall t, {\mathit{disk}},{\mathit{peg}} \\ & \sum_{\mathit{disk},\mathit{peg_2}}\color{darkblue}{\mathit{size}}_{\mathit{disk}}\cdot \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} \le \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_1}}} && \forall t,{\mathit{peg}}_1 \\ & \sum_{\mathit{disk},\mathit{peg_1}}\color{darkblue}{\mathit{size}}_{\mathit{disk}}\cdot \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} \le \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_2}}} && \forall t,{\mathit{peg}}_2 \\ & \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}}} \in \{0,1\} \\ & \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\in \{0,1\} \\ & \color{darkred}{\mathit{done}}_t \in \{0,1\}\\ & \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} \in [1,\dots,n] \end{align}$

Inventory Model

$\begin{align}\max&\sum_t \color{darkred}{\mathit{done}}_t \\ & \sum_{{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}=1-\color{darkred}{\mathit{done}}_t&& \forall t\\& \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}},{\mathit{peg}}} = 0 && \forall t,{\mathit{disk}},{\mathit{peg}} \\ & \sum_{\mathit{disk}} \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{pegB}}} \ge n\cdot \color{darkred}{\mathit{done}}_t && \forall t\\ & \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}_1}} = \color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}_1}} - \sum_{\mathit{peg}_2} \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} + \sum_{\mathit{peg}_2} \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_2},{\mathit{peg}_1}} && \forall t,\mathit{disk},\mathit{peg}_1 \\ & \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} \le \color{darkblue}{\mathit{size}}_{\mathit{disk}} \cdot \color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}} + \color{darkblue}n \cdot (1-\color{darkred}{\mathit{inv}}_{t-1,{\mathit{disk}},{\mathit{peg}}}) && \forall t, {\mathit{disk}},{\mathit{peg}} \\ & \sum_{\mathit{disk},\mathit{peg_2}}\color{darkblue}{\mathit{size}}_{\mathit{disk}}\cdot \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} \le \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_1}}} && \forall t,{\mathit{peg}}_1 \\ & \sum_{\mathit{disk},\mathit{peg_1}}\color{darkblue}{\mathit{size}}_{\mathit{disk}}\cdot \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}} \le \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg_2}}} && \forall t,{\mathit{peg}}_2 \\ & \color{darkred}{\mathit{inv}}_{t,{\mathit{disk}},{\mathit{peg}}} \in \{0,1\} \\ & \color{darkred}{\mathit{move}}_{t,{\mathit{disk}},{\mathit{peg}_1},{\mathit{peg}_2}}\in \{0,1\} \\ & \color{darkred}{\mathit{done}}_t \in \{0,1\}\\ & \color{darkred}{\mathit{smallestPrev}}_{t,{\mathit{peg}}} \in [1,\dots,n] \end{align}$

The results for a

$n=3$ disk problem are:

----    157 VARIABLE move.L  disk is moved from one peg to another

           pegA.pegB   pegA.pegC   pegB.pegC   pegC.pegA   pegC.pegB

t1.disk1           1
t2.disk2                       1
t3.disk1                                   1
t4.disk3           1
t5.disk1                                               1
t6.disk2                                                           1
t7.disk1           1


----    157 VARIABLE inv.L  inventory after move

                pegA        pegB        pegC

t1.disk1                   1.000
t1.disk2       1.000
t1.disk3       1.000
t2.disk1                   1.000
t2.disk2                               1.000
t2.disk3       1.000
t3.disk1                               1.000
t3.disk2                               1.000
t3.disk3       1.000
t4.disk1                               1.000
t4.disk2                               1.000
t4.disk3                   1.000
t5.disk1       1.000
t5.disk2                               1.000
t5.disk3                   1.000
t6.disk1       1.000
t6.disk2                   1.000
t6.disk3                   1.000
t7.disk1                   1.000
t7.disk2                   1.000
t7.disk3                   1.000

The corresponding plot generated from this solution is:

GAMS Model

$onText

   Towers of Hanoi
   
$offText


*-----------------------------------------------------------------------------------
* size of problem
*-----------------------------------------------------------------------------------

$set n 3
$set makeplot 1

*-----------------------------------------------------------------------------------
* recursion in Python
*-----------------------------------------------------------------------------------


$onEmbeddedCode Python:

# returns number of moves
def hanoi(n,A,B,C):
   if n==0: return 0
   n1 = hanoi(n-1,A,C,B)
   print(f"move disk from peg {A} to {B}")
   n2 = hanoi(n-1,C,B,A)
   return n1+1+n2
   
N = %n%
print(f"N={N}")
cnt = hanoi(N,'A','B','C')
print(f"{cnt} moves")

$offEmbeddedCode


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

* max number of moves
* we can guess (overestimate) or use the known value 2^n-1
*$set tmax 25
$eval tmax 2**%n%-1

sets
   dummy 'for ordering' /initial/
   peg /pegA,pegB,pegC/
   disk  /disk1*disk%n%/
   t 'moves or timesteps' /t1*t%tmax%/
;
display peg,disk,t;

alias (peg,peg1,peg2);

Parameters
    n 'number of disks' /%n%/
    size(disk) 'size of each disk: 1..n'
    maxsize 'maximum size of disk'
    initial(t,disk,peg) 'initial inventory (state)'
    final(disk,peg) 'final inventory (state)'
;

abort$(card(t)<2**n-1) "increase size of set t";

size(disk) = ord(disk);
maxsize = smax(disk,size(disk));

* initial state: all disks on peg A
* has t index to make it easier to use in the constraints
initial('t1',disk,'pegA') = 1;

* final state: all disks on peg B
final(disk,'pegB') = 1;

display n,size,maxsize,initial,final;


*-----------------------------------------------------------------------------------
* MIP model
*-----------------------------------------------------------------------------------


binary variable
   move(t,disk,peg,peg)   'disk is moved from one peg to another'
   inv(t,disk,peg)        'inventory after move'
   done(t)                'we are done'   
;
variable smallestPrev(t,peg) 'smallest disk in previous iteration';
smallestPrev.lo(t,peg) = 1;
smallestPrev.up(t,peg) = maxsize;

variable z 'obj';

equations  
   move1disk(t)           'move 1 disk (or 0 if done)'
   bndsmall(t,disk,peg)   'bound on smallest size of disk'
   smallfrom(t,peg)       'need to select smallest disk from source peg' 
   smallto(t,peg)         'size restriction on destination peg' 
   inventory(t,disk,peg)  'inventory balance' 
   isdone(t)              'we are done'
   ordering(t)            'optional'
   assignment(t,disk)     'optional assignment constraint'
   obj                    'objective'  
;


* if done=1: move zero disks
* if done=0: move exactly one disk
move1disk(t).. sum((disk,peg1,peg2),move(t,disk,peg1,peg2)) =e= 1-done(t-1);

* don't move a disk from and to the same peg
move.fx(t,disk,peg,peg) = 0;

* we can only move the smallest disk
parameter initsmallest(t,peg);
initsmallest('t1',peg) = maxsize;
initsmallest('t1','pegA') = smin(disk,size(disk));
display initsmallest;

smallestPrev.fx('t1',peg) = initsmallest('t1',peg);
bndsmall(t,disk,peg)$(ord(t)>1).. smallestPrev(t,peg) =l= size(disk)*inv(t-1,disk,peg) + n*(1-inv(t-1,disk,peg));
smallfrom(t,peg)..     sum((disk,peg2),size(disk)*move(t,disk,peg,peg2)) =l= smallestPrev(t,peg);
smallto(t,peg)..       sum((disk,peg1),size(disk)*move(t,disk,peg1,peg)) =l= smallestPrev(t,peg);

* inventory balance
inventory(t,disk,peg)..
  inv(t,disk,peg) =e= inv(t-1,disk,peg) - sum(peg2,move(t,disk,peg,peg2)) + sum(peg1,move(t,disk,peg1,peg)) + initial(t,disk,peg);
  
* done if all disks are at peg B
isdone(t).. sum(disk,inv(t,disk,'pegB')) =g= n*done(t);

obj.. z =e= sum(t,done(t));

* optional constraints
ordering(t-1).. done(t) =g= done(t-1);
assignment(t,disk).. sum(peg,inv(t,disk,peg)) =e= 1;

* optional fixes
set lastt(t);
lastt(t) = ord(t)=card(t);
done.fx(lastt) = 1;
inv.fx(lastt,disk,peg) = final(disk,peg); 

model m /all/;

*-----------------------------------------------------------------------------------
* Solve
*-----------------------------------------------------------------------------------

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

* raw results
option move:0:2:2;
display smallestPrev.l,move.l,inv.l,done.l;

*-----------------------------------------------------------------------------------
* Reporting
*-----------------------------------------------------------------------------------

set sinv(*,disk,peg) 'include initial, and stop after done';
sinv('initial',disk,peg) = initial('t1',disk,peg);
loop(t,
   sinv(t,disk,peg) = inv.l(t,disk,peg)>0.5;
   break$(done.l(t)>0.5);
);
display sinv;


*-----------------------------------------------------------------------------------
* Visualization
*-----------------------------------------------------------------------------------

abort.noError$(n>6 or %makeplot%=0) "skipping plot";

$set svg hanoiplots.html

file f /%svg%/;
put f;

put '<style>table,th,td {border-collapse: collapse;}</style>'/;
put '<h2>Towers of Hanoi Results</h2>'/;
put 'Number of disks: %n%'/;

put '<table border="1">'/;
put '<tr>'/;

$eval n2 %n%+2

alias (t1,*);
parameter ndisks(*,peg) 'number of disks on each peg';
ndisks(t1,peg) = sum(sinv(t1,disk,peg),1);
display ndisks;

scalar nd,dw,x,y,w,k /0/;


loop(t1$sum(sinv(t1,disk,peg),1),
   if (k=4,
       put "</tr><tr>"/;
       k = 0;
    );
   k = k + 1;

   put '<td style="text-align: center">'/;
   put '<svg height="100" width="300" viewBox="0 0 12 %n2%">'/;
* draw pegs   
   put '<line x1="3" y1="1" x2="3" y2="%n2%" style="stroke:brown;stroke-width:0.1"/>'/;
   put '<line x1="6" y1="1" x2="6" y2="%n2%" style="stroke:brown;stroke-width:0.1"/>'/;
   put '<line x1="9" y1="1" x2="9" y2="%n2%" style="stroke:brown;stroke-width:0.1"/>'/;

   loop(peg,
      nd = ndisks(t1,peg);
      loop(sinv(t1,disk,peg),
         dw = size(disk);
         y = n+1-nd+1;
         w = size(disk)*3/n;
         x = 3*ord(peg)-0.5*w;
*         display x,w,y;
         put '<rect x="',x:4:2,'" y="',y:4:2,'" height="1" width="',w:4:2,'" fill="lightblue"/>'/;
         put '<text x="',(3*ord(peg)):0:0,'" y="',(y+0.6):3:1,'" dominant-baseline="middle" text-anchor="middle" font-size="0.6">',(ord(disk)):0:0,'</text>'/;
         nd = nd-1;
      );
   );
   put '</svg><br>',t1.tl:0/;
   put '</td>'/;
);

put '</tr></table>';
putclose;
executetool 'win32.ShellExecute "%svg%"';

Shortest path network model

An alternative is to use a shortest path network model. The nodes are the different states. Each disk can be placed on one of three pegs, so the number of nodes is equal to

$3^n$ where

$n$ is the number of disks. For

$n=3$ we have 27 nodes:

----     30 SET disk  

disk1,    disk2,    disk3


----     30 SET peg  

pegA,    pegB,    PegC


----     30 SET node  

node1 ,    node2 ,    node3 ,    node4 ,    node5 ,    node6 ,    node7 ,    node8 ,    node9 ,    node10,    node11,    node12
node13,    node14,    node15,    node16,    node17,    node18,    node19,    node20,    node21,    node22,    node23,    node24
node25,    node26,    node27


----     30 SET nodes  

                    pegA        pegB        PegC

node1 .disk1         YES
node1 .disk2         YES
node1 .disk3         YES
node2 .disk1         YES
node2 .disk2         YES
node2 .disk3                     YES
node3 .disk1         YES
node3 .disk2         YES
node3 .disk3                                 YES
node4 .disk1         YES
node4 .disk2                     YES
node4 .disk3         YES
node5 .disk1         YES
node5 .disk2                     YES
node5 .disk3                     YES
node6 .disk1         YES
node6 .disk2                     YES
node6 .disk3                                 YES
node7 .disk1         YES
node7 .disk2                                 YES
node7 .disk3         YES
node8 .disk1         YES
node8 .disk2                                 YES
node8 .disk3                     YES
node9 .disk1         YES
node9 .disk2                                 YES
node9 .disk3                                 YES
node10.disk1                     YES
node10.disk2         YES
node10.disk3         YES
node11.disk1                     YES
node11.disk2         YES
node11.disk3                     YES
node12.disk1                     YES
node12.disk2         YES
node12.disk3                                 YES
node13.disk1                     YES
node13.disk2                     YES
node13.disk3         YES
node14.disk1                     YES
node14.disk2                     YES
node14.disk3                     YES
node15.disk1                     YES
node15.disk2                     YES
node15.disk3                                 YES
node16.disk1                     YES
node16.disk2                                 YES
node16.disk3         YES
node17.disk1                     YES
node17.disk2                                 YES
node17.disk3                     YES
node18.disk1                     YES
node18.disk2                                 YES
node18.disk3                                 YES
node19.disk1                                 YES
node19.disk2         YES
node19.disk3         YES
node20.disk1                                 YES
node20.disk2         YES
node20.disk3                     YES
node21.disk1                                 YES
node21.disk2         YES
node21.disk3                                 YES
node22.disk1                                 YES
node22.disk2                     YES
node22.disk3         YES
node23.disk1                                 YES
node23.disk2                     YES
node23.disk3                     YES
node24.disk1                                 YES
node24.disk2                     YES
node24.disk3                                 YES
node25.disk1                                 YES
node25.disk2                                 YES
node25.disk3         YES
node26.disk1                                 YES
node26.disk2                                 YES
node26.disk3                     YES
node27.disk1                                 YES
node27.disk2                                 YES
node27.disk3                                 YES

The directed arcs $\mathit{node}_1 \rightarrow \mathit{node}_2$ have the following properties:

The difference between $\mathit{node}_1$ and $\mathit{node}_2$ is just one disk being moved from one peg to another. E.g., there is no arc between node 1 and node 5, as the difference between them is 2.
The disk being moved is the smallest on its peg, both at $\mathit{node}_1$ and $\mathit{node}_2$ . E.g., these is no arc between node 2 and node 5, as disk 2 is not the smallest on peg A in node 2.

If we check this systematically and carefully, we get the following set of feasible moves (and thus arcs):

----     63 SET diskMoved  feasible moves

                    disk1       disk2       disk3

node1 .node10         YES
node1 .node19         YES
node2 .node3                                  YES
node2 .node11         YES
node2 .node20         YES
node3 .node2                                  YES
node3 .node12         YES
node3 .node21         YES
node4 .node7                      YES
node4 .node13         YES
node4 .node22         YES
node5 .node8                      YES
node5 .node14         YES
node5 .node23         YES
node6 .node9                      YES
node6 .node15         YES
node6 .node24         YES
node7 .node4                      YES
node7 .node16         YES
node7 .node25         YES
node8 .node5                      YES
node8 .node17         YES
node8 .node26         YES
node9 .node6                      YES
node9 .node18         YES
node9 .node27         YES
node10.node1          YES
node10.node16                     YES
node10.node19         YES
node11.node2          YES
node11.node17                     YES
node11.node20         YES
node12.node3          YES
node12.node18                     YES
node12.node21         YES
node13.node4          YES
node13.node15                                 YES
node13.node22         YES
node14.node5          YES
node14.node23         YES
node15.node6          YES
node15.node13                                 YES
node15.node24         YES
node16.node7          YES
node16.node10                     YES
node16.node25         YES
node17.node8          YES
node17.node11                     YES
node17.node26         YES
node18.node9          YES
node18.node12                     YES
node18.node27         YES
node19.node1          YES
node19.node10         YES
node19.node22                     YES
node20.node2          YES
node20.node11         YES
node20.node23                     YES
node21.node3          YES
node21.node12         YES
node21.node24                     YES
node22.node4          YES
node22.node13         YES
node22.node19                     YES
node23.node5          YES
node23.node14         YES
node23.node20                     YES
node24.node6          YES
node24.node15         YES
node24.node21                     YES
node25.node7          YES
node25.node16         YES
node25.node26                                 YES
node26.node8          YES
node26.node17         YES
node26.node25                                 YES
node27.node9          YES
node27.node18         YES

The number of directed arcs is known [3]: $3\cdot(3^n-1)$ , so we have an easy check if our procedure to generate arcs is correct. For $n=3$ the number of arcs is 78.

To solve this model, we use a standard LP formulation:

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:

----     93 VARIABLE f.L  flow

             node5       node8      node10      node14      node16      node25      node26

node1                                1.000
node5                                            1.000
node8        1.000
node10                                                       1.000
node16                                                                   1.000
node25                                                                               1.000
node26                   1.000

A more meaningful report created from this solution is:

----    125 SET trace  

move1.node1 .node10.disk1.pegA.pegB
move2.node10.node16.disk2.pegA.PegC
move3.node16.node25.disk1.pegB.PegC
move4.node25.node26.disk3.pegA.pegB
move5.node26.node8 .disk1.PegC.pegA
move6.node8 .node5 .disk2.PegC.pegB
move7.node5 .node14.disk1.pegA.pegB

This is the same solution as from the previous model.

GAMS Model

$onText

  Towers of Hanoi, network formulation

$offText


*-----------------------------------------------------------------------------------
* size of problem
*-----------------------------------------------------------------------------------

$set n 3

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

$eval nn 3**%n%

set
  node /node1*node%nn%/
  disk /disk1*disk%n%/
  peg /pegA,pegB,PegC/
  nodes(node,disk,peg) / system.powersetRight /
;

scalars
   n  'number of disks' /%n%/
   nn 'number of nodes' /%nn%/
;

display disk,peg,node,nodes;

*-----------------------------------------------------------------------------------
* arcs
*-----------------------------------------------------------------------------------

set arc(node,node);

alias (node,node1,node2), (disk,disk1,disk2), (peg,peg1,peg2), (nodes,nodes1,nodes2);

set smallest(node,disk,peg) 'smallest disk on peg';
smallest(node,disk,peg) = smin(disk2$nodes(node,disk2,peg),ord(disk2)) = ord(disk);
display smallest;

parameter nodeVal(node,disk) 'peg number of (node,disk)';
nodeVal(node,disk) = sum(nodes(node,disk,peg),ord(peg));
display nodeVal;

set nodediff1(node,node) '(node1,node2) combos with one difference';
nodediff1(node1,node2) = sum(disk$(abs(nodeVal(node1,disk)-nodeVal(node2,disk))>0.5),1)=1;
display nodediff1;

set nodeDiff1disk(node,node,disk) 'augment nodediff1 with disk id';
nodeDiff1Disk(nodeDiff1(node1,node2),disk) = abs(nodeVal(node1,disk)-nodeVal(node2,disk))>0.5;
display nodeDiff1Disk;

set diskMoved(node,node,disk) 'feasible moves';
singleton set p1(peg),p2(peg);
loop(nodeDiff1disk(node1,node2,disk),
   p1(peg) = nodeval(node1,disk)=ord(peg);
   p2(peg) = nodeval(node2,disk)=ord(peg);
   diskMoved(node1,node2,disk) = smallest(node1,disk,p1) and smallest(node2,disk,p2);  
);
display diskMoved;

arc(node1,node2) = sum(diskMoved(node1,node2,disk),1);

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

Parameters
   supply(node)
   demand(node)  
;
supply(node)$(sum(nodes(node,disk,'pegA'),1) = n) = 1;
demand(node)$(sum(nodes(node,disk,'pegB'),1) = n) = 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(node1)..
   sum(arc(node2,node1),f(node2,node1)) + supply(node1) =e= sum(arc(node1,node2),f(node1,node2)) + demand(node1);

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

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

set
    move/move1*move100/
    trace(move,node,node,disk,peg,peg)    
;

singleton sets
  cur(node) 'current node' /node1/
  next(node) 'next node'
;

cur(node) = supply(node)>0;

loop(move$card(cur),

    next(node) = f.l(cur,node)>0.5;

    loop(diskMoved(cur,next,disk),
        p1(peg) = nodeval(cur,disk)=ord(peg);
        p2(peg) = nodeval(next,disk)=ord(peg);
        trace(move,cur,next,disk,p1,p2) = yes;
    );

    cur(node) = next(node);  
);

option trace:0:0:1;
display trace;

Conclusion

Interestingly, we could use the same strategies we previously applied for the Wolf, Goat, and Cabbage problem [2]: a MIP model where the dynamics are handled by an inventory balance equation and a shortest path model.

From a performance point of view, the network approach is much better: shortest path problems are easy LPs, and we can even use specialized algorithms if we want. The MIP model is not so easy when we increase the size of the problem a bit. Interestingly, the open source solver HiGHS seems to do better than Cplex on this particular model (caveat: I only ran a few tests). It is also noted that the implementations of both models (MIP and network model) are not carefully optimized for speed. I just wanted to demonstrate the feasibility of these approaches for this problem. The shortest path approach has been suggested by others, so that is not my idea.

I believe these models are a bit more general than the recursive algorithm: we can start from any given configuration, while the recursion can only start with all disks on one peg. (Similar for the final configuration: we are flexible about this.)

References

Tower of Hanoi, https://en.wikipedia.org/wiki/Tower_of_Hanoi
Wolf, goat, and cabbage problem: MIP and Network Model, https://yetanothermathprogrammingconsultant.blogspot.com/2025/03/wolf-goat-and-cabbage-problem-mip-and.html
Hanoi graph, https://en.wikipedia.org/wiki/Hanoi_graph

1 comment:

Jacob UrickApril 5, 2025 at 10:57 AM
HiGHS is quite formidable. While I would be shocked if we see it outperform commercial solvers on standard benchmark sets like MIPLIB20XX, it is absolutely worth running a comparison to see how it performs on a specific problem vs a commercial solver - you may find that the free solution is faster and better.

Tuesday, April 1, 2025

Towers of Hanoi: inventory and network formulation