###### 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

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
/dev/stdin:86: warning: final NL missing before end of file
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
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 [ ]: