RALP

library(lpSolveAPI)

# Input
n <- 3   # n is cross-section
d <- 5    # number of worst case losses (based on the marginal; to be matched - dependence is unknown )

# simple example
# Create matrix with identical rows, filled with inverse probabilities
data <- matrix(1 : d, 
               ncol = n, 
               nrow = d, 
               byrow = FALSE)
x <- as.numeric(data)

# we define the d*(d*n) decision variable of 0/1 such that sum( solution[1:(d*n)]*x ) = sum for the first period

# Objective function
# S1 
# Repeat the matrix d times
# Dimension [1 : d * (n * d) ]
obj <- c(x , rep(0, (d - 1) * n * d)) 

# Type of objective variables
types <- rep("B", length(obj))


# Constrain coefficients
# n * d rows to select one element per column 
#  -eps <= S1-S2, S2-S3, S3-S4,.. <= + eps : 2*(d-1) comparisons
mat <- matrix(0,  
              ncol = length(obj),
              nrow = n * d +  (d - 1),
              byrow = TRUE)

# Select each element exactly once
for (row in 1 : (n * d)) {
  for( i in 1 : d){
    mat[row,(i - 1) * (n * d) + row] = 1   
  }
}

for (row in 1 : (d - 1)) {
  mat[n * d + row, (row - 1) * n * d + c(1 : (2 * (n * d)))] = c(x, -x)
}

# Dir
dir <- c(rep("==", n * d), rep(c("<="), d - 1) )

# Right hand side
eps = 0
rhs <- c(rep(1, n * d), rep(eps, d - 1))

lp <- make.lp(n * d + d - 1, n * d ^ 2)

# Set objective function
set.objfn(lp, obj)
# 1:(d*n) komt overeen met S1

# Loop over all n * d * d columns
for(j in 1 : ncol(mat)) {
  set.column(lp, j, mat[, j])
}

# Set direction
set.constr.type(lp, dir)

# Set variable type to binary
set.type(lp, 1:length(obj), "binary")

# Set right hand side
set.rhs(lp, rhs)

lp.control(lp, sense = "max")$sense

solve(lp)
get.primal.solution(lp)

write.lp(lp, "lpfilename.lp", "lp")