title: How to set MIP start (initial solution) with Gurobi solver from PuLP?
body: <p>I'm using the <code>PuLP</code> module in Python to formulate a mixed integer program. I am trying to work out how to set a <code>MIP start</code> (i.e. a feasible solution for the program to start from) via the <code>PuLP</code> interface.</p>
  
  <p>Details on how to set <code>MIP start</code> are given <a href="https://www.gurobi.com/documentation/6.5/examples/mip_starts.html" rel="noreferrer">here</a></p>
  
  <p>And the developer of the <code>PuLP</code> package claims that you can access the full Gurobi model via the <code>PuLP</code> interface <a href="https://groups.google.com/forum/#!topic/pulp-or-discuss/-aiL7jQyD-k" rel="noreferrer">here</a></p>
  
  <p>Pasted below are two complete models. I have made these as small as possible whilst preventing the gurobi solver from finding the optimal value using a heuristic.</p>
  
  <p>I have attempted to set an initial solution (to the optimal values) in both models, but in the <code>PuLP</code> model it is ignored, but in the <code>gurobipy</code> model it works as expected.</p>
  
  <p>How do you set an initial solution for the Gurobi solve via the PuLP interface?</p>
  
  <pre><code>from pulp import *
  
  prob = LpProblem("min example",LpMinimize)
  
  x1=LpVariable("x1",0,None,LpInteger)
  x2=LpVariable("x2",0,None,LpInteger)
  x3=LpVariable("x3",0,None,LpInteger)
  x4=LpVariable("x4",0,None,LpInteger)
  
  # Objective function
  prob += 3*x1 + 5*x2 + 6*x3 + 9*x4
  
  # A constraint
  prob += -2*x1 + 6*x2 -3*x3 + 4*x4 &gt;= 2, "Con1"
  prob += -5*x1 + 3*x2 + x3 + 3*x4 &gt;= -2, "Con2"
  prob += 5*x1 - x2 + 4*x3 - 2*x4 &gt;= 3, "Con3"
  
  # Choose solver, and set it to problem, and build the Gurobi model
  solver = pulp.GUROBI()
  prob.setSolver(solver)
  prob.solver.buildSolverModel(prob)
  
  # Attempt to set an initial feasible solution (in this case to an optimal solution)
  prob.solverModel.getVars()[0].start = 1
  prob.solverModel.getVars()[1].start = 1
  prob.solverModel.getVars()[2].start = 0
  prob.solverModel.getVars()[3].start = 0
  
  # Solve model
  prob.solve()
  
  # Status of the solution is printed to the screen
  print "Status:", LpStatus[prob.status]
  
  # Each of the variables is printed with it's resolved optimum value
  for v in prob.variables():
      print v.name, "=", v.varValue
  
  # Optimised objective function value is printed to the screen
  print "OF = ", value(prob.objective)
  </code></pre>
  
  <p>Which returns:</p>
  
  <pre><code>Optimize a model with 3 rows, 4 columns and 12 nonzeros
  Coefficient statistics:
    Matrix range    [1e+00, 6e+00]
    Objective range [3e+00, 9e+00]
    Bounds range    [0e+00, 0e+00]
    RHS range       [2e+00, 3e+00]
  Found heuristic solution: objective 12
  Presolve removed 0 rows and 1 columns
  Presolve time: 0.00s
  Presolved: 3 rows, 3 columns, 9 nonzeros
  Variable types: 0 continuous, 3 integer (0 binary)
  
  Root relaxation: objective 7.400000e+00, 1 iterations, 0.00 seconds
  
      Nodes    |    Current Node    |     Objective Bounds      |     Work
   Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
  
       0     0    7.40000    0    1   12.00000    7.40000  38.3%     -    0s
  H    0     0                       8.0000000    7.40000  7.50%     -    0s
  
  Explored 0 nodes (1 simplex iterations) in 0.00 seconds
  Thread count was 8 (of 8 available processors)
  
  Optimal solution found (tolerance 1.00e-04)
  Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0%
  ('Gurobi status=', 2)
  Status: Optimal
  x1 = 1.0
  x2 = 1.0
  x3 = -0.0
  x4 = -0.0
  OF =  8.0
  </code></pre>
  
  <p>Secondly I can implement the same model using the <code>gurobipy</code> module, but in this case the MIP start is actually used:</p>
  
  <pre><code>from gurobipy import *
  
  m = Model("min example")
  m.modelSense = GRB.MINIMIZE
  
  objFcnCoeffs = [3, 5, 6, 9]
  xVars = []
  for i in range(4):
      xVars.append(m.addVar(vtype=GRB.INTEGER, obj=objFcnCoeffs[i], name="Open%d" % i))
  
  # Update model to integrate new variables
  m.update()
  
  # Constraints
  m.addConstr(-2*xVars[0] + 6*xVars[1] -3*xVars[2] + 4*xVars[3] &gt;= 2, "Con1")
  m.addConstr(-5*xVars[0] + 3*xVars[1] + xVars[2] + 3*xVars[3] &gt;= -2, "Con2")
  m.addConstr(5*xVars[0] - xVars[1] + 4*xVars[2] - 2*xVars[3] &gt;= 3, "Con3")
  
  
  # Attempt to set an initial feasible solution (in this case to an optimal solution)
  startValues = [1, 1, 0, 0]
  for i in range(4):
      xVars[i].start = startValues[i]
  
  # Solve model
  m.optimize()
  
  # Print solution
  print('\nTOTAL COSTS: %g' % m.objVal)
  for i in range(4):
      print('\n xVar[%s] = %g' % i, xVars[i])
  </code></pre>
  
  <p>Which returns:</p>
  
  <pre><code>Optimize a model with 3 rows, 4 columns and 12 nonzeros
  Coefficient statistics:
    Matrix range    [1e+00, 6e+00]
    Objective range [3e+00, 9e+00]
    Bounds range    [0e+00, 0e+00]
    RHS range       [2e+00, 3e+00]
  Found heuristic solution: objective 12
  Presolve removed 0 rows and 1 columns
  Presolve time: 0.00s
  Presolved: 3 rows, 3 columns, 9 nonzeros
  
  Loaded MIP start with objective 8
  
  Variable types: 0 continuous, 3 integer (0 binary)
  
  Root relaxation: infeasible, 0 iterations, 0.00 seconds
  
  Explored 0 nodes (0 simplex iterations) in 0.00 seconds
  Thread count was 8 (of 8 available processors)
  
  Optimal solution found (tolerance 1.00e-04)
  Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0%
  
  TOTAL COSTS: 8
  
   xVar[0] = 1
  
   xVar[1] = 1
  
   xVar[2] = 0
  
   xVar[3] = 0
  </code></pre>