Inventory.R

#=============================================================
#     Inventory simulation model under uncertainty 
#         based on Montecarlo simulation·
#=============================================================
# 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
q<-2000    # Quantity of product order per cycle of product
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)

fil <- tempfile("data")

##############################################################
# "R" simulation replications of "n0" days per cycle
##############################################################
InvNegTime=0
Tinv=0
TNegInv=0

r<-1

for(r in 1:R){
  ##############################################################
  #        Starting simulation replications
  ##############################################################

  totalcost[r]=0

# Daily inventory at the end of the day 1
  inventory[1,r] <- q - 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 the day
    #inventory[i,r] <- inventory[i-1,r]-rtriangle(1,a,b,md)
    inventory[i,r] <- inventory[i-1,r]-demand(id)
    
    # Daily cost  
    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
      
  
}#end for r
#============================================================
#        Statistical report of simulation
#===========================================================



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

  Probab = InvNegTime/R

  print(Probab)

## [1] 0.187

  summary(totalcost)

##        V1        
##  Min.   :104183  
##  1st Qu.:134072  
##  Median :148892  
##  Mean   :149890  
##  3rd Qu.:163500  
##  Max.   :214326

 min(inventory)

## [1] -529.523

 mean(inventory)

## [1] 1062.443

 max(inventory)

## [1] 1999.774

 quantile(inventory,0.25)

##      25% 
## 641.3822

 quantile(inventory,0.50)

##      50% 
## 1083.294

 quantile(inventory,0.75)

##      75% 
## 1518.576

t.test(totalcost)

## 
##  One Sample t-test
## 
## data:  totalcost
## t = 231.37, df = 999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  148618.8 151161.4
## sample estimates:
## mean of x 
##  149890.1

t.test(inventory)

## 
##  One Sample t-test
## 
## data:  inventory
## t = 232.63, df = 13999, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1053.491 1071.395
## sample estimates:
## mean of x 
##  1062.443

hist(inventory,col="gold",breaks=nclass.Sturges,freq =TRUE,main="Inventory",xlab="inventory",ylab="Frequency")

hist(inventory,col="green",breaks=nclass.scott,freq =TRUE,main="Total cost of inventory",xlab="Cost",ylab="Frequency")

hist(inventory,col="lightyellow",breaks=20,freq =TRUE,main="Inventory",xlab="inventory",ylab="Frequency")

hist(inventory,col="lightgreen",breaks=nclass.scott,freq =TRUE,main="Inventory",xlab="inventory",ylab="Frequency")

plot.ecdf(inventory,main="Cumulative Distribution Function",xlab="inventory",ylab=TeX('$F_{I}(i)$'), col = "red")
abline(v = 0, lty = 2, col = "green")

hist(totalcost,col="gold",breaks=nclass.Sturges,freq =TRUE,main="Total cost of inventory",xlab="Total Cost",ylab="Frequency")

hist(totalcost,col="green",breaks=nclass.scott,freq =TRUE,main="Total cost of inventory",xlab="Cost",ylab="Frequency")

hist(totalcost,col="lightblue",breaks=10,freq =TRUE,main="Inventory cost",xlab="Cost",ylab="Frequency")