J.C. Kantor (Kantor.1@nd.edu)
This notebook demonstrates the use of linear programming to analyze the case study
Owen Hall and Kenneth Ko. "Landhills Winery: Developing an Optimal Blending Plan," 3 pages, HBS Case W14167-PDF-ENG, published May 21, 2014.
The notebook uses MathProg/GLPK to represent the model and calculate solutions.
%matplotlib inline
from pylab import *
from IPython.core.display import HTML
HTML(open("styles/custom.css", "r").read())
%%script glpsol -m /dev/stdin -o /dev/stdout -y display.txt --out output
set GRAPES;
set WINES;
param gAcid{GRAPES};
param gSugr{GRAPES};
param gAlch{GRAPES};
param gQnty{GRAPES};
param gCost{GRAPES};
param wAcid{WINES};
param wSugr{WINES};
param wAlch{WINES};
param wPric{WINES};
var x{WINES,GRAPES}, integer >= 0;
var wProd{WINES} >= 0;
var gCons{GRAPES} >= 0;
maximize Profit :
sum{w in WINES} wPric[w]*wProd[w] - sum{g in GRAPES} gCost[g]*gCons[g];
# balances
subject to grapes {g in GRAPES} :
sum{w in WINES} x[w,g] == gCons[g];
subject to wines {w in WINES} :
sum{g in GRAPES} x[w,g] == wProd[w];
# bound on grapes
subject to capacity {g in GRAPES} :
gCons[g] <= gQnty[g];
# acidity
subject to acidity {w in WINES} :
sum{g in GRAPES} x[w,g]*gAcid[g] <= wProd[w]*wAcid[w];
# sugar
subject to sugar {w in WINES} :
sum{g in GRAPES} x[w,g]*gSugr[g] <= wProd[w]*wSugr[w];
# Varietal content
subject to v1 {w in {'wSB11','wSL10','wSL11','wCabS'}} :
sum {g in {'gSB11','gSL10','gSL11'}} x[w,g] >= 0.75*wProd[w];
subject to v2 {w in {'wMerl'}} :
sum {g in {'gMerl'}} x[w,g] >= 0.75*wProd[w];
# alchohol content
subject to a1 {w in WINES} :
sum{g in GRAPES} x[w,g]*gAlch[g] >= 10.0*wProd[w];
subject to a2 {w in WINES} :
sum{g in GRAPES} x[w,g]*gAlch[g] <= 15.0*wProd[w];
# vintage
subject to vn1 :
x['wSB11','gSB11'] >= 0.95*wProd['wSB11'];
subject to vn2 :
x['wSL10','gSL10'] >= 0.95*wProd['wSL10'];
subject to vn3 :
x['wSL11','gSL11'] >= 0.95*wProd['wSL11'];
solve;
printf "Profit = %9.2f\n\n", Profit;
for {w in WINES} {
printf "%6s ", w;
printf "%8.1f ", wProd[w];
for {g in GRAPES} printf "%8.1f ", x[w,g];
printf "\n";
}
data;
param : GRAPES : gAcid, gSugr, gAlch, gQnty, gCost :=
gSB11 0.35 0.12 13.5 50000 2.35
gSL10 0.75 0.25 15.3 60000 2.60
gSL11 0.55 0.30 11.5 30000 2.10
gMerl 0.25 0.08 15.7 200000 1.55;
param : WINES : wAcid, wSugr, wPric :=
wSB11 0.7 0.2 9.00
wSL10 0.7 0.2 9.00
wSL11 0.7 0.2 9.00
wCabS 0.7 0.3 5.50
wMerl 0.3 1.0 2.95;
end;
print(open('display.txt').read())
Profit = 809905.50 wSB11 52631.0 50000.0 0.0 0.0 2631.0 wSL10 0.0 0.0 0.0 0.0 0.0 wSL11 0.0 0.0 0.0 0.0 0.0 wCabS 93061.0 0.0 60000.0 9796.0 23265.0 wMerl 121224.0 0.0 0.0 20204.0 101020.0
print output
GLPSOL: GLPK LP/MIP Solver, v4.52 Parameter(s) specified in the command line: -m /dev/stdin -o /dev/stdout -y display.txt Reading model section from /dev/stdin... Reading data section from /dev/stdin... 73 lines were read Generating Profit... Generating grapes... Generating wines... Generating capacity... Generating acidity... Generating sugar... Generating v1... Generating v2... Generating a1... Generating a2... Generating vn1... Generating vn2... Generating vn3... Model has been successfully generated GLPK Integer Optimizer, v4.52 42 rows, 29 columns, 186 non-zeros 20 integer variables, none of which are binary Preprocessing... 37 rows, 25 columns, 169 non-zeros 20 integer variables, none of which are binary Scaling... A: min|aij| = 8.000e-02 max|aij| = 1.570e+01 ratio = 1.962e+02 GM: min|aij| = 4.751e-01 max|aij| = 2.105e+00 ratio = 4.430e+00 EQ: min|aij| = 2.257e-01 max|aij| = 1.000e+00 ratio = 4.430e+00 2N: min|aij| = 1.562e-01 max|aij| = 1.500e+00 ratio = 9.600e+00 Constructing initial basis... Size of triangular part is 37 Solving LP relaxation... GLPK Simplex Optimizer, v4.52 37 rows, 25 columns, 169 non-zeros * 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0) * 14: obj = 8.099113856e+05 infeas = 0.000e+00 (0) OPTIMAL LP SOLUTION FOUND Integer optimization begins... + 14: mip = not found yet <= +inf (1; 0) + 63921: mip = not found yet <= 8.099075263e+05 (15571; 0) + 89096: mip = not found yet <= 8.099075263e+05 (21702; 0) +106858: mip = not found yet <= 8.099075263e+05 (26028; 0) +121867: mip = not found yet <= 8.099075263e+05 (29681; 0) +135030: mip = not found yet <= 8.099075263e+05 (32887; 0) +147134: mip = not found yet <= 8.099075263e+05 (35835; 0) +157857: mip = not found yet <= 8.099075263e+05 (38447; 0) +167910: mip = not found yet <= 8.099075263e+05 (40893; 0) +177522: mip = not found yet <= 8.099075263e+05 (43230; 0) +186554: mip = not found yet <= 8.099075263e+05 (45426; 0) +194800: mip = not found yet <= 8.099075263e+05 (47428; 0) Time used: 60.0 secs. Memory used: 20.6 Mb. +202782: mip = not found yet <= 8.099075263e+05 (49361; 0) +205452: >>>>> 8.099055000e+05 <= 8.099058632e+05 < 0.1% (50008; 1) +205453: mip = 8.099055000e+05 <= tree is empty 0.0% (0; 100017) INTEGER OPTIMAL SOLUTION FOUND Time used: 61.7 secs Memory used: 20.8 Mb (21850447 bytes) Model has been successfully processed Writing MIP solution to `/dev/stdout'... Problem: stdin Rows: 42 Columns: 29 (20 integer, 0 binary) Non-zeros: 186 Status: INTEGER OPTIMAL Objective: Profit = 809905.5 (MAXimum) No. Row name Activity Lower bound Upper bound ------ ------------ ------------- ------------- ------------- 1 Profit 809906 2 grapes[gSB11] 0 -0 = 3 grapes[gSL10] 0 -0 = 4 grapes[gSL11] 0 -0 = 5 grapes[gMerl] 0 -0 = 6 wines[wSB11] 0 -0 = 7 wines[wSL10] 0 -0 = 8 wines[wSL11] 0 -0 = 9 wines[wCabS] 0 -0 = 10 wines[wMerl] 0 -0 = 11 capacity[gSB11] 50000 50000 12 capacity[gSL10] 60000 60000 13 capacity[gSL11] 30000 30000 14 capacity[gMerl] 126916 200000 15 acidity[wSB11] -18683.9 -0 16 acidity[wSL10] 0 -0 17 acidity[wSL11] 0 -0 18 acidity[wCabS] -8938.65 -0 19 acidity[wMerl] 0 -0 20 sugar[wSB11] -4315.72 -0 21 sugar[wSL10] 0 -0 22 sugar[wSL11] 0 -0 23 sugar[wCabS] -8118.3 -0 24 sugar[wMerl] -107081 -0 25 v1[wSB11] 10526.8 -0 26 v1[wSL10] 0 -0 27 v1[wSL11] 0 -0 28 v1[wCabS] 0.25 -0 29 v2[wMerl] 10102 -0 30 a1[wSB11] 189997 -0 31 a1[wSL10] 0 -0 32 a1[wSL11] 0 -0 33 a1[wCabS] 465304 -0 34 a1[wMerl] 606120 -0 35 a2[wSB11] -73158.3 -0 36 a2[wSL10] 0 -0 37 a2[wSL11] 0 -0 38 a2[wCabS] -0.5 -0 39 a2[wMerl] 0 -0 40 vn1 0.55 -0 41 vn2 0 -0 42 vn3 0 -0 No. Column name Activity Lower bound Upper bound ------ ------------ ------------- ------------- ------------- 1 x[wSB11,gSB11] * 50000 0 2 x[wSL10,gSB11] * 0 0 3 x[wSL11,gSB11] * 0 0 4 x[wCabS,gSB11] * 0 0 5 x[wMerl,gSB11] * 0 0 6 x[wSB11,gSL10] * 0 0 7 x[wSL10,gSL10] * 0 0 8 x[wSL11,gSL10] * 0 0 9 x[wCabS,gSL10] * 60000 0 10 x[wMerl,gSL10] * 0 0 11 x[wSB11,gSL11] * 0 0 12 x[wSL10,gSL11] * 0 0 13 x[wSL11,gSL11] * 0 0 14 x[wCabS,gSL11] * 9796 0 15 x[wMerl,gSL11] * 20204 0 16 x[wSB11,gMerl] * 2631 0 17 x[wSL10,gMerl] * 0 0 18 x[wSL11,gMerl] * 0 0 19 x[wCabS,gMerl] * 23265 0 20 x[wMerl,gMerl] * 101020 0 21 wProd[wSB11] 52631 0 22 wProd[wSL10] 0 0 23 wProd[wSL11] 0 0 24 wProd[wCabS] 93061 0 25 wProd[wMerl] 121224 0 26 gCons[gSB11] 50000 0 27 gCons[gSL10] 60000 0 28 gCons[gSL11] 30000 0 29 gCons[gMerl] 126916 0 Integer feasibility conditions: KKT.PE: max.abs.err = 0.00e+00 on row 0 max.rel.err = 0.00e+00 on row 0 High quality KKT.PB: max.abs.err = 0.00e+00 on row 0 max.rel.err = 0.00e+00 on row 0 High quality End of output
%%script glpsol -m /dev/stdin -o /dev/stdout -y display.txt --out output
set VARIETALS := {'SB2011','SLO2010','SLO2011'};
set NONVARIETALS := {'Cab','Merlot'};
set CABERNETs := {'SB2011','SLO2010','SLO2011','Cab'};
param price{WINES};
param acidity{GRAPES};
param sugar{GRAPES};
param alcohol{GRAPES};
param quantity{GRAPES};
param cost{GRAPES};
var x{WINES,GRAPES} >= 0;
var wine{WINES} >= 0;
var grape{GRAPES} >= 0;
maximize PROFIT: sum{w in WINES} wine[w]*price[w] - sum{g in GRAPES} grape[g]*cost[g];
subject to g_cons{g in GRAPES}: sum{w in WINES} x[w,g] = grape[g];
subject to w_prod{w in WINES}: sum{g in GRAPES} x[w,g] = wine[w];
subject to acid{w in WINES}
subject to qty{g in GRAPES}: sum{w in WINES} x[w,g] <= quantity[g];
solve;
data;
param : WINES : price :=
V1 9.00
V2 9.00
V3 9.00
NV 5.50
M 2.95;
param : GRAPES : acidity, sugar, alcohol, quantity, cost :=
G1 0.35 0.12 13.5 50000 2.35
G2 0.75 0.25 15.3 60000 2.60
G3 0.55 0.30 11.5 30000 2.10
G4 0.25 0.08 15.7 200000 1.55;
end;
Looking at the model output
print output
GLPSOL: GLPK LP/MIP Solver, v4.52 Parameter(s) specified in the command line: -m /dev/stdin -o /dev/stdout -y display.txt Reading model section from /dev/stdin... Reading data section from /dev/stdin... 41 lines were read Generating PROFIT... Generating g_cons... Generating w_prod... Generating qty... Model has been successfully generated GLPK Simplex Optimizer, v4.52 14 rows, 29 columns, 78 non-zeros Preprocessing... 4 rows, 20 columns, 20 non-zeros Scaling... A: min|aij| = 1.000e+00 max|aij| = 1.000e+00 ratio = 1.000e+00 Problem data seem to be well scaled Constructing initial basis... Size of triangular part is 4 * 0: obj = 0.000000000e+00 infeas = 0.000e+00 (0) * 4: obj = 2.413500000e+06 infeas = 0.000e+00 (0) OPTIMAL LP SOLUTION FOUND Time used: 0.0 secs Memory used: 0.1 Mb (142786 bytes) Model has been successfully processed Writing basic solution to `/dev/stdout'... Problem: stdin Rows: 14 Columns: 29 Non-zeros: 78 Status: OPTIMAL Objective: PROFIT = 2413500 (MAXimum) No. Row name St Activity Lower bound Upper bound Marginal ------ ------------ -- ------------- ------------- ------------- ------------- 1 PROFIT B 2.4135e+06 2 g_cons[G1] NS 0 -0 = 2.35 3 g_cons[G2] NS 0 -0 = 2.6 4 g_cons[G3] NS 0 -0 = 2.1 5 g_cons[G4] NS 0 -0 = 1.55 6 w_prod[V1] NS 0 -0 = -9 7 w_prod[V2] NS 0 -0 = -9 8 w_prod[V3] NS 0 -0 = -9 9 w_prod[NV] NS 0 -0 = -5.5 10 w_prod[M] NS 0 -0 = -2.95 11 qty[G1] NU 50000 50000 6.65 12 qty[G2] NU 60000 60000 6.4 13 qty[G3] NU 30000 30000 6.9 14 qty[G4] NU 200000 200000 7.45 No. Column name St Activity Lower bound Upper bound Marginal ------ ------------ -- ------------- ------------- ------------- ------------- 1 x[V1,G1] B 50000 0 2 x[V2,G1] NL 0 0 < eps 3 x[V3,G1] NL 0 0 < eps 4 x[NV,G1] NL 0 0 -3.5 5 x[M,G1] NL 0 0 -6.05 6 x[V1,G2] B 60000 0 7 x[V2,G2] NL 0 0 < eps 8 x[V3,G2] NL 0 0 < eps 9 x[NV,G2] NL 0 0 -3.5 10 x[M,G2] NL 0 0 -6.05 11 x[V1,G3] B 30000 0 12 x[V2,G3] NL 0 0 < eps 13 x[V3,G3] NL 0 0 < eps 14 x[NV,G3] NL 0 0 -3.5 15 x[M,G3] NL 0 0 -6.05 16 x[V1,G4] B 200000 0 17 x[V2,G4] NL 0 0 < eps 18 x[V3,G4] NL 0 0 < eps 19 x[NV,G4] NL 0 0 -3.5 20 x[M,G4] NL 0 0 -6.05 21 wine[V1] B 340000 0 22 wine[V2] B 0 0 23 wine[V3] B 0 0 24 wine[NV] B 0 0 25 wine[M] B 0 0 26 grape[G1] B 50000 0 27 grape[G2] B 60000 0 28 grape[G3] B 30000 0 29 grape[G4] B 200000 0 Karush-Kuhn-Tucker optimality conditions: KKT.PE: max.abs.err = 0.00e+00 on row 0 max.rel.err = 0.00e+00 on row 0 High quality KKT.PB: max.abs.err = 0.00e+00 on row 0 max.rel.err = 0.00e+00 on row 0 High quality KKT.DE: max.abs.err = 0.00e+00 on column 0 max.rel.err = 0.00e+00 on column 0 High quality KKT.DB: max.abs.err = 0.00e+00 on row 0 max.rel.err = 0.00e+00 on row 0 High quality End of output