Run multiple linear optimisation in parallel using R software

Question

I set up a small linear programming for optimizing crop yields given a yield response curve to Nitrogen (with a linear response with a plateau corresponding to the maximum reachable yield). I wrote this LP program using the OPEN-OR CLP solver in R software, by using the ROI package.

Here the structure of my simplified problem, given only two crops:

crop yields: Y1 and Y2 crops area: A1 and 12 crops nitrogen uptake per unit of area: N1 and N2 total nitrogen available: N_total yield response curve: Y1 = k1 * N1, Y2 = k2 * N2

In matrix form:

         variables      Y1      Y2       N1       N2
objective function      A1      A2        0        0
constraint matrix       1        0       -k1       0     =   0
                        0        1        0       -k2    =   0
                        1        0        0        0     <=  Y1_max
                        0        1        0        0     <=  Y2_max
                        0        0        A1       A2    <=  N_total  

Now, I need to run this optimization several time (and given 61 crops), where each time the values of k, Y_max, and N_total changes. I built a vector for each of these parameter, and each vector contains about 9 billion values (The original data were stored in raster maps of the entire globe with a 5 arc min resolution, that I transformed in vectors). Since the computational time is huge, I would like to parallelise my code, and I thought to use a foreach loop using the foreach package and the parallel operator %dopar%. Nevertheless, the parallelisation do not work.

Here my code in details. I first build up the constraints matrix starting from my arrays (each array has dimension dim=c(1,9231201,61) containing the values of each variable , and then I compute the optimization:

i <- 1


sparse_matrix_1 <- simple_triplet_diag_matrix(rep(1, 61), nrow=61)
sparse_matrix_2 <- simple_triplet_diag_matrix(c(matrix_yieldmax_Nmax[i,j,1],-matrix_yieldmax_Nmax[i,j,2],-matrix_yieldmax_Nmax[i,j,3],-matrix_yieldmax_Nmax[i,j,4],-matrix_yieldmax_Nmax[i,j,5],-matrix_yieldmax_Nmax[i,j,6],-matrix_yieldmax_Nmax[i,j,7],-matrix_yieldmax_Nmax[i,j,8],-matrix_yieldmax_Nmax[i,j,9],-matrix_yieldmax_Nmax[i,j,10],-matrix_yieldmax_Nmax[i,j,11],-matrix_yieldmax_Nmax[i,j,12],-matrix_yieldmax_Nmax[i,j,13],-matrix_yieldmax_Nmax[i,j,14],-matrix_yieldmax_Nmax[i,j,15],-matrix_yieldmax_Nmax[i,j,16],-matrix_yieldmax_Nmax[i,j,17],-matrix_yieldmax_Nmax[i,j,18],-matrix_yieldmax_Nmax[i,j,19],-matrix_yieldmax_Nmax[i,j,20],-matrix_yieldmax_Nmax[i,j,21],-matrix_yieldmax_Nmax[i,j,22],-matrix_yieldmax_Nmax[i,j,23],-matrix_yieldmax_Nmax[i,j,24],-matrix_yieldmax_Nmax[i,j,25],-matrix_yieldmax_Nmax[i,j,26],-matrix_yieldmax_Nmax[i,j,27],-matrix_yieldmax_Nmax[i,j,28],-matrix_yieldmax_Nmax[i,j,29],-matrix_yieldmax_Nmax[i,j,30],-matrix_yieldmax_Nmax[i,j,31],-matrix_yieldmax_Nmax[i,j,32],-matrix_yieldmax_Nmax[i,j,33],-matrix_yieldmax_Nmax[i,j,34],-matrix_yieldmax_Nmax[i,j,35],-matrix_yieldmax_Nmax[i,j,36],-matrix_yieldmax_Nmax[i,j,37],-matrix_yieldmax_Nmax[i,j,38],-matrix_yieldmax_Nmax[i,j,39],-matrix_yieldmax_Nmax[i,j,40],-matrix_yieldmax_Nmax[i,j,41],-matrix_yieldmax_Nmax[i,j,42],-matrix_yieldmax_Nmax[i,j,43],-matrix_yieldmax_Nmax[i,j,44],-matrix_yieldmax_Nmax[i,j,46],-matrix_yieldmax_Nmax[i,j,47],-matrix_yieldmax_Nmax[i,j,48],-matrix_yieldmax_Nmax[i,j,49],-matrix_yieldmax_Nmax[i,j,50],-matrix_yieldmax_Nmax[i,j,51],-matrix_yieldmax_Nmax[i,j,52],-matrix_yieldmax_Nmax[i,j,53],-matrix_yieldmax_Nmax[i,j,54],-matrix_yieldmax_Nmax[i,j,55],-matrix_yieldmax_Nmax[i,j,56],-matrix_yieldmax_Nmax[i,j,57],-matrix_yieldmax_Nmax[i,j,58],-matrix_yieldmax_Nmax[i,j,59],-matrix_yieldmax_Nmax[i,j,60],-matrix_yieldmax_Nmax[i,j,61]),nrow=61)
sparse_matrix_3 <- simple_triplet_diag_matrix(rep(1, 61), nrow=61)
sparse_matrix_4 <- simple_triplet_zero_matrix(61)
sparse_matrix_5 <- simple_triplet_matrix(rep(1, 61), seq(1:61), rep(0,61))
sparse_matrix_6 <- simple_triplet_matrix(rep(1, 61), seq(1:61), c(matrix_area_vec[i,j,1],matrix_area_vec[i,j,2],matrix_area_vec[i,j,3],matrix_area_vec[i,j,4],matrix_area_vec[i,j,5],matrix_area_vec[i,j,6],matrix_area_vec[i,j,7],matrix_area_vec[i,j,8],matrix_area_vec[i,j,9],matrix_area_vec[i,j,10],matrix_area_vec[i,j,11],matrix_area_vec[i,j,12],matrix_area_vec[i,j,13],matrix_area_vec[i,j,14],matrix_area_vec[i,j,15],matrix_area_vec[i,j,16],matrix_area_vec[i,j,17],matrix_area_vec[i,j,18],matrix_area_vec[i,j,19],matrix_area_vec[i,j,20],matrix_area_vec[i,j,21],matrix_area_vec[i,j,22],matrix_area_vec[i,j,23],matrix_area_vec[i,j,24],matrix_area_vec[i,j,25],matrix_area_vec[i,j,26],matrix_area_vec[i,j,27],matrix_area_vec[i,j,28],matrix_area_vec[i,j,29],matrix_area_vec[i,j,30],matrix_area_vec[i,j,31],matrix_area_vec[i,j,32],matrix_area_vec[i,j,33],matrix_area_vec[i,j,34],matrix_area_vec[i,j,35],matrix_area_vec[i,j,36],matrix_area_vec[i,j,37],matrix_area_vec[i,j,38],matrix_area_vec[i,j,39],matrix_area_vec[i,j,40],matrix_area_vec[i,j,41],matrix_area_vec[i,j,42],matrix_area_vec[i,j,43],matrix_area_vec[i,j,44],matrix_area_vec[i,j,45],matrix_area_vec[i,j,46],matrix_area_vec[i,j,47],matrix_area_vec[i,j,48],matrix_area_vec[i,j,49],matrix_area_vec[i,j,50],matrix_area_vec[i,j,51],matrix_area_vec[i,j,52],matrix_area_vec[i,j,53],matrix_area_vec[i,j,54],matrix_area_vec[i,j,55],matrix_area_vec[i,j,56],matrix_area_vec[i,j,57],matrix_area_vec[i,j,58],matrix_area_vec[i,j,59],matrix_area_vec[i,j,60],matrix_area_vec[i,j,61]))

top_sparse_matrix <- cbind(sparse_matrix_1, sparse_matrix_2)
middle_sparse_matrix <- cbind(sparse_matrix_3,sparse_matrix_4)
bottom_sparse_matrix <- cbind(sparse_matrix_5, sparse_matrix_6)
final_sparse_matrix <- rbind(top_sparse_matrix, middle_sparse_matrix, bottom_sparse_matrix)

dir <- c("==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" ,"==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "==" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" ,"<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=" , "<=")
rhs <-c(rep(0,61),matrix_yield_max_vec[i,j,1],matrix_yield_max_vec[i,j,2],matrix_yield_max_vec[i,j,3],matrix_yield_max_vec[i,j,4],matrix_yield_max_vec[i,j,5],matrix_yield_max_vec[i,j,6],matrix_yield_max_vec[i,j,7],matrix_yield_max_vec[i,j,8],matrix_yield_max_vec[i,j,9],matrix_yield_max_vec[i,j,10],matrix_yield_max_vec[i,j,11],matrix_yield_max_vec[i,j,12],matrix_yield_max_vec[i,j,13],matrix_yield_max_vec[i,j,14],matrix_yield_max_vec[i,j,15],matrix_yield_max_vec[i,j,16],matrix_yield_max_vec[i,j,17],matrix_yield_max_vec[i,j,18],matrix_yield_max_vec[i,j,19],matrix_yield_max_vec[i,j,20],matrix_yield_max_vec[i,j,21],matrix_yield_max_vec[i,j,22],matrix_yield_max_vec[i,j,23],matrix_yield_max_vec[i,j,24],matrix_yield_max_vec[i,j,25],matrix_yield_max_vec[i,j,26],matrix_yield_max_vec[i,j,27],matrix_yield_max_vec[i,j,28],matrix_yield_max_vec[i,j,29],matrix_yield_max_vec[i,j,30],matrix_yield_max_vec[i,j,31],matrix_yield_max_vec[i,j,32],matrix_yield_max_vec[i,j,33],matrix_yield_max_vec[i,j,34],matrix_yield_max_vec[i,j,35],matrix_yield_max_vec[i,j,36],matrix_yield_max_vec[i,j,37],matrix_yield_max_vec[i,j,38],matrix_yield_max_vec[i,j,39],matrix_yield_max_vec[i,j,40],matrix_yield_max_vec[i,j,41],matrix_yield_max_vec[i,j,42],matrix_yield_max_vec[i,j,43],matrix_yield_max_vec[i,j,44],matrix_yield_max_vec[i,j,45],matrix_yield_max_vec[i,j,46],matrix_yield_max_vec[i,j,47],matrix_yield_max_vec[i,j,48],matrix_yield_max_vec[i,j,49],matrix_yield_max_vec[i,j,50],matrix_yield_max_vec[i,j,51],matrix_yield_max_vec[i,j,52],matrix_yield_max_vec[i,j,53],matrix_yield_max_vec[i,j,54],matrix_yield_max_vec[i,j,55],matrix_yield_max_vec[i,j,56],matrix_yield_max_vec[i,j,57],matrix_yield_max_vec[i,j,58],matrix_yield_max_vec[i,j,59],matrix_yield_max_vec[i,j,60],matrix_yield_max_vec[i,j,61],matrix_N_tot[j])



foreach (j=1:9231201, .packages=c("ROI.plugin.clp", "ROI")) %dopar% {

    contst1 <- L_constraint(final_sparse_matrix, dir, rhs)   
    obj <- matrix_area_vec[1,j,]
    obj<- c(obj, rep(0,61))
    obj <- L_objective(obj, names = c("Y1" , "Y2" , "Y3" , "Y4" , "Y5" , "Y6" , "Y7" , "Y8" , "Y9" , "Y10" , "Y11" , "Y12" , "Y13" , "Y14" , "Y15" , "Y16" , "Y17" , "Y18" , "Y19" , "Y20" , "Y21" , "Y22" , "Y23" , "Y24" , "Y25" , "Y26" , "Y27" , "Y28" , "Y29" , "Y30" , "Y31" , "Y32" , "Y33" , "Y34" , "Y35" , "Y36" , "Y37" , "Y38" , "Y39" , "Y40" , "Y41" , "Y42" , "Y43" , "Y44" , "Y45" , "Y46" , "Y47" , "Y48" , "Y49" , "Y50" , "Y51" , "Y52" , "Y53" , "Y54" , "Y55" , "Y56" , "Y57" , "Y58" , "Y59" , "Y60" , "Y61" , "N1" , "N2" , "N3" , "N4" , "N5" , "N6" , "N7" , "N8" , "N9" , "N10" , "N11" , "N12" , "N13" , "N14" , "N15" , "N16" , "N17" , "N18" , "N19" , "N20" , "N21" , "N22" , "N23" , "N24" , "N25" , "N26" , "N27" , "N28" , "N29" , "N30" , "N31" , "N32" , "N33" , "N34" , "N35" , "N36" , "N37" , "N38" , "N39" , "N40" , "N41" , "N42" , "N43" , "N44" , "N45" , "N46" , "N47" , "N48" , "N49" , "N50" , "N51" , "N52" , "N53" , "N54" , "N55" , "N56" , "N57" , "N58" , "N59" , "N60" , "N61"))
    types <- rep("C", 61)

    prob <- OP(obj, final_const, types = types, maximum= TRUE)
    ROI_solve(prob, solver = "clp")

    matrix_final_yield[1,j,] <- solution$solution[1:3]
}

Does anybody know why the for each loop does not work (i.e. it does not produce any results not errors, and when I stop R from running the counter j did not change its value)? Any advice on how I could run this optimization in parallel in a different way?

Answer 1

As far as I can see you never assign the solution, so maybe change

ROI_solve(prob, solver = "clp")

to

solution <- ROI_solve(prob, solver = "clp")