RandomLPVer4.r

#===================================================================%
# A linear programming problem under uncertainty using Monte Carlo  %
#===================================================================%
#            Modeled & coded by Ebert Brea (c)2021                  %
#===================================================================%
require(lpSolve)

## Loading required package: lpSolve

library(latex2exp)
rm(list=ls()) 
if(file.exists("filltxt")){close(filtxt)}
#=================================================%
#         Parameters of the algorithm             %
#=================================================%

all.int=TRUE

data<- 
  data.frame("Products" = c("Cables(Und)", "Engines(Und)", "Pipelines(Und)", "Aluminum(Und)"), 
             "w_k" = c(3, 2, 0.5, 1), 
             "r_k" = c(8, 5, 1, 2), 
             "M_k"   = c(6, 4, 10, 10))

#set decision variables
objective.in <- c(8,5,1,2)

# constraint matrix
const.mat<- 
  matrix(c(3, 2, 0.5, 1,
           1, 0, 0, 0, 
           0, 1, 0, 0, 
           0, 0, 1, 0, 
           0, 0, 0, 1), 
         nrow = 5, byrow = TRUE)

#Direction for constraints
const.dir<- c("<=", "<=", "<=", "<=", "<=")

##############################################################
#        Simulation of "R" replications
##############################################################
R=2000 # Number of simulation replications %

# For defining matrix according to number of simulation replications

profit <-matrix( nrow=R, ncol = 1)
solution <- matrix( nrow=R, ncol = 9, byrow=TRUE)


# defining constraints
weight<-30


for(r in 1:R){
  ##############################################################
  # Starting the simulation replications
  ##############################################################
  # defining constraints
  
  solution[r,5]<-ceiling(runif(1,3,10)) #Cables
  solution[r,6]<-ceiling(runif(1,1,6))  #Engines
  solution[r,7]<-ceiling(runif(1,4,12)) #Pipelines
  solution[r,8]<-ceiling(runif(1,5,12)) #Aluminum
  
  # RHS for constraints
  
  const.rhs<- c(weight, solution[r,5], solution[r,6], solution[r,7], solution[r,8])
  
  opt=lp(direction = "max", objective.in , const.mat, const.dir, const.rhs, int.vec=1:4)$solution#
  
  solution[r,9]=0
  
  for(j in 1:4){solution[r,j]=opt[j]}
  for(j in 1:4){solution[r,9]=solution[r,9]+opt[j]*objective.in[j]}
  profit[r]=solution[r,9]
} 


summary(profit)

##        V1      
##  Min.   :59.0  
##  1st Qu.:75.0  
##  Median :78.0  
##  Mean   :76.7  
##  3rd Qu.:79.0  
##  Max.   :80.0

summary(solution[,1])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.000   5.000   7.000   6.918   8.000  10.000

summary(solution[,2])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   3.000   3.003   4.000   6.000

summary(solution[,3])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    2.00    3.54    6.00   12.00

summary(solution[,4])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   1.402   2.000  11.000

hist(profit,col="blue",breaks=nclass.Sturges,freq=TRUE,main=" ",xlab="profit",ylab="Frequency")#%

hist(solution[,1],col="green",breaks=nclass.Sturges,freq=TRUE,main=" ",xlab="Cables",ylab="Frequency")#%

hist(solution[,2],col="green",breaks=nclass.Sturges,freq=TRUE,main=" ",xlab="Engines",ylab="Frequency")#%

hist(solution[,3],col="green",breaks=nclass.Sturges,freq=TRUE,main=" ",xlab="Pipelines",ylab="Frequency")#%

hist(solution[,4],col="green",breaks=nclass.Sturges,freq=TRUE,main=" ",xlab="Aluminum",ylab="Frequency")#%

hist(solution[,9],col="green",breaks=20,freq=TRUE,main="Histogram",xlab="profit",ylab="Frequency")#%

plot.ecdf(profit,main="Cumulative Distribution Function",xlab="profit",ylab="F(i)", col = "red")#%