Solve Hashiwokakero puzzle using Pyomo

This saturday I decided to solve an interesting puzzle with the follwoing simple rules.

Rules:

1. Each circle with number must have as many lines from it as the number in the circle.
2. It is not allowed to have more than two lines between two circles.
3. The lines can not cross eachother.
4. All the lines must either go horizontal or vertical. It is not allowed to draw lines diagonally.
5. All lines must be connected.

Some solved examples:

But is it possible to formulate it as a MILP or constraint programming ?

The answer is yes.

In this episod, I do formulate it as a MILP model and solve it using Pyomo and CBC solver.

MILP Model :

Let's explain the model.

It looks trivial and easy, however, it is anything but easy.

First we need to define the decision variables:

• Xij shows the number of links from node i to node j (in order to reduce the computation burden , we just let i>j take non-zero values), it takes value from {0,1,2} (integer variable)
• Uij shows if the node i to j is connected or not. (binary variable)
• flowij is a real variable in [0,1] interval for ensuring the connectivity of nodes
• source_i is the source node (it is only available in node 1)
• Allowed Set A: is the set of allowed links . it contains (i,j) if node i can be connected to link j (the line between them is vertical or horizontal)
• Overlap Set O: this set is a very important one in this formulation. If link (i,j) and link (m,n) can not simultaneously exst then these two links will be stored in O.

I think this specific set requires more explanation

Here are some examples of links that they do overlap (and should not exist in the final solution simultaneously)

• Case: A both links can take values
• Case B: Not allowed
• Case C: Not allowed
• Case D: Allowed

Now let's explain the model :

Objective function is defined as the total number of used links

1. The value of node n (Vn) is equal to the outgoing/incoming links to this node
2. If link (i,j) and link (n,m) are in non-overlap set then they can't be both on the final solution
3. if Xij takes a non-zero value then Uij =1
4. To ensure the connectivity of the final graph, a single source at node 1 is considered (source_n is only non-zero at node 1) and a nodal balance constraint is considered at each node.
5. floweij can happen only if the link between i and j is on

Pyomo code

model = AbstractModel()
model.i = RangeSet(Nc)
model.j = Set(initialize=model.i)

model.X = Var(model.i,model.j,bounds=(0,2),initialize=0, within=Integers)
model.U = Var(model.i,model.j,bounds=(0,1),initialize=0, within=Binary)
model.flow = Var(model.i,model.j,bounds=(0,1),initialize=0, within=Reals)
model.s = Var(bounds=(0,1),initialize=0, within=Reals)

def rule_C1(model,c): 
    return   data[c,'v'] == sum(model.X[c,j] for j in model.j if (c,j) in lines) +\
                            sum(model.X[j,c] for j in model.j if (j,c) in lines)
model.C1 = Constraint(model.i,rule=rule_C1)

def rule_C2(model,m,n,i,j): 
    L1 = [(data[i,'x'],data[i,'y']), (data[j,'x'],data[j,'y'])]
    L2 = [(data[m,'x'],data[m,'y']), (data[n,'x'],data[n,'y'])]
    line = LineString([L1[0], L1[1]])
    other = LineString([L2[0], L2[1]])
    a = line.intersection(other)
    a = line.intersection(other)
    cancel_con = len(set([i,j,m,n])) <4
    con0 = (m,n) in lines and (i,j) in lines and (m,n)!=(i,j) and (m*n>i*j)
    con1 = con0 and (not cancel_con) and line.intersects(other)
    con2 = con0 and (cancel_con and a.length> 1) and line.intersects(other) and (m,n)!=(i,j) and (m*n>i*j)

    if con1 or con2:
        return   model.U[i,j]+model.U[m,n] <= 1
    else:
        return Constraint.Skip
model.C2 = Constraint(model.i,model.i,model.i,model.i, rule=rule_C2)

def rule_C3(model,i,j):
    if (i,j) in lines:
        return   model.X[i,j]<= 2*model.U[i,j]
    else:
        return Constraint.Skip
model.C3 = Constraint(model.i,model.j, rule=rule_C3)

def rule_C3B(model,i,j):
    if (i,j) in lines:
        return   model.X[i,j]>= model.U[i,j]
    else:
        return Constraint.Skip
model.C3B = Constraint(model.i,model.j, rule=rule_C3B)

def rule_C4A(model,i):
    if i ==1 :
        return model.s - 0.001 == sum( model.flow[i,j]-model.flow[j,i] for j in model.j if (i,j) in ALL_lines ) 
    else:
        return         - 0.001 == sum( model.flow[i,j]-model.flow[j,i] for j in model.j if (i,j) in ALL_lines ) 
model.C4A = Constraint(model.i, rule=rule_C4A)

def rule_C4B(model,i,j):
    if (i,j) in lines:
        return   model.flow[i,j]<= model.U[i,j]
    elif (j,i) in lines:
        return   model.flow[i,j]<= model.U[j,i]
    else:
        return Constraint.Skip
model.C4B = Constraint(model.i,model.j, rule=rule_C4B)

def rule_OF(model):
    return  sum(model.U[i,j] for (i,j) in lines) 
model.obj1 = Objective(rule=rule_OF, sense=minimize)
instance = model.create_instance()          

Results:

The thick/thin lines represent two parallel /one single lines.

Some Real world applications :

• Project/task scheduling with special competency requirements
• Strategic allocation of assets

Sometimes I have been asked "what is the benefit of solving these toy puzzles" ?

Solving puzzle might seem to be waste of time but by solving these puzzles you obtain the following skills:

• Effitiently transforming the oral problem description into the math model (critical thinking)
• Data preparation
• Model debugging
• Visualizing the output

The Pyomo code is available on Github

Github: https://lnkd.in/eJWCJcPC