J.C. Kantor (Kantor.1@nd.edu)
This Jupyter notebook demonstrates the formulation and solution of the well-known "Newsvendor Problem" using GLPK/Mathprog and Gurobi Python.
The newsvendor problem is a two stage decision problem with recourse. The vendor needs to decide how much inventory to order today to fulfill an uncertain demand. The data includes the unit cost, price, and salvage value of the product being sold, and a probabilistic forecast of demand. The objective is to maximize expected profit.
As shown in lecture, this problem can be solved with a plot, and the solution interpreted in terms of a cumulative probability distribution. The advantage of a MathProg model is that additional constraints or other criteria may be considered, such as risk aversion.
There is an extensive literature on the newsvendor problem which has been studied since at least 1888. See here for a thorough discussion.
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from gurobipy import *
%%script glpsol -m /dev/stdin
# Example: Newsvendor.mod
/* Unit Price Data */
param r >= 0; # Price
param c >= 0; # Cost
param w >= 0; # Salvage value
/* Price data makes sense only if Price > Cost > Salvage */
check: c <= r;
check: w <= c;
/* Probabilistic Demand Forecast */
set SCENS; # Scenarios
param D{SCENS} >= 0; # Demand
param Pr{SCENS} >= 0; # Probability
/* Probabilities must sum to one. */
check: sum{k in SCENS} Pr[k] = 1;
/* Expected Demand */
param ExD := sum{k in SCENS} Pr[k]*D[k];
/* Lower Bound on Profit: Expected Value of the Mean Solution */
param EVM := -c*ExD + sum{k in SCENS} Pr[k]*(r*min(ExD,D[k])+w*max(ExD-D[k],0));
/* Upper Bound on Profit: Expected Value with Perfect Information */
param EVPI := sum{k in SCENS} Pr[k]*(r-c)*D[k];
/* Two Stage Stochastic Programming */
var x >= 0; # Stage 1 (Here and Now): Order Quqntity
var y{SCENS}>= 0; # Stage 2 (Scenario Dep): Actual Sales
var ExProfit; # Expected Profit
/* Maximize Expected Profit */
maximize OBJ: ExProfit;
/* Goods sold are limited by the order quantities and the demand */
s.t. PRFT: ExProfit = -c*x + sum{k in SCENS} Pr[k]*(r*y[k] + w*(x-y[k]));
s.t. SUPL {k in SCENS}: y[k] <= x;
s.t. DMND {k in SCENS}: y[k] <= D[k];
solve;
table Table_EVM {k in SCENS} OUT "CSV" "evm.csv" "Table":
k~Scenario,
Pr[k]~Probability,
D[k]~Demand,
ExD~Order,
min(ExD,D[k])~Sold,
max(ExD-D[k],0)~Salvage,
-c*ExD + r*min(ExD,D[k]) + w*max(ExD-D[k],0)~Profit;
table Table_EVPI {k in SCENS} OUT "CSV" "evpi.csv" "Table":
k~Scenario,
Pr[k]~Probability,
D[k]~Demand,
D[k]~Order,
D[k]~Sold,
0~Salvage,
-c*D[k] + r*D[k]~Profit;
table Table_SP {k in SCENS} OUT "CSV" "evsp.csv" "Table":
k~Scenario,
Pr[k]~Probability,
D[k]~Demand,
x~Order,
y[k]~Sold,
x-y[k]~Salvage,
-c*x + r*y[k] + w*(x-y[k])~Profit;
data;
/* Problem Data corresponds to a hypothetical case of selling programs prior
to a home football game. */
param r := 10.00; # Unit Price
param c := 6.00; # Unit Cost
param w := 2.00; # Unit Salvage Value
param: SCENS: Pr D :=
HiDmd 0.25 250
MiDmd 0.50 125
LoDmd 0.25 75 ;
end;
GLPSOL: GLPK LP/MIP Solver, v4.60 Parameter(s) specified in the command line: -m /dev/stdin Reading model section from /dev/stdin... Reading data section from /dev/stdin... 86 lines were read Checking (line 10)... Checking (line 11)... Checking (line 19)... Generating OBJ... Generating PRFT... Generating SUPL... Generating DMND... Model has been successfully generated GLPK Simplex Optimizer, v4.60 8 rows, 5 columns, 15 non-zeros Preprocessing... 3 rows, 4 columns, 6 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 3 * 0: obj = -0.000000000e+00 inf = 0.000e+00 (3) * 5: obj = 4.000000000e+02 inf = 0.000e+00 (0) OPTIMAL LP SOLUTION FOUND Time used: 0.0 secs Memory used: 0.1 Mb (142447 bytes) Writing Table_EVM... Writing Table_EVPI... Writing Table_SP... Model has been successfully processed
evm = pd.read_csv("evm.csv")
print(evm)
ev_evm = sum(evm['Probability']*evm['Profit'])
print("Expected Value for the Mean Scenario = {:6.2f}".format(ev_evm))
Scenario Probability Demand Order Sold Salvage Profit 0 HiDmd 0.25 250 143.75 143.75 0.00 575 1 MiDmd 0.50 125 143.75 125.00 18.75 425 2 LoDmd 0.25 75 143.75 75.00 68.75 25 Expected Value for the Mean Scenario = 362.50
evpi = pd.read_csv("evpi.csv")
print(evpi)
ev_evpi = sum(evpi['Probability']*evpi['Profit'])
print("Expected Value with Perfect Information = {:6.2f}".format(ev_evpi))
Scenario Probability Demand Order Sold Salvage Profit 0 HiDmd 0.25 250 250 250 0 1000 1 MiDmd 0.50 125 125 125 0 500 2 LoDmd 0.25 75 75 75 0 300 Expected Value with Perfect Information = 575.00
evsp = pd.read_csv("evsp.csv")
print(evsp)
ev_evsp = sum(evsp['Probability']*evsp['Profit'])
print("Expected Value by Stochastic Programming = {:6.2f}".format(ev_evsp))
Scenario Probability Demand Order Sold Salvage Profit 0 HiDmd 0.25 250 125 125 0 500 1 MiDmd 0.50 125 125 125 0 500 2 LoDmd 0.25 75 125 75 50 100 Expected Value by Stochastic Programming = 400.00
print("Value of Perfect Information = {:6.2f}".format(ev_evpi-ev_evsp))
Value of Perfect Information = 175.00
print("Value of the Stochastic Solution = {:6.2f}".format(ev_evsp-ev_evm))
Value of the Stochastic Solution = 37.50
r = 1.00
c = 0.60
w = 0.25
scenarios = [['Low Demand',75,.25],['High Demand',200,.75]]
def profit(D,x):
return r*min([D,x]) + w*max([0,x-D]) - c*x
def exprofit(x):
v = 0
for s in scenarios:
v += s[2]*profit(s[1],x)
return profit(04)
x = np.linspace(0,400,400)
plt.plot(x,[exprofit(xx) for xx in x])
plt.xlabel('Order size')
plt.ylabel('Expected Profit')
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-37-b0cfaf956ad0> in <module>() 17 18 x = np.linspace(0,400,400) ---> 19 plt.plot(x,[exprofit(xx) for xx in x]) 20 plt.xlabel('Order size') 21 plt.ylabel('Expected Profit') /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/pyplot.py in plot(*args, **kwargs) 3316 mplDeprecation) 3317 try: -> 3318 ret = ax.plot(*args, **kwargs) 3319 finally: 3320 ax._hold = washold /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/__init__.py in inner(ax, *args, **kwargs) 1890 warnings.warn(msg % (label_namer, func.__name__), 1891 RuntimeWarning, stacklevel=2) -> 1892 return func(ax, *args, **kwargs) 1893 pre_doc = inner.__doc__ 1894 if pre_doc is None: /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/axes/_axes.py in plot(self, *args, **kwargs) 1405 1406 for line in self._get_lines(*args, **kwargs): -> 1407 self.add_line(line) 1408 lines.append(line) 1409 /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/axes/_base.py in add_line(self, line) 1785 line.set_clip_path(self.patch) 1786 -> 1787 self._update_line_limits(line) 1788 if not line.get_label(): 1789 line.set_label('_line%d' % len(self.lines)) /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/axes/_base.py in _update_line_limits(self, line) 1807 Figures out the data limit of the given line, updating self.dataLim. 1808 """ -> 1809 path = line.get_path() 1810 if path.vertices.size == 0: 1811 return /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/lines.py in get_path(self) 987 """ 988 if self._invalidy or self._invalidx: --> 989 self.recache() 990 return self._path 991 /Users/jeff/anaconda/lib/python3.5/site-packages/matplotlib/lines.py in recache(self, always) 683 y = ma.asarray(yconv, np.float_).filled(np.nan) 684 else: --> 685 y = np.asarray(yconv, np.float_) 686 y = y.ravel() 687 else: /Users/jeff/anaconda/lib/python3.5/site-packages/numpy/core/numeric.py in asarray(a, dtype, order) 480 481 """ --> 482 return array(a, dtype, copy=False, order=order) 483 484 def asanyarray(a, dtype=None, order=None): TypeError: float() argument must be a string or a number, not 'function'
[exprofit(x) for x in x]
[<function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>, <function __main__.profit>]
map(exprofit,x)
<map at 0x113bbd5f8>
from gurobipy import *
import pandas as pd
# problem parameters
# Problem Data corresponds to a hypothetical case of selling programs prior
# to a home football game
r = 10.00 # Unit Price
c = 6.00 # Unit Cost
w = 2.00 # Unit Salvage Value
scenarios = pd.DataFrame([
['Hi Demand', 0.25, 250],
['Med Demand', 0.50, 125],
['Lo Demand', 0.25, 75]
])
scenarios.columns = ['Scenario','Probability','Demand']
scenarios
Scenario | Probability | Demand | |
---|---|---|---|
0 | Hi Demand | 0.25 | 250 |
1 | Med Demand | 0.50 | 125 |
2 | Lo Demand | 0.25 | 75 |
m = Model('Newsvendor Problem')
x = m.addVar()
y = m.addVars(len(scenarios.index))
m.update()
p = sum([scenarios['Probability'][s]*(r*y[s] + w*(x-y[s]) - c*x) for s in scenarios.index])
for s in scenarios.index:
m.addConstr(y[s] <= scenarios['Demand'][s])
m.addConstr(y[s] <= x)
m.setObjective(p,GRB.MAXIMIZE)
m.optimize()
x.X
Optimize a model with 6 rows, 4 columns and 9 nonzeros Coefficient statistics: Matrix range [1e+00, 1e+00] Objective range [2e+00, 4e+00] Bounds range [0e+00, 0e+00] RHS range [8e+01, 2e+02] Presolve removed 3 rows and 0 columns Presolve time: 0.01s Presolved: 3 rows, 4 columns, 6 nonzeros Iteration Objective Primal Inf. Dual Inf. Time 0 1.1500000e+03 4.500000e+02 0.000000e+00 0s 2 4.0000000e+02 0.000000e+00 0.000000e+00 0s Solved in 2 iterations and 0.01 seconds Optimal objective 4.000000000e+02
125.0
p
<gurobi.LinExpr: 2.5 C1 + 0.5 C0 + -0.5 C1 + -1.5 C0 + 5.0 C2 + C0 + -1.0 C2 + -3.0 C0 + 2.5 C3 + 0.5 C0 + -0.5 C3 + -1.5 C0>