InventoryByQ.R

#================================================================
#     Inventory simulation model under uncertainty 
#         based on Montecarlo simulation
#
# This model estimates the total cost and inventory for
# different value of quantity of product order per cycle, 
# namely, given by Nq*StepQ.
#
# The model also shows a figure of unsatisfied demand probability.
#================================================================
#  Coded and developed by Ebert Brea (c)2021
#===============================================================
rm(list=ls())

require(graphics)
require(latex2exp)

## Loading required package: latex2exp

require(stringi)

## Loading required package: stringi

require(stringr)

## Loading required package: stringr

require(lhs)

## Loading required package: lhs

library(triangle)

if(file.exists("filltxt")){close(filtxt)}
#==============================================================================
#           Parameters of the simulation algorithm
#==============================================================================
R=1000     # Number of simulation replications 
n0<-14     # Number of days per replacement cycle of product
Nq<-12     # Number of product quantity sampling
StepQ<-200 # Step of product quantity sampling
a<-0       # Minimum Demand of product
b<-250     # Maximum Demand of product 
md<- 100   # mode of triangular distribution
hi<-10     # Inventory cost per day and per product unit 
ho<-20     # Opportunity cost due to missing product, per day and product unit 
###############################################################################
#               Demand function according to id indicator
###############################################################################
#  id=1   Uniform distribution demand
#  id=2   Triangular distribution demand

id<-1

demand <-function(id){
  if(id==1){return(runif(1,a,b))}
  if(id==2){ return(rtriangle(1,a,b,md))}
}
########################################

inventory <-matrix( nrow=n0, ncol = R, byrow=TRUE)

totalcost <-matrix( nrow=R, ncol = 1)

risk <-matrix( nrow=Nq, ncol = 2, byrow=TRUE)

TC <- matrix( nrow=Nq, ncol = R)

TCSummary  <- matrix( nrow=Nq, ncol = 6)

INV <- matrix( nrow=Nq, ncol = R)

fil <- tempfile("data")

##############################################################
# "R" simulation replications of "n0" days per cycle
##############################################################
for(iq in 1:Nq){

InvNegTime=0
SamplingTime=0

Tinv=0
TNegInv=0

risk[iq,1] = iq*StepQ

for(r in 1:R){
  ##############################################################
  #        Starting simulation replications
  ##############################################################
  
  totalcost[r]=0
  
  # Daily inventory at the end of the day 1
  inventory[1,r] <- risk[iq,1] - demand(id) 
  
  # Daily cost  
  if(inventory[1,r] >= 0){
    totalcost[r] <- totalcost[r] + hi*inventory[1,r]}
  else{
    totalcost[r] <- totalcost[r] - 1*ho*inventory[1,r]
    NegInv <- 1
    TNegInv <- TNegInv + inventory[1,r]
  }  
  
  for(i in 2:n0){
    NegInv = 0
    # Daily inventory at the end of each ith day
    inventory[i,r] <- inventory[i-1,r]-demand(id)
    
    # Daily cost at the end of each ith day  
    if(inventory[i,r] >= 0){
      totalcost[r] <- totalcost[r] + hi*inventory[i,r]}
    else{
      totalcost[r] <- totalcost[r] - 1*ho*inventory[i,r]
      NegInv <- 1
      TNegInv <- TNegInv + inventory[i,r]
    }
    
  }#end for i
  
  InvNegTime <- InvNegTime + NegInv
  
#============================================================
#        Collect Total Inventory Cost by StepQ & r
#============================================================

    TC[iq,r] <- totalcost[r]

#============================================================
#        Collect Inventory Cost by StepQ & r
#============================================================
  
  INV[iq,r] <- inventory[n0,r]
  
}#end for r

#============================================================
#        Estimation of Pr{inventory < 0} by StepQ
#===========================================================

    risk[iq,2] = InvNegTime/R

#============================================================
#        Statistic od Total cost by StepQ
#===========================================================

  TCSummary[iq,1] <- min(TC[iq,])
  TCSummary[iq,2] <- mean(TC[iq,])
  TCSummary[iq,3] <- max(TC[iq,])

  TCSummary[iq,4] <- quantile(TC[iq,],0.25)
  TCSummary[iq,5] <- quantile(TC[iq,],0.50)
  TCSummary[iq,6] <- quantile(TC[iq,],0.75)

}#end iq

#============================================================
#        Graphic of Pr{inventory < 0} vs. q 
#===========================================================

plot(ylim=c(0, 1),risk[,1],risk[,2],main=TeX('$p\\{\\sum_{k=1}^{n_{0}} D_{k}>q\\}$'),xlab=TeX('$q$'),ylab="Probability",type="l",lwd = 2,col="blue")
abline(h = 1, lty = 2, col = "red")
abline(h = 0, lty = 2, col = "blue")

print(risk)

##       [,1]  [,2]
##  [1,]  200 1.000
##  [2,]  400 1.000
##  [3,]  600 1.000
##  [4,]  800 1.000
##  [5,] 1000 0.998
##  [6,] 1200 0.968
##  [7,] 1400 0.894
##  [8,] 1600 0.731
##  [9,] 1800 0.409
## [10,] 2000 0.185
## [11,] 2200 0.048
## [12,] 2400 0.004

#============================================================
#        Graphic of Total Inventory Cost vs. q 
#===========================================================

dfCost <- as.data.frame(TCSummary)

MaxLimit <-ceiling(max(TCSummary[,3])/1000)

plot(ylim=c(0, MaxLimit+1),risk[,1],dfCost$V1/1000,main="Total cost",xlab=TeX('$q$'),ylab="MUS$",type="l",lwd = 2,col="blue")
lines(risk[,1],dfCost$V2/1000,col="darkgreen",type="l",lwd = 2)
lines(risk[,1],dfCost$V3/1000,col="red",type="l",lwd = 2)
abline(h = MaxLimit, lty = 2, col = "black")
abline(h = 0, lty = 2, col = "black")
legend("topright",c("Maximum total cost","Average total cost","Minimum total cost"), fill=c("red","darkgreen","blue"))