Refinery Example

This model is an example of a production problem. In production planning problems, choices must be made regarding what resources to use to produce what products. These problems are common across a broad range of manufacturing and mining industries.

In this example we’ll model and solve a production planning problem. In this case the application is to optimize the output of a refinery. The company can buy an restricted amount of two different kinds of crude oils, which can be processed to gasoline (motor fuel), jet fuel or oil. The aim is to create an optimal production plan that maximizes the total profit, while considering production capacity and other restrictions.

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.


H. Paul Williams, Model Building in Mathematical Programming, fifth edition (Page 258-260, 356-357)


Problem Description

An oil refinery purchases two types of crude oil (Crude 1 and Crude 2) and refines them through a four-step process (distillation, reforming, cracking and blending) to produce gas and fuels that are sold.

Step One: Distillation

The distillation process separates each crude oil into six fractions according to their boiling points. These fractions are: light naphtha, medium naphtha, light oil, heavy oil and residuum. The octane numbers for the light, medium and heavy naphthas are, respectively, 90, 80 and 70.

The fractions into which one barrel of each type of crude separates into are as follows:

  Light Naphtha Medium Naphtha Heavy Naphtha Light Oil Heavy Oil Residuum
Crude 1 0.1 0.2 0.2 0.12 0.2 0.13
Crude 2 0.15 0.25 0.18 0.08 0.19 0.12





There is a small amount of waste in distillation.


Step Two: Reforming

After distillation, the resulting naphthas can be blended into different grades of gas, or they can go through a process called reforming. The output of the reforming process is a product known as reformed gasoline with an octane number of 115.

The yields of reformed gasoline from each barrel of the different naphthas are given as follows:

  • One barrel of light naphtha yields 0.6 barrels of reformed gasoline
  • One barrel of medium naphtha yields 0.52 barrels of reformed gasoline
  • One barrel of heavy naphtha yields 0.45 barrels of reformed gasoline

Step Three: Cracking

Light and heavy oils can be blended into jet fuel or be put though a process known as catalytic cracking. The catalytic cracker produces cracked oil and cracked gasoline.

Cracked gasoline has an octane number of 105 with the following yields:

  • One barrel of light oil yields 0.68 barrels of cracked oil and 0.28 barrels of cracked gasoline
  • One barrel of heavy oil yields 0.75 barrels of cracked oil and 0.2 barrels of cracked gasoline

Cracked oil is used for blending fuel oil and jet fuel; cracked gasoline is used for blending gasoline. Residuum can be used for either producing lube-oil or blending into jet fuel and fuel oil:

  • One barrel of residuum yields 0.5 barrels of lube-oil

Step Four: Blending

Gasoline

There are two kinds of gasoline, regular and premium, made by blending naphtha, reformed gasoline and cracked gasoline. The only requirement is that regular gasoline must have an octane of at least 84 and premium gasoline must have an octane number of at least 94. It is assumed that octane numbers blend linearly by volume.

Jet Fuel

Jet fuel must have a vapor pressure that does not exceed 1.0 kg cm2. The vapor pressures for light, heavy cracked oils and residuum are, respectively, 1.0, 0.6, 1.5 and 0.05 kg cm2. It is assumed that vapor pressures blend linearly by volume.

Fuel Oil

To produce fuel oil, light oil, cracked oil, heavy oil and residuum must be blending in the ratio of 10:4:3:1. The availability and capacity limitations are as follows:

  • The daily availability of crude 1 is 20 000 barrels
  • The daily availability of crude 2 is 30 000 barrels
  • At most 45 000 barrels of crude can be distilled per day
  • At most 10 000 barrels of naphtha can be reformed per day
  • At most 8000 barrels of oil can be cracked per day
  • The daily production of lube oil must be between 500 and 1000 barrels
  • Premium gasoline production must be at least 40% of regular gasoline production

Each final product has the following profit contribution per barrel (pennies per barrel):

  Profit Contribution
Premium Gasoline 700
Regular Gasoline 600
Jet Fuel 400
Fuel Oil 350
Lube Oil 150

How should the operations of the refinery be planned in order to maximize total profit?

Model Formulation

We define continous non-negative variables:

\[
\begin{array}{rl}
CR1 &\text{crude 1}\\
CR2 &\text{crude 2}\\
LN &\text{light naphtha}\\
MN &\text{medium naphtha}\\
HN &\text{heavy naphtha }\\
LO &\text{light oil}\\
HO &\text{heavy oil}\\
R &\text{residuum}\\
LNRG &\text{light naphtha used to produce reformed gasoline }\\
MNRG &\text{medium naphtha used to produce reformed gasoline }\\
HNRG &\text{heavy naphtha used to produce reformed gasoline}\\
RG &\text{reformed gasoline}\\
LOCGO &\text{light oil used to produce cracked oil and cracked gasoline }\\
HOCGO &\text{heavy oil used to produce cracked oil and cracked gasoline }\\
CG &\text{cracked gasoline}\\
CO &\text{cracked oil}\\
LNPMF &\text{light naphtha used to produce premium motor fuel }\\
LNRMF &\text{light naphtha used to produce regular motor fuel }\\
MNPMF &\text{medium naphtha used to produce premium motor fuel }\\
MNRMF &\text{medium naphtha used to produce regular motor fuel }\\
HNPMF &\text{heavy naphtha used to produce premium motor fuel }\\
HNRMF &\text{heavy naphtha used to produce regular motor fuel }\\
RGPMF &\text{reformed gasoline used to produce premium motor fuel }\\
RGRMF &\text{reformed gasoline used to produce regular motor fuel }\\
CGPMF &\text{cracked gasoline used to produce premium motor fuel }\\
CGRMF &\text{cracked gasoline used to produce regular motor fuel}\\
LOJF &\text{light oil used to produce jet fuel }\\
HOJF &\text{heavy oil used to produce jet fuel }\\
RJF &\text{residuum used to produce jet fuel }\\
COJF &\text{cracked oil used to produce jet fuel}\\
RLBO &\text{residuum used to produce lube-oil}\\
PMF &\text{premium motor fuel}\\
RMF &\text{regular motor fuel}\\
JF &\text{jet fuel}\\
FO &\text{fuel oil}\\
LBO &\text{lube-oil}\\
\end{array}
\]

The variables can be formulated as:
\begin{equation} cr1, cr2, ln, mn, hn, lo, ho, r, lnrg, mnrg, hnrg, rg, locgo, hocgo \geq 0
\end{equation} \begin{equation} cg, co, lnpmf, lnrmf, mnpmf, mnrmf, hnpmf, hnrmf, rgpmf, rgrmf \geq 0
\end{equation} \begin{equation} cgpmf, cgrmf, lojf, hojf, rjf, cojf, lofo, hofo, rfo, cofo \geq 0
\end{equation} \begin{equation} rlbo, pmf, rmf, jf, fo, lbo \geq 0 \end{equation}


The limited availability of the crude oils gives simple upper-bounding constraints:
\begin{equation} cr1 \leq 20000 \end{equation} \begin{equation} cr2 \leq 30000 \end{equation}


The distillation capacity constraint is:
\begin{equation} cr1+cr2 \leq 45000 \end{equation}


The reforming capacity constraint is:
\begin{equation} lnrg+mnrg+hnrg \leq 10000 \end{equation}

The cracking capacity constraint is:
\begin{equation} locgo+hocgo \leq 8000 \end{equation}


The stipulation concerning production of lube-oil gives the following lower- and upper-bounding constraints:
\begin{equation} 500 \leq lbo \leq 1000 \end{equation}


The quantity of light naphtha produced depends on the quantities of crude oil used, taking into account the way in which each crude splits under distillation. This gives:
\begin{equation} 0.1 *cr1 +0.15*cr2 = ln \end{equation}


Similar constraints exist for MN, HN, LO, HO and R:
\begin{equation} 0.2 *cr1 +0.25*cr2 = mn \end{equation} \begin{equation} 0.2 *cr1 +0.18*cr2 = hn \end{equation} \begin{equation} 0.12*cr1 +0.08*cr2 = lo \end{equation} \begin{equation} 0.2 *cr1 +0.19*cr2 = ho \end{equation} \begin{equation} 0.13*cr1 +0.12*cr2 = r \end{equation}


The quantity of reformed gasoline produced depends on the quantities of the naphthas used in the reforming process. This gives the constraint:
\begin{equation} 0.6*lnrg+0.52*mnrg+0.45*hnrg = rg \end{equation}


The quantities of cracked oil and cracked gasoline produced depend on the quantities of light and heavy oil used. This gives the constraints:
\begin{equation} 0.68*locgo+0.75*hocgo = co \end{equation} \begin{equation} 0.28*locgo + 0.2*hocgo = cg \end{equation}


The quantities of light naphtha used for reforming and blending are equal to the quantities available. This gives:
\begin{equation} lnrg+lnpmf+lnrmf = ln \end{equation}


Similar constraints exist for MN and HN:
\begin{equation} mnrg+mnpmf+mnrmf = mn \end{equation} \begin{equation} hnrg+hnpmf+hnrmf = hn \end{equation}


The quantities of light oil used for cracking and blending are equal to the quantities available.

For the blending of fuel oil, the proportion of light oil is fixed at 10/18. Therefore separate variables have not been introduced for this proportion as it is determined by the variable LO. This gives:
\begin{equation} locgo+lojf+0.55*fo = lo \end{equation}


Similar constraints exist for HO, CO and R, also involving fixed proportions of FO, and for CG and RG:
\begin{equation} hocgo+hojf+0.17*fo = ho \end{equation} \begin{equation} cojf+0.22*fo = co \end{equation} \begin{equation} rlbo+rjf+0.055*fo = r \end{equation} \begin{equation} cgpmf+cgrmf = cg \end{equation} \begin{equation} rgrmf+rgpmf = rg \end{equation}


The quantity of premium gasoline produced is equal to the total quantity of its ingredients. This gives:
\begin{equation} lnpmf+mnpmf+hnpmf+rgpmf+cgpmf = pmf \end{equation}


Similar constraints exist for RMF and JF:
\begin{equation} lnrmf+mnrmf+hnrmf+rgrmf+cgrmf = rmf \end{equation} \begin{equation} lojf+hojf+cojf+rjf = jf \end{equation}


The quantity of lube-oil produced (and sold) is 0.5 times the quantity of residuum used. This gives:
\begin{equation} 0.5*rlbo = lbo \end{equation}


Premium motor fuel production must be at least 40% of regular motor fuel production, giving:
\begin{equation} pmf \geq 0.4*rmf \end{equation}


It is necessary to stipulate that the octane number of premium gasoline does not drop below 94. This is done by the constraint:
\begin{equation} 90*lnpmf+80*mnpmf+70*hnpmf+115*rgpmf+105*cgpmf \geq 94*pmf \end{equation}


There is a similar constraint for RMF:
\begin{equation} 90*lnrmf+80*mnrmf+70*hnrmf+115*rgrmf+105*cgrmf \geq 84*rmf \end{equation}


For jet fuel we have the constraint imposed by vapor pressure:
\begin{equation} lojf+0.6*hojf+1.5*cojf+0.05*rjf \leq jf \end{equation}


The only variables involving a profit (or cost) are the final products. This gives an objective (in dollars) to be maximized of:
\begin{equation} \max 7*pmf+6*rmf+4*jf+3.5*fo+1.5*lbo \end{equation}


The full model (with the second objective) can be stated as:
\begin{equation} \max 7*pmf+6*rmf+4*jf+3.5*fo+1.5*lbo \end{equation} \begin{equation} cr1 \leq 20000 \end{equation} \begin{equation} cr2 \leq 30000 \end{equation} \begin{equation} cr1+cr2 \leq 45000 \end{equation} \begin{equation} lnrg+mnrg+hnrg \leq 10000 \end{equation} \begin{equation} locgo+hocgo \leq 8000 \end{equation} \begin{equation} 500 \leq lbo \leq 1000 \end{equation} \begin{equation} 0.1 *cr1 +0.15*cr2 = ln \end{equation} \begin{equation} 0.2 *cr1 +0.25*cr2 = mn \end{equation} \begin{equation} 0.2 *cr1 +0.18*cr2 = hn \end{equation} \begin{equation} 0.12*cr1 +0.08*cr2 = lo \end{equation} \begin{equation} 0.2 *cr1 +0.19*cr2 = ho \end{equation} \begin{equation} 0.13*cr1 +0.12*cr2 = r \end{equation} \begin{equation} 0.6*lnrg+0.52*mnrg+0.45*hnrg = rg \end{equation} \begin{equation}
0.68*locgo+0.75*hocgo = co \end{equation} \begin{equation} 0.28*locgo + 0.2*hocgo = cg \end{equation} \begin{equation} lnrg+lnpmf+lnrmf = ln \end{equation} \begin{equation} mnrg+mnpmf+mnrmf = mn \end{equation} \begin{equation} hnrg+hnpmf+hnrmf = hn \end{equation} \begin{equation} locgo+lojf+0.55*fo = lo \end{equation} \begin{equation} hocgo+hojf+0.17*fo = ho \end{equation} \begin{equation} cojf+0.22*fo = co \end{equation} \begin{equation} rlbo+rjf+0.055*fo = r \end{equation} \begin{equation} cgpmf+cgrmf = cg \end{equation} \begin{equation} rgrmf+rgpmf = rg \end{equation} \begin{equation} lnpmf+mnpmf+hnpmf+rgpmf+cgpmf = pmf \end{equation} \begin{equation} lnrmf+mnrmf+hnrmf+rgrmf+cgrmf = rmf \end{equation} \begin{equation} lojf+hojf+cojf+rjf = jf \end{equation} \begin{equation} 0.5*rlbo = lbo \end{equation} \begin{equation} pmf \geq 0.4*rmf \end{equation} \begin{equation} 90*lnpmf+80*mnpmf+70*hnpmf+115*rgpmf+105*cgpmf \geq 94*pmf \end{equation} \begin{equation} 90*lnrmf+80*mnrmf+70*hnrmf+115*rgrmf+105*cgrmf \geq 84*rmf \end{equation} \begin{equation} lojf+0.6*hojf+1.5*cojf+0.05*rjf \leq jf \end{equation} \begin{equation} cr1, cr2, ln, mn, hn, lo, ho, r, lnrg, mnrg, hnrg, rg, locgo, hocgo \geq 0 \end{equation} \begin{equation} cg, co, lnpmf, lnrmf, mnpmf, mnrmf, hnpmf, hnrmf, rgpmf, rgrmf \geq 0 \end{equation} \begin{equation} cgpmf, cgrmf, lojf, hojf, rjf, cojf, lofo, hofo, rfo, cofo \geq 0 \end{equation} \begin{equation} rlbo, pmf, rmf, jf, fo, lbo \geq 0 \end{equation}

Implementation with comments

First, we import the Gurobi Python Modue and initialize the data structures.

from gurobipy import *

# tested with Python 3.5.2 & Gurobi 7.0.1

crude_numbers = range(1,2+1)
petrols = ["Premium_fuel", "Regular_fuel"]
end_product_names = ["Premium_fuel", "Regular_fuel", "Jet_fuel", "Fuel_oil", "Lube_oil"]
distillation_products_names = ["Light_naphtha", "Medium_naphtha", "Heavy_naphtha",
                               "Light_oil", "Heavy_oil", "Residuum"]
naphthas = ["Light_naphtha", "Medium_naphtha", "Heavy_naphtha"]
intermediate_oils = ["Light_oil", "Heavy_oil"]
cracking_products_names = ["Cracked_gasoline", "Cracked_oil"]
used_for_motor_fuel_names = ["Light_naphtha", "Medium_naphtha", "Heavy_naphtha",
                             "Reformed_gasoline", "Cracked_gasoline"]
used_for_jet_fuel_names = ["Light_oil", "Heavy_oil", "Residuum", "Cracked_oil"]

crude_bounds = {1:20000, 2:30000}
lb_lube_oil = 500
ub_lube_oil = 1000

max_crude = 45000
max_reform = 10000
max_cracking = 8000

distillation_splitting_coefficients = {"Light_naphtha": (0.1, 0.15),
                          "Medium_naphtha": (0.2, 0.25),
                         "Heavy_naphtha": (0.2, 0.18),
                         "Light_oil": (0.12, 0.08),
                         "Heavy_oil": (0.2, 0.19),
                         "Residuum": (0.13, 0.12)}

cracking_splitting_coefficients = {("Light_oil","Cracked_oil"): 0.68,
                                   ("Heavy_oil","Cracked_oil"): 0.75,
                                   ("Light_oil","Cracked_gasoline"): 0.28,
                                   ("Heavy_oil","Cracked_gasoline"): 0.2}

reforming_splitting_coefficients = {"Light_naphtha": 0.6, "Medium_naphtha":0.52, "Heavy_naphtha":0.45}
end_product_profit = {"Premium_fuel":7, "Regular_fuel":6, "Jet_fuel":4, "Fuel_oil":3.5, "Lube_oil":1.5}
blending_coefficients = {"Light_oil": 0.55, "Heavy_oil": 0.17, "Cracked_oil": 0.22, "Residuum": 0.055}

lube_oil_factor = 0.5
pmf_rmf_ratio = 0.4

octance_number_coefficients = {
    "Light_naphtha":90,
    "Medium_naphtha":80,
    "Heavy_naphtha":70,
    "Reformed_gasoline":115,
    "Cracked_gasoline":105,
}
octance_number_fuel = {"Premium_fuel": 94,"Regular_fuel": 84}

vapor_pressure_constants = [0.6, 1.5, 0.05]

Next, we create a model and the variables.

model = Model('Refinery_Optimization')

# Variables
crudes = model.addVars(crude_numbers, ub=crude_bounds, name="cr")    
end_products = model.addVars(end_product_names, name="end_prod")
end_products["Lube_oil"].lb= lb_lube_oil
end_products["Lube_oil"].ub= ub_lube_oil
distillation_products = model.addVars(distillation_products_names, name="dist_prod")
reform_usage = model.addVars(naphthas, name="napthas_to_reformed_gasoline")
reformed_gasoline = model.addVar(name="reformed_gasoline")
cracking_usage = model.addVars(intermediate_oils,name="intermediate_oils_to_cracked_gasoline")
cracking_products = model.addVars(cracking_products_names,  name="cracking_prods")
used_for_regular_motor_fuel = model.addVars(used_for_motor_fuel_names, name="motor_fuel_to_regular_motor_fuel")
used_for_premium_motor_fuel = model.addVars(used_for_motor_fuel_names, name="motot_fuel_to_premium_motor_fuel")
used_for_jet_fuel = model.addVars(used_for_jet_fuel_names, name="jet_fuel")
used_for_lube_oil = model.addVar(vtype=GRB.CONTINUOUS,name="residuum_used_for_lube_oil")

Next, we insert the constraints. The distillation capacity constraint is:

# Max Crude
model.addConstr(crudes.sum() <= max_crude, "max_crude")

The reforming capacity constraint is:

# Reforming Capacity
model.addConstr(reform_usage.sum() <= max_reform, "max_reform")

The cracking capacity constraint is:

model.addConstr(cracking_usage.sum() <= max_cracking, "max_cracking")

The quantity of distillation products produced depends on the quantity of crude oil used, taking into account the way in which each crude splits under distillation. This gives:

# Splitting
model.addConstrs((quicksum(distillation_splitting_coefficients[dpn][crude-1]*crudes[crude] 
	for crude in crudes) == distillation_products[dpn]  
	for dpn in distillation_products_names), "splitting_distillation")

The quantity of reformed gasoline produced depends on the quantities of naphthas used in the reforming process. This gives the constraint:

# Reforming
model.addConstr(reform_usage.prod(reforming_splitting_coefficients) == reformed_gasoline,
                    "splitting_reforming")

The quantities of cracked oil and cracked gasoline produced depend on the quantities of light and heavy oil used. This gives the constraints:

# Cracking
model.addConstrs((quicksum(cracking_splitting_coefficients[oil, crack_prod]*cracking_usage[oil]
                           for oil in intermediate_oils) == cracking_products[crack_prod]
                  for crack_prod in cracking_products_names),
                 name="splitting_cracking")

The quantities of naphthas used for reforming and blending are equal to the quantities available. This gives:

# Continuity
model.addConstrs((reform_usage[naphtha] +
                    used_for_regular_motor_fuel[naphtha] +
                    used_for_premium_motor_fuel[naphtha] ==
                    distillation_products[naphtha] for naphtha in naphthas), "continuity")


model.addConstr(used_for_regular_motor_fuel["Cracked_gasoline"] +
                used_for_premium_motor_fuel["Cracked_gasoline"] ==
                cracking_products["Cracked_gasoline"], "continuity_cracked_gasoline")

model.addConstr(used_for_regular_motor_fuel["Reformed_gasoline"] +
                used_for_premium_motor_fuel["Reformed_gasoline"] ==
                reformed_gasoline, "continuity_reformed_gasoline")

model.addConstr(used_for_premium_motor_fuel.sum() == end_products["Premium_fuel"], "continuity_premium_fuel")

model.addConstr(used_for_regular_motor_fuel.sum() == end_products["Regular_fuel"], "continuity_regular_fuel")

model.addConstr(used_for_jet_fuel.sum() == end_products["Jet_fuel"], "continuity_jet_fuel")

For the blending of fuel oil, the proportion of light oil/heavy oil/ cracked oil/ residuum is fixed. Therefore separate variables have not been introduced for this proportion as it is determined by the variables. This gives:

# Blending
model.addConstr(cracking_usage["Light_oil"]+
                used_for_jet_fuel["Light_oil"]+
                blending_coefficients["Light_oil"]*end_products["Fuel_oil"] ==
                distillation_products["Light_oil"], "fixed_proportion_light_oil_for_blending")
model.addConstr(cracking_usage["Heavy_oil"]+
                used_for_jet_fuel["Heavy_oil"]+
                blending_coefficients["Heavy_oil"]*end_products["Fuel_oil"] ==
                distillation_products["Heavy_oil"], "fixed_proportion_heavy_oil_for_blending")
model.addConstr(used_for_jet_fuel["Cracked_oil"]+
                blending_coefficients["Cracked_oil"]*end_products["Fuel_oil"] ==
                cracking_products["Cracked_oil"], "fixed_proportion_cracked_oil_for_blending")
model.addConstr(used_for_lube_oil +
                used_for_jet_fuel["Residuum"]+
                blending_coefficients["Residuum"]*end_products["Fuel_oil"] ==
                distillation_products["Residuum"], "fixed_proportion_residuum_for_blending")
model.addConstr(used_for_regular_motor_fuel["Cracked_gasoline"] +
                used_for_premium_motor_fuel["Cracked_gasoline"] ==
                cracking_products["Cracked_gasoline"], "continuity_cracked_gasoline")
model.addConstr(used_for_regular_motor_fuel["Reformed_gasoline"] +
                used_for_premium_motor_fuel["Reformed_gasoline"] ==
                reformed_gasoline, "continuity_reformed_gasoline")
The quantity of premium motor fuel produced is equal to the total quantity of its ingredients. This gives

model.addConstr(quicksum(used_for_premium_motor_fuel[name] for name in used_for_motor_fuel_names) ==
                end_products["Premium_fuel"],
                "continuity_premium_fuel")

Similar constraints exist for regular fuel and jet fuel.

model.addConstr(quicksum(used_for_regular_motor_fuel[name] for name in used_for_motor_fuel_names) ==
                end_products["Regular_fuel"],
                "continuity_regular_fuel")
model.addConstr(quicksum(used_for_jet_fuel[name] for name in used_for_jet_fuel_names)
                == end_products["Jet_fuel"], "continuity_jet_fuel")

The quantity of lube-oil produced (and sold) is 0.5 times the quantity of residuum used. This gives:

#  lube-oil is 0.5 of residuum used
model.addConstr(lube_oil_factor*used_for_lube_oil == end_products["Lube_oil"],
                "lube-oil_is_0.5_of_residuum_used")

Premium motor fuel production must be at least 40% of regular motor fuel production, giving:

# pmf/rmf must be 40%
model.addConstr(end_products["Premium_fuel"] >= pmf_rmf_ratio*end_products["Regular_fuel"],
                "pmf/rmf_must_be_40%")

It is necessary to stipulate that the octane number of premium motor (regular motor) fuel does not drop below 94 (84). This is done by the constraint:

# octane numbers
model.addConstr(used_for_regular_motor_fuel.prod(octance_number_coefficients) >=
                octance_number_fuel["Regular_fuel"] * end_products["Regular_fuel"],
                "Octane_number_regular_fuel")

model.addConstr(used_for_premium_motor_fuel.prod(octance_number_coefficients) >=
                octance_number_fuel["Premium_fuel"] * end_products["Premium_fuel"],
                "Octane_number_premium_fuel")

For jet fuel we have the constraint imposed by vapor pressure:

# Vapour pressure
model.addConstr(used_for_jet_fuel["Light_oil"] +
                vapor_pressure_constants[0]*used_for_jet_fuel["Heavy_oil"] +
                vapor_pressure_constants[1]*used_for_jet_fuel["Cracked_oil"] +
                vapor_pressure_constants[2]*used_for_jet_fuel["Residuum"] <= end_products["Jet_fuel"],
                "vapour_pressure")

This model has 29 constraints together with simple bounds on three variables.

A note should be made concerning the blending of fuel oil where the ingredients (light, heavy, and cracked oil and residuum) are used in fixed proportions. It might be preferable to think of the production of FO as an activity. It is common in the oil industry to think in terms of activities rather than quantities. In Section 3.4 of the H.P. Williams book, model formulations are discussed where activities represent the extreme modes of operation of a process. In this example we have a special case of a process with one mode of operation. The level of this activity then fixes the ratios of the ingredients, in a case such as this, automatically.

The only variables involving a profit (or cost) are the final products. This gives an objective (in dollars) to be maximized of:

# Profit
model.setObjective(end_products.prod(end_product_profit), GRB.MAXIMIZE)

Next, this starts 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("refinery-output.sol")

Gurobi Output and Analysis

Optimize a model with 29 rows, 36 columns and 106 nonzeros
Coefficient statistics:
  Matrix range     [5e-02, 1e+02]
  Objective range  [2e+00, 7e+00]
  Bounds range     [5e+02, 3e+04]
  RHS range        [8e+03, 4e+04]
Presolve removed 13 rows and 14 columns
Presolve time: 0.00s
Presolved: 16 rows, 22 columns, 72 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0    7.5812816e+05   2.301700e+05   0.000000e+00      0s
      15    2.1136513e+05   0.000000e+00   0.000000e+00      0s

Solved in 15 iterations and 0.00 seconds
Optimal objective  2.113651348e+05

Analyzing the Output

The optimal solution results in a profit of $211 365.

The optimal values of the variables are given below:

Variable Value Variable Value
CR1 15 000 MNPMF 3537
CR2 30 000 MNRMF 6962
LN 6000 HNPMF 0
MN 10 500 HNRMF 2993
HN 8400 RGPMF 1344
LO 4200 RGRMF 1089
HO 8700 CGPMF 1936
R 5550 CGRMF 0
LNRG 0 LOJF 0
MNRG 0 HOJF 4900
HNRG 5407 RJF 4550
Variable Value Variable Value
RG 2433 COJF 5706
LOCGO 4200 RLBO 1000
HOCGO 3800 PMF 6818
CG 1936 RMF 17 044
CO 5706 JF 15 156
LNPMF 0 FO 0
LNRMF 6000 LBO 500