HW3_Group2

組員：唐思琪、林書霆、柯炯名、簡佑臻、王婷、胡祐銓

Q1

(1)

library(pracma) # Hooke-Jeeves algorithm 套件

# 需求隨價格進行變化，輸入價格傳回對應需求量。
D=function(price, a=200, b=35){
  D.val=a-b*price+rnorm(1,0,20) # noise 為常態分配
  return(max(D.val,0))
}

profit1 = function(x){
  S = 5000 # 模擬五千次
  revenue = rep(0,S)   # 初始化
  price=floor(x[2])    # 取整數
  quantity=floor(x[1]) # 取整數
  for(i in 1:S){
    demand = D(price)  # 從demand分佈產生需求數量
    # 計算營收(套用公式)
    revenue[i] = price*min(demand, quantity) - 
      quantity*1 + 
      0.5*max(quantity-demand, 0) - 
      max(demand-quantity, 0)*1 
  }
  mean(revenue) # 營收期望值
}

# Hooke-Jeeves 最佳化
result = fminsearch(profit1,
               c(10,3),
               lower=c(10,1), 
               upper=c(200,5.5),
               method=c("Hooke-Jeeves"),
               minimize=FALSE)
result

$xmin
[1] 83.617156  4.212381

$fmin
[1] 165.4497

$count
[1] 493

$convergence
[1] 0

$info
$info$solver
[1] "Hooke-Jeeves"

$info$iterations
[1] 26

(2)

# 輸入價格回傳對應需求，noise 改成指數分佈
expD=function(price,a=200,b=35){
  D.val=a-b*price+rexp(1, rate=1/10) # noise 改為指數分佈
  return(max(D.val, 0))
}

profit2 = function(x){
  S = 5000
  revenue = rep(0, S)
  price=floor(x[2])
  quantity=floor(x[1])
  for(i in 1:S){
    demand = expD(price)
    revenue[i] = price*min(demand,quantity) - quantity*1 + 0.5*max(quantity-demand,0) - max(demand-quantity,0)*1
  }
  mean(revenue)
}

result2 = fminsearch(profit2,
                     c(10,3),
                     lower=c(10,1), 
                     upper=c(200,5.5),
                     method=c("Hooke-Jeeves"),
                     minimize=FALSE)
result2

$xmin
[1] 81.745067  4.709776

$fmin
[1] 199.3369

$count
[1] 479

$convergence
[1] 0

$info
$info$solver
[1] "Hooke-Jeeves"

$info$iterations
[1] 26

Q2

(1)

S = 10000 # 模擬次數
before_demand = round(rnorm(S,9000,2000),0) # 模擬比賽前的需求
price = 25
qty = seq(500,21000,1000) # 訂購數量區間

# 初始化賽前/賽後營收
early_earning = matrix(nrow = S, ncol = length(qty)) 
late_earning = matrix(nrow = S, ncol = length(qty))
# 初始化比賽前的剩餘銷量與最終營收
remain_qty = matrix(nrow=S, ncol=length(qty))
final_earning = matrix(nrow=S, ncol=length(qty))
# 初始化營收期望值與CVAR
avg_revenue = c()
cvar = c()

# 模擬比賽後的需求
win = sample(c(0,1), S, prob = c(0.6,0.4), replace=TRUE) # 模擬球賽的輸贏
after_demand = matrix(nrow=S, ncol=2)
for(i in 1:S){
  after_demand[i,1] = win[i] # column1存放球賽的結果
  if(win[i] == 0){ # 輸球，賽後結果服從gamma分配(mu=2000,sigma=1000) -> ~gamma(scale=500,shape=4)
    after_demand[i,2] = rgamma(1,shape=4,scale=500) # 賽後輸球的需求存入column2
  } else{
    after_demand[i,2] = rnorm(1,6000,2000) # 賽後贏球的需求存入column2
  }
}

# 模擬營收
for(q in 1:length(qty)){
  for(i in 1:S){
    early_earning[i,q] = min(qty[q], before_demand[i])*15 # 賽前獲利 = 賣出min(進貨量,比賽前需求)*淨利
    remain_qty[i,q] = qty[q]-before_demand[i] # 剩下的商品數量
    if(win[i]==0){ # 輸球
      salvage = 15-12.5 # 殘值 每件只賺2.5元
      # 賽後獲利 = min(剩餘數量,比賽後需求)*殘值利潤-剩下賣不出去的成本(虧損10元)
      late_earning[i,q] = min(max(remain_qty[i,q], 0), round(after_demand[i], 0))*salvage - 
        10*max((remain_qty[i,q] - after_demand[i]), 0)
    } else{ # 贏球
      salvage = 25-10
      # 賽後獲利 = min(剩餘數量,比賽後需求)*殘值利潤-剩下賣不出去的成本(虧損10元)
      late_earning[i,q] = min(max(remain_qty[i,q],0), round(after_demand[i],0))*salvage -
        10*max((remain_qty[i,q] - after_demand[i]), 0)
    }
    final_earning[i,q] = early_earning[i,q] + late_earning[i,q]
  }
}

for(i in 1:length(qty)){
  avg_revenue[i] = mean(final_earning[,i])   # 進貨數量所對應的平均收益
  cvar[i] = quantile(final_earning[,i], 0.1) # 進貨數量所對應的最差收益值
}

# 畫出訂貨數量的區間
plot(qty, avg_revenue, type='l',xlab="production quantity",lwd=3)
points(qty[which.max(avg_revenue)], max(avg_revenue), col="red",pch=1,cex=1)

# 計算最佳訂購量與其獲利
qty[which.max(avg_revenue)] # 最佳訂購量

[1] 9500

max(avg_revenue) # 最佳訂購量下的獲利

[1] 115441

(2)

qty[which.max(cvar)] # 在最大化CVAR(q=10%)條件下最佳訂購量

[1] 6500

max(cvar) # 在最大化CVAR(q=10%)條件下最佳訂購量下的獲利

[1] 95872.5

(3)

# 計算最佳service level
cu1=25-10
cu2=15-12.5
co=10
normal_service_level = cu1/(cu1+co)    # 贏球的service_level
disscount_service_level = cu2/(cu2+co) # 輸球的service_level

opt_before_demand = quantile(before_demand,normal_service_level) # 取得球賽前滿足服務水準的銷量
opt_lose_demand = quantile(after_demand[after_demand[,1]==0,],disscount_service_level) # 輸球的滿足服務水準的銷量
opt_win_demand = quantile(after_demand[after_demand[,1]==1,],normal_service_level) # 贏球的滿足服務水準的銷量

# 初始化
early_earning2 = matrix(nrow = S, ncol = 1)
remain_qty2 = matrix(nrow = S, ncol = 1)
late_earning2 = matrix(nrow = S, ncol = 1)
final_earning2 = matrix(nrow = S, ncol = 1)

# 模擬獲利
for(i in 1:S){
  early_earning2[i] = min(before_demand[i],opt_before_demand)*15
  remain_qty2[i] = before_demand[i]-opt_before_demand
  if(win[i]==0){
    late_earning2[i] = min(max(remain_qty2[i],0), opt_lose_demand)*cu2 -
      10*max((remain_qty2[i]-opt_lose_demand),0)
  }else{
    late_earning2[i] = min(max(remain_qty2[i],0), opt_win_demand)*cu1 -
      10*max((remain_qty2[i]-opt_win_demand),0)
  }
  final_earning2[i] = early_earning2[i]+late_earning2[i]
}

mean(final_earning2)-max(avg_revenue)

[1] 10848.02

Q3

S = 1000
dict = c('jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec')
c1 = c('jan','sep')
c2 = c('may','jun','jul','aug')
c3 = c('feb','mar','apr','oct','nov','dec')
history_quantity = c(80,90,120,105,95,75,70,70,110,105,90,65) #from Jun
p = matrix(nrow = 12, ncol = S)
p[4,] = 63 #四月初剩餘的分析師人數
Q = seq(11,110,1)
econmy.simulation = rnorm(S,0,0.05)#模擬每年的經濟變化
#個月的分析師需求受到經濟與noise影響
Month_Demand = function(month,n){
  his = history_quantity[month]
  #econmy = rnorm(1,0,0.05)
  econmy = econmy.simulation[n]#對應年的經濟狀況
  noise = rnorm(1,0,0.1) #noise的影響
  month_demand = round(his*(1+econmy)*(1+noise))
  return(month_demand)
}
#初始化每年每月的分析師需求
d = matrix(nrow = 12,ncol=1000)
for (s in 1:S) {
  for (i in 1:12) {
    d[i,s] = Month_Demand(i,s)
  }
}

revenue = matrix(nrow=12, ncol=S)
Q_Revenue_table = matrix(nrow = 2,ncol = 100)
avg_rev = rep(0, 100)
h = 1

(a)

simulation = function(Q){
  for (q in Q){ # 所需要聘的分析師數量
    for (s in c(1:S)){    # 模擬到第幾年
      for (i in c(4:11)){ # 初始化該年有的分析師數目
        if (dict[i] %in% c1){ # 1、9月留存率0.8~1
          retention_rate = runif(1, 0.8, 1)
        } else if (dict[i] %in% c2){ # 5、6月、7、8月留存率0.95~1
          retention_rate = runif(1, 0.95, 1)
        } else { # 2、3、4月、10、11、12月留存率0.9~1
          retention_rate = runif(1, 0.9, 1)
        }

        # 6月進行招募7月的分析師數為6月的留存人數+接受發出offer的人數；其他時間為上個月的分析師數＊留存率 
        if (i == 6){
          p[i+1,s] = round(p[i,s]*retention_rate + rbinom(1, q, 0.7))
        } else {
          p[i+1,s] = round(p[i,s]*retention_rate)
        }
      p[1,s] = round(p[12,s]*retention_rate)
      p[2,s] = round(p[1,s]*retention_rate)
      p[3,s] = round(p[2,s]*retention_rate)
      p[4,s] = round(p[3,s]*retention_rate) # 新一年4月的人數
      
      }
      # 比較模擬的需求數量與模擬的現存分析師算出獲利
      for (j in c(1:12)) {
        if (p[j,s] > d[j,s]) {     # 分析師數>需求數
          revenue[j,s] = 10000*d[j,s] - 6000*p[j,s]
        }else if(p[j,s] < d[j,s]){ # 分析師數<需求數
          revenue[j,s] = 4000*p[j,s] + 400*(d[j,s]-p[j,s])
        } else { # 分析師數=需求數
          revenue[j,s] = 4000*d[j,s]
        }
      }
    }
    Q_Revenue_table[1,h] = q # 紀錄數量
    Q_Revenue_table[2,h] = mean(colMeans(revenue)) # 數量所模擬出的年平均（1000年）
    h = h+1
  }
  Q_Revenue_table
}

result = simulation(Q)
result[1,which.max(result[2,])]

[1] 61

result[2,which.max(result[2,])]

[1] 248249.1

plot(x=Q, y=result[2,],
     xlab = "Q7-Flexible",
     ylab = "E-Flexible")
    #  ylim = c(15000, 350000))
points(x = result[1,which.max(result[2,])], y = max(result[2,]), col = "red",pch=22)

(b)

revenue = matrix(nrow=12, ncol=S)
Q_Revenue_table = matrix(nrow = 2,ncol = 100)
avg_rev = rep(0, 100)
h = 1

simulation2 = function(Q){
  for (q in Q){
    half_q = round(q/2)
    for (s in c(1:S)){
      for (i in c(4:11)){
        if (dict[i] %in% c1){
          retention_rate = runif(1, 0.8, 1)
        } else if (dict[i] %in% c2){
          retention_rate = runif(1, 0.95, 1)
        } else {
          retention_rate = runif(1, 0.9, 1)
        }

        if (i == 6 | i ==8){ # 在6月與8月都會進行招募
          p[i+1,s] = round(p[i,s]*retention_rate + rbinom(1, half_q, 0.7))
        } else {
          p[i+1,s] = round(p[i,s]*retention_rate)
        }

      p[1,s] = round(p[12,s]*retention_rate)
      p[2,s] = round(p[1,s]*retention_rate)
      p[3,s] = round(p[2,s]*retention_rate)
      p[4,s] = round(p[3,s]*retention_rate) # 新一年4月的人數
      
      }

      for (j in c(1:12)) {
        if (p[j,s] > d[j,s]){
          revenue[j,s] = 10000*d[j,s] - 6000*p[j,s]
        }else if(p[j,s] < d[j,s]){
          revenue[j,s] = 4000*p[j,s] + 400*(d[j,s]-p[j,s])
        } else { # equal
          revenue[j,s] = 4000*d[j,s]
        }
      }

    }
    Q_Revenue_table[1,h] = q
    Q_Revenue_table[2,h] = mean(colMeans(revenue))
    h = h+1
  }
  Q_Revenue_table
}

result = simulation2(Q)
result[1,which.max(result[2,])] # 數量

[1] 63

result[2,which.max(result[2,])] # 獲利

[1] 266721.9

plot(x=Q, y=result[2,],
     xlab = "Q7-Flexible",
     ylab = "E-Flexible")
points(x = result[1, which.max(result[2,])], y = max(result[2,]), col = "red", pch=22)

(c)

(a) Fixed start strategy

# Fixed start strategy + 12月進行招募
Fixed_Start = function(x){
  Q7 = floor(x[1])
  Q11 = floor(x[2])
  for (s in c(1:S)){
    for (i in c(4:11)){
      if (dict[i] %in% c1){
        retention_rate = runif(1, 0.8, 1)
      } else if (dict[i] %in% c2){
        retention_rate = runif(1, 0.95, 1)
      } else {
        retention_rate = runif(1, 0.9, 1)
      }

      if (i == 6){ # 6月招募7月上線
        p[i+1,s] = round(p[i,s]*retention_rate + rbinom(1, Q7, 0.7))
      }else if(i ==11){
        p[i+1,s] = round(p[i,s]*retention_rate + rbinom(1, Q11, 0.7))
      } 
      else {
        p[i+1,s] = round(p[i,s]*retention_rate)
      }
    p[1,s] = round(p[12,s]*retention_rate)
    p[2,s] = round(p[1,s]*retention_rate)
    p[3,s] = round(p[2,s]*retention_rate)
    p[4,s] = round(p[3,s]*retention_rate) # 新一年4月的人數
  }

    for (j in c(1:12)) {
      if (p[j,s] > d[j,s]) {
        revenue[j,s] = 10000*d[j,s] - 6000*p[j,s]
      }else if(p[j,s] < d[j,s]){
        revenue[j,s] = 4000*p[j,s] + 400*(d[j,s]-p[j,s])
      } else { # equal
        revenue[j,s] = 4000*d[j,s]
      }
    }
  }
  fix_start_revenue = mean(colMeans(revenue))
}

FsRevenue = fminsearch(Fixed_Start,
               c(11,55), # 起始點
               lower=c(11,11), 
               upper=c(110,110),
               method=c("Hooke-Jeeves"),
               minimize=FALSE)

FsRevenue # 獲利

$xmin
[1] 31.74335 41.76246

$fmin
[1] 260567.4

$count
[1] 446

$convergence
[1] 0

$info
$info$solver
[1] "Hooke-Jeeves"

$info$iterations
[1] 26

(b) Flexible start strategy

# Flexible start strategy + 12月進行招募
Flexible_Start = function(x){
  half_Q7 = floor(x[1]/2)
  Q11 = floor(x[2])
  for (s in c(1:S)){
    for (i in c(4:11)){
      if (dict[i] %in% c1){
        retention_rate = runif(1, 0.8, 1)
      } else if (dict[i] %in% c2){
        retention_rate = runif(1, 0.95, 1)
      } else {
        retention_rate = runif(1, 0.9, 1)
      }

      if (i == 6 || i ==8){  # 6月招募七月上線一半、8月招募9月上線一半
        p[i+1,s] = round(p[i,s]*retention_rate + rbinom(1, half_Q7, 0.7))
      }else if(i ==11){
        p[i+1,s] = round(p[i,s]*retention_rate + rbinom(1, Q11, 0.7))
      } 
      else {
        p[i+1,s] = round(p[i,s]*retention_rate)
      }
      
    p[1,s] = round(p[12,s]*retention_rate)
    p[2,s] = round(p[1,s]*retention_rate)
    p[3,s] = round(p[2,s]*retention_rate)
    p[4,s] = round(p[3,s]*retention_rate) # 新一年4月的人數
  }

    for (j in c(1:12)) {
      if (p[j,s] > d[j,s]) {
        revenue[j,s] = 10000*d[j,s] - 6000*p[j,s]
      }else if(p[j,s] < d[j,s]){
        revenue[j,s] = 4000*p[j,s] + 400*(d[j,s]-p[j,s])
      } else { # equal
        revenue[j,s] = 4000*d[j,s]
      }
    }
  }
  flex_start_revenue = mean(colMeans(revenue))
}

FlexRevenue = fminsearch(Flexible_Start,
               c(11,55), # 起始點
               lower=c(11,11), 
               upper=c(110,110),
               method=c("Hooke-Jeeves"),
               minimize=FALSE)
FlexRevenue

$xmin
[1] 67.86424 11.02790

$fmin
[1] 275723.8

$count
[1] 460

$convergence
[1] 0

$info
$info$solver
[1] "Hooke-Jeeves"

$info$iterations
[1] 26