Blending Craft Beer

Introduction

In this notebook we will use:

We wil learn how to model problems in a way that scales to hundreds or thousands of variables and how to get all valuable information from the results generated by PuLP.

You cannot test this notebook unless you have a Gurobi license in your local environment. You also need to install the package gurobipy. If you do not have a Gurobi license, you can request a license for free at Gurobi Academic License. Run the following cell to install gurobipy.

[ ]:

!pip install gurobipy

Problem Definition

Duff Beer is an american brewerie that sells alcoholic beverages worldwide. Duff beer has a very large market share in the US. One of the new business lines of Duff Beer is in the craft beer sector. Duff beer buys 3 types of wort from local suppliers. The availability of wort (in liters) and the price for each of them is specified in Table 1.

Table 1: daily availability and cost of worts

Wort

Availability (liters)

Cost (€cents/liter)

Type 1

2500

120

Type 2

1200

095

Type 3

2000

150

The blend of the three types of worts must comply with the following quality requirements:

  • For Beer Angels, no less than 60% of Type 3 and no more of 20% of Type 2. The price of one bottle of Beer Angel is 4.10€s/liter.

  • For Beer Beast, no less than 15% of Type 3 and no more than 60% of Type 2. The price of one bottle of Beer Beast is 2.80€s/liter.

  • For Beer Cactus, no more than 50% of type 2. The selling price is 2.45€/liter.

Write a linear program to find the most profitable blend of worts to elaborate the three kinds of beers

Problem Model

Indices

We can define the following indices to express the problem in a compact form:

  • i: wort type \(i \in [1, 2, 3]\)

  • j: Beer brand \(j \in [A, B, C]\) where A notes Angels, B notes Beast, and C notes Cactus

Decision variables

The decision variables are:

\(x_{ij}\): Amount of wort of type i (1, 2, 3) in Beer j: A (Angels), B (Beast), C (Cactus), (ie i=[1,2,3], j=[A,B,C])

Objective Function

The objective function is:

\[\max z = \sum_{i=1}^{3}\sum_{j=A}^{C}(d_j - c_i)*x_{ij}\]

where \(d_j\) is hte selling price of a bottle of beer j and \(c_i\) is the cost per liter of wort type i.

\[\max z=410(x_{1A}+x_{2A}+x_{3A})+280(x_{1B}+x_{2B}+x_{3B})+245(x_{1C}+x_{2C}+x_{3C})-120(x_{1A}+x_{1B}+x_{1C})-95(x_{2A}+x_{2B}+x_{2C})-150(x_{3A}+x_{3B}+x_{3C})\]
\[\max z=290x_{1A}+315x_{2A}+260x_{3A}+160x_{1B}+185x_{2B}+130x_{3B}+125x_{1C}+150x_{2C}+95x_{3C}\]

Constraints

Subject to the following constraints:

Availability constraints The availability constraints can be expressed in a compact form as:

\[\sum_jx_{ij} \leq a_i \forall i\]
\[x_{1A}+x_{1B}+x_{1C} \leq 2500\]
\[x_{2A}+x_{2B}+x_{2C} \leq 1200\]
\[x_{3A}+x_{3B}+x_{3C} \leq 2000\]

Angels quality requirement constraints

\[x_{3A} \geq 0.6(x_{1A}+x_{2A}+x_{3A})\]
\[x_{2A} \leq 0.2(x_{1A}+x_{2A}+x_{3A})\]

Beast quality requirement constraints

\[x_{3B} \geq 0.15(x_{1B}+x_{2B}+x_{3B})\]
\[x_{2B} \leq 0.6(x_{1B}+x_{2B}+x_{3B})\]

Cactus quality requirement constraints

\[x_{2C} \leq 0.5(x_{1C}+x_{2C}+x_{3C})\]
[ ]:



[ ]:

import pandas as pd
import numpy as np
from IPython.display import display, Markdown
import pulp
# Uncomment to check if Gurobi is installed
# print(pulp.list_solvers(onlyAvailable=True))

# Instantiate our problem class
model = pulp.LpProblem("Maximising profits for Duff Beer", pulp.LpMaximize)
worts = ['Type 1', 'Type 2', 'Type 3']
beers = ['Angel', 'Beast', 'Cactus']

variables = pulp.LpVariable.dicts("Quantity in liters",
                                     ((i, j) for i in worts for j in beers),
                                     lowBound=0,
                                     cat='Continuous')
price = [410, 280, 245]
cost = [120, 95, 150]

coefficients = [price[j]-cost[i] for j in range(len(price)) for i in range(len(cost))]

# Objective Function
model += (
    pulp.lpSum([
        ((price[j]-cost[i]) * variables[(worts[i], beers[j])]
        for i in range(len(worts)) for j in range(len(beers)))])
)

# Constraints
model += pulp.lpSum([variables['Type 1', j] for j in beers]) <= 2500, "Type 1 availability"
model += pulp.lpSum([variables['Type 2', j] for j in beers]) <= 1200, "Type 2 availability"
model += pulp.lpSum([variables['Type 3', j] for j in beers]) <= 2000, "Type 3 availability"

model += variables['Type 3', 'Angel']>=0.6*pulp.lpSum([variables[i,'Angel'] for i in worts]), "Angel >=60% Type 3 requirement"
model += variables['Type 2', 'Angel']<=0.2*pulp.lpSum([variables[i,'Angel'] for i in worts]), "Angel <=20% Type 2 requirement"

model += variables['Type 3', 'Beast']>=0.15*pulp.lpSum([variables[i,'Beast'] for i in worts]), "Beast >=15% Type 3 requirement"
model += variables['Type 2', 'Beast']<=0.6*pulp.lpSum([variables[i,'Beast'] for i in worts]), "Beast <=60% Type 2 requirement"

model += variables['Type 2', 'Cactus']<=0.5*pulp.lpSum([variables[i,'Angel'] for i in worts]), "Cactus <=50% Type 2 requirement"


model.solve(solver=pulp.GUROBI(msg = 0))
pulp.LpStatus[model.status]

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: "{:.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["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)