title: PuLP integer programming: solver returns 'Undefined' for anything beyond toy example
body: <p>I am trying to find an integer programming formulation for the following problem: </p>
  
  <p>Sort a list of items according to 2 criteria, represented by the following cost: </p>
  
  <ul>
  <li>'Positional cost': dependent on the position of the item in the list</li>
  <li>'Neighbour cost': dependent on two items being positioned next to eachother in the list (ie the items are immediate neighbours) </li>
  </ul>
  
  <p>The following apply:</p>
  
  <ul>
  <li>Every item can only be placed once </li>
  <li>The list has a beginning position P_1 (-> the item in that position does not have a preceeding list item) </li>
  <li>The list has an end position P_X (-> the item in that position does not have a subsequent item) </li>
  </ul>
  
  <p>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?<br>
  Any input and comments are highly appreciated. </p>
  
  <pre><code>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 : 
  #     -&gt; 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): 
  #     -&gt; 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"  
  #     -&gt; (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 &gt;= 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])
  </code></pre>
  
  <p><strong>Update:</strong>
  Thanks for the comments. After removing constraints this works now.  </p>