在方程1-4中，Y_i的是二元变量，应该如何将联动约束纳入优化问题的公式中？
minimize
22X11 + 20.5X12 + 21X13 + 19X14 + 21X21 + 21X22 + 18X23 + 17X24 +
22X31 + 24X32 + 22X33 + 23X34 + 22X41 + 16X42 + 17X43 + 17.5X44 +
1500 Y1 + 1700 Y2 + 1500 Y3 + 1400 Y4
subject to constraints
2.00 X11 + 1.80 X21 + 2.30 X31 + 2.10 X41 >= 750
2.80 X12 + 2.30 X22 + 2.20 X32 + 2.60 X42 >= 800
1.70 X13 + 1.75 X23 + 1.60 X33 + 1.90 X43 >= 1000
2.40 X14 + 1.90 X24 + 2.60 X34 + 2.40 X44 >= 300
X11 - 1500 Y1 <= 0
X21 - 2000 Y2 <= 0
X31 - 1500 Y3 <= 0
X41 - 1800 Y4 <= 0
X12 - 500 Y1 >= 0
X22 - 500 Y2 >= 0
X32 - 500 Y3 >= 0
X42 - 500 Y4 >= 0
Y1 = 1 if X11>0 and 0 otherwise
Y1 = 1 if X12>0 and 0 otherwise
Y2 = 1 if X21>0 and 0 otherwise
Y2 = 1 if X22>0 and 0 otherwise
Y3 = 1 if X31>0 and 0 otherwise
Y3 = 1 if X32>0 and 0 otherwise
Y4 = 1 if X41>0 and 0 otherwise
Y4 = 1 if X42>0 and 0 otherwise
在没有这些链接变量的情况下，用R代码找到最优解。
library(lpSolve)
f.obj <- c(22,20.5,21,19,21,21,18,17,22,24,22,23,22,16,17,17.5,1500,1700,1500,1400)
const1 = c(2.00,1.80,2.30,2.10)
const2 = c(2.80,2.30,2.20,2.60)
const3 = c(1.70,1.75,1.60,1.90)
const4 = c(2.40,1.90,2.60,2.40)
f.con <- rbind(const1,const2,const3,const4)
f.dir <- c(rep(">=",2), ">=", ">=")
f.rhs <- c(750,800,1000,300)
res <- lp("min", f.obj, f.con,f.dir,f.rhs, binary.vec=c(1,2))
res$solution
res$objective
res$objval
我相信你的意思是 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.
然而，正如我在回答你的 前一个问题 当我们添加伪参数的时候，我们需要用到
现在这是简单的部分。只需在现有方程的末尾添加参数即可。我总是建议给目标函数起个名字，如果有人觉得很难跟踪的话。
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].
现在这又是棘手的部分。首先要注意的是，在我的 回答 对上一个问题的回答是，你的约束条件的尺寸必须与你的目标函数相匹配。在这个问题中，你有 length(f.obj) = 28（包括新的伪参数），所以 都 限制条件需要有 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中，我将代码添加到了 本库，可以加载并执行。
我在这里说明的解决方法，是一个很有说明性的步骤方法。这样一来，就很简单了，但在现实中，人们会用 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 ]
}
而在这种情况下，标准的优化问题不一定有用。
这里有一个选择。
将X_ij的转换为小数。我只是用比这里看到的数值大一个数量级。
library(lpSolve)
n <- 20
denom <- 10e3
f.obj <- c(c(22,20.5,21,19,
21,21,18,17,
22,24,22,23,
22,16,17,17.5) / denom,
1500,1700,1500,1400)
#identify location of variables for easy identification
vars <- c(paste0("X", c(outer(1:4, 1:4, paste0))), paste0("Y", 1:4))
vars <- setNames(seq_along(vars), vars)
#X11 X21 X31 X41 X12 X22 X32 X42 X13 X23 X33 X43 X14 X24 X34 X44 Y1 Y2 Y3 Y4
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
按照OP的约束，同时适当的收缩系数。
const01 <- replace(rep(0, n), 1:4, c(2.00,1.80,2.30,2.10))
const02 <- replace(rep(0, n), 5:8, c(2.80,2.30,2.20,2.60))
const03 <- replace(rep(0, n), 9:12, c(1.70,1.75,1.60,1.90))
const04 <- replace(rep(0, n), 13:16, c(2.40,1.90,2.60,2.40))
#X11, X12, Y1
const05 <- replace(rep(0, n), c(1,17), c(-1, 1500 / denom))
const06 <- replace(rep(0, n), c(5,17), c(1, -500 / denom))
#X21, X22, Y2
const07 <- replace(rep(0, n), c(2,18), c(-1, 2000 / denom))
const08 <- replace(rep(0, n), c(6,18), c(1, -500 / denom))
#X31, X32, Y3
const09 <- replace(rep(0, n), c(3,19), c(-1, 1500 / denom))
const10 <- replace(rep(0, n), c(7,19), c(1, -500 / denom))
#X41, X42, Y4
const11 <- replace(rep(0, n), c(4,20), c(-1, 1800 / denom))
const12 <- replace(rep(0, n), c(8,20), c(1, -500 / denom))
对于连接变量，使用X <=Y，所以当X为非零时，Y必须为1。
#X11, X12, Y1
const11 <- replace(rep(0, n), c(1,17), c(-1, 1))
const12 <- replace(rep(0, n), c(5,17), c(-1, 1))
#X21, X22, Y2
const13 <- replace(rep(0, n), c(2,18), c(-1, 1))
const14 <- replace(rep(0, n), c(6,18), c(-1, 1))
#X31, X32, Y3
const15 <- replace(rep(0, n), c(3,19), c(-1, 1))
const16 <- replace(rep(0, n), c(7,19), c(-1, 1))
#X41, X42, Y4
const17 <- replace(rep(0, n), c(4,20), c(-1, 1))
const18 <- replace(rep(0, n), c(8,20), c(-1, 1))
调用求解器，同时适当收缩RHS。
f.con <- do.call(rbind, mget(ls(pattern="^const")))
f.dir <- rep(">=", nrow(f.con))
f.rhs <- c(c(750,800,1000,300) / denom, rep(0, nrow(f.con) - 4))
res <- lp("min", f.obj, f.con, f.dir, f.rhs, binary.vec=17:20)
res$solution
# [1] 0.00000000 0.00000000 0.00000000 0.03571429 0.00000000 0.00000000 0.00000000 0.03076923
# [9] 0.00000000 0.00000000 0.00000000 0.05263158 0.00000000 0.00000000 0.01153846 0.00000000
#[17] 0.00000000 0.00000000 0.00000000 1.00000000
反过来收缩。
(optval <- setNames(replace(res$solution, 1:16, denom * res$solution[1:16]), names(vars)))
和最佳值。
X11 X21 X31 X41 X12 X22 X32 X42 X13 X23
0.0000 0.0000 0.0000 357.1429 0.0000 0.0000 0.0000 307.6923 0.0000 0.0000
X33 X43 X14 X24 X34 X44 Y1 Y2 Y3 Y4
0.0000 526.3158 0.0000 0.0000 115.3846 0.0000 0.0000 0.0000 0.0000 1.0000
Obj值：
sum(optval * replace(f.obj, 1:16, denom * f.obj[1:16]))
#[1] 27483.29