Maximize daily profit problem using Python and Pulp

I am having trouble while trying to maximize the daily profit of a chocolate factory.

Problem: Maximize daily profit of a factory that produces individual chocolates(each one having a profit and a maximum daily production capacity), and special packs (each of them containing 3 chocolates and a profit). Example:

# Limit Variables
# n - number of chocolates
# p - the number of special packs
# max_chocolates - the total maximum chocolate production
# capacity in a day 
n = 5
p = 2
max_chocolates = 150

# For the chocolates: Each chocolate (1,2,3,4,5) has a profit and a maximum daily production capacity
profit_chocolates = [50, 30, 45, 40, 35]
capacity_chocolates = [27, 33, 30, 37, 35]

# For the special packs: Each pack is made of 3 chocolates and a profit for said pack
pack_components = [(1, 3, 5), (2, 3, 4)]
profit_packs = [130, 130]

Next, I present the code to get the maximum daily profit:

# Problem initialization
prob = LpProblem("Maximize_chocolates_and_packs_profit", LpMaximize)

# Decision Variables:
chocolate_vars = [LpVariable(f"x{i}", lowBound=0, cat="Integer") for i in range(n)]
packs_vars = [LpVariable(f"p{i}", lowBound=0, cat="Integer") for i in range(p)]


# Objective Function
prob += lpSum(profit_chocolates[i] * chocolate_vars[i] for i in range(n)) + \
        lpSum(profit_packs[i] * packs_vars[i] for i in range(p))

# Constraints:
# For the maximum number of chocolates we are allowed to produce in a day
prob += lpSum(chocolate_vars[i] for i in range(n)) <= max_chocolates

# For the maximum daily production of each chocolate
for i in range(n):
    prob += chocolate_vars[i] <= capacity_chocolates[i]

# For the special packs:
# Profit Constraint: The profit of each pack should not be lower or equal to the
# sum of the profits of the chocolates said pack contains 
for i in range(p):
    prob += lpSum(profit_chocolates[j - 1] for j in pack_components[i]) >= profit_packs[i]

# Capacity Constraint for the packs: Each pack should take into consideration the
# capacity of the chocolates it contains - if the capacity of single chocolate of the pack
# exceeds the daily production limit of that chocolate then we are not able to produce the pack 
for i in range(p):
    for j in pack_components[i]:
        prob += (chocolate_vars[j - 1]) <= capacity_chocolates[j - 1]

# Each decision variable must be greater or equal to 0
for i in range(n):
    prob += chocolate_vars[i] >= 0

for i in range(p):
    prob += packs_vars[i] >= 0

prob.solve()

for i in range(len(chocolate_vars)):
    print(chocolate_vars[i].varValue)

for i in range(len(packs_vars)):
    print(packs_vars[i].varValue)

print(f"Maximum daily profit: {int(prob.objective.value())}")

For the above input data, the expected result (maximum profit) should be 6440. However, the result I get is 6035.

I think the difference in results has to do with the special packs, since they depend on each chocolate.

Could you please help me find out what I'm missing / doing wrong?

Thanks

Solution

You are on the right track. The core piece you are missing is how to get the correct total of each type of chocolate produced.

What is really needed is some mechanism to interrogate the packs to get the quantity of individual chocolates needed by production of that pack because what you'd really like to do is say:

production[choc] = individual_pieces[choc] + quantity_in_packs[choc]

That way when you do your "for each chocolate" type of constraint, you get the right total. The solution below does that with a dictionary. Actually a dictionary of dictionaries so that you can ask it:

quantity_in_pack[pack][choc]

to get that number, then when you do "for each chocolate" you need to sum across the packs to get the grand total (shown). There are a couple other ways to do this probably... You could probably flip the dictionary around and index by chocolate, then pack, but this seems logical.

Anyhow, you really only have 2 constraints:

total number produced
quantity of each chocolate produced

Notes:

I left some comments in and commented out some unnecessary code
I renamed the sets. Easier on the eyes and avoids the j-1 index pain, etc.
Making lists or sets of indices is easier to read than the range(n) stuff, so I like doing that for readability and then you can do math-like expressions fairly easy: for car in cars etc.
I have zipped a lot of things into dictionaries here, which is easy, or you could just create them that way from scratch.
I hope I didn't just do your winter Linear Programming project for you.

Code:

from pulp import *

# Limit Variables
# n - number of chocolates
# p - the number of special packs
# max_chocolates - the total maximum chocolate production
# capacity in a day 
n = 5
chocolates = ['extra dark', 'milk', 'raspberry', 'caramel', 'sprinkles']
p = 2
packs = ['mega deluxe', 'mammas pride']
max_chocolates = 150

# For the chocolates: Each chocolate (1,2,3,4,5) has a profit and a maximum daily production capacity
profit_chocolates = [50, 30, 45, 40, 35]
profit_chocolates = dict(zip(chocolates, profit_chocolates))
capacity_chocolates = [27, 33, 30, 37, 35]
capacity_chocolates = dict(zip(chocolates, capacity_chocolates))

# For the special packs: Each pack is made of 3 chocolates and a profit for said pack
pack_components = [(1, 3, 5), (2, 3, 4)]
# this will be very helpful...  break down the pack contents and quantities:
pack_breakdown = {  
                    'mega deluxe'  : {chocolates[i-1]: 1 for i in pack_components[0]},
                    'mammas pride' : {chocolates[i-1]: 1 for i in pack_components[1]}
                 }
profit_packs = [130, 130]
profit_packs = dict(zip(packs, profit_packs))


# Problem initialization
prob = LpProblem("Maximize_chocolates_and_packs_profit", LpMaximize)

# Decision Variables:
# Let the make_chocolate be JUST for individual sale, not the aggregate number.
make_chocolate = LpVariable.dicts("make_chocolate", chocolates, lowBound=0, cat="Integer") 
make_pack = LpVariable.dicts("make_pack", packs, lowBound=0, cat="Integer")


# This is GOOD / Correct
# Objective Function
prob += lpSum(profit_chocolates[c] * make_chocolate[c] for c in chocolates) + \
        lpSum(profit_packs[p] * make_pack[p] for p in packs)

# Constraints:
# Need to capture the individual + the total qty from packs
# For the maximum number of chocolates we are allowed to produce in a day
# .values() is used here to sum up the quantities of each type, which 
# *could* be more than 1
prob += lpSum(make_chocolate[c] for c in chocolates) \
        + lpSum(make_pack[p] * sum(pack_breakdown[p].values()) for p in packs) \
        <= max_chocolates

# For the maximum daily production of each chocolate
# we need to add the individuals + the qty from produced packs...
for c in chocolates:
    # we can make a convenient expression to sum up the qty of c in all packs...
    qty_in_packs = lpSum(pack_breakdown[p].get(c, 0) * make_pack[p] for p in packs)

    prob += make_chocolate[c] + qty_in_packs <= capacity_chocolates[c]

# For the special packs:
# THIS SHOULD NOT BE NEEDED.  Your profit expression is good.
# Profit Constraint: The profit of each pack should not be lower or equal to the
# sum of the profits of the chocolates said pack contains 
# for i in range(p):
#     prob += lpSum(profit_chocolates[j - 1] for j in pack_components[i]) >= profit_packs[i]

# Capacity Constraint for the packs: Each pack should take into consideration the
# capacity of the chocolates it contains - if the capacity of single chocolate of the pack
# exceeds the daily production limit of that chocolate then we are not able to produce the pack 
# for i in range(p):
#     for j in pack_components[i]:
#         prob += (make_chocolate[j - 1]) <= capacity_chocolates[j - 1]

# # Each decision variable must be greater or equal to 0
# for i in range(n):
#     prob += make_chocolate[i] >= 0

# for i in range(p):
#     prob += make_pack[i] >= 0

prob.solve()

print('Individual chocolate plan:')
for c in chocolates:
    print(f'  {c: >10}:  {make_chocolate[c].varValue}')

print('\nPack plan:')
for p in packs:
    print(f'  {p: >15}:  {make_pack[p].varValue}')

print('\nProduction Limit check (total/limit)')
for c in chocolates:
    amount_produced = make_chocolate[c].varValue
    amount_produced += sum(pack_breakdown[p].get(c, 0) * make_pack[p].varValue for p in packs)
    print(f'  {c: >10}: {amount_produced}/{capacity_chocolates[c]}')

print(f"\nMaximum daily profit: {int(prob.objective.value())}")

Output:

Result - Optimal solution found

Objective value:                6440.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

Individual chocolate plan:
  extra dark:  27.0
        milk:  0.0
   raspberry:  0.0
     caramel:  7.0
   sprinkles:  26.0

Pack plan:
      mega deluxe:  0.0
     mammas pride:  30.0

Production Limit check (total/limit)
  extra dark: 27.0/27
        milk: 30.0/33
   raspberry: 30.0/30
     caramel: 37.0/37
   sprinkles: 26.0/35

Maximum daily profit: 6440