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

【问题标题】：Is there a way to define a complex objective function in an R optimizer?有没有办法在 R 优化器中定义复杂的目标函数？
【发布时间】：2021-10-29 08:37:34
【问题描述】：

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

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

简单示例：抓取下面 4 行的组合，其 col 总和比任何其他 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 和一些大 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 列时，上述结果返回不可行。鉴于我在较小的示例中测试了各种场景，我的直觉是我的模型定义没问题...

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

【问题讨论】：

有趣的问题！一个最小的例子会很有帮助。
添加了一个例子。让我知道这是否有帮助或不清楚。
也许我理解错了，但在示例中，似乎没有列总和可以超过其目标？
哈哈，这可能是真的！我完全随机生成数据并且没有看。编辑以改变这一点
抱歉，这不是真的。我认为我的解释清楚地说明了我要解决的问题，您列出的内容不是找到最佳行组合的解决方案。你的建议，如果我解释正确的话，将意味着生成所有行组合并从那里开始。虽然对于只有 10 行的示例来说这是可能的，但对于我的实际问题，大约 500 行（即 500*499*498*497*496 可能的组合）是不切实际的。

标签： r optimization ompr

【解决方案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 是默认值，但你肯定会想玩弄这个号码。
在某些情况下（即数据集）可能需要目标函数不能提供良好的指导：如果选择不同的行不会改变数字在达到目标的列中，算法无法判断新解决方案是否比现有解决方案更好前一个。因此，添加更连续的指导（例如，通过一些与目标的距离）可能会有所帮助。
更新：上面的计算实际上做了很多那没有必要。当一个新的候选解决方案被评估，实际上没有必要重新计算全列总和。相反，只调整先前解决方案的总和由更改行。（对于小型数据集，这无关紧要。）

【讨论】：

【解决方案2】：

这并不是您在python 中提出的要求，但也许它会向您展示使用整数编程执行此操作的方法。您应该能够在 R 中复制这一点，因为 R 中有多个求解器的绑定，包括 CBC，这是我在下面使用的，它适用于整数程序。

我还使用pyomo 来构建求解器的数学模型。我认为通过一些研究，您可以在 R 中找到等效的方法。一开始的语法只是摄取数据（我只是将其粘贴到 .csv 文件中）。其余的应该是可读的。

好/坏...

这几乎可以立即解决您的玩具问题。可以证明5行可以超过所有列的总数。

对于更多的列，它可能会大大陷入困境。我用大量随机数矩阵进行了几次测试......这对求解器来说非常具有挑战性，因为它无法轻松识别“好”行。通过放宽解决方案的容差，我可以让它在合理的时间内用随机值（以及随机的总行并乘以 5（选择的数量......只是为了使其具有挑战性）来解决 500x100。

如果您真的有 10K 列，那么只有几种方法可以工作... 1. 您有几行可以覆盖所有列总数（求解器应该很快发现这一点）或 2. 有一些模式（除了随机噪声）到可以指导求解器的数据/总数，以及 3. 使用基于大比率的间隙（或时间限制）

import pyomo.environ as pyo
import pandas as pd
import numpy as np

df = pd.read_csv("data.csv", header=None)  # this is the data from the post

# uncomment this below for a randomized set of data
# df = pd.DataFrame(
#     data = np.random.random(size=(500,100)))
# df.iloc[-1] = df.iloc[-1]*5

# convert to dictionary
data = df.iloc[:len(df)-1].stack().to_dict()
col_sums = df.iloc[len(df)-1].to_dict()

limit = 5  # max number or rows selected

m = pyo.ConcreteModel('row picker')

### SETS
m.R = pyo.Set(initialize=range(len(df)-1))
m.C = pyo.Set(initialize=range(len(df.columns)))

### Params
m.val = pyo.Param(m.R, m.C, initialize=data)
m.tots = pyo.Param(m.C, initialize=col_sums)

### Variables
m.sel = pyo.Var(m.R, domain=pyo.Binary)  # indicator for which rows are selected
m.abv = pyo.Var(m.C, domain=pyo.Binary)  # indicator for which column is above total

### OBJECTIVE
m.obj = pyo.Objective(expr=sum(m.abv[c] for c in m.C), sense=pyo.maximize)

### CONSTRAINTS
# limit the total number of selections...
m.sel_limit = pyo.Constraint(expr=sum(m.sel[r] for r in m.R) <= limit)

# link the indicator variable to the column sum 
def c_sum(m, c):
    return sum(m.val[r, c] * m.sel[r] for r in m.R) >= m.tots[c] * m.abv[c]
m.col_sum = pyo.Constraint(m.C, rule=c_sum)

### SOLVE
print("...built... solving...")
solver = pyo.SolverFactory('cbc', options={'ratio': 0.05})
result = solver.solve(m)
print(result)

### Inspect answer ...
print("rows to select: ")
for r in m.R:
    if m.sel[r]:
        print(r, end=', ')

print("\ncolumn sums from those rows")
tots = [sum(m.val[r,c]*m.sel[r].value for r in m.R) for c in m.C]
print(tots)
print(f'percentage of column totals exceeded:  {len([1 for c in m.C if m.abv[c]])/len(m.C)*100:0.2f}%')

产量：

Problem: 
- Name: unknown
  Lower bound: -10.0
  Upper bound: -10.0
  Number of objectives: 1
  Number of constraints: 11
  Number of variables: 20
  Number of binary variables: 20
  Number of integer variables: 20
  Number of nonzeros: 10
  Sense: maximize
Solver: 
- Status: ok
  User time: -1.0
  System time: 0.0
  Wallclock time: 0.0
  Termination condition: optimal
  Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
    Black box: 
      Number of iterations: 0
  Error rc: 0
  Time: 0.013128995895385742
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

rows to select: 
0, 2, 3, 8, 9, 
column sums from those rows
[261.0, 203.0, 153.0, 174.0, 312.0, 194.0, 287.0, 246.0, 214.0, 321.0]
percentage of column totals exceeded:  100.00%
[Finished in 845ms]

编辑：

我看到您的编辑遵循与上述解决方案类似的模式。

对于较大的实例化，您获得“不可行”的原因是，当值更大并且相加更多时，您的 Big-M 不再足够大。您应该预先分析您的矩阵并将BIG_M 设置为目标行中的最大值，这将足以覆盖任何间隙（通过检查）。这将使您在BIG_M 上不会出现大量超调，这也会产生后果。

我对您的 r 模型进行了一些调整。我的r 语法很糟糕，但试试这个：

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] <= BIG_M*A[j], j =1:nt)%>%
  # NOT NEEDED:  add_constraint(s[j] >= A[j], j =1:nt)%>%
  # DON'T include A[j]:  add_constraint(sum_expr(vals[i,j]*x[i], i = 1:nr) + A[j] + s[j] >= targets[j], j=1:nt) %>%   
  add_constraint(sum_expr(vals[i,j]*x[i], i = 1:nr) + s[j] >= targets[j], j=1:nt) %>%  
  # REMOVE unneded indexing for i:  set_objective(sum_expr(A[j], i = 1:nr, j = 1:nt), "min")
  # and just minimize.  No need to multiply by a large constant here.
  set_objective(sum_expr(A[j], j = 1:nt), "min")

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

【讨论】：

【解决方案3】：

恕我直言，这是一个线性规划建模问题：我们能否将问题表述为“归一化”线性问题，例如可以通过 ompr 或 ROI 解决（我会添加 lpSolveAPI）？

我相信这是可能的，尽管我没有时间提供完整的表述。以下是一些想法：

作为参数，即固定值，我们有

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

我们感兴趣的决策变量是 x[1...nr] 作为二进制变量（如果选择行，则为 1，否则为 0）。

显然，一个约束是sum(x[i],i)==4——我们选择的行数。

为了目标，我会引入辅助变量，例如

y[j] = 1, if sum_{i=1..nr} x[i]*vals[i,j]>= targets[j]

（否则为 0）用于j=1...nt。现在y这个定义不兼容线性规划，需要线性化。如果我们可以假设val[i,j] 和targets[j] 大于或等于零，那么我们可以将y[j] 定义为二进制变量，如下所示：

x'vals[,j]-t[j]*y[j] >= 0

(x'y 表示内积，即sum(x[i]*y[i], i)。) 在x'vals[,j]>=t[j] 的情况下，值y[j]==1 是有效的。在x'vals[,j]<t[j] 的情况下，强制执行y[j]==0。

对于目标max sum(y[j],j)，我们应该得到问题的正确表述。不需要大M。但是引入了关于非负性的额外假设。

【讨论】：

您好，您对 Big-M 的建议非常有帮助。我实现了一些对我的小例子很有效的东西，但在应用于更大的集合时却不可行。我正在用我的代码编辑我的初始帖子。
编辑了更多想法的答案。

【解决方案4】：

您在这里要解决的问题称为“混合整数程序”，并且围绕它设计了很多（主要是商业）软件。

您的典型 R 函数（例如 optim）由于这种限制几乎没有任何好处，但您可以使用专门的软件（例如 CBC），只要您能够在标准 MIP 中构建问题结构（在这种情况下，要优化的变量是数据中每一行的二进制变量）。

作为替代方案，您还可以查看包 nloptr 及其全局无衍生黑盒优化器，您可以在其中输入这样的函数（设置变量边界）并让它优化它一些通用的启发式方法。

【讨论】：

感谢您的回复。看着 CBC 和 nloptr，我仍在努力定义我的目标函数。我知道我的决策变量是代表每行包含的一组二进制变量。但是，我正在寻求有关如何定义“最大化 countif ( sum(value column) > target value for column ) * (vector of binary decision variables)”的指导。正如最初提到的，我可以定义“最大化列总和”，但我不知道如何将其更改为“将列总和计数为 1，如果它 > 那列的目标 val，0 如果
是的，这很难以标准形式进行编码，并且您可能必须根据每列的行值总和为列/阈值加上不等式约束引入额外的二进制变量（在顶部在我的脑海中，我不能 100% 确定是否可以将其构图为简单的 MIP）。您还可以查看or-tools for python，它允许以不太标准的方式构建问题并在内部或在 R 的 ROI 上进行翻译。