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

Question

在 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 列时，上述返回不可行。鉴于我在较小的示例中测试了各种场景，我的直觉是我的模型定义是好的......

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

Answer 1

这并不是您在 中所要求的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"))

Answer 2

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

这是一个草图。从任意有效的解决方案开始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 并返回一个稍微改变的版本（原始 的“邻居” 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 是默认值，但你肯定想玩弄这个数字。

• 可能存在目标函数不能提供良好指导的情况（即数据集）：如果选择不同的行不会改变满足目标的列数，则算法无法判断新的解决方案是否优于先前的解决方案一。因此，添加更连续的引导（例如，通过一些到目标的距离）可能会有所帮助。

• 更新：上面的计算实际上做了很多不必要的事情。当评估一个新的候选解决方案时，实际上不需要重新计算整个列的总和。相反，仅通过更改的行调整先前解决方案的总和。（对于小型数据集，这无关紧要。）

Answer 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]andtargets[j]大于或等于 0，那么我们可以定义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。但是引入了关于非负性的额外假设。

Answer 4

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

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

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