我想使用具有分段常数目标函数的R来解决最小化问题。其思想是，对于我的(整数)决策变量x的较低值，会产生比较高值更高的惩罚成本。考虑到一些限制，我想尽量减少总罚款成本。
因此，我的程序如下：
min P(x)
s.t. A x <= b
x >= 0
我在编程分段常数目标函数p(x)时遇到了一个问题，其中p是向量x的所有元素的和。我知道它不能与linprog库中的函数lp()组合使用。然而，如果不指定大量额外的变量，我就无法找到实现这一点的方法。此外，广泛的互联网搜索并没有提供任何有用的想法。
让我举一个例子，说明这个函数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找到分段常数目标函数的最小值？
我对你的问题有两种方法。您的问题似乎是一个整数编程，其中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=4656个候选者。因此，用蛮力解决这个问题不会花那么长时间；对于每个组合，检查它是否满足约束和值。
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,2]，值为117。
LP配方
蛮力是无效的，你显然是在寻找一种使用LP解算器解决问题的方法。为了为您的问题制定LP，我将重新定义一组二进制变量如下(不过，这可能是您所说的"大量额外变量")：
y_{ij} = 1(x_i = j) for each i = 1,...,6 and j = 1,...,6
这里，1()是一个接受1的函数，当且仅当其中的语句为true时。
然后将这些变量矢量化为
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
克罗内克的产品便于这种转换。
最后，对于每个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