A Price-Setting Newsvendor Problem (Zhan & Shen 2005)

Assume ε follows a normal distribution with mean 0 and standard deviation 20, write a simulation program, which simulates S=5000 demand points and returns the expected profit given your decisions on Q and p. Here Q and p are continuous real numbers (正實數).


import packages
library(pracma)
library(genalg)
set variables
sim_times <- 5000
v <- 0.5
s <- 1
c <- 1
a <- 200
b <- 35
e_norm <- rnorm(n=sim_times, mean=0, sd=20)
e_exp <- rexp(n=sim_times, rate = 1/10)
define functions
cal_demand <- function(a, b, p, e) {
  res <- a - b*p + e
  return(res)
}
cal_profit <- function(a, b, p, e, Q, D, c, v, s) {
  D <- cal_demand(a=a, b=b, p=p, e=e)
  res <- p*min(Q, D) - c*Q + v*max((Q - D), 0) - s*max((D - Q), 0)
  return(res)
}
sim_func_HJ <- function(pQ, e) {
  p <- pQ[1]
  Q <- pQ[2]
  sim_res <- c()
  for (sim_time in 1:sim_times) {
    sim_res[sim_time] <- cal_profit(a, b, p, e[sim_time], Q, D, c, v, s)
  }
  res <- mean(sim_res)
  return(res)
}


  1. Use the Hooke-Jeeves algorithm to find the best Q, p, and expected profit:
# Hooke-Jeeves algorithm
optres.HJ <- fminsearch(sim_func_HJ,   # object function
                   c(5, 180),                # initial values
                   e=e_norm,
                   lower=c(1, 10),          # lower bound
                   upper=c(5.5, 200),         # upper bound
                   method=c("Hooke-Jeeves"), # function method
                   minimize=FALSE)           # maximize objective function

optres.HJ
## $xmin
## [1]   3.343994 105.660781
## 
## $fmin
## [1] 179.1865
## 
## $count
## [1] 474
## 
## $convergence
## [1] 0
## 
## $info
## $info$solver
## [1] "Hooke-Jeeves"
## 
## $info$iterations
## [1] 26
cat('maximum expected profit: ', optres.HJ$fmin)
## maximum expected profit:  179.1865
cat('the best p: ', optres.HJ$xmin[1], '\nthe best Q: ', optres.HJ$xmin[2])
## the best p:  3.343994 
## the best Q:  105.6608


  1. Use the genetic algorithm to find the best Q, p, and expected profit:
# Genetic algorithm
sim_func_GA <- function(pQ) {
  p <- pQ[1]
  Q <- pQ[2]
  sim_res <- c()
  for (sim_time in 1:sim_times) {
    sim_res[sim_time] <- cal_profit(a, b, p, e_norm[sim_time], Q, D, c, v, s)
  }
  res <- mean(sim_res)
  return(-res)
}
optres.GA <- rbga(stringMin=c(1, 10),
                  stringMax=c(5.5, 200),
                  iters=50,
                  popSize=50,
                  evalFunc=sim_func_GA,
                  showSettings=TRUE,
                  verbose=FALSE)
## GA Settings
##   Type                  = 
##   Population size       = 50
##   Number of Generations = 50
##   Elitism               = 10
##   Mutation Chance       = 0.333333333333333
## 
## Search Domain
##   Var 1 = [1,5.5]
##   Var 2 = [10,200]
plot(optres.GA)

optres.GA
## $type
## [1] "floats chromosome"
## 
## $stringMin
## [1]  1 10
## 
## $stringMax
## [1]   5.5 200.0
## 
## $popSize
## [1] 50
## 
## $iters
## [1] 50
## 
## $suggestions
## NULL
## 
## $population
##          [,1]      [,2]
##  [1,] 3.38489 104.25053
##  [2,] 3.38489 104.25053
##  [3,] 3.38489 104.25053
##  [4,] 3.38489 104.25053
##  [5,] 3.38489 104.25053
##  [6,] 3.38489 104.25053
##  [7,] 3.38489 104.25053
##  [8,] 3.38489 104.25053
##  [9,] 3.38489 104.25053
## [10,] 3.38489 104.25053
## [11,] 3.38489 100.38453
## [12,] 3.48149 104.25053
## [13,] 3.38489 104.25053
## [14,] 3.38489 111.98253
## [15,] 3.38489 104.25053
## [16,] 3.38489 104.25053
## [17,] 3.48149 111.98253
## [18,] 3.38489 104.25053
## [19,] 3.38489 100.38453
## [20,] 3.38489 104.25053
## [21,] 3.38489 111.98253
## [22,] 3.38489 100.38453
## [23,] 3.38489 104.25053
## [24,] 3.38489 111.98253
## [25,] 3.28829 108.11653
## [26,] 3.19169 100.38453
## [27,] 3.38489 115.84853
## [28,] 3.38489 104.25053
## [29,] 3.28829 108.11653
## [30,] 3.38489 100.38453
## [31,] 3.38489 104.25053
## [32,] 3.48149 104.25053
## [33,] 3.38489 104.25053
## [34,] 3.19169 104.25053
## [35,] 3.38489 100.38453
## [36,] 3.28829 100.38453
## [37,] 3.38489 104.25053
## [38,] 3.28829 104.25053
## [39,] 3.28829 104.25053
## [40,] 3.48149 104.25053
## [41,] 3.09509 108.11653
## [42,] 3.28829 104.25053
## [43,] 3.38489 104.25053
## [44,] 3.38489 108.11653
## [45,] 3.48149 104.25053
## [46,] 3.38489 104.25053
## [47,] 3.38489 111.98253
## [48,] 3.38489  96.51853
## [49,] 3.28829 115.84853
## [50,] 3.38489 115.84853
## 
## $elitism
## [1] 10
## 
## $mutationChance
## [1] 0.3333333
## 
## $evaluations
##  [1] -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282
##  [8] -179.1282 -179.1282 -179.1282 -178.7449 -178.3598 -179.1282 -178.1333
## [15] -179.1282 -179.1282 -176.6743 -179.1282 -178.7449 -179.1282 -178.1333
## [22] -178.7449 -179.1282 -178.1333 -179.0699 -175.8353 -177.0283 -179.1282
## [29] -179.0699 -178.7449 -179.1282 -178.3598 -179.1282 -177.4482 -178.7449
## [36] -177.8768 -179.1282 -178.8283 -178.8283 -178.3598 -176.3339 -178.8283
## [43] -179.1282 -178.8638 -178.3598 -179.1282 -178.1333 -177.6139 -177.8888
## [50] -177.0283
## 
## $best
##  [1] -177.6347 -177.6347 -178.5242 -178.6062 -178.6062 -178.7421 -179.0595
##  [8] -179.0595 -179.0595 -179.0595 -179.0595 -179.0595 -179.0595 -179.0595
## [15] -179.0595 -179.0595 -179.1234 -179.1234 -179.1234 -179.1234 -179.1234
## [22] -179.1234 -179.1234 -179.1234 -179.1234 -179.1282 -179.1282 -179.1282
## [29] -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282
## [36] -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282
## [43] -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282 -179.1282
## [50] -179.1282
## 
## $mean
##  [1]  -85.42583 -124.75285 -131.61069 -127.86690 -145.75143 -138.31263
##  [7] -142.78742 -147.74326 -143.45019 -144.26812 -149.71915 -137.94208
## [13] -135.21243 -135.22271 -137.57390 -152.22573 -146.49746 -151.88366
## [19] -147.08984 -148.24344 -148.18720 -152.57597 -152.00061 -151.28364
## [25] -137.20517 -133.26581 -120.68820 -121.04340 -124.86587 -110.72414
## [31]  -96.43473 -127.30789 -115.92518 -133.13749 -129.14900 -133.41012
## [37] -127.59095 -136.57347 -148.37359 -160.92084 -162.30709 -165.73745
## [43] -162.37215 -163.48453 -166.38559 -174.46995 -174.76103 -177.26523
## [49] -177.43218 -178.55466
## 
## attr(,"class")
## [1] "rbga"
cat('maximum expected profit: ', -min(optres.GA$evaluations))
## maximum expected profit:  179.1282
cat('the best p: ', optres.GA$population[which.min(optres.GA$evaluations), 1], '\nthe best Q: ', optres.GA$population[which.min(optres.GA$evaluations), 2])
## the best p:  3.38489 
## the best Q:  104.2505


  1. Which algorithm gives you the optimal decisions?

Hooke-Jeeves algorithm


  1. Repeat (a) - (c) by assuming ε follows an exponential distribution with mean 10:
  2. Hooke-Jeeves algorithm
# Hooke-Jeeves algorithm
optres.HJ <- fminsearch(sim_func_HJ,   # object function
                   c(5, 180),                # initial values
                   e=e_exp,
                   lower=c(1, 10),          # lower bound
                   upper=c(5.5, 200),         # upper bound
                   method=c("Hooke-Jeeves"), # function method
                   minimize=FALSE)           # maximize objective function

optres.HJ
## $xmin
## [1]  3.484375 99.058764
## 
## $fmin
## [1] 208.6432
## 
## $count
## [1] 402
## 
## $convergence
## [1] 0
## 
## $info
## $info$solver
## [1] "Hooke-Jeeves"
## 
## $info$iterations
## [1] 26
cat('maximum expected profit: ', optres.HJ$fmin)
## maximum expected profit:  208.6432
cat('the best p: ', optres.HJ$xmin[1], '\nthe best Q: ', optres.HJ$xmin[2])
## the best p:  3.484375 
## the best Q:  99.05876
  1. Genetic algorithm
# Genetic algorithm
sim_func_GA <- function(pQ) {
  p <- pQ[1]
  Q <- pQ[2]
  sim_res <- c()
  for (sim_time in 1:sim_times) {
    sim_res[sim_time] <- cal_profit(a, b, p, e_exp[sim_time], Q, D, c, v, s)
  }
  res <- mean(sim_res)
  return(-res)
}
optres.GA <- rbga(stringMin=c(1, 10),
                  stringMax=c(5.5, 200),
                  iters=50,
                  popSize=50,
                  evalFunc=sim_func_GA,
                  showSettings=TRUE,
                  verbose=FALSE)
## GA Settings
##   Type                  = 
##   Population size       = 50
##   Number of Generations = 50
##   Elitism               = 10
##   Mutation Chance       = 0.333333333333333
## 
## Search Domain
##   Var 1 = [1,5.5]
##   Var 2 = [10,200]
plot(optres.GA)

optres.GA
## $type
## [1] "floats chromosome"
## 
## $stringMin
## [1]  1 10
## 
## $stringMax
## [1]   5.5 200.0
## 
## $popSize
## [1] 50
## 
## $iters
## [1] 50
## 
## $suggestions
## NULL
## 
## $population
##           [,1]      [,2]
##  [1,] 3.462577  98.22277
##  [2,] 3.462577  98.22277
##  [3,] 3.462577  98.22277
##  [4,] 3.462577  98.22277
##  [5,] 3.462577  98.22277
##  [6,] 3.462577  98.22277
##  [7,] 3.462577  98.22277
##  [8,] 3.462577  98.22277
##  [9,] 3.462577  98.22277
## [10,] 3.462577  98.22277
## [11,] 3.269377  98.22277
## [12,] 3.365977 109.82077
## [13,] 3.365977  98.22277
## [14,] 3.462577 109.82077
## [15,] 3.462577  94.35677
## [16,] 3.365977 102.08877
## [17,] 3.462577  98.22277
## [18,] 3.559177  98.22277
## [19,] 3.559177 105.95477
## [20,] 3.559177  98.22277
## [21,] 3.462577 105.95477
## [22,] 3.462577  98.22277
## [23,] 3.559177  98.22277
## [24,] 3.559177 102.08877
## [25,] 3.462577 102.08877
## [26,] 3.462577  98.22277
## [27,] 3.462577  98.22277
## [28,] 3.559177  98.22277
## [29,] 3.462577  94.35677
## [30,] 3.462577  98.22277
## [31,] 3.462577  98.22277
## [32,] 3.462577  98.22277
## [33,] 3.559177  98.22277
## [34,] 3.365977 102.08877
## [35,] 3.655777  98.22277
## [36,] 3.462577 102.08877
## [37,] 3.462577  98.22277
## [38,] 3.559177  98.22277
## [39,] 3.462577  98.22277
## [40,] 3.462577  90.49077
## [41,] 3.462577  94.35677
## [42,] 3.462577 102.08877
## [43,] 3.365977  94.35677
## [44,] 3.365977  98.22277
## [45,] 3.559177  98.22277
## [46,] 3.559177  86.62477
## [47,] 3.462577  98.22277
## [48,] 3.462577  98.22277
## [49,] 3.559177  98.22277
## [50,] 3.462577  98.22277
## 
## $elitism
## [1] 10
## 
## $mutationChance
## [1] 0.3333333
## 
## $evaluations
##  [1] -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633
##  [8] -208.5633 -208.5633 -208.5633 -205.0669 -207.2485 -207.5272 -206.8525
## [15] -207.7658 -208.1369 -208.5633 -208.3977 -206.8945 -208.3977 -207.8926
## [22] -208.5633 -208.3977 -207.8700 -208.5195 -208.5633 -208.5633 -208.3977
## [29] -207.7658 -208.5633 -208.5633 -208.5633 -208.3977 -208.1369 -207.1801
## [36] -208.5195 -208.5633 -208.3977 -208.5633 -205.6865 -207.7658 -208.5195
## [43] -205.7327 -207.5272 -208.3977 -204.9262 -208.5633 -208.5633 -208.3977
## [50] -208.5633
## 
## $best
##  [1] -206.7579 -206.7579 -206.7579 -206.7579 -208.5034 -208.5106 -208.5106
##  [8] -208.5106 -208.5106 -208.5106 -208.5106 -208.5633 -208.5633 -208.5633
## [15] -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633
## [22] -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633
## [29] -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633
## [36] -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633
## [43] -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633 -208.5633
## [50] -208.5633
## 
## $mean
##  [1]  -73.85725 -116.59258 -142.87133 -173.13148 -183.99725 -171.18074
##  [7] -177.18945 -179.25865 -161.85591 -186.19269 -185.84855 -175.72294
## [13] -168.17985 -154.82002 -163.61336 -178.01475 -171.64396 -166.93193
## [19] -153.29073 -167.66363 -165.06316 -160.84203 -184.16084 -170.02121
## [25] -181.10034 -145.25868 -142.88455 -145.05200 -145.67947 -118.26703
## [31] -128.45080 -144.54103 -155.76226 -141.03649 -136.39536 -157.69379
## [37] -156.46570 -165.67504 -164.29991 -179.78745 -183.44745 -184.80815
## [43] -178.71342 -195.33059 -199.00036 -202.49240 -204.44184 -205.00994
## [49] -206.47907 -208.02220
## 
## attr(,"class")
## [1] "rbga"
cat('maximum expected profit: ', -min(optres.GA$evaluations))
## maximum expected profit:  208.5633
cat('the best p: ', optres.GA$population[which.min(optres.GA$evaluations), 1], '\nthe best Q: ', optres.GA$population[which.min(optres.GA$evaluations), 2])
## the best p:  3.462577 
## the best Q:  98.22277
  1. Which algorithm gives you the optimal decisions?

Hooke-Jeeves algorithm