*This notebook contains material from ND-Pyomo-Cookbook;
content is available on Github.*
*The text is released under the CC-BY-NC-ND-4.0 license,
and code is released under the MIT license.*

Keywords: blending, cbc usage

In [4]:

```
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import shutil
import sys
import os.path
if not shutil.which("pyomo"):
!pip install -q pyomo
assert(shutil.which("pyomo"))
if not (shutil.which("cbc") or os.path.isfile("cbc")):
if "google.colab" in sys.modules:
!apt-get install -y -qq coinor-cbc
else:
try:
!conda install -c conda-forge coincbc
except:
pass
assert(shutil.which("cbc") or os.path.isfile("cbc"))
import pyomo.environ as pyomo
```

A brewery receives an order for 100 gallons of 4% ABV (alchohol by volume) beer. The brewery has on hand beer A that is 4.5% ABV that cost USD 0.32 per gallon to make, and beer B that is 3.7% ABV and cost USD 0.25 per gallon. Water could also be used as a blending agent at a cost of USD 0.05 per gallon. Find the minimum cost blend that meets the customer requirements.

We will use this problem as an opportunity to write a Python function that accepts data on raw materials and customer specifications to produce the lowest cost blend.

The first step is to represent the problem data in a generic manner that could, if needed, be extended to include additional blending components. Here we use a dictionary of materials, each key denoting a blending agent. For each key there is a sub-dictionary containing attributes of each blending component.

In [5]:

```
data = {
'A': {'abv': 0.045, 'cost': 0.32},
'B': {'abv': 0.037, 'cost': 0.25},
'W': {'abv': 0.000, 'cost': 0.05},
}
```

If we let subscript $c$ denote a blending component from the set of blending components $C$, and denote the volume of $c$ used in the blend as $x_c$, the cost of the blend is

\begin{align} \mbox{cost} & = \sum_{c\in C} x_c P_c \end{align}where $P_c$ is the price per unit volume of $c$. Using the Python data dictionary defined above, the price $P_c$ is given by `data[c]['cost']`

.

The customer requirement is produce a total volume $V$. Assuming ideal solutions, the constraint is given by

\begin{align} V & = \sum_{c\in C} x_c \end{align}where $x_c$ denotes the volume of component $c$ used in the blend.

The product composition is specified as 4% alchohol by volume. Denoting this as $\bar{A}$, the constraint may be written as

\begin{align} \bar{A} & = \frac{\sum_{c\in C}x_c A_c}{\sum_{c\in C} x_c} \end{align}where $A_c$ is the alcohol by volume for component $c$. As written, this is a nonlinear constraint. Multiplying both sides of the equation by the denominator yields a linear constraint

\begin{align} \bar{A}\sum_{c\in C} x_c & = \sum_{c\in C}x_c A_c \end{align}A final form for this constraint can be given in either of two versions. In the first version we subtract the left-hand side from the right to give

\begin{align} 0 & = \sum_{c\in C}x_c \left(A_c - \bar{A}\right) & \mbox{ Version 1 of the linear blending constraint} \end{align}Alternatively, the summation on the left-hand side corresponds to total volume. Since that is known as part of the problem specification, the blending constraint could also be written as

\begin{align} \bar{A}V & = \sum_{c\in C}x_c A_c & \mbox{ Version 2 of the linear blending constraint} \end{align}Which should you use? Either will generally work well. The advantage of version 1 is that it is fully specified by a product requirement $\bar{A}$, which is sometimes helpful in writing elegant Python code.

A Pyomo implementation of this blending model is shown in the next cell. The model is contained within a Python function so that it can be more easily reused for additional calculations, or eventually for use by the process operator.

Note that the pyomo library has been imported with the prefix `pyomo`

. This is good programming practive to avoid namespace collisions with problem data.

In [6]:

```
vol = 100
abv = 0.040
def beer_blend(vol, abv, data):
C = data.keys()
model = pyomo.ConcreteModel()
model.x = pyomo.Var(C, domain=pyomo.NonNegativeReals)
model.cost = pyomo.Objective(expr = sum(model.x[c]*data[c]['cost'] for c in C))
model.vol = pyomo.Constraint(expr = vol == sum(model.x[c] for c in C))
model.abv = pyomo.Constraint(expr = 0 == sum(model.x[c]*(data[c]['abv'] - abv) for c in C))
solver = pyomo.SolverFactory('cbc')
solver.solve(model)
print('Optimal Blend')
for c in data.keys():
print(' ', c, ':', model.x[c](), 'gallons')
print()
print('Volume = ', model.vol(), 'gallons')
print('Cost = $', model.cost())
beer_blend(vol, abv, data)
```

In [ ]:

```
```