The latest version of this IPython notebook is available at http://github.com/jckantor/ESTM60203 for noncommercial use under terms of the Creative Commons Attribution Noncommericial ShareAlike License (CC BY-NC-SA 4.0).

J.C. Kantor ([email protected])

Newsvendor Problem

This IPython notebook demonstrates the formulation and solution of the well-known "Newsvendor Problem" using GLPK/Mathprog.

Initializations

In [3]:
from IPython.core.display import HTML
HTML(open("styles/custom.css", "r").read())
Out[3]:

Background

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.

Solution

In [1]:
%%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.52
Parameter(s) specified in the command line:
 -m /dev/stdin
Reading model section from /dev/stdin...
Reading data section from /dev/stdin...
/dev/stdin:86: warning: final NL missing before end of file
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.52
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  infeas =  0.000e+00 (0)
*     5: obj =   4.000000000e+02  infeas =  0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.1 Mb (142439 bytes)
Writing Table_EVM...
Writing Table_EVPI...
Writing Table_SP...
Model has been successfully processed

Expected Value for the Mean Scenario (EVM)

In [2]:
import pandas
evm = pandas.read_csv("evm.csv")
display(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

3 rows × 7 columns

Expected Value for the Mean Scenario = 362.50

Expected Value with Perfect Information (EVPI)

In [3]:
evpi = pandas.read_csv("evpi.csv")
display(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

3 rows × 7 columns

Expected Value with Perfect Information = 575.00

Expected Value by Stochastic Programming

In [4]:
evsp = pandas.read_csv("evsp.csv")
display(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

3 rows × 7 columns

Expected Value by Stochastic Programming = 400.00

Value of Perfect Information

In [5]:
print "Value of Perfect Information = {:6.2f}".format(ev_evpi-ev_evsp)
Value of Perfect Information = 175.00

Value of the Stochastic Solution

In [6]:
print "Value of the Stochastic Solution = {:6.2f}".format(ev_evsp-ev_evm)
Value of the Stochastic Solution =  37.50
In [6]:
 
In [30]:
r = 1.00
c = 0.60
w = 0.25

def profit(D,x):
    return r*min([D,x]) + w*max([0,x-D]) - c*x
In [31]:
scenarios = [['Low Demand',75,.25],['High Demand',200,.75]]
In [33]:
def exprofit(x):
    v = 0
    for s in scenarios:
        v += s[2]*profit(s[1],x)
    return profit

x = linspace(0,400,400)
exprofit(100)
Out[33]:
<function __main__.profit>
In [23]:
x = linspace(0,400,400)
plot(x,map(exprofit,x))
xlabel('Order size')
ylabel('Expected Profit')
Out[23]:
<matplotlib.text.Text at 0x108b82790>
In [ ]: