Finding the longest simple path on the divisor graph using linear programming

Question

I've been working on this problem posed by Jane Street. The same problem has been asked by NRICH and 538 Riddler

"Write down a chain of integers between 1 and 100, with no repetition, such that if x and y are consecutive numbers in the chain, then x evenly divides y or y evenly divides x. Here is an example of such a chain, with length 12: 37, 74, 2, 8, 4, 16, 48, 6, 3, 9, 27, 81 What is the longest chain you can find?"

After some searching I found he list of longest paths for a given integer on OEIS with the following comment

Can be formulated as an optimal subtour problem by introducing a depot node 0 that is adjacent to all other nodes.
An integer linear programming formulation is as follows.
For {i,j} in E, let binary decision variable x_{i,j} indicate whether edge {i,j} is traversed, and
for i in N let binary decision variable y_i indicate whether node i is visited.
The objective is to maximize Sum_{i in N \ {0}} y_i.
The constraints are Sum_{{i,j} in E: k in {i,j}} x_{i,j} = 2 y_k for all k in N, y_0 = 1,
as well as (dynamically generated) generalized subtour elimination constraints
Sum_{i in S, j in S: {i,j} in E} x_{i,j} <= Sum_{i in S \ {k}} y_i for all S subset N \ {0} and k in S.
- Rob Pratt, Dec 28 2020

The page lists the first 200 terms.

I've attempted to code this up in python using pulp but the code is very slow and not reliable. Looking for some help to correct the code below and some insight on how someone managed expand this to calculate 200 terms.

import pulp,networkx as nx,time,matplotlib.pyplot as plt
from itertools import chain,combinations

# implementing the function above but starts to get very slow as the size of the powerset goes up at 2^n
# the subtour elimination isn't working well. Think should be including self in the powerset and the 
# comparison should be < not <=

def powerset(iterable):
    "powerset excluding null and self"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(1,len(s)))

########################################################################

start = time.time()

length = 16

model=pulp.LpProblem('tsp',pulp.LpMaximize)

###############
#  Variables  #
###############
N = pulp.LpVariable.dicts("N",[i for i in range(length+1)],upBound=1,lowBound=0,cat='Binary')
E = pulp.LpVariable.dicts("E",[(i,j) for i in range(1,length) for j in range(i+1,length+1) if j % i==0],upBound=1,lowBound=0,cat='Binary')
for i in range(1,length+1):
    E[(0,i)] = pulp.LpVariable("E_(%s,_%s)" % (0,i),upBound=1,lowBound=0,cat='Binary')
    
model += N[0] == 1   

###################
#    Objective    #
###################
model += pulp.lpSum([N[i] for i in N ])    

###############
# Constraints #
###############

# if visited has 2 edges otherwise zero
for i in range(length+1) :
    model += pulp.lpSum([E[i,j] if (i,j) in E else E[(j,i)] if (j,i) in E else 0 for j in range(length+1)]) == 2*N[i]   

# check every subset (excluding zero) has fewer edges than nodes
for P in powerset(range(1,length+1)):
    model += pulp.lpSum([E[(i,j)] for i in P for j in P if (i,j) in E]) <= pulp.lpSum([N[i] for i in P])
  
print("Setup took {:,.4f} seconds\n ".format(time.time()-start))   

###################
# Solve and print #
###################
    
model.solve()

e = {k: v.varValue  for k,v in E.items()}
n = {k: v.varValue  for k,v in N.items()}

print("Took {:,.4f} seconds\n ".format(time.time()-start))   

g = nx.Graph()
print("Path Length",len([g.add_edge(i,j) for i,j in e if e[(i,j)]==1 and i*j !=0])+1)
nodes = sorted(g.degree, key=lambda x: x[1])
print([*nx.all_simple_paths(g,nodes[0][0],nodes[1][0])])
print("\nMissing",[k for k in n if n[k] ==0])

fig, ax = plt.subplots(figsize=(10,5))
fig.suptitle("Graph for n ={}".format(length),fontsize=20)
pos = nx.spring_layout(g, seed=10)  
options = { "font_size": 8,"node_size": 200,"node_color": "white","edgecolors": "black","linewidths": 1,"width": 1,}
nx.draw_networkx(g, pos, **options)
plt.show()