HW4_Group2

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

Q1

##Multi-Period: A Simple Inventory Model

#Holding cost per item per day
h.val=0.6

unit.price=400
unit.cost=250
#Profit margin per TV sold
unit.profit=unit.price-unit.cost

#Parameters of (r, q) inventory policy
r=30
q=130

#cost per delivery from supplier to the dealer
shipping.cost=500

#The number of simulation runs
M=5000
#The number of days
T_value=360

#The experiment will be repeated for 
# m in 1:M repetitions
# t in 1:T duration

#Vector initialization
#SL: service level
SL=rep(NA,M)
daily.profit=rep(NA,M)
set.seed(5566)

start.time <- Sys.time()
for(m in 1:M){
  #Initialization: Set inventory on-hand=q in day 0
  inv.onhand=c(q,rep(NA,T_value-1))
  sales=rep(NA,T_value)
  orderQ=rep(0,T_value)
  #Keep track of inventory on-order but not arrived yet
  inv.onorder=c()
  inv.onorder.stamp=c()
  #The number of deliveries from the supplier
  delivery=0
  #The total amount of lost sales 
  loss=0
  #The number of stockouts in a repetition
  stockout=0
  #Simulate Normal demand for T periods
  d=rnorm(T_value, mean = 12, sd = 2)
  #
  for(t in 1:T_value){
     if(t==1){
       #Selling process
       sales[t]=min(inv.onhand[t],d[t])
       #Check stockout
       if(d[t]>sales[t]){
         stockout=stockout+1
         loss=loss+(d[t]-sales[t])
       }
       #Compute inventory position
       inv.position=inv.onhand[1]
       #Ordering mechanism of (r, q)
       if(inv.position<=r){
         #Compute order quantity
         Q=q-inv.position
         #Update inventory on-order
         inv.onorder=c(inv.onorder,Q)
         #Simulate stochastic lead time
         time.to.arrive=t+sample(c(4,5),1,prob=c(0.5,0.5))+1
         #Update time stamp for invenory on-order
         inv.onorder.stamp=c(inv.onorder.stamp,time.to.arrive)
         #Record order quantity
         orderQ[t]=Q
       }
       #Update inventory on-hand in the end of each day
       inv.onhand[t+1]=inv.onhand[t]-sales[t]
     }
     #
     if(t>1){
        #Check if any inventory on-order should arrive
        if(any(inv.onorder.stamp==t)){
           #Update the number of deliveries
           delivery=delivery+1
           #Compute the total of arrived inventories
           index=which(inv.onorder.stamp==t)
           arrival=sum(inv.onorder[index])
           #Update inventory on-hand before starting the day
           inv.onhand[t]=inv.onhand[t]+arrival
           #Removed those just arrived from inventory on-order
           inv.onorder=inv.onorder[-index]
           inv.onorder.stamp=inv.onorder.stamp[-index]
           #cat("arrival time:",t)
        } #end t==1
       #Record sales
       sales[t]=min(inv.onhand[t],d[t])
       #Check if any stockout takes place
       if(d[t]>sales[t]){
         stockout=stockout+1
         loss=loss+d[t]-sales[t]
       }
       #Update inventory position
       inv.position=inv.onhand[t]+sum(inv.onorder)
       #Ordering mechanism of (r, q) inventory policy
       if(inv.position<=r){
         Q=q-inv.position
         inv.onorder=c(inv.onorder,Q)
         time.to.arrive=t+sample(c(4,5),1,prob=c(0.5,0.5))+1
         inv.onorder.stamp=c(inv.onorder.stamp,time.to.arrive)
         orderQ[t]=Q
       }
       #Update inventory on-hand in the end of the day
       inv.onhand[t+1]=inv.onhand[t]-sales[t]
     } #end t>1
     #cat("day:",t,";sales:",sales[t],";onhand:",inv.onhand[t],";onorder:",
     #     inv.onorder,";order:",orderQ[t],"\n") 
  } #end for t
  
  SL[m]=(T_value-stockout)/T_value  #Calculate the service rate; 
    
  #Calculate the average daily profit
  daily.profit[m]=(unit.profit*sum(sales)-sum(h.val*inv.onhand)-
                   shipping.cost*delivery-unit.profit*loss)/T_value
  #Simulation progress marker
  if(m%%100==0){cat("Repetitions:",m,"\n")} 

} #end for m

Repetitions: 100 
Repetitions: 200 
Repetitions: 300 
Repetitions: 400 
Repetitions: 500 
Repetitions: 600 
Repetitions: 700 
Repetitions: 800 
Repetitions: 900 
Repetitions: 1000 
Repetitions: 1100 
Repetitions: 1200 
Repetitions: 1300 
Repetitions: 1400 
Repetitions: 1500 
Repetitions: 1600 
Repetitions: 1700 
Repetitions: 1800 
Repetitions: 1900 
Repetitions: 2000 
Repetitions: 2100 
Repetitions: 2200 
Repetitions: 2300 
Repetitions: 2400 
Repetitions: 2500 
Repetitions: 2600 
Repetitions: 2700 
Repetitions: 2800 
Repetitions: 2900 
Repetitions: 3000 
Repetitions: 3100 
Repetitions: 3200 
Repetitions: 3300 
Repetitions: 3400 
Repetitions: 3500 
Repetitions: 3600 
Repetitions: 3700 
Repetitions: 3800 
Repetitions: 3900 
Repetitions: 4000 
Repetitions: 4100 
Repetitions: 4200 
Repetitions: 4300 
Repetitions: 4400 
Repetitions: 4500 
Repetitions: 4600 
Repetitions: 4700 
Repetitions: 4800 
Repetitions: 4900 
Repetitions: 5000 

end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken

Time difference of 4.741276 secs

summary(daily.profit)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  647.4   714.7   728.8   729.4   743.2   818.7 

summary(SL)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6528  0.6778  0.6806  0.6818  0.6861  0.7139 

##Search for optimal (R*, Q*)
R=seq(0,60,5)
Q=seq(90,360,10)

search.val=as.matrix(expand.grid(R,Q))
dim(search.val)

[1] 364   2

#The (r, q) model as a function for simulation
InvModel=function(i){
  r=search.val[i,1]
  q=search.val[i,2]
  #Vector initialization
  SL=rep(NA,M)
  daily.profit=rep(NA,M)
  for(m in 1:M){
    #Initialization: Set inventory on-hand=q in day 0
    inv.onhand=c(q,rep(NA,T_value-1))
    sales=rep(NA,T_value)
    orderQ=rep(0,T_value)
    #Keep track of inventory on-order but not arrived yet
    inv.onorder=c()
    inv.onorder.stamp=c()
    #The number of deliveries from the supplier
    delivery=0
    #The total amount of lost sales 
    loss=0
    #The number of stockouts in a repetition
    stockout=0
    #Simulate Normal demand for T periods
    d=rnorm(T_value, mean = 12, sd = 2)
    #
    for(t in 1:T_value){
      if(t==1){
        #Selling process
        sales[t]=min(inv.onhand[t],d[t])
        #Check stockout
        if(d[t]>sales[t]){
          stockout=stockout+1
          loss=loss+(d[t]-sales[t])
        }
        #Compute inventory position
        inv.position=inv.onhand[1]
        #Ordering mechanism of (r, q)
        if(inv.position<=r){
          #Compute order quantity
          Q=q-inv.position
          #Update inventory on-order
          inv.onorder=c(inv.onorder,Q)
          #Simulate stochastic lead time
          time.to.arrive=t+sample(c(4,5),1,prob=c(0.5,0.5))+1
          #Update time stamp for invenory on-order
          inv.onorder.stamp=c(inv.onorder.stamp,time.to.arrive)
          #Record order quantity
          orderQ[t]=Q
        }
        #Update inventory on-hand in the end of each day
        inv.onhand[t+1]=inv.onhand[t]-sales[t]
      }
      #
      if(t>1){
        #Check if any inventory on-order should arrive
        if(any(inv.onorder.stamp==t)){
          #Update the number of deliveries
          delivery=delivery+1
          #Compute the total of arrived inventories
          index=which(inv.onorder.stamp==t)
          arrival=sum(inv.onorder[index])
          #Update inventory on-hand before starting the day
          inv.onhand[t]=inv.onhand[t]+arrival
          #Removed those just arrived from inventory on-order
          inv.onorder=inv.onorder[-index]
          inv.onorder.stamp=inv.onorder.stamp[-index]
          #cat("arrival time:",t)
        } #end t==1
        #Record sales
        sales[t]=min(inv.onhand[t],d[t])
        #Check if any stockout takes place
        if(d[t]>sales[t]){
          stockout=stockout+1
          loss=loss+d[t]-sales[t]
        }
        #Update inventory position
        inv.position=inv.onhand[t]+sum(inv.onorder)
        #Ordering mechanism of (r, q) inventory policy
        if(inv.position<=r){
          Q=q-inv.position
          inv.onorder=c(inv.onorder,Q)
          time.to.arrive=t+sample(c(4,5),1,prob=c(0.5,0.5))+1
          inv.onorder.stamp=c(inv.onorder.stamp,time.to.arrive)
          orderQ[t]=Q
        }
        #Update inventory on-hand in the end of the day
        inv.onhand[t+1]=inv.onhand[t]-sales[t]
      } #end t>1
      #cat("day:",t,";sales:",sales[t],";onhand:",inv.onhand[t],";onorder:",
      #     inv.onorder,";order:",orderQ[t],"\n") 
    } #end for t
    
    SL[m]=(T_value-stockout)/T_value  #Calculate the service rate; 
    
    #Calculate the average daily profit
    daily.profit[m]=(unit.profit*sum(sales)-sum(h.val*inv.onhand)-
                       shipping.cost*delivery-unit.profit*loss)/T_value
    
  } #end for m
  c(mean(SL),mean(daily.profit),sd(daily.profit))
}


InvModel(1)

##Use parallel CPU computing for optimal search
library(foreach)
library(doParallel)
registerDoParallel(cores=6)
getDoParWorkers()

[1] 6

start.time <- Sys.time()
sim.results=foreach(i=1:nrow(search.val),
             .combine=rbind,.verbose=F) %dopar% InvModel(i)

##Apply the Hooke-Jeeves search algorithm
InvModel=function(par){
  r=par[1]
  q=par[2]
  #Vector initialization
  SL=rep(NA,M)
  daily.profit=rep(NA,M)
  for(m in 1:M){
    #Initialization: Set inventory on-hand=q in day 0
    inv.onhand=c(q,rep(NA,T_value-1))
    sales=rep(NA,T_value)
    orderQ=rep(0,T_value)
    #Keep track of inventory on-order but not arrived yet
    inv.onorder=c()
    inv.onorder.stamp=c()
    #The number of deliveries from the supplier
    delivery=0
    #The total amount of lost sales 
    loss=0
    #The number of stockouts in a repetition
    stockout=0
    #Simulate Normal demand for T periods
    d=rnorm(T_value, mean = 12, sd = 2)
    #
    for(t in 1:T_value){
      if(t==1){
        #Selling process
        sales[t]=min(inv.onhand[t],d[t])
        #Check stockout
        if(d[t]>sales[t]){
          stockout=stockout+1
          loss=loss+(d[t]-sales[t])
        }
        #Compute inventory position
        inv.position=inv.onhand[1]
        #Ordering mechanism of (r, q)
        if(inv.position<=r){
          #Compute order quantity
          Q=q-inv.position
          #Update inventory on-order
          inv.onorder=c(inv.onorder,Q)
          #Simulate stochastic lead time
          time.to.arrive=t+sample(c(4,5),1,prob=c(0.5,0.5))+1
          #Update time stamp for invenory on-order
          inv.onorder.stamp=c(inv.onorder.stamp,time.to.arrive)
          #Record order quantity
          orderQ[t]=Q
        }
        #Update inventory on-hand in the end of each day
        inv.onhand[t+1]=inv.onhand[t]-sales[t]
      }
      #
      if(t>1){
        #Check if any inventory on-order should arrive
        if(any(inv.onorder.stamp==t)){
          #Update the number of deliveries
          delivery=delivery+1
          #Compute the total of arrived inventories
          index=which(inv.onorder.stamp==t)
          arrival=sum(inv.onorder[index])
          #Update inventory on-hand before starting the day
          inv.onhand[t]=inv.onhand[t]+arrival
          #Removed those just arrived from inventory on-order
          inv.onorder=inv.onorder[-index]
          inv.onorder.stamp=inv.onorder.stamp[-index]
          #cat("arrival time:",t)
        } #end t==1
        #Record sales
        sales[t]=min(inv.onhand[t],d[t])
        #Check if any stockout takes place
        if(d[t]>sales[t]){
          stockout=stockout+1
          loss=loss+d[t]-sales[t]
        }
        #Update inventory position
        inv.position=inv.onhand[t]+sum(inv.onorder)
        #Ordering mechanism of (r, q) inventory policy
        if(inv.position<=r){
          Q=q-inv.position
          inv.onorder=c(inv.onorder,Q)
          time.to.arrive=t+sample(c(4,5),1,prob=c(0.5,0.5))+1
          inv.onorder.stamp=c(inv.onorder.stamp,time.to.arrive)
          orderQ[t]=Q
        }
        #Update inventory on-hand in the end of the day
        inv.onhand[t+1]=inv.onhand[t]-sales[t]
      } #end t>1
      #cat("day:",t,";sales:",sales[t],";onhand:",inv.onhand[t],";onorder:",
      #     inv.onorder,";order:",orderQ[t],"\n") 
    } #end for t
    
    SL[m]=(T_value-stockout)/T_value  #Calculate the service rate; 
    
    #Calculate the average daily profit
    daily.profit[m]=(unit.profit*sum(sales)-sum(h.val*inv.onhand)-
                       shipping.cost*delivery-unit.profit*loss)/T_value
    
  } #end for m
  #c(mean(SL),mean(daily.profit),sd(daily.profit))
  mean(daily.profit)
}

library(pracma)
res=fminsearch(InvModel,c(30,130),
         lower=c(20,80), 
         upper=c(150,250),
         method=c("Hooke-Jeeves"),
         minimize=FALSE)

res

$xmin
[1]  87.63159 215.07648

$fmin
[1] 1700.547

$count
[1] 479

$convergence
[1] 0

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

$info$iterations
[1] 26