Blending Craft Beer

Try me

Open In ColabBinder

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.

[ ]:

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

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 the selling price of a bottle of beer j and \(c_i\) is the cost per liter of wort type i.

\(\max z =4 10(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\)

Where \(a_i\) is the availability of wort type 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\)

Quality constraints The proportion of a given type of wort i’ in a type of beer j’ can be expressed as:

\(\frac{x_{i'j'}}{\sum_i(x_{ij'})} \forall i', j'\)

If we compare this to the minimum proportion \(rmin_{i'j'}\) we obtain:

\(\frac{x_{i'j'}}{\sum_i(x_{ij'})} \geq rmin_{i'j'} \forall i', j'\)

Or

\(x_{i'j'} \geq rmin_{i'j'}*\sum_i(x_{ij'}) \forall i', j'\)

If we take the denominator to the RHS. We can obtain a similar expression for the maximum proportion \(rmax_{i'j'}\):

\(x_{i'j'} \geq rmax_{i'j'}*\sum_i(x_{ij'}) \forall i', j'\)

From these expressions, we can derive the constraints for every type of beer:

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
# 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()
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"])

# First we add the solution. We apply a lambda function to get only two decimals:
var_df["Solution"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.varValue))
# We do the same for the reduced cost:
var_df["Reduced cost"] = var_df["Variables"].apply(lambda item: "{:.2f}".format(item.dj))


# We use the display function to represent the results:
display(var_df)


const_dict = dict(model.constraints)
con_df = pd.DataFrame.from_records(list(const_dict.items()), exclude=["Expression"], columns=["Constraint", "Expression"])

#We create a list of records from the dictionary and exclude the Expression to have a more compact solution.
con_df = pd.DataFrame.from_records(list(const_dict.items()), exclude=["Expression"], columns=["Constraint", "Expression"])

#Now we add columns for the solution, the slack and shadow price

con_df["Right Hand Side"] = con_df["Constraint"].apply(lambda item: "{:.2f}".format(-const_dict[item].constant))
con_df["Slack"] = con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].slack))
con_df["Shadow Price"] = con_df["Constraint"].apply(lambda item: "{:.2f}".format(const_dict[item].pi))

# And we display the results
display(con_df)

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

[ ]: