Blending Problem

Problem definition

The Kahuna company manufactures sausages using three kinds of meat. The relevant information about the ingredients is provided in the table below:

Ingredient

Cost (€/kg)

Availability (kg)

Pork

4.32

30

Wheat

2.46

20

Starch

1.86

17

The company makes two types of sausages:

  • Economy (>=40% Pork)

  • Premium (>=60% Pork)

One sausage is 50 grams (0.05 kg)

According to government regulations, the most starch we can use in our sausages is 25%

We have a contract with a butcher, and have already purchased 23 kg pork, that will go bad if it’s not used.

We have a demand for 350 economy sausages and 500 premium sausages.

Write a linear program to figure out how to most cost effectively blend our sausages.

Let’s model our problem

pe = Pork in the economy sausages (kg)
we = Wheat in the economy sausages (kg)
se = Starch in the economy sausages (kg)
pp = Pork in the premium sausages (kg)
wp = Wheat in the premium sausages (kg)
sp = Starch in the premium sausages (kg)

We want to minimise costs such that:

Cost = 4.32(pe + pp) + 2.46(we + wp) + 1.86(se + sp)

With the following constraints:

pe + we + se = 350 * 0.05
pp + wp + sp = 500 * 0.05
pe ≥ 0.4(pe + we + se)
pp ≥ 0.6(pp + wp + sp)
se ≤ 0.25(pe + we + se)
sp ≤ 0.25(pp + wp + sp)
pe + pp ≤ 30
we + wp ≤ 20
se + sp ≤ 17
pe + pp ≥ 23
[ ]:

!pip install pandas
!pip install numpy
!pip install pulp

[1]:

import pandas as pd
import numpy as np
from IPython.display import display, Markdown
import pulp

[2]:

# Instantiate our problem class
model = pulp.LpProblem("Cost minimising blending problem", pulp.LpMinimize)

Here we have 6 decision variables, we could name them individually but this wouldn’t scale up if we had hundreds/thousands of variables (you don’t want to be entering all of these by hand multiple times).

We’ll create a couple of lists from which we can create tuple indices.

[3]:

# Construct our decision variable lists
sausage_types = ['economy', 'premium']
ingredients = ['pork', 'wheat', 'starch']

Each of these decision variables will have similar characteristics (lower bound of 0, continuous variables). Therefore we can use PuLP’s LpVariable object’s dict functionality, we can provide our tuple indices.

These tuples will be keys for the ing_weight dict of decision variables

[4]:

variables = pulp.LpVariable.dicts("weight kg",
                                     ((i, j) for i in sausage_types for j in ingredients),
                                     lowBound=0,
                                     cat='Continuous')

PuLP provides an lpSum vector calculation for the sum of a list of linear expressions.

Whilst we only have 6 decision variables, I will demonstrate how the problem would be constructed in a way that could be scaled up to many variables using list comprehensions.

[5]:

#coefficients
coefficients = [4.32, 2.46, 1.86]

# Objective Function
model += (
    pulp.lpSum([
        coefficients[j] * variables[(i, ingredients[j])]
        for i in sausage_types for j in range(len(ingredients))])
)

Now we add our constraints, bear in mind again here how the use of list comprehensions allows for scaling up to many ingredients or sausage types

[6]:

# Constraints
# 350 economy and 500 premium sausages at 0.05 kg
model += pulp.lpSum([variables['economy', j] for j in ingredients]) == 350 * 0.05, "Economy demand"
model += pulp.lpSum([variables['premium', j] for j in ingredients]) == 500 * 0.05, "Premium demand"

# Economy has >= 40% pork, premium >= 60% pork
model += variables['economy', 'pork'] >= (
    0.4 * pulp.lpSum([variables['economy', j] for j in ingredients])), "40% Pork in Economy"

model += variables['premium', 'pork'] >= (
    0.6 * pulp.lpSum([variables['premium', j] for j in ingredients])), "60$ Pork in Premium"

# Sausages must be <= 25% starch
model += variables['economy', 'starch'] <= (
    0.25 * pulp.lpSum([variables['economy', j] for j in ingredients])), "25% Starch in Economy"

model += variables['premium', 'starch'] <= (
    0.25 * pulp.lpSum([variables['premium', j] for j in ingredients])), "25% Starch in Premium"

# We have at most 30 kg of pork, 20 kg of wheat and 17 kg of starch available
model += pulp.lpSum([variables[i, 'pork'] for i in sausage_types]) <= 30, "Pork Availability"
model += pulp.lpSum([variables[i, 'wheat'] for i in sausage_types]) <= 20, "Wheat Availability"
model += pulp.lpSum([variables[i, 'starch'] for i in sausage_types]) <= 17, "Starch Availability"

# We have at least 23 kg of pork to use up
model += pulp.lpSum([variables[i, 'pork'] for i in sausage_types]) >= 23, "Pork Stock"

[7]:

# Solve our problem
model.solve(solver=pulp.solvers.GUROBI(msg = 0))
pulp.LpStatus[model.status]

Academic license - for non-commercial use only

[7]:

'Optimal'

[8]:

total_cost = pulp.value(model.objective)
display(Markdown("Total cost is %0.2f €"%total_cost))

display(Markdown("The following table shows the decision variables: "))
var_df = pd.DataFrame.from_dict(variables, orient="index",
                                columns = ["Variables"])
var_df["Solution (GRB)"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.solverVar.X))
var_df["Reduced cost (GRB)"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.solverVar.RC))
var_df["Objective Coefficient (GRB)"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.solverVar.Obj))
var_df["Objective Lower bound (GRB)"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.solverVar.SAObjLow) if item.solverVar.SAObjLow > -0.1 else "-Inf" )
var_df["Objective Upper bound (GRB)"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.solverVar.SAObjUp) if item.solverVar.SAObjUp != item.solverVar.UB else "Inf")


display(var_df)


const_dict = dict(model.constraints)
con_df = pd.DataFrame.from_records(list(const_dict.items()), exclude=["Expression"], columns=["Constraint", "Expression"])
con_df["Right Hand Side"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.RHS))
con_df["Slack"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.Slack))
con_df["Shadow Price"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.Pi))
con_df["Min RHS"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.SARHSLow) if const_dict[item].solverConstraint.SARHSLow > -1e10 else "-Inf")
con_df["Max RHS"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.SARHSUp) if const_dict[item].solverConstraint.SARHSUp < 1e10 else "Inf" )


display(Markdown("The following table shows the constraints: "))
display(con_df)

Total cost is 140.96 €

The following table shows the decision variables:

Variables Solution (GRB) Reduced cost (GRB) Objective Coefficient (GRB) Objective Lower bound (GRB) Objective Upper bound (GRB)
(economy, pork) weight_kg_('economy',_'pork') 7.00 0.00 4.32 4.32 Inf
(economy, wheat) weight_kg_('economy',_'wheat') 6.12 0.00 2.46 1.86 2.46
(economy, starch) weight_kg_('economy',_'starch') 4.38 0.00 1.86 -Inf 2.46
(premium, pork) weight_kg_('premium',_'pork') 16.00 0.00 4.32 2.46 4.32
(premium, wheat) weight_kg_('premium',_'wheat') 2.75 0.00 2.46 2.46 4.32
(premium, starch) weight_kg_('premium',_'starch') 6.25 0.00 1.86 -Inf 2.46

The following table shows the constraints:

Constraint Right Hand Side Slack Shadow Price Min RHS Max RHS
0 Economy_demand 17.50 0.00 2.31 10.62 20.00
1 Premium_demand 25.00 0.00 2.31 21.33 26.67
2 40%_Pork_in_Economy 0.00 0.00 0.00 -2.75 1.00
3 60$_Pork_in_Premium 0.00 -1.00 0.00 -Inf 1.00
4 25%_Starch_in_Economy 0.00 0.00 -0.60 -4.38 6.12
5 25%_Starch_in_Premium 0.00 0.00 -0.60 -6.25 2.75
6 Pork_Availability 30.00 7.00 0.00 23.00 Inf
7 Wheat_Availability 20.00 11.12 0.00 8.88 Inf
8 Starch_Availability 17.00 6.38 0.00 10.62 Inf
9 Pork_Stock 23.00 0.00 1.86 22.00 25.75