Pulp problem gets the right objective value but fails to find the minimal solution

Question

I converted an excel solver problem to Pulp. However, while the pulp program calculates the objective function value correctly, it doesn't arrive at the right answer. I can't figure out what is wrong in the code. Can anyone identify what I may be doing wrong here?

*edit: Have a look at the problem in excel here. Solver finds a much better solution.

I know the objective functions in pulp and excel are both the same, since if I plug in the pulp solution in the excel model, excel gives me the same objective function value as pulp. But for some reason, pulp cannot find the low cost that excel does.

Thank you again for your help!

from pulp import *
import pandas as pd

init_inventory = 500
init_workers = 100

reg_work_hours = 160
max_ot_hours = 20

cost_hiring = 1600
cost_firing = 2000

monthly_wage = 1500
ot_wage_rate = 13

labor_hours_per_shoe = 4
mat_cost_per_shoe = 15
hold_cost_per_month = 3

Months = ["Month 1", "Month 2", "Month 3", "Month 4"]

ShoesNeeded = dict(zip(Months, [3000, 5000, 2000, 1000]))


prob = LpProblem("ProductionPlanning", LpMinimize)

hire = LpVariable.dicts("hire", Months, 0, None, LpInteger)
fire = LpVariable.dicts("fire", Months, 0, None, LpInteger)

avail_workers = {Months[0] : init_workers + hire[Months[0]] - fire[Months[0]]}

for m in range(1,4):
  avail_workers[Months[m]] = avail_workers[Months[m-1]] + hire[Months[m]] - fire[Months[m]]

reg_hours = {m: avail_workers[m] * reg_work_hours for m in Months}
max_ot_avail = {m: avail_workers[m] * max_ot_hours for m in Months}

ot_used = LpVariable.dicts("ot_used", Months, 0, None, LpInteger)

total_labor_hours = {m: reg_hours[m] + ot_used[m] for m in Months}

prod_capacity = {m: total_labor_hours[m] / labor_hours_per_shoe for m in Months}

shoes_produced = LpVariable.dicts("Shoes_produced", Months, 0, None, LpInteger)

invent_ending = {Months[0] : init_inventory + shoes_produced[Months[0]] - ShoesNeeded[Months[0]]}

for m in range(1,4):
  invent_ending[Months[m]] = invent_ending[Months[m-1]] + shoes_produced[Months[m]] - ShoesNeeded[Months[m]]

prob += lpSum(hire[m] * cost_hiring +
              fire[m] * cost_firing + 
              avail_workers[m] * monthly_wage +
              ot_used[m] * ot_wage_rate + 
              shoes_produced[m] * mat_cost_per_shoe +
              invent_ending[m] * hold_cost_per_month for m in Months), "Total Cost"

for m in Months:
  prob += ot_used[m] <= max_ot_avail[m], f"Max OT for {m}"
  prob += shoes_produced[m] <= prod_capacity[m], f"Production capacity for {m}"
  prob += invent_ending[m] >= ShoesNeeded[m], f"Minimum quantity needed for {m}"

prob.solve()

for v in prob.variables():
    print(v.name, "=", v.varValue)

print("Total Cost = ", value(prob.objective))

Answer 1

How do you define the "right" answer?

If the objective function value is computed correctly and reaches minimum value, then it is possible that the optimal solution is not unique, and that the solver you have used with pulp (CBC) simply found another optimal solution. There is no reason to say it is not "the right answer", as long as the minimal value is reached and that all constraints are satisfied.

If you wish to have a different optimal value, you can add no good cuts for integer variables to exclude the current solution.

You could also add some additional constraints to enforce part of the solution you found in Excel, and check if that guides the solver in pulp (CBC) to the same solution.

Answer 2

Avail_workers seems to be an array of variable, indexed by month. It's like a derived decision variable as it depends upon two other decision variables- hire and fire
Add a constraint with $ availworker_m$ being a variable (no need to declare it integer)
Add this constraint
$ availWorker_m = availWorker_{m-1}+hire_m-fire_m \ \ \forall m \in$ Months[1:]
Also, as already there the initial month[0]:
$ availworker_0 = initWorker+hire_0$

Modify the below constraints
$shoesproduced_m\times labourHoursPerShoe \le availWorker_m\times regWorkHours +otUsed_m$