Food Manufacture II Example

Both this model and Food Manufacture I are examples of a blending problem. In blending optimization problems, multiple raw materials are combined in a way the meets the stated constraints for the lowest cost. These problems are common in numerous industries including, for example, the oil industry (blending different types of crude oil at a refinery) and agriculture (manufacturing feed that meets the different nutritional requirements of different animals).

In this particular example, we’ll model and solve a production planning problem where we must create a final product from a number of ingredients, each of which has different costs, restrictions and features. The aim is to create an optimal multi-period production plan that maximizes profit. More details can be found on the Problem Description and Model Formulation tabs below.

Compared to Food Manufacturing I, this example has additional constraints that change the problem type from a linear program (LP) to a mixed integer program (MIP), which makes it harder to solve. You can see these additional constraints at the end of the Problem Description section below.

Note: you can download the model, implemented in Python, here. More information on this type of model can be found in the fifth edition of Model Building in Mathematical Programming, by H. Paul Williams.


Reference: H. Paul Williams, Model Building in Mathematical Programming, fifth edition (Pages 255, 298-300, 349-350)


Problem Description

A manufacturer needs to refine several raw oils and blend them together to produce a given food product that can be sold. The raw oils needed can be divided into two categories:

Category Oil
Vegetable oils: VEG 1
  VEG 2
Non-vegetable oils: OIL 1
  OIL2
  OIL 3

The manufacturer can choose to buy raw oils for the current month and/or buy them on the futures market for delivery in a subsequent month. Prices for immediate deliver and in the futures market are given below in $/ton:

Month VEG 1 VEG 2 OIL 1 OIL 2 OIL 3
January 110 120 130 110 115
February 130 130 110 90 115
March 110 140 130 100 95
April 120 110 120 120 125
May 100 120 150 110 105
June 90 100 140 80 135

There are a number of additional factors that must be taken into account. These include:

  1. The final food product sells for $150 per ton.
  2. Each category of oil (vegetable and non-vegetable) needs to be refined on a different production line.
  3. There is limited refinement capacity such that in any given month a maximum of 200 tons of vegetable oils and 250 tons of non-vegetable oils can be refined.
  4. Also, there is no waste in the refinement process, so the sum of the raw oils refined will equal the amount of refined oils available.
  5. The cost of refining the oils may be ignored.

In addition to the refining limits above, there are limits to the amount of raw oils that can be stored for future use, and there is a cost associated with each ton of oil stored. The limit is 1000 tons of each raw oil and the storage cost is $5 per ton per month. The manufacturer can not store the produced food product or the refined oils.

The final food product must have a hardness between three and six on a given hardness scale. For the purposes of the model, hardness blends linearly and the hardness of each raw oil is:

Oils Hardness
VEG 1 8.8
VEG 2 6.1
OIL 1 2.0
OIL2 4.2
OIL 3 5.0

At the start of January there are 500 tons of each type of raw oil in storage. For the purpose of the model, this should also be the level of raw oils in storage at the end of June.

This version of the Food Manufacturing problem adds the following additional constraints to the first version:

  • The maximum number of oils used in a month is three.
  • If an oil is used during a month, the minimum quantity used must be 20 tons.
  • The use of VEG1 or VEG2 in a given month requires the use of OIL3 in that same month.

Given the above information, what monthly buying and manufacturing decisions should be made in order to maximize profit?

This problem is based on a larger model built for the margarine producer Van den Bergs and Jurgens and discussed in Williams and Redwood (1974).

Model Formulation

Let \(T\) be a set of time periods, where \(t_0 \in T\) is the first period and \(t_e \in T\) the last period. Let \(V\) be the set of vegetable oils and \(N\) be the set of non-vegetable oils. We set \(O\) as the union of \(V\) and \(N\). For each time period \(t \in T\) and each oil \(o \in O\) we introduce continuous non-negative variables \(b_{t,o}\), \(u_{t,o}\), \(s_{t,o}\). \(b_{t,o}\) describes how much oil we buy of the type o in the time period \(t\). \(u_{t,o}\) describes how much oil we use of the type o in the time period \(t\). \(s_{t,o}\) describes how much oil we store of the type o in the time period \(t\). For each time period \(t \in T\) we introduce a continuous non-negative variable \(f_t\), that describes how much food is produced in the period \(t\). Furthermore, for each oil \(o \in O\) we are given the positive number \(h_o\), which indicates the hardness of the oil. This is a property of each oil. For each time period \(t \in T\) and each oil \(o \in O\), we are given the estimated purchase price \(g_{t,o}\) for the oil in that period. We are also given the sale price \(p\) for the final product, which is a mix of different oils. Let i be the initial storage amount. Let \(K=\{(a,b,c) \in O^3 : a \neq b ,b \neq c,a \neq c\}\) be a triple of oils, for which we have a relation, which says if we choose \(a\) or \(b\), then we also have to choose \(c\) in our product mix in the current time period. The variables can be formulated as: \begin{equation} b_{t,o},u_{t,o},s_{t,o} \geq 0 \quad \forall t \in T, \ \forall o \in O \end{equation} \begin{equation} f_t \geq 0 \quad \forall t \in T \end{equation} \begin{equation} d_{t,o} \in \{0,1\} \quad \forall t \in T, \ \forall o \in O. \end{equation} The balance constraint states that the amount of this kind of oil in the last period and that what is bought equals the amount of what is used and stored. This can be formulated as: \begin{equation} s_{t-1,o} + b_{t,o} = u_{t,o}+ s_{t,o} \quad \forall t \in T \setminus t_0, \ \forall o \in O \end{equation} \begin{equation} i + b_{t_0,o} = u_{t_0,o} + s_{t_0,o} \quad \forall o \in O. \end{equation} We set \(s_{t-1,o}=i\) for the first time period \(t_0 \in T\), since there is no period before the first one. The endstore constraints force that at the end of the last period the storage contains the initial amount of each product. The problem description demands that the storage is as full as it was in the beginning. This can be formulated as: \begin{equation} s_{t_e,o} = i \quad \forall o \in O. \end{equation} The capacity constraints restrict the amount of veg and non-veg oils which can be processed per period. Per month, only 200 tons of vegetable oil and 250 tons of non-vegetable oil can be processed to the final product. This can be formulated as: \begin{equation} \sum_{o \in V} u_{t,o} \leq 200 \quad \forall t \in T \end{equation} \begin{equation} \sum_{o \in N} u_{t,o} \leq 250 \quad \forall t \in T. \end{equation} The hardness constraints limit the hardness of the final product, which needs to remain between 3 and 6. Each oil has a certain hardness. The final product may be made up of different oils. The hardness of the final product is measured by the hardness of each ingredient multiplied by its share of the final product. It is assumed that the hardness blends linearily. This can be stated as: \begin{equation} 3 f_t \leq \sum_{o \in V} h_o u_{t,o} \leq 6f_t \quad \forall t \in T. \end{equation} The conserve constraints ensure that the amount of products used in each period equals the food produced in that period. This makes sure that all oil that is used is also processed to the final product (food). This can be formulated as: \begin{equation} \sum_{o \in V} u_{t,o} = f_t \quad \forall t \in T. \end{equation} The Condition 1 constraints und Condition 2 constraints force that the binary variable for each product and each period is set to one if and only if the continuous variable used for the same product and the same period is non zero. The binary variable is called an indicator variable since it is linked to a continuous variable and indicates if it is non-zero. Condition 1 constraints ensure also that each final product is only made up by at most three ingredients. Condition 2 constraint also force that if any product is used in any period then at least 20 tons is used. This can be formulated as: \begin{equation} u_{t,o}-20d_{t,o} \geq 0 \quad \forall t \in T,\ \forall o \in O \end{equation} \begin{equation} u_{t,o}-200d_{t,o} \leq 0 \quad \forall t \in T, \ o \in V \end{equation} \begin{equation} u_{t,o}-250d_{t,o} \leq 0 \quad \forall t \in T, \forall o \in N. \end{equation} Condition 3 constraints ensure that if either \(o_1 \in O \) or \(o_2 \in O\) is used in any period then also \(o_3 \in O\) has to be used. This can be stated as: \begin{equation} d_{t,o_1}+ d_{t,o_2}-2 d_{t,o_3} \quad \forall t \in T, \ \forall (o_1,o_2,o_3 ) \in K. \end{equation} The objective is to maximize the profit of the company. This is calculated at therevenue minus costs for buying and storing purchased products (ingredients). This can be stated as: \begin{equation} \max \sum_{t \in T} (p f_t) - \sum_{t \in T} \sum_{o \in O} (g_{t,o} u_{t,o}+5 s_{t,o}). \end{equation} The full model can be stated as: \begin{equation} \max \sum_{t \in T} (p f_t) - \sum_{t \in T} \sum_{o \in O} (g_{t,o} u_{t,o}+5 s_{t,o}) \end{equation} \begin{equation} s_{t-1,o} + b_{t,o} = u_{t,o}+ s_{t,o} \quad \forall t \in T \setminus t_0, \ \forall o \in O \end{equation} \begin{equation} i + b_{t_0,o} = u_{t_0,o} + s_{t_0,o} \quad \forall o \in O \end{equation} \begin{equation} s_{t_e,o} = i \quad \forall o \in O \end{equation} \begin{equation} \sum_{o \in V} u_{t,o} \leq 200 \quad \forall t \in T \end{equation} \begin{equation} \sum_{o \in N} u_{t,o} \leq 250 \quad \forall t \in T \end{equation} \begin{equation} 3 f_t \leq \sum_{o \in V} h_o u_{t,o} \leq 6f_t \quad \forall t \in T \end{equation} \begin{equation} \sum_{o \in V} u_{t,o} = f_t \quad \forall t \in T \end{equation} \begin{equation} u_{t,o}-20 d_{t,o} \geq 0 \quad \forall t \in T,\ \forall o \in O \end{equation} \begin{equation} u_{t,o}-200 d_{t,o} \leq 0 \quad \forall t \in T, \ o \in V \end{equation} \begin{equation} u_{t,o}-250 d_{t,o} \leq 0 \quad \forall t \in T, \forall o \in N \end{equation} \begin{equation} d_{t,o_1}+ d_{t,o_2}-2 d_{t,o_3} \quad \forall t \in T, \ \forall (o_1,o_2,o_3) \in K \end{equation} \begin{equation} b_{t,o},u_{t,o},s_{t,o} \geq 0 \quad \forall t \in T, \ \forall o \in O \end{equation} \begin{equation} f_t \geq 0 \quad \forall t \in T \end{equation} \begin{equation} d_{t,o} \in \{0,1\} \quad \forall t \in T, \ \forall o \in O. \end{equation}

Implementation with comments

First, we import the Gurobi Python Module and initialize the data structures with the given data.

from gurobipy import *

# tested with Python 3.5.2 & Gurobi 7.0.1

time_periods = ["January", "February", "March", "April", "May", "June"]

oils = ["VEG1", "VEG2", "OIL1", "OIL2", "OIL3"]

prices = tupledict({
	('January', 'VEG1'): 110,
	('January', 'VEG2'): 120,
	('January', 'OIL1'): 130,
	('January', 'OIL2'): 110,
	('January', 'OIL3'): 115,
	('February', 'VEG1'): 130,
	('February', 'VEG2'): 130,
	('February', 'OIL1'): 110,
	('February', 'OIL2'): 90,
	('February', 'OIL3'): 115,
	('March', 'VEG1'): 110,
	('March', 'VEG2'): 140,
	('March', 'OIL1'): 130,
	('March', 'OIL2'): 100,
	('March', 'OIL3'): 95,
	('April', 'VEG1'): 120,
	('April', 'VEG2'): 110,
	('April', 'OIL1'): 120,
	('April', 'OIL2'): 120,
	('April', 'OIL3'): 125,
	('May', 'VEG1'): 100,
	('May', 'VEG2'): 120,
	('May', 'OIL1'): 150,
	('May', 'OIL2'): 110,
	('May', 'OIL3'): 105,
	('June', 'VEG1'): 90,
	('June', 'VEG2'): 100,
	('June', 'OIL1'): 140,
	('June', 'OIL2'): 80,
	('June', 'OIL3'): 135
})


hardness = {"VEG1": 8.8, "VEG2": 6.1, "OIL1": 2.0, "OIL2": 4.2, "OIL3": 5.0}

price = 150
IStore = 500
vegCapa = 200
oilCapa = 250

hardness_lb = 3
hardness_ub = 6
store_pricing = 5
min_oil_used = 20

Next, we create a model and the variables. We set the UpdateMode parameter to 1 (which simplifies the code - see the documentation for more details). For each period we create a variable which will take the value of the food produced. For each product (five kinds of oils) and each period we will create variables for the amount, which get purchased, used and stored.

For each period and each product we need a binary variable, which indicates if this product is used in the current period.

model = Model('Food Manufacture II')
# Quantity of food produced in each period
food = model.addVars(time_periods, name="Food")
# Quantity bought of each product in each period
buy = model.addVars(time_periods, oils, name = "Buy")
# Quantity used of each product  in each period
use = model.addVars(time_periods, oils, name = "Use")
# Quantity stored of each product  in each period
store = model.addVars(time_periods, oils, name = "Store")
# binary variables =1, if use>0
d = model.addVars(time_periods, oils, vtype=GRB.BINARY, name = "d")

Next, we insert the constraints. The balance constraints ensure that the amount of oil that is in the storage in the last period and the amount that gets purchased equals the amount that is used and stored for each oil in the current period. This makes sure that all oil in the model was either in the initial storage or is purchased in some month.

#Initial Balance
model.addConstrs((IStore + buy[time_periods[0], oil] ==
     use[time_periods[0], oil] + store[time_periods[0], oil] for oil in oils), "Initial_Balance")

# Balance
model.addConstrs(
    (store[time_periods[time_periods.index(time_period)-1], oil] + buy[time_period, oil] ==
     	use[time_period, oil] + store[time_period, oil] for oil in oils for time_period in time_periods 
     	if time_period != time_periods[0]), "Balance")

The end balance constraints force that at the end of the last period the storage contains the initial amount of each oil. The problem description demands that the storage is as full as in the beginning.

#End Balance
model.addConstrs((store[time_periods[-1], oil] == IStore for oil in oils),
                 "End_Balance")

The capacity constraints restrict the amount of veg and non-veg oils which can be processed per period. Per month only 200 tons of vegetable oil and 250 tons of non-vegetable oil can be processed to the final product.

# Capacity1 & Capacity2
model.addConstrs(
    (quicksum(use[time_period, oil] for oil in oils if "VEG" in oil) <= vegCapa
     	for time_period in time_periods), "Capacity_Veg")
model.addConstrs(
    (quicksum(use[time_period, oil] for oil in oils if "OIL" in oil) <= oilCapa
     	for time_period in time_periods), "Capacity_Oil")

The hardness constraints limit the hardness of the final product, which needs to remain between 3 and 6. Each oil has a certain hardness. The final product may be made up of different oils. The hardness of the final product is measured by the hardness of each ingredient multiplied by its share of the final product. It is assumed that the hardness blends linearily.

# Hardness
model.addConstrs(
    (quicksum(hardness[oil] * use[time_period, oil] for oil in oils)
     	>= hardness_lb * food[time_period] for time_period in time_periods), "Hardness_lower")
model.addConstrs(
    (quicksum(hardness[oil] * use[time_period, oil] for oil in oils)
     	<= hardness_ub * food[time_period] for time_period in time_periods), "Hardness_upper")

The conserve constraints ensure that the amount of products used in each period equals the amount of food produced in that period. This makes sure that all oil that is used is also processed to the final product (food).

# Conserve
model.addConstrs((use.sum(time_period) == food[time_period]
                  for time_period in time_periods), "Conserve")

Condition 1 constraints ensure that each final product is only made up of at most three ingredients.

# Cond 1 (force condition: d=1 -> use>0)
model.addConstrs((d.sum(time_period) <= 3 for time_period in time_periods),"Cond1")

Condition 2 constraints force that if any product is used in any period then at least 20 tons is used.

Condition 2 constraints also force that the binary variable for each product and each period is set to one if and only if the continuous variable used for the same product and the same period is non-zero. The binary variable is called an indicator variable since it is linked to a continuous variable and indicates if it is non-zero.

# Cond 2
##### Pure MIP formulation
model.addConstrs(
	(use[time_period, oil] - min_oil_used * d[time_period, oil] >= 0
	 	for time_period in time_periods for oil in oils), "Cond2")
model.addConstrs(
	(use[time_period, oil] - vegCapa * d[time_period, oil] <= 0
	 	for time_period in time_periods
	 	for oil in oils if "VEG" in oil), "Cond2a")
model.addConstrs(
	(use[time_period, oil] - oilCapa * d[time_period, oil] <= 0
	 	for time_period in time_periods
	 	for oil in oils if "OIL" in oil), "Cond2a")

Condition 2 can also be modeled using Gurobi's general constraints (from version 7.0 onwards). The alternative formulation of condition 2 is given below

# Alternative formulation - Using Gurobi's General Constraints
for time_period in time_periods:
    for oil in oils:
        model. addGenConstrIndicator(d[time_period, oil], 0,  use[time_period, oil] == 0,
                        name = "Cond2a_{}_{}".format(oil, time_period))
        model. addGenConstrIndicator(d[time_period, oil], 1,  use[time_period, oil] >=  min_oil_used,
                        name = "Cond2b_{}_{}".format(oil, time_period))

Condition 3 constraints ensure that if vegetable one or vegetable two are used, then oil three must also be used.

# Cond 3 If Veg1 or Veg2 are used then Oil3 must also be used in this month
##### Pure MIP formulation
model.addConstrs((d[time_period, "VEG1"] <= d[time_period, "OIL3"] for time_period in time_periods),"Cond3a")
model.addConstrs((d[time_period, "VEG2"] <= d[time_period, "OIL3"] for time_period in time_periods), "Cond3b")

The alternative formulation of Condition 3 using Gurobi General Constraints
# Alternative formulation - Using Gurobi's General Constraints
for time_period in time_periods:
    model. addGenConstrIndicator(d[time_period, "VEG1"], 1,  d[time_period, "OIL3"] == 1, name = "Cond3a_{}".format(time_period))
    model. addGenConstrIndicator(d[time_period, "VEG2"], 1,  d[time_period, "OIL3"] == 1, name = "Cond3b_{}".format(time_period))

The objective is to maximize the profit of the company. This is calculated as revenue minus costs for buying and storing of the purchased products (ingredients).

# Objective
obj = price*food.sum() - buy.prod(prices) - store_pricing*store.sum()
model.setObjective(obj, GRB.MAXIMIZE) # maximize profit

Next, we start the optimization and Gurobi tries to find the optimal solution.

model.optimize()
for v in model.getVars():
    if v.X != 0:
        print("%s %f" % (v.Varname, v.X))

Note: If you want to write your solution to a file, rather than print it to the terminal, you can use the model.write() command. An example implementation is:

  model.write("food-manufacture-ii-output.sol")

Gurobi Output and Analysis

 Optimize a model with 71 rows, 126 columns and 288 nonzeros
Model has 72 general constraints
Variable types: 96 continuous, 30 integer (30 binary)
Coefficient statistics:
  Matrix range     [1e+00, 9e+00]
  Objective range  [5e+00, 2e+02]
  Bounds range     [1e+00, 1e+00]
  RHS range        [3e+00, 5e+02]
Presolve removed 60 rows and 85 columns
Presolve time: 0.01s
Presolved: 173 rows, 173 columns, 442 nonzeros
Presolved model has 42 SOS constraint(s)
Variable types: 71 continuous, 102 integer (90 binary)
Found heuristic solution: objective -75000.00000

Root relaxation: objective 1.078426e+05, 140 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 107842.593    0    6 -75000.000 107842.593   244%     -    0s
H    0     0                    68805.555556 107842.593  56.7%     -    0s
H    0     0                    78205.555556 107842.593  37.9%     -    0s
     0     0 105650.000    0    5 78205.5556 105650.000  35.1%     -    0s
     0     2 105650.000    0    5 78205.5556 105650.000  35.1%     -    0s
H   16    15                    88828.703704 105650.000  18.9%   7.0    0s
H   16    15                    89428.703704 105650.000  18.1%   7.0    0s
H   16    15                    99078.703704 105650.000  6.63%   7.0    0s
H   32    27                    100278.70370 105650.000  5.36%   6.7    0s

Cutting planes:
  Implied bound: 5
  Inf proof: 2

Explored 841 nodes (4417 simplex iterations) in 0.10 seconds
Thread count was 4 (of 4 available processors)

Solution count 7: 100279 99078.7 89428.7 ... -75000
Pool objective bound 100279

Optimal solution found (tolerance 1.00e-04)
Best objective 1.002787037037e+05, best bound 1.002787037037e+05, gap 0.0000%

Analyzing the Output

When originally designed, this model proved comparatively hard to solve (see Food manufacture I). The profit (revenue from sales − cost of raw oils) resulting from this plan is $100 279. There are alternative, equally good, solutions.  

Buy Use Store
January Nothing 85.2 tons VEG 1
114.8 tons VEG 2
250 tons OIL 3
414.8 tons VEG 1
385.2 tons VEG 2
500 tons OIL 1
500 tons OIL 2
250 tons OIL 3
February 190 tons OIL 2

85.2 tons VEG 1
114.8 tons VEG 2
250 tons OIL 3

329.6 tons VEG 1
270.4 tons VEG 2
500 tons OIL 1
690 tons OIL 2
0 tons OIL 3
March 540 tons OIL 3 200 tons VEG 2
20 tons OIL 3
329.6 tons VEG 1
70.4 tons VEG 2
500 tons OIL 1
460 tons OIL 2
520 tons OIL 3
April Nothing 155 tons VEG 1
230 tons OIL 2
20 tons OIL 3
174.6 tons VEG 1
70.4 tons VEG 2
500 tons OIL 1
230 tons OIL 2
500 tons OIL 3
May 40 tons OIL 3 155 tons VEG 1
230 tons OIL 2
20 tons OIL 3

19.6 tons VEG 1
70.4 tons VEG 2
500 tons OIL 1
0 tons OIL 2
520 tons OIL 3

June

480.4 tons VEG 1
629.6 tons VEG 2
730 tons OIL 2
200 tons VEG 2
230 tons OIL 2
20 tons OIL 3
500 tons each oil (stipulated)