Petroleum Blending (Solved)¶
Try me¶
Problem Definition¶
Harry Stamper is an engineer working on a petroleum company. The company produces three grades of motor oil – super, premium, and extra. The three grades are made from the same 3 ingredients, named component 1, component 2 and component 3. The company wants to determine the optimal mix of the 3 components in each grade of motor oil that will maximize profit.
The maximum availability of each component and their cost per barrel are as follows:
Component |
Max barrels available / day |
Cost / barrel |
|---|---|---|
Component 1 |
4.500 |
12€ |
Component 2 |
2.700 |
10€ |
Component 3 |
3.500 |
14€ |
To ensure the appropriate blend, each grade must have a minimum amount of component 1 plus a combination of other components as follows:
Grade |
Component 1 |
Component 2 |
Component 3 |
Selling price/barrel |
|---|---|---|---|---|
Super |
At least 50% |
No more than 30% |
23€ |
|
Premium |
At least 40% |
No more than 25% |
20€ |
|
Extra |
At least 60% |
At least 10% |
18€ |
The company wants to make at least 3000 barrels of each blend. Formulate a LP to help Harry Stamper find the optimal blend that maximises profit.
Problem Model¶
The decision variables are:
\(x_{ij}\) = amount in barrels of component i in grade j
Objective function
max \(z = 23*(x_{11} + x_{21} + x_{31}) + 20*(x_{12} + x_{22} + x_{32})+18*(x_{13} + x_{23} + x_{33})-12*(x_{11} + x_{12} + x_{13}) - 10*(x_{21} + x_{22} + x_{23})-14*(x_{31} + x_{32} + x_{33})\)
Constraints
\(x_{11} \geq 0.5*(x_{11} + x_{21} + x_{31})\) Component 1 requirement in Super
\(x_{21} \leq 0.3*(x_{11} + x_{21} + x_{31})\) Component 2 requirement in Super
\(x_{12} \geq 0.4*(x_{12} + x_{22} + x_{32})\) Component 1 requirement in Premium
\(x_{32} \leq 0.25*(x_{12} + x_{22} + x_{32})\) Component 3 requirement in Premium
\(x_{13} \geq 0.6*(x_{13} + x_{23} + x_{33})\) Component 1 requirement in Extra
\(x_{23} \geq 0.1*(x_{13} + x_{23} + x_{33})\) Component 2 requirement in Extra
\(x_{11}+x_{12}+x_{13} \leq 4500\) Component 1 availability
\(x_{21}+x_{22}+x_{23} \leq 2700\) Component 1 availability
\(x_{31}+x_{32}+x_{33} \leq 3500\) Component 1 availability
\(x_{11}+x_{21}+x_{31} \geq 3000\) Super minimum production
\(x_{12}+x_{22}+x_{32} \geq 3000\) Premium minimum production
\(x_{13}+x_{23}+x_{33} \geq 3000\) Extra minimum production
[ ]:
!pip install pandas
!pip install numpy
!pip install pulp
[50]:
import pandas as pd
import numpy as np
from IPython.display import display, Markdown
import pulp
[51]:
# Instantiate our problem class
model = pulp.LpProblem("Maximising profits for Harry Stamper", pulp.LpMaximize)
[52]:
# Construct our decision variable lists
components = ['Component 1', 'Component 2', 'Component 3']
blends = ['Super', 'Premium', 'Extra']
#decision variables
variables = pulp.LpVariable.dicts("Quantity in barrels",
((i, j) for i in blends for j in components),
lowBound=0,
cat='Continuous')
# barrels is Super Component 1, super component 2, ..., extra component 3: x11, x21, x31, ...
cost = [12, 10, 14]
price = [23, 20, 18]
coefficients =[price[i] - cost[j] for i in range(len(price)) for j in range(len(cost))]
# Objective Function
model += (
pulp.lpSum([
((price[i]-cost[j]) * variables[(blends[i], components[j])]
for i in range(len(blends)) for j in range(len(components)))])
)
# Requirements
model += variables['Super','Component 1']>= 0.5*pulp.lpSum([variables['Super',components[j]] for j in range(len(components))]), "Component 1 requirement in super"
model += variables['Super','Component 2']<= 0.3*pulp.lpSum([variables['Super',components[j]] for j in range(len(components))]), "Component 2 requirement in super"
model += variables['Premium','Component 1']>= 0.4*pulp.lpSum([variables['Premium',components[j]] for j in range(len(components))]), "Component 1 requirement in premium"
model += variables['Premium','Component 3']<= 0.25*pulp.lpSum([variables['Premium',components[j]] for j in range(len(components))]), "Component 3 requirement in premium"
model += variables['Extra','Component 1']>= 0.6*pulp.lpSum([variables['Extra',components[j]] for j in range(len(components))]), "Component 1 requirement in Extra"
model += variables['Extra','Component 2']>= 0.1*pulp.lpSum([variables['Extra',components[j]] for j in range(len(components))]), "Component 2 requirement in Extra"
model += pulp.lpSum([variables[blends[i], 'Component 1'] for i in range(len(blends))]) <= 4500, "Component 1 availability"
model += pulp.lpSum([variables[blends[i], 'Component 2'] for i in range(len(blends))]) <= 2700, "Component 2 availability"
model += pulp.lpSum([variables[blends[i], 'Component 3'] for i in range(len(blends))]) <= 3500, "Component 3 availability"
model += pulp.lpSum([variables['Super', components[j]] for j in range(len(components))]) >= 3000, "Super minimum production"
model += pulp.lpSum([variables['Premium', components[j]] for j in range(len(components))]) >= 3000, "Premium minimum production"
model += pulp.lpSum([variables['Extra', components[j]] for j in range(len(components))]) >= 3000, "Extra minimum production"
[53]:
# Solve our problem
model.solve(solver=pulp.solvers.GUROBI(msg = 0))
pulp.LpStatus[model.status]
[53]:
'Optimal'
[54]:
total_profit = pulp.value(model.objective)
display(Markdown("Total profit is %0.2f €"%total_profit))
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: "{:.3f}".format(item.solverVar.X))
var_df["Reduced cost (GRB)"] = var_df["Variables"].apply(lambda item: "{:.3f}".format(item.solverVar.RC))
var_df["Objective Coefficient (GRB)"] = var_df["Variables"].apply(lambda item: "{:.3f}".format(item.solverVar.Obj))
var_df["Objective Lower bound (GRB)"] = var_df["Variables"].apply(lambda item: "{:.3f}".format(item.solverVar.SAObjLow) if item.solverVar.SAObjLow > -0.1 else "-Inf" )
var_df["Objective Upper bound (GRB)"] = var_df["Variables"].apply(lambda item: "{:.3f}".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["Shadow Price"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.Pi))
con_df["Slack"]=con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].solverConstraint.Slack))
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 profit is 76800.00 €
The following table shows the decision variables:
| Variables | Solution (GRB) | Reduced cost (GRB) | Objective Coefficient (GRB) | Objective Lower bound (GRB) | Objective Upper bound (GRB) | |
|---|---|---|---|---|---|---|
| (Super, Component 1) | Quantity_in_barrels_('Super',_'Component_1') | 1500.000 | 0.000 | 11.000 | 8.000 | Inf |
| (Super, Component 2) | Quantity_in_barrels_('Super',_'Component_2') | 450.000 | 0.000 | 13.000 | 13.000 | 13.000 |
| (Super, Component 3) | Quantity_in_barrels_('Super',_'Component_3') | 1050.000 | 0.000 | 9.000 | 9.000 | 9.000 |
| (Premium, Component 1) | Quantity_in_barrels_('Premium',_'Component_1') | 1200.000 | 0.000 | 8.000 | -Inf | 11.000 |
| (Premium, Component 2) | Quantity_in_barrels_('Premium',_'Component_2') | 1050.000 | 0.000 | 10.000 | -Inf | 10.000 |
| (Premium, Component 3) | Quantity_in_barrels_('Premium',_'Component_3') | 750.000 | 0.000 | 6.000 | 6.000 | 10.800 |
| (Extra, Component 1) | Quantity_in_barrels_('Extra',_'Component_1') | 1800.000 | 0.000 | 6.000 | -Inf | 17.333 |
| (Extra, Component 2) | Quantity_in_barrels_('Extra',_'Component_2') | 1200.000 | 0.000 | 8.000 | 8.000 | 25.000 |
| (Extra, Component 3) | Quantity_in_barrels_('Extra',_'Component_3') | 0.000 | -0.000 | 4.000 | -Inf | 4.000 |
The following table shows the constraints:
| Constraint | Right Hand Side | Shadow Price | Slack | Min RHS | Max RHS | |
|---|---|---|---|---|---|---|
| 0 | Component_1_requirement_in_super | 0.00 | -18.00 | 0.00 | -850.00 | 0.00 |
| 1 | Component_2_requirement_in_super | 0.00 | 0.00 | 450.00 | -450.00 | Inf |
| 2 | Component_1_requirement_in_premium | 0.00 | -18.00 | 0.00 | -450.00 | 0.00 |
| 3 | Component_3_requirement_in_premium | 0.00 | 0.00 | 0.00 | -450.00 | 450.00 |
| 4 | Component_1_requirement_in_Extra | 0.00 | -18.00 | 0.00 | -450.00 | 0.00 |
| 5 | Component_2_requirement_in_Extra | 0.00 | 0.00 | -900.00 | -Inf | 900.00 |
| 6 | Component_1_availability | 4500.00 | 20.00 | 0.00 | 4500.00 | 6200.00 |
| 7 | Component_2_availability | 2700.00 | 4.00 | 0.00 | 2250.00 | 3150.00 |
| 8 | Component_3_availability | 3500.00 | 0.00 | 1700.00 | 1800.00 | Inf |
| 9 | Super_minimum_production | 3000.00 | 0.00 | -0.00 | -Inf | 3000.00 |
| 10 | Premium_minimum_production | 3000.00 | -1.20 | 0.00 | -0.00 | 3000.00 |
| 11 | Extra_minimum_production | 3000.00 | -6.80 | 0.00 | 0.00 | 3000.00 |