我建立了一个小的线性规划，用于优化作物产量，给出氮的产量响应曲线（具有对应于最大可达产量的平台的线性响应）。我使用ROI软件包在R软件中使用OPEN-OR CLP求解器编写了这个LP程序。
这是我简化问题的结构，仅给出了两种作物：
作物产量：Y1和Y2作物面积：单位面积A1和12作物氮吸收量：N1和N2总氮量：N_total产量响应曲线：Y1 = k1 * N1，Y2 = k2 * N2
以矩阵形式：
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
现在，我需要多次运行此优化（并给出61个作物），其中每次k，Y_max和N_total的值都会发生变化。我为每个参数构建了一个向量，每个向量包含大约90亿个值（原始数据存储在整个地球的栅格地图中，以5弧分的分辨率，我在向量中转换）。由于计算时间很长，我想将我的代码并行化，并且我想使用foreach包和并行运算符foreach来使用%dopar%循环。然而，并行化不起作用。
这里我的代码详细说明。我首先从我的数组开始构建约束矩阵（每个数组的维度dim = c（1,9231201,61）包含每个变量的值，然后我计算优化：
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]
}
有没有人知道为什么每个循环都不起作用（即它不产生任何结果而不是错误，当我停止R运行计数器时j没有改变它的值）？关于如何以不同方式并行运行此优化的任何建议？
据我所见，你永远不会分配解决方案，所以可能会改变
ROI_solve(prob, solver = "clp")
至
solution <- ROI_solve(prob, solver = "clp")