在R中的ompr包中，我如何重新表述我的目标/约束/变量，以避免问题太大？ [英] In the ompr package in R, how can I rephrase my objective/constraints/variables so as to avoid the "problem too large" error?

问题描述

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

以下是来自服务器的故障报告的代码和输出：

#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.8 GB的文件？我遇到的R是否对内存或对象大小设置了上限？难道这些函数不能考虑有这么多决策点的目标吗？

我可以确认，我已经能够在boats.n.ij的一个非常小的子集上运行模型的缩小版本，它优化得很好，所以我不认为这是我的模型规范的问题，但我可能是错的……我还应该明确地指出，我对不涉及用R求解该模型的解决方案不感兴趣，因为这是这里明确的目标。但是，如果有更强大的包可用，我对使用其他包持开放态度(尽管我不喜欢这个包)。

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

编辑：请注意，@Nicola对目标的改革将设置并解决问题，但原始的约束和/或变量将不再与其具有相同的关系，因此它将适合不同于我想要适应的模型。在最初的构造中，由于涉及预算参数B的限制，x[i]中最多10个值，因此决策变量u[i，j]中最多10个唯一的i值将被允许为1。在@Nicola的版本中，在u[i，j]中允许超过10个唯一的i值为1。实际上，至少对我来说，最初编写的约束与@Nicola的目标是如何相互作用的，如果有的话，我是不清楚的。然而，我怀疑像@Nicola这样的目标肯定可以用来利用我的boats.n.ij矩阵的稀疏性，以避免问题太大&错误，但它需要相应地修改变量和/约束。我更改了问题的标题，以便更清楚地说明我要寻找的是什么--我希望避免错误，但在其他方面适合等价的模型。

第二次编辑：@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值的列表。

这篇关于在R中的ompr包中，我如何重新表述我的目标/约束/变量，以避免问题太大？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！