r - 有没有办法在 R 优化器中定义复杂的目标函数？

Question

在 R 中，我正在尝试优化以下内容:选择行来最大化其总和超过特定值的列数，该值因列而异+行选择的一些其他基本约束。

R 中是否有任何东西可以让您将逻辑合并到目标函数中？即在约 10k 列上最大化 countif( sum(value column) > 列的目标值)，选择 5 行，并选择约 500 行。

简单示例:抓取下面 4 行的组合，其列总和超过目标的频率比任何其他 4 行组合更频繁。

  +--------+------+------+------+------+------+------+------+------+------+-------+
    |   x    | col1 | col2 | col3 | col4 | col5 | col6 | col7 | col8 | col9 | col10 |
    +--------+------+------+------+------+------+------+------+------+------+-------+
    | row1   |   82 |   73 |   50 |   11 |   76 |   12 |   46 |   64 |    5 |    44 |
    | row2   |    2 |   33 |   35 |   55 |   52 |   18 |   13 |   86 |   72 |    39 |
    | row3   |   94 |    5 |   10 |   21 |   90 |   62 |   54 |   54 |    7 |    17 |
    | row4   |   27 |   10 |   28 |   87 |   27 |   83 |   62 |   56 |   54 |    86 |
    | row5   |   17 |   50 |   34 |   30 |   80 |    7 |   96 |   91 |   32 |    21 |
    | row6   |   73 |   75 |   32 |   71 |   37 |    1 |   13 |   76 |   10 |    34 |
    | row7   |   98 |   13 |   87 |   49 |   27 |   90 |   28 |   75 |   55 |    21 |
    | row8   |   45 |   54 |   25 |    1 |    3 |   75 |   84 |   76 |    9 |    87 |
    | row9   |   40 |   87 |   44 |   20 |   97 |   28 |   88 |   14 |   66 |    77 |
    | row10  |   18 |   28 |   21 |   35 |   22 |    9 |   37 |   58 |   82 |    97 |
    | target |  200 |  100 |  125 |  135|  250 |  89 |  109 |  210|  184 |   178 |
    +--------+------+------+------+------+------+------+------+------+------+-------+

编辑 + 更新:我使用 ompr、ROI 和一些 Big M 逻辑实现了以下内容。

nr <- 10 # number of rows
nt <- 15 # number of target columns
vals <- matrix(sample.int(nr*nt, nr*nt), nrow=nr, ncol=nt)

targets <- vector(length=nt)
targets[1:nt] <- 4*mean(vals)


model <- MIPModel() %>%
  add_variable(x[i], i = 1:nr, type = "binary") %>%
  add_constraint(sum_expr(x[i], i = 1:nr)==4)%>%
  add_variable(A[j], j = 1:nt, type = "binary") %>%
  add_variable(s[j], j = 1:nt, type = "continuous",lb=0) %>%
  add_constraint(s[j] <= 9999999*A[j], j =1:nt)%>%
  add_constraint(s[j] >= A[j], j =1:nt)%>%
  add_constraint(sum_expr(vals[i,j]*x[i], i = 1:nr) + A[j] + s[j] >= targets[j], j=1:nt) %>%    
    set_objective(sum_expr(-9999999*A[j], i = 1:nr, j = 1:nt), "max")

model <- solve_model(model,with_ROI(solver = "glpk"))

该模型非常适合解决小问题，包括那些不存在超出每列目标的解决方案的问题。

但是，当我将列数更改为 150 时，上面的结果会返回 Infeasible。鉴于我在较小的示例上测试了各种场景，我的预感是我的模型定义没问题...

关于为什么这是不可行的任何建议？或者也许是定义模型的更优化方法？

Answer 1

您可以尝试本地搜索算法。它可能只给你一个“好的”解决方案；但作为交换，它非常灵活。

这是一个草图。例如，从任意有效的解决方案 x 开始对于您的示例数据

x <- c(rep(TRUE, 4), rep(FALSE, 6))
## [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE

定义目标函数:

obj_fun <- function(x, table, target, ...) {
    -sum(colSums(table[x, ]) >= target)
}

给定一个表和一个目标向量，它选择行在x中定义并计算行数达到或超过目标。我写-sum因为我将使用一个最小化的实现目标函数。

-obj_fun(x, table, target)
## [1] 7

因此，对于所选的初始解决方案，7 列总和等于或大于目标。

那么你需要一个邻域函数。这需要一个解决方案 x 并返回稍有更改的版本(a原始x的“邻居”)。这是一个邻居函数更改 x 中的一行。

nb <- function(x, ...) {
    true  <- which( x)
    false <- which(!x)
  
    i <-  true[sample.int(length( true), size = 1)]
    j <- false[sample.int(length(false), size = 1)]
    x[i] <- FALSE
    x[j] <- TRUE
    x
}


x
## [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE

nb(x)
## [1] FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE
##     ^^^^^                                      ^^^^

这是您的数据:

library("orgutils")
tt <- readOrg(text = "
    |   x    | col1 | col2 | col3 | col4 | col5 | col6 | col7 | col8 | col9 | col10 |
    |--------+------+------+------+------+------+------+------+------+------+-------+
    | row1   |   82 |   73 |   50 |   11 |   76 |   12 |   46 |   64 |    5 |    44 |
    | row2   |    2 |   33 |   35 |   55 |   52 |   18 |   13 |   86 |   72 |    39 |
    | row3   |   94 |    5 |   10 |   21 |   90 |   62 |   54 |   54 |    7 |    17 |
    | row4   |   27 |   10 |   28 |   87 |   27 |   83 |   62 |   56 |   54 |    86 |
    | row5   |   17 |   50 |   34 |   30 |   80 |    7 |   96 |   91 |   32 |    21 |
    | row6   |   73 |   75 |   32 |   71 |   37 |    1 |   13 |   76 |   10 |    34 |
    | row7   |   98 |   13 |   87 |   49 |   27 |   90 |   28 |   75 |   55 |    21 |
    | row8   |   45 |   54 |   25 |    1 |    3 |   75 |   84 |   76 |    9 |    87 |
    | row9   |   40 |   87 |   44 |   20 |   97 |   28 |   88 |   14 |   66 |    77 |
    | row10  |   18 |   28 |   21 |   35 |   22 |    9 |   37 |   58 |   82 |    97 |
    | target |  200 |  100 |  125 |   135|  250  |  89 |  109 |   210|  184 |   178 |
")


table  <- tt[1:10, -1]
target <- tt[11,   -1]

运行搜索；在这种情况下，使用称为“阈值接受”。我使用包 NMOF 中的实现(我坚持这一点)。

library("NMOF")
x0 <- c(rep(TRUE, 4), rep(FALSE, 6))
sol <- TAopt(obj_fun,
             list(neighbour = nb,     ## neighbourhood fun
          x0 = sample(x0),    ## initial solution
          nI = 1000,          ## iterations
                  OF.target = -ncol(target)  ## when to stop
                 ),
             target = target,
             table = as.matrix(table))

rbind(Sums = colSums(table[sol$xbest, ]), Target = target)       
##        col1 col2 col3 col4 col5 col6 col7 col8 col9 col10
## Sums    222  206  216  135  252  148  175  239  198   181
## Target  200  100  125  135  250   89  109  210  184   178

正如我所说，这只是一个草图，具体取决于如何你的实际问题是大而重要的，有很多需要考虑的要点:

• 最重要的是:nI设置搜索次数迭代。 1000 是默认值，但您肯定会想玩弄这个数字。

• 在某些情况下(即数据集)，目标函数没有提供良好的指导:如果选择不同的行不会改变数字满足目标的列数，算法无法判断新的解决方案是否比现有的更好前一个。因此，增加更连续的指导(例如，通过到目标的一定距离)可能会有所帮助。

• 更新:上面的计算实际上做了很多事情那没有必要。当新的候选解出现时进行了评估，实际上没有必要重新计算全列总和。相反，只需调整先前解决方案的总和已更改行。 (对于小型数据集，这并不重要。)