我应该如何合并链接约束以满足使用 Binary.vec 的 Y1 = 1、Y2 = 2、Y3 = 3 和 Y4 = 4 的约束？
我相信你的意思是Y[i] = 1 [ for i in 1:4 ]。
binary.vec 参数是一个
像 int.vec 这样的数值向量，给出了二进制变量的索引。
现在这里的索引指的是优化问题。因此我们需要索引来匹配Y[i] [for i in 1:4 ] 的位置。在这种情况下，它是您的目标函数中的后 4 个，因此我们可以将它们的索引设为
binary.vec <- seq(length(f.obj), length(f.obj) - 3)
但是请注意，您当前的问题不是在寻找严格的二进制Y[i]，而是在给定其他参数的值的情况下进一步寻找Y[i]。正如它所写的那样，它不是线性约束，因此我们需要做一些诡计。
这里的问题是您需要Y[i] = 1 if X[i,1] > 0 and X[i, 2] > 0, 0 otherwise，而简单地添加Y[i] %in% 0:1 不一定会支持约束。有几种方法可以解决这个问题，但我发现最直观的一种方法是使用 pseudo-parameters 重新表述问题
C * z[i, j] > x[i, j] [ for all i, j ]
z[i, j] %in% 0:1 [ for all i, j ]
C = Inf
z[i, j] = 1 if X[i, j] > 0, 0 otherwise [ for i in 1:4, j in 1:2 ]
Y[i] <= 0.5 * z[i, 1] + 0.5 * z[i, 2]
这里的诀窍是现在我们只有线性约束。请注意，如果只有z[i, 1] = 1 这些约束会强制Y[i] = 0。
我们增加了相当多的参数数量，但这个问题是线性，因为所有方程都以y = a * x + b 的形式编写。
但是，正如我在对您的previous question 的回答中所描述的，当我们添加我们需要的伪参数时
- 扩展我们的目标函数，使其包含这些伪参数（让求解器了解优化函数中的参数）
- 扩展我们的约束，以保持尺寸正确（让求解器了解这些对剩余约束的影响）
扩展优化问题
现在这是简单的部分。只需在现有方程的末尾添加参数即可。如果发现难以跟踪，我总是建议命名目标函数
n <- length(f.obj)
names(f.obj)[seq(n - 4)] <- paste0('X[', rep(1:4, 2), ', ', rep(1:2, each = 4), ']')
names(f.obj)[-seq(n - 4)] <- paste0('Y[', 1:4, ']')
f.obj <- c(f.obj, numeric(8))
n <- length(f.obj)
names(f.obj)[seq(n - 7, n)] <- paste0('z[', rep(1:4, 2), ', ', rep(1:2, each = 4), ']')
tail(f.obj, 12)
Y[1] Y[2] Y[3] Y[4] z[1, 1] z[2, 1] z[3, 1] z[4, 1] z[1, 2] z[2, 2] z[3, 2] z[4, 2]
1500 1700 1500 1400 0 0 0 0 0 0 0 0
请注意，它们的权重都为 0，因为它们的效果被 X[i, j] 和 Y[i] 的效果所捕获。
扩展你的约束
现在这又是棘手的部分。首先请注意，在我的answer 中对上一个问题，您的约束的尺寸必须与您的目标函数相匹配。在这个问题中，你有 length(f.obj) = 28（包括新的伪参数），所以 all 约束需要有 28 列。这里的一般想法是我们（形象地）将sum_[i, j] 0 * X[i, j] + sum_[i] 0 * Y[i] 添加到每个约束。因此你的约束应该是
cons1 <- numeric(n)
names(cons1) <- names(f.obj)
cons12 <- cons11 <- cons10 <- cons9 <- cons8 <- cons7 <- cons6 <- cons5 <- cons4 <- cons3 <- cons2 <- cons1
cons1[1:4] <- c(2, 1.8, 2.3, 2.1)
cons2[5:8] <- c(2.8, 2.3, 2.2, 2.6)
cons3[9:12] <- c(1.7, 1.75, 1.6, 1.9)
cons4[13:16] <- c( 2.4, 1.9, 2.6, 2.4)
#X[i, j] - Cy[i] * Y[i] <= 0
cons5[c('X[1, 1]', 'Y[1]')] <- c(1, -1500)
cons6[c('X[2, 1]', 'Y[2]')] <- c(1, -2000)
cons7[c('X[3, 1]', 'Y[3]')] <- c(1, -1500)
cons8[c('X[4, 1]', 'Y[4]')] <- c(1, -1800)
cons9[c('X[1, 2]', 'Y[1]')] <- c(1, -500)
cons10[c('X[2, 2]', 'Y[2]')] <- c(1, -500)
cons11[c('X[3, 2]', 'Y[3]')] <- c(1, -500)
cons12[c('X[4, 2]', 'Y[4]')] <- c(1, -500)
#Pseudo-constraints
## Start with C * z[i, j] > x[i, j] for some very large C (here i choose 1e12, but could be any number large enough so X[i, j] does not exceed it!)
pseudo_z_1 <- numeric(n)
names(pseudo_z_1) <- names(f.obj)
pseudo_z_8 <- pseudo_z_7 <- pseudo_z_6 <- pseudo_z_5 <- pseudo_z_4 <- pseudo_z_3 <- pseudo_z_2 <- pseudo_z_1
C <- 1e12
pseudo_z_1[c('X[1, 1]', 'z[1, 1]')] <- c(-1, C)
pseudo_z_2[c('X[2, 1]', 'z[2, 1]')] <- c(-1, C)
pseudo_z_3[c('X[3, 1]', 'z[3, 1]')] <- c(-1, C)
pseudo_z_4[c('X[4, 1]', 'z[4, 1]')] <- c(-1, C)
pseudo_z_5[c('X[1, 2]', 'z[1, 2]')] <- c(-1, C)
pseudo_z_6[c('X[2, 2]', 'z[2, 2]')] <- c(-1, C)
pseudo_z_7[c('X[3, 2]', 'z[3, 2]')] <- c(-1, C)
pseudo_z_8[c('X[4, 2]', 'z[4, 2]')] <- c(-1, C)
## Next create pseudo constraints for Y[i] <= 0.5 * Z[i, 1] + 0.5 * Z[i, 2]
pseudo_y_1 <- numeric(n)
names(pseudo_y_1) <- names(f.obj)
pseudo_y_4 <- pseudo_y_3 <- pseudo_y_2 <- pseudo_y_1
pseudo_y_1[c('Y[1]', 'z[1, 1]', 'z[1, 2]')] <- c(1, -0.5, -0.5)
pseudo_y_2[c('Y[2]', 'z[2, 1]', 'z[2, 2]')] <- c(1, -0.5, -0.5)
pseudo_y_3[c('Y[3]', 'z[3, 1]', 'z[3, 2]')] <- c(1, -0.5, -0.5)
pseudo_y_4[c('Y[4]', 'z[4, 1]', 'z[4, 2]')] <- c(1, -0.5, -0.5)
#Combine them all together
cons <- c(paste0('cons', 1:12), paste0('pseudo_z_', 1:8), paste0('pseudo_y_', 1:4))
cons.mat <- do.call(rbind, mget(cons)) #mget finds all variables given in cons.
呼...然后我们必须再次创建我们的direction 和rhs。我们必须确保binary.vec 将z[i, j] 指定为二进制。让我们按顺序执行此操作
cons.dir <- c(rep('>=', 4), rep('<=', 8), rep('>=' 8), rep('<=', 4))
names(cons.dir) <- cons
cons.rhs <- c(750, 800, 1000, 300, rep(0, 8), rep(0, 8), rep(0, 4))
#Update binary.vec to include z[i, j]
binary.vec <- c(binary.vec, seq(n, n - 7))
现在我们应该已经准备好一切准备好了
lpSolve::lp('min', objective.in = f.obj,
const.mat = cons.mat,
const.dir = cons.dir,
const.rhs = cons.rhs,
binary.vec = binary.vec)
res$objval
[1] 27483.29
names(res$solution) <- names(f.obj)
options(scipen = 999) #avoid scientific notation. Just to make it readable for us.
round(res$solution[res$solution != 0], 3)
X[4, 1] X[4, 2] X[4, 3] X[3, 4] Y[4] z[4, 1] z[4, 2]
357.143 307.692 526.316 115.385 1.000 1.000 1.000
这似乎符合我们的限制条件。由于复制粘贴到 R 中比较麻烦，我已将代码添加到 this repository，可以加载和执行。
旁注
我在这里说明解决方案的方式是一个非常说明性的分步方法。这使得跟随变得简单，但实际上人们会使用for-loops 或允许更自然地编写约束的软件。如果问题变得比这更复杂，我强烈建议使用替代包，它也可能实现更快的算法。
在满足链接变量之前是否应该先生成最优解？
在您当前的问题中，您没有任何链接变量。链接变量通常是从子优化问题获得的变量，例如
Max sum_i a[i] * y[i] * x[i] [for i in 1:10]
st.
y[i] %in% 0:1
x[i] = argmin{
3 x[i]^2 + 4
st.
x[i] != x[j] [ for all i, j ]
}
在这种情况下，标准优化问题不一定有用。