PuLP integer programming: solver returns 'Undefined' for anything beyond toy example

Question

I am trying to find an integer programming formulation for the following problem:

Sort a list of items according to 2 criteria, represented by the following cost:

• 'Positional cost': dependent on the position of the item in the list
• 'Neighbour cost': dependent on two items being positioned next to eachother in the list (ie the items are immediate neighbours)

The following apply:

• Every item can only be placed once
• The list has a beginning position P_1 (-> the item in that position does not have a preceeding list item)
• The list has an end position P_X (-> the item in that position does not have a subsequent item)

I used PuLP to formulate as described below and tried solving with GLPK_CMD and COIN_CMD. It works for 3 items, but returns 'Undefined' for 4 or more items. As far as I can tell, an order such as the one described in DataFrame 'feasible_solution_example' will not violate the constraints. can anyone venture a guess why the solver does not find a solution?
Any input and comments are highly appreciated.

import pulp 
import numpy as np
import pandas as pd 


############# Size of problem: number of items in the list to be sorted
length_of_list=12
list_of_items=['I'+str(x).zfill(3)  for x in np.arange(1, length_of_list+1)]
list_of_neighbours=['N_I'+str(x).zfill(3)  for x in np.arange(1,length_of_list+1)]
list_of_positions=['P'+str(x).zfill(3) for x in np.arange(1,length_of_list+1)]


############# Cost matrices
random_position_cost=np.random.randint(-length_of_list, length_of_list, [length_of_list,length_of_list])
positional_cost=pd.DataFrame(random_position_cost, index=list_of_items, columns=list_of_positions).astype(int)
random_neighbour_cost=np.random.randint(-length_of_list, length_of_list, [length_of_list,length_of_list])
neigbour_cost=pd.DataFrame(data=random_neighbour_cost, index=list_of_items, columns=list_of_neighbours)

############# This is/should be an example of a feasible solution 
feasible_solution_example=pd.DataFrame(index=list_of_items, columns=(list_of_items + list_of_neighbours))
feasible_solution_example.loc[:, list_of_items]=np.identity(length_of_list)
feasible_solution_example.loc[:, list_of_neighbours[0]]=0
feasible_solution_example.loc[:, list_of_neighbours[1:]]=np.identity(length_of_list)[:,:-1]
feasible_solution_example=feasible_solution_example.astype(int)

############# Model definition 
model = pulp.LpProblem("Optimal order problem ", pulp.LpMinimize)

# Decision variables: 
item_in_position = pulp.LpVariable.dicts("item_in_position ",
                                     ((a, b) for a in feasible_solution_example.index for b in list_of_positions),
                                     cat='Binary')

item_neighbours = pulp.LpVariable.dicts("item_neighbours",
                                     ((a, b) for a in feasible_solution_example.index for b in list_of_neighbours if str('N_' + a)!=b),
                                     cat='Binary')

# Objective Function
model += (
    pulp.lpSum([
        (positional_cost.loc[i, j] * item_in_position[(i, j)])
        for i in feasible_solution_example.index for j in list_of_positions] + 
        [(neigbour_cost.loc[i, j] * item_neighbours[(i, j)])
        for i in feasible_solution_example.index for j in list_of_neighbours if str('N_' + i)!=j]
        )
)    


## Constraints: 
# 1- Every item can take only one position...    
for cur_item in feasible_solution_example.index: 
    model += pulp.lpSum([item_in_position[cur_item, j] for j in list_of_positions]) == 1

# 2-...and every position can only be taken once
for cur_pos in list_of_positions: 
    model += pulp.lpSum([item_in_position[i, cur_pos] for i in feasible_solution_example.index]) == 1

# 3-...item in position 1 can not be the neighbour (ie the preceeding item) to any other item : 
#     -> all neighbour values for any item x + the value for position 1 for this item must add up to 1 
for cur_neighbour in list_of_neighbours: 
    model += pulp.lpSum([item_neighbours[cur_item, cur_neighbour] for cur_item in list_of_items if cur_neighbour!=str('N_' + cur_item)] + 
                         item_in_position[cur_neighbour.split('_')[1], list_of_positions[0]] ) == 1
# 4-...item in the last position can not have any neighbours (ie any preceeding items): 
#     -> all neighbour values of all items + value for the last position must sum up to 1: 
for cur_item in feasible_solution_example.index: 
    model += pulp.lpSum(([item_neighbours[cur_item, cur_neighbour] for cur_neighbour in list_of_neighbours if cur_neighbour!=str('N_' + cur_item)] + 
                         [item_in_position[cur_item,list_of_positions[-1]]])) ==1


# 5-... When two items are neighbours the cost for that combination needs to be "switched on"  
#     -> (e.g. I001 in position P001 and I002 in P002 then I001 and N_I002 must be eqaul to 1) 
for item_x in np.arange(0,length_of_list-1): 
    for neighbour_y in np.arange(0,length_of_list-1): 
        if item_x!=neighbour_y: 
            # Item x would be neighbour with lot y IF:
            for cur_position in np.arange(0,length_of_list-1):
                #...they would have subsequent positions...    
                model += item_neighbours[list_of_items[item_x], list_of_neighbours[neighbour_y]] +1 >= item_in_position[list_of_items[item_x],list_of_positions[cur_position]]  + item_in_position[list_of_items[neighbour_y],list_of_positions[cur_position+1]]


############# Solve 
model.solve()
model.solve(pulp.GLPK_CMD())
model.solve(pulp.COIN_CMD())
print(pulp.LpStatus[model.status])

Update: Thanks for the comments. After removing constraints this works now.