Blending Whisky with Z3

April 29, 2021

A whisky producer uses three kind of licor to make their whisky: \(A\), \(B\) and \(C\).

Three kind of whisky can be made using the licor in the following proportions:

The producer has the following stock and price list from its licor supplier:

Note: the prices are in $ per liter and the stock are in liters.

The goal is to maximize the profit.

Bad plan: non-linear system

Let’s begin with the bounds of each licor:

We introduced a new variable per combination of licor and whisky: \(A_E\) stands for the ratio of licor \(A\) used in the whisky B, \(B_F\) stands for the ratio of licor \(B\) used in the whisky \(F\) and so on.

These variables are in liters of licor per liter of whisky (or they are unit-less if you prefer).

And we are fuck.

These three equations are not longer linear because we have the product of two variable (\(A_E E\) for example).

We need to rethink our strategy.

The blending

Think in \(A_E\) not as the ratio but as the amount of liters itself of licor \(A\) used in \(E\).

This is a blending problem where some variables represent a part, fraction or subcomponent of another.

Now the system is linear again.

What about the restrictions of each recipe?

It would be wonderful to put something like this that it is a literal translation of the problem into inequalities:

But again, that makes the system non-linear.

Instead we do the following:

We complete the recipe of each whisky with the remaining licor:

So what did we do? We split each licor (input \(I\)) into \(N\) blending variables (\(I_o\)), one for each whisky type (output \(O\)).

• We ensured that the sum of the blending variables for the same licor (input) summed up the total amount of that licor (\(A_E + A_F + A_G = A\)) – we did a partition.

• We restricted each blending variable based on a proportion of each amount of whisky (\(0.6 E <= A_E <= 1 E\)) – these are the blending rules.

• We ensured that the sum of the blending variable across all the licor (input) for the same whisky (output) summed up the total amount of that whisky.

  \(A_E\)	\(A_F\)	\(A_G\)	→ \(A\)
  \(B_E\)	\(B_F\)	\(B_G\)	→ \(B\)
  \(C_E\)	\(C_F\)	\(C_G\)	→ \(C\)
  ↓	↓	↓
  \(E\)	\(F\)	\(G\)

The goal

And finally, this is the goal to maximize:

Find the optimum with Z3

Setup the engine, create the variables and ensure that them are non-negative.

>>> from z3 import Reals, Optimize
>>> s = Optimize()

>>> A, B, C, E, F, G = Reals('A B C E F G')
>>> A_E, A_F, A_G = Reals('A_E A_F A_G')
>>> B_E, B_F, B_G = Reals('B_E B_F B_G')
>>> C_E, C_F, C_G = Reals('C_E C_F C_G')

>>> s.add(*[x >= 0 for x in [A, B, C, E, F, G]])
>>> s.add(*[x >= 0 for x in [A_E, A_F, A_G]])
>>> s.add(*[x >= 0 for x in [B_E, B_F, B_G]])
>>> s.add(*[x >= 0 for x in [C_E, C_F, C_G]])

Limit the amount of supplies

>>> s.add(
...     A <= 2000,
...     B <= 2500,
...     C <= 1200
... )

Split the licor for each whisky

>>> s.add(
...     A_E + A_F + A_G == A,
...     B_E + B_F + B_G == B,
...     C_E + C_F + C_G == C
... )

Blending rules:

>>> s.add(
...     # minimum amount of licor for whisky E
...     0.6 * E <= A_E,
...     0 <= C_E,
...     # maximum amount of licor for whisky E
...     A_E <= E,
...     C_E <= 0.2 * E,
...
...     # minimum amount of licor for whisky F
...     0.15 * F <= A_F,
...     0 <= C_F,
...     # maximum amount of licor for whisky F
...     A_F <= F,
...     C_F <= 0.8 * F,
...
...     # the same for whisky G
...     0 <= C_G,
...     C_G <= 0.5 * G
... )

>>> s.add(
...     A_E + B_E + C_E == E,
...     A_F + B_F + C_F == F,
...     A_G + B_G + C_G == G
... )

Profit!

>>> costs = A * 700 + B * 500 + C * 400
>>> income = E * 680 + F * 570 + G * 450
>>> profit = s.maximize(income - costs)

>>> s.check()
sat

>>> profit.value()
3590000/9

>>> m = s.model()
>>> print("Licor:", "A =", m[A], "B =", m[B], "C =", m[C])
Licor: A = 2000 B = 2500 C = 1200

>>> # We use m[.] to retrieve simple variables
>>> print("Whisky:", "E =", m[E], "F =", m[F], "G =", m[G])
Whisky: E = 22900/9 F = 28400/9 G = 0

>>> # We use m.eval() to evaluate complex expressions in the model context
>>> print("Costs =", m.eval(costs), "Income =", m.eval(income), "Profit =", profit.value())
Costs = 3130000 Income = 31760000/9 Profit = 3590000/9

Related tags: z3, smt, sat, solver, linear optimization, blending