在 R 的 ompr 包中，如何重新表述我的 objective/constraints/variables 以避免 "problem too large" 错误？
In the ompr package in R, how can I rephrase my objective/constraints/variables so as to avoid the "problem too large" error?
我正在尝试学习使用 ompr 软件包在 R 中拟合线性整数规划优化模型，该模型是一位同事之前使用 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
为什么我无法将 objective 添加到我的模型中？是否只是我在暗示存在三个非常大的矩阵（决策矩阵 u、boats.n.ij 矩阵及其乘积矩阵）？是不是因为模型已经一个6.8gb左右的文件？我正在 运行 加入的 R 是否对内存或 object 大小设置了上限？这些功能是否不能考虑具有这么多决策点的 objective？
我可以确认我已经能够 运行 模型的 scaled-down 版本在 boats.n.ij 的一个非常小的子集上优化得很好，所以我不认为这是我的模型规范的问题，但我可能是错的......我还应该明确声明我对不涉及在 R 中解决此模型的解决方案不感兴趣，因为这里是表达式 objective。但是，如果有更强大的可用软件包，我愿意使用其他软件包（尽管我喜欢其他软件包）。
注意：与我引用的论文不同，我已经消除了对我同事使用的名为 b.ij 的向量的需求，所以这不是这里的问题。
编辑：请注意，@nicola 对 objective 的改造将设置并解决，但原始约束 and/or 变量将不再与它具有相同的关系，因此它是合适的与我想要安装的型号不同的型号。在最初的构造中，x[i] 中最多只有 10 个值，因此决策变量 u[i,j] 中 i 的最多 10 个唯一值将被允许为 1s，这要归功于涉及我们的约束预算参数 B。在@nicola 的版本中，u[i,j] 中允许 i 的 10 个以上的唯一值是 1。实际上，至少我不清楚最初编写的约束如何与@nicola 的 objective 交互，如果有的话。但是，我怀疑像@nicola's 这样的 objective 肯定可以用来利用我的 boats.n.ij 矩阵的稀疏性以避免“问题太大”错误，但它需要变量 and/constraints 进行相应修改。我更改了问题的标题，以便更清楚地了解我在寻找什么——我想避免错误 ，但在其他方面适合等效模型。
第二次编辑：@nicola 的解决方案毕竟有效！但是，自从我发布此问题以来，由于 ompr 的更新，变量和约束需要进行一些修改。请参阅以下玩具示例：
library(Matrix)
library(slam)
library(dplyr)
library(tidyr)
library(ROI)
library(ompr)
library(ompr.roi)
library(Rglpk)
library(ROI.plugin.glpk)
library(lattice)
set.seed(101)
N = 500
boats = rpois(N*N, 2)
keep = sample(c(0,1), N*N, replace=T, prob = c(0.97, 0.03))
boat.dat = boats*keep
boats.n.ij = Matrix(boat.dat, nrow=N, ncol=N, sparse =T)
diag(boats.n.ij) = 0
boats.n.ij[1:10, 1:10]
n.ij = N
B = 5
mod1 = MIPModel() %>%
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(y[j], type = "binary", j = 1:n.ij) %>%
add_constraint(x[i] == y[j], i = 1:n.ij, j = 1:n.ij, i == j) %>%
add_constraint(sum_over(x[i], i = 1:n.ij) <= B) %>%
add_constraint(u[i,j] <= x[i] + y[j], i = 1:n.ij, j = 1:n.ij)
boatsSTM = as.simple_triplet_matrix(boats.n.ij)
#setting the objective function
mod.2nd = mod1 %>% set_objective(sum_over(u[boatsSTM$i[k], boatsSTM$j[k]] * boatsSTM$v[k], k = 1:length(boatsSTM$i)))
mod.2nd.solved = mod.2nd %>%
solve_model(with_ROI("glpk", verbose=TRUE))
testB = get_solution(mod.2nd.solved, u[i,j])
test2B = pivot_wider(testB, names_from = j, values_from = value) %>% dplyr::select(-variable, -i)
test3B = as.matrix(test2B, nrow=100)
levelplot(test3B)
一次尝试：
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 值的列表。
我正在尝试学习使用 ompr 软件包在 R 中拟合线性整数规划优化模型，该模型是一位同事之前使用 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
为什么我无法将 objective 添加到我的模型中？是否只是我在暗示存在三个非常大的矩阵（决策矩阵 u、boats.n.ij 矩阵及其乘积矩阵）？是不是因为模型已经一个6.8gb左右的文件？我正在 运行 加入的 R 是否对内存或 object 大小设置了上限？这些功能是否不能考虑具有这么多决策点的 objective？
我可以确认我已经能够 运行 模型的 scaled-down 版本在 boats.n.ij 的一个非常小的子集上优化得很好，所以我不认为这是我的模型规范的问题，但我可能是错的......我还应该明确声明我对不涉及在 R 中解决此模型的解决方案不感兴趣，因为这里是表达式 objective。但是，如果有更强大的可用软件包，我愿意使用其他软件包（尽管我喜欢其他软件包）。
注意：与我引用的论文不同，我已经消除了对我同事使用的名为 b.ij 的向量的需求，所以这不是这里的问题。
编辑：请注意，@nicola 对 objective 的改造将设置并解决，但原始约束 and/or 变量将不再与它具有相同的关系，因此它是合适的与我想要安装的型号不同的型号。在最初的构造中，x[i] 中最多只有 10 个值，因此决策变量 u[i,j] 中 i 的最多 10 个唯一值将被允许为 1s，这要归功于涉及我们的约束预算参数 B。在@nicola 的版本中，u[i,j] 中允许 i 的 10 个以上的唯一值是 1。实际上，至少我不清楚最初编写的约束如何与@nicola 的 objective 交互，如果有的话。但是，我怀疑像@nicola's 这样的 objective 肯定可以用来利用我的 boats.n.ij 矩阵的稀疏性以避免“问题太大”错误，但它需要变量 and/constraints 进行相应修改。我更改了问题的标题，以便更清楚地了解我在寻找什么——我想避免错误 ，但在其他方面适合等效模型。
第二次编辑：@nicola 的解决方案毕竟有效！但是，自从我发布此问题以来，由于 ompr 的更新，变量和约束需要进行一些修改。请参阅以下玩具示例：
library(Matrix)
library(slam)
library(dplyr)
library(tidyr)
library(ROI)
library(ompr)
library(ompr.roi)
library(Rglpk)
library(ROI.plugin.glpk)
library(lattice)
set.seed(101)
N = 500
boats = rpois(N*N, 2)
keep = sample(c(0,1), N*N, replace=T, prob = c(0.97, 0.03))
boat.dat = boats*keep
boats.n.ij = Matrix(boat.dat, nrow=N, ncol=N, sparse =T)
diag(boats.n.ij) = 0
boats.n.ij[1:10, 1:10]
n.ij = N
B = 5
mod1 = MIPModel() %>%
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(y[j], type = "binary", j = 1:n.ij) %>%
add_constraint(x[i] == y[j], i = 1:n.ij, j = 1:n.ij, i == j) %>%
add_constraint(sum_over(x[i], i = 1:n.ij) <= B) %>%
add_constraint(u[i,j] <= x[i] + y[j], i = 1:n.ij, j = 1:n.ij)
boatsSTM = as.simple_triplet_matrix(boats.n.ij)
#setting the objective function
mod.2nd = mod1 %>% set_objective(sum_over(u[boatsSTM$i[k], boatsSTM$j[k]] * boatsSTM$v[k], k = 1:length(boatsSTM$i)))
mod.2nd.solved = mod.2nd %>%
solve_model(with_ROI("glpk", verbose=TRUE))
testB = get_solution(mod.2nd.solved, u[i,j])
test2B = pivot_wider(testB, names_from = j, values_from = value) %>% dplyr::select(-variable, -i)
test3B = as.matrix(test2B, nrow=100)
levelplot(test3B)
一次尝试：
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 值的列表。