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)
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:
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:
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).
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")
# 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")
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%
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 |
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 |
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) |