opt-sim-exalpe.R

library(lpSolve)
#Defining parameters
q <- c(30,20,40,25)
w <- c(11,34,23)

muq <- c(0.15,0.55,0.3)
mur <- c(0.25,0.35,0.4)

yqtpr <- c(0.45,0.5,0.55,0.45)
yqtnr <- c(0.25,0.15, 0.35,0.5)

ywtpr <- c(0.25,0.35,0.45)
ywtnr <- c(0.05,0.35,0.05)

sq1 <- c(0.33, 0.66, 0.65, 0.49)
sq2 <- c(0.49, 0.36, 0.68, 0.58)
sq3 <- c(0.27, 0.57, 0.37, 0.48)

rw1 <- c(0.32, 0.28, 0.21)
rw2 <- c(0.29, 0.31, 0.53)
rw3 <- c(0.38, 0.48, 0.23)

#Costs
cq <- c(1,1,1,1)
cr <- c(1,1,1)

#Defining coefficients for objective function
obcof <- matrix(data = 0, nrow = 2, ncol = 21)
concof <- matrix(data = 0, nrow = 7, ncol = 7)
#obfpe <- w*(1-rw1)*(1-ywtnr) + w*(1-rw2)*(1-ywtnr) + w*(1-rw3)*(1-ywtnr)
obcof[1, c(1:3)] <- w*(1-rw1)*(1-ywtnr)*mur[1]
obcof[1,c(4:6)] <- w*(1-rw2)*(1-ywtnr)*mur[2]
obcof[1,c(7:9)] <- w*(1-rw3)*(1-ywtnr)*mur[3]
obcof[1,c(10:13)] <- q*(1-sq1)*(1-yqtnr)*muq[1]
obcof[1,c(14:17)] <- q*(1-sq2)*(1-yqtnr)*muq[2]
obcof[1,c(18:21)] <- q*(1-sq3)*(1-yqtnr)*muq[3]

obcof[2, c(1:3)] <- w*(1-rw1)*(1-ywtpr)*mur[1]
obcof[2,c(4:6)] <- w*(1-rw2)*(1-ywtpr)*mur[2]
obcof[2,c(7:9)] <- w*(1-rw3)*(1-ywtpr)*mur[3]
obcof[2,c(10:13)] <- q*(1-sq1)*(1-yqtpr)*muq[1]
obcof[2,c(14:17)] <- q*(1-sq2)*(1-yqtpr)*muq[2]
obcof[2,c(18:21)] <- q*(1-sq3)*(1-yqtpr)*muq[3]

#Defining constraints for the constraints
concof[1,c(1:3)] <- cr
concof[1,c(4:7)] <- cq

concof[2,c(1:3)] <- w*(rw1)*(ywtpr)*mur[1]
concof[2,c(4:7)] <- q*(sq1)*(yqtpr)*muq[1]

concof[3,c(1:3)] <- w*(rw2)*(ywtpr)*mur[2]
concof[3,c(4:7)] <- q*(sq2)*(yqtpr)*muq[2]

concof[4,c(1:3)] <- w*(rw3)*(ywtpr)*mur[3]
concof[4,c(4:7)] <- q*(sq3)*(yqtpr)*muq[3]

concof[5,c(1:3)] <- w*(1-rw1)*(ywtnr)*mur[1]
concof[5,c(4:7)] <- q*(1-sq1)*(yqtnr)*muq[1]

concof[6,c(1:3)] <- w*(1-rw2)*(ywtnr)*mur[2]
concof[6,c(4:7)] <- q*(1-sq2)*(yqtnr)*muq[2]

concof[7,c(1:3)] <- w*(1-rw3)*(ywtnr)*mur[3]
concof[7,c(4:7)] <- q*(1-sq3)*(yqtnr)*muq[3]

#Defining RHS
condir <- c("<=",">=",">=",">=","<=","<=","<=")
conrhs <- c(1000,2,2,2,300,300,300)

#Defining overall objective coefficient
ob.in1 <- c(sum(obcof[1,1], obcof[1,4], obcof[1,7]),sum(obcof[1,2], obcof[1,5], obcof[1,8]),sum(obcof[1,3], obcof[1,6], obcof[1,9]),
           sum(obcof[1,10], obcof[1,14], obcof[1,18]),sum(obcof[1,11], obcof[1,15], obcof[1,19]),sum(obcof[1,12], obcof[1,16], obcof[1,20]),
           sum(obcof[1,13], obcof[1,17], obcof[1,21]))
ob.in2 <- c(sum(obcof[2,1], obcof[2,4], obcof[2,7]),sum(obcof[2,2], obcof[2,5], obcof[2,8]),sum(obcof[2,3], obcof[2,6], obcof[2,9]),
            sum(obcof[2,10], obcof[2,14], obcof[2,18]),sum(obcof[2,11], obcof[2,15], obcof[2,19]),sum(obcof[2,12], obcof[2,16], obcof[2,20]),
            sum(obcof[2,13], obcof[2,17], obcof[2,21]))

ob.in <- ob.in1 + ob.in2

optimum <- lp(direction = "min", ob.in, concof, condir, conrhs, all.bin = TRUE)
optimum$solution

## [1] 0 0 0 0 0 1 0