0

我想使用具有分段常数目标函数的 R 来解决最小化问题。这个想法是，对于我的（整数）决策变量x的较低值，会产生比较高值更高的惩罚成本。考虑到一些限制，我想尽量减少总罚款成本。

所以，我的程序如下所示：

min   P(x)
s.t.  A x <= b
        x >= 0

我在编程分段常数目标函数P(x)时遇到问题，其中P(x)是向量x的所有元素的总和。我知道它不能与库中的函数结合lp()使用linprog。但是，如果不指定大量额外的变量，我无法找到一种方法。此外，广泛的互联网搜索并没有提供任何有用的想法。

让我举一个例子来说明这个函数P的样子

[,1] [,2] [,3] [,4] [,5] [,6]
  21   11    9    9    0    0
  45   28   17   17    6    0
  14    0    0    0    0    0
  17   11   11    5    5    0
  26   11    0    0    0    0
  38   18   18   13    0    0

这应该按以下方式阅读：如果x1=2，则会产生 11 的惩罚成本。如果x6=4，则会产生 13 的惩罚成本。也就是说，对于x=c(2, ..., 4)，我们有 ，P=c(11, ..., 13)总惩罚成本（目标值）是sum(11, ..., 13)。

我的矩阵A（它完全是单模的）和向量b如下所示。

A <- matrix(c(1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1),nrow=6) b <- c(4,5,1,5,2,4).

所以，我的问题是：

如何使用 R 找到分段常数目标函数的最小值？

4

1 回答 1

0

对于您的问题，我有两种方法。您的问题似乎是整数编程，其中 x1 ... x6 可能取 {1, 2, 3, 4, 5, 6} 中的值。

Pmat <- matrix(c(25, 11, 9, 9, 0, 0,
                 45, 28, 17, 17, 6, 0,
                 14, 0, 0, 0, 0, 0,
                 17, 11, 11, 5, 5, 0,
                 26, 11, 0, 0, 0, 0,
                 38, 18, 18, 13, 0, 0),
                byrow = TRUE, nrow = 6, ncol = 6)
A <- matrix(c(1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1),nrow=6)
b <- c(4,5,1,5,2,4)

蛮力

由于您有 6 个这样的变量，因此只有 6^6 = 46656 个候选者。因此，用蛮力解决这个问题不会花费那么长时间；对于每个组合，检查它是否满足约束和值。

d <- expand.grid(x1 = 1:6, x2 = 1:6, x3 = 1:6, x4 = 1:6, x5 =1:6, x6 = 1:6)
condition <- logical(nrow(d))
value <- numeric(nrow(d))
for (i in seq(nrow(d))) {
  x <- as.matrix(as.numeric(d[i, 1:6]), ncol = 1)
  d$condition[i] <- all(A %*% x <= b)
  d$value[i] <- sum(Pmat[cbind(1:6, x)])
}
d <- d[order(d$value, decreasing = FALSE), ]
d <- d[order(d$condition, decreasing = TRUE), ]
print(head(d))

这得到

#     x1 x2 x3 x4 x5 x6 condition value
#7789  1  3  1  1  1  2      TRUE   117
#7783  1  2  1  1  1  2      TRUE   128
#229   1  3  1  2  1  1      TRUE   131
#13    1  3  1  1  1  1      TRUE   137
#223   1  2  1  2  1  1      TRUE   142
#7777  1  1  1  1  1  2      TRUE   145

所以解决方案是 [1, 3, 1, 1, 1, 2]，值为 117。

LP配方

蛮力效率不高，您显然正在寻找一种使用 LP 求解器来求解的方法。为了为您的问题制定 LP，我将重新定义一组二进制变量，如下所示（不过，这可能是您所说的“大量额外变量”）：

y_{ij} = 1(x_i = j)  for each i = 1,...,6 and j = 1,...,6

这里， 1() 是一个函数，当且仅当内部语句为真时才取 1。

然后将这些变量向量化为

y = [y_{11}, y_{12}, ..., y_{16}, y_{21}, ..., y_{66}]

这是 LP 的新控制变量。请注意，大小y为 36。

然后，您需要将P矩阵向量化为 的权重向量y。您还需要转换矩阵A。例如，如果原始约束是

x_1 + x_2 <= b

那么这相当于

  y_{11} + 2*y_{12} + 3*y_{13} + 4*y_{14} + 5*y_{15} + 6*y_{16} 
+ y_{21} + 2*y_{22} + 3*y_{23} + 4*y_{24} + 5*y_{25} + 6*y_{26}  <= b

Kronecker 产品便于这种转换。

最后，对于每个i，必须只有一个j这样的y_{ij} = 1。因此你需要

y_{i1} + y_{i2} + ... + y_{i6} = 1  for each i.

Kronecker 产品对此也有帮助。

library(lpSolve)
Pvec <- as.numeric(t(Pmat))
A1 <- kronecker(A, matrix(1:6, nrow = 1))
b1 <- b
A2 <- kronecker(diag(6), matrix(rep(1, 6), nrow = 1))
b2 <- rep(1, 6)
o <- lp(direction = "min",
        objective.in = Pvec, const.mat = rbind(A1, A2),
        const.dir = c(rep("<=", 6), rep("=", 6)), const.rhs = c(b1, b2),
        all.bin = TRUE)
# show the variable names
tmp <- expand.grid(1:6, 1:6)
vname <- paste("x", tmp[,2], tmp[,1], sep = "")
names(o$solution) <- vname
print(o)
print(o$solution)

这会产生与蛮力相同的解决方案。

Success: the objective function is 117 

x11 x12 x13 x14 x15 x16 x21 x22 x23 x24 x25 x26 x31 x32 x33 x34 x35 x36 
  1   0   0   0   0   0   0   0   1   0   0   0   1   0   0   0   0   0 
x41 x42 x43 x44 x45 x46 x51 x52 x53 x54 x55 x56 x61 x62 x63 x64 x65 x66 
  1   0   0   0   0   0   1   0   0   0   0   0   0   1   0   0   0   0
于 2015-12-15T06:13:00.660 回答