在 R 中的 ompr 包中，是什么导致了“问题太大”的错误？

Question

我正在尝试学习在 R 中使用 ompr 包来拟合线性整数规划优化模型，这是一位同事之前使用 CPLEX/GAMS 拟合的（具体来说，这里描述的那个：Haight et al. 2021）。我正在我大学的 Linux 超级计算服务器上运行我的实现，它有 248gb 的内存，我认为这足以完成这项工作。

这是我的代码和服务器故障报告的输出：

#Read in the necessary pre-generated data and packages

library(pacman); library(dplyr); library(ROI); library(ompr); library(ompr.roi)
n.ij = readRDS(file="nij1.rds") #An indexing vector.
B = 10 #Budget constraint--inspect only 10 lakes maximum

#Initialize model prior to setting the objective.
mod1 = MILPModel() %>% 
add_variable(u[i, j], type = "binary", i = 1:n.ij, j = 1:n.ij) %>%
add_variable(x[i], type = "binary", i = 1:n.ij) %>% 
add_variable(x[j], type = "binary", j = 1:n.ij) %>% 
add_constraint(x[i] + x[j] >= u[i,j], i = 1:n.ij, j = 1:n.ij) %>% 
add_constraint(sum_expr(x[i], i = 1:n.ij) <= B)
 
#Read in the relevant adjacency matrix of boat movements between every pair of lakes.
boats.n.ij = readRDS(file="boatsnij1.rds")

#Some system and object size info.
system(paste0("cat /proc/",Sys.getpid(),"/status | grep VmSize"))
VmSize: 13017708 kB
object.size(mod1)
6798778288 bytes

#Now, set objective with this specific boats.n.ij file.
mod1.full = mod1 %>% 
set_objective(sum_expr(u[i,j] * boats.n.ij[i, j], i = 1:n.ij, j = 1:n.ij))

Error in subCsp_ij(x, i, j, drop = drop) : 
  Cholmod error 'problem too large' at file ../Core/cholmod_sparse.c, line 89
Calls: %>% ... [ -> callGeneric -> eval -> eval -> [ -> [ -> subCsp_ij
Execution halted

为了创建可重现的示例，n.ij 和 boats.n.ij 的模拟版本可以按如下方式生成：

library(Matrix)

boats = rpois(7940*7940, 2)
keep = sample(c(0,1), 7940*7940, replace=T, prob = c(0.8, 0.2))
boat.dat = boats*keep

boats.n.ij = matrix(boat.dat, nrow=7940, ncol=7940)
diag(boats.n.ij) = 0
boats.n.ij = Matrix(boats.n.ij, sparse = T)

boats.n.ij[1:10, 1:10]

n.ij = 1:7940

为什么我无法将目标添加到我的模型中？只是我在暗示存在三个非常大的矩阵（决策矩阵u、boats.n.ij 矩阵和它们的乘积矩阵）？是因为模型已经是一个大约 6.8gb 的文件吗？我遇到的 R 是否对内存或对象大小施加了上限？这些功能是否无法考虑具有这么多决策点的目标？

我可以确认我已经能够在 boats.n.ij 的一个很小的子集上运行模型的缩小版本，优化得很好，所以我认为这不是我的模型规范的问题，但是我可能是错的......我还应该明确说明我对不涉及在 R 中解决这个模型的解决方案不感兴趣，因为这是这里的明确目标。但是，如果有更强大的包可用，我愿意使用其他包（尽管我喜欢这个包）。

注意：与我引用的论文不同，我不再需要我的同事使用的名为 b.ij 的向量，所以这不是这里的问题。

Answer 1

一次尝试：

require(slam)
boatsSTM<-as.simple_triplet_matrix(boats.n.ij)
...

#setting the objective function
set_objective(sum_expr(u[boatsSTM$i[k], boatsSTM$j[k]] * boatsSTM$v[k], k = 1:length(boatsSTM$i)))

我们利用矩阵的稀疏性。在一个简单的三元组矩阵中，您只需列出不为零的值，这意味着如果未列出元素则等于零。这些值用(i, j, v) 三元组表示，其中i 表示行索引，j 表示列索引，v 表示值。因此，例如，(2, 4, 10.32) 三元组表示m[2, 4] = 10.32。

在您的sum_expr 行中，我们利用了这一点，只添加不为零的元素。我们不会将u 的每个元素与boats 的每个元素相乘，因为大多数都是零并且与总和无关；相反，我们只是对重要的元素执行上述操作。

slam 包实现了简单的三元组矩阵，它的根是i、j 和v 值的列表。