R 中具有条件约束的线性规划
Linear programming with conditional constraints in R
我有一个线性规划问题,我试图从一些二进制资源中 select 来优化价值,基本上是一个背包问题。我遇到的问题是不同的资源具有共同的特征,我想确保我的最终解决方案具有 0 或 2 个具有特定特征的资源。有什么办法可以做到这一点?尽管进行了广泛的搜索,但我一直无法想到或找到一个。在我的数据中,决策变量是资源,约束是这些资源的特征。考虑以下代码:
library(lpSolve)
const_mat_so<-matrix(c(
c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,1,0,0,1,0,1)
,c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,1,1,0,0,1,1)
,c(0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,0,1,0,1,0,0)
,c(1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0,0,0,0,0,0,0)
,c(8800, 8500, 7600, 8600, 8400, 7500, 7000, 8500, 8800, 7700, 6700,5500,1200,6700,9500,8700,6500)
,c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,0,0,1,0,1,0)
,c(0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,1,1,0,0,0,0)
,c(0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,1,1,1,0,1,0)
,c(0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,1,0)
,c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,0,0,0,1,0,0)
,c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,0,0,0,0,0,0)
),nrow=15,byrow = TRUE)
const_dir_so<-c("=","=","=","=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=")
max_cost_so = 25000
objective_so = c(21.0, 19.3, 19.2, 18.8, 18.5, 16.6, 16.4, 16.4, 16.0, 16.0, 14.9, 14.6, 14.0, 13.9,12.0,5.5,24.6)
const_rhs_so<-c(1,1,1,1,25000,3,3,3,2,2,2,2,2,2,2)
x = lp ("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin=TRUE, all.int=TRUE
)
> x
Success: the objective function is 68.1
> x$solution
[1] 1 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0
虽然上面产生了一个解决方案,但这不是我想要的解决方案,因为我实际上希望最后七个约束为 >=2 或 0。我不知道如何编写代码或是否可能。任何帮助,将不胜感激。我不是线性规划高手,所以请原谅对这种方法的任何误解。
我的理解是,最后 7 个约束中的每一个都大于 2 或等于 0,即不是 1。
1) 只有 7 个这样的约束,所以有 2^7 = 128 种可能性,这足够小,我们可以 运行 每个使用公式的可能性他在没有过多 运行 时间的情况下提出问题,然后找到其中的最大值。
dec2bin
以 10 为基数(即十进制)并将其转换为 0 和 1 的二进制向量。 运行 它在 0 到 127 之间的每个数字上给出二进制数,使得 1 对应于 >= 2 的约束(其余等于 0)。
dec2bin <- function(dec, digits = 7) {
# see
tail(rev(as.integer(intToBits(dec))), digits)
}
runLP <- function(i) {
bin <- dec2bin(i)
n <- length(const_rhs_so) # 15
ix <- seq(to = n, length = length(bin)) # indexes of last 7 constraints, i.e. 9:15
const_dir_so[ix] <- ifelse(bin, ">=", "=")
const_rhs_so[ix] <- 2*bin
lp("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin = TRUE)
}
lpout <- lapply(0:127, runLP)
ixmax <- which.max(sapply(lpout, "[[", "objval"))
ans <- lpout[[ixmax]]
ans
ans$solution
tail(c(const_mat_so %*% ans$solution), 7)
给予:
> ans
Success: the objective function is 62
> ans$solution
[1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
> tail(c(const_mat_so %*% ans$solution), 7) # last 7 constraint values
[1] 0 0 0 0 0 0 0
2) 在@Erwin Kalvelagen 的第二个替代方案中,它指的是约束变量,但我认为他的回答中的 x
是 LHS 的值最后 7 个约束之一。也就是说,如果 C
是最初的最后 7 个约束的矩阵,则用这 14 个约束替换原来的 7 个约束:
Cx + D1 y <= 0
Cx + D2 y >= 0
其中D1是对角元素为任意足够大的负数的对角矩阵,D2是对角元素全为-2的对角矩阵。在这里,我们正在优化二进制变量的 x
和 y
向量。 x
变量与问题中的一样,并且有 7 个新的 y
二进制变量,使得 y[i] 为 0 以将最后 7 个原始约束中的第 i 个约束为 0 或 1 以将其约束为2个或更多。 y
变量在 (1) 中称为 bin
。 objective中y
个变量的系数全为零
lpSolve R 代码方面:
objective_so2 <- c(objective_so, numeric(7))
const_mat_so2 <- cbind(rbind(const_mat_so, const_mat_so[9:15, ]),
rbind(matrix(0, 8, 7), diag(-100, 7), diag(-2, 7)))
const_dir_so2 <- c(const_dir_so, rep(">=", 7))
const_rhs_so2 <- c(const_rhs_so[1:8], numeric(14))
x2 = lp ("max", objective_so2, const_mat_so2, const_dir_so2, const_rhs_so2, all.bin = TRUE)
给出与 (1) 相同的值 62。 y
变量(最后 7 个)均为 0,也对应于 (1)。这也提供了双重检查,因为两种方法现在给出了一致的答案。
> x2
Success: the objective function is 62
> x2$solution
[1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
我相信 LpSolve 支持半连续变量。具有下限 L 和上限 U 的半连续变量可以假定值为 0 或介于 L 和 U 之间。我不确定 R 包 lpSolve 是否支持此变量类型。
但是我们可以用一个额外的二进制变量 y
和额外的约束来模拟它。所以你需要让你的 x
变量连续(或者整数,如果你只想要整数值)并添加约束:
2*y <= x <= U*y
其中 U
是 x
的上限。
lpSolveAPI
软件包为 "lp_solve" 提供了更高级的接口。正如@Erwin Kalvelagen 提到的,"lp_solve" 和 lpSolveAPI
支持半连续变量(半连续决策变量可以在其上限和下限之间取允许值,也可以取零)。约束矩阵使您能够将第 9-15 个约束公式的输出转换为第 18-24 个变量。例如(关于第9个约束),当x6 + x11 + x14 + x16 - x18 = 0
,x6 + x11 + x14 + x16 = x18
。所以我认为你可以通过半连续变量 x18
.
来控制 x6 + x11 + x14 + x16
library(lpSolveAPI)
## add 18-24th cols to define the 18-24th variables
const_mat_so2 <- cbind(const_mat_so, rbind(matrix(0, nrow = 8, ncol = 7), diag(-1, 7)))
## [EDITED] make a model and set a constraint matrix and objective coefs
model <- make.lp(nrow(const_mat_so2), 0)
for(i in 1:ncol(const_mat_so2)) add.column(model, const_mat_so2[,i])
set.constr.type(model, c(const_dir_so[-c(9:15)], rep("=", 7)))
set.rhs(model, c(const_rhs_so[-c(9:15)], rep(0, 7))) # each original output - 18-24th = 0
set.objfn(model, c(objective_so, rep(0, 7))) # 18-24th are 0
## define semi-continuous and bounds.
set.semicont(model, col = 18:24)
set.bounds(model, lower = rep(1.9, 7), col = 18:24) # default upper is Inf.
## define other things
set.type(model, col = 1:17, type = "binary") # original variable
set.type(model, col = 18:24, type = "integer") # outputs of original constraint formulas
lp.control(model, sense = "max") # do maximize
# write.lp(model, "filename.lp", "lp") # if you want to watch the whole model
solve(model)
get.variables(model)
# [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 [18] 0 0 0 0 0 0 0
get.objective(model)
# [1] 62
t(const_mat_so %*% res[1:17])
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
# [1,] 1 1 1 1 22300 1 0 0 0 0 0 0 0 0 0
我有一个线性规划问题,我试图从一些二进制资源中 select 来优化价值,基本上是一个背包问题。我遇到的问题是不同的资源具有共同的特征,我想确保我的最终解决方案具有 0 或 2 个具有特定特征的资源。有什么办法可以做到这一点?尽管进行了广泛的搜索,但我一直无法想到或找到一个。在我的数据中,决策变量是资源,约束是这些资源的特征。考虑以下代码:
library(lpSolve)
const_mat_so<-matrix(c(
c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,1,0,0,1,0,1)
,c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,1,1,0,0,1,1)
,c(0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1,0,1,0,1,0,0)
,c(1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0,0,0,0,0,0,0)
,c(8800, 8500, 7600, 8600, 8400, 7500, 7000, 8500, 8800, 7700, 6700,5500,1200,6700,9500,8700,6500)
,c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,0,0,1,0,1,0)
,c(0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,1,1,0,0,0,0)
,c(0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,0,0)
,c(0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,1,1,1,0,1,0)
,c(0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,0,0,0,0,1,0)
,c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,0,0,0,1,0,0)
,c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,0,0,0,0,0,0)
),nrow=15,byrow = TRUE)
const_dir_so<-c("=","=","=","=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=","<=")
max_cost_so = 25000
objective_so = c(21.0, 19.3, 19.2, 18.8, 18.5, 16.6, 16.4, 16.4, 16.0, 16.0, 14.9, 14.6, 14.0, 13.9,12.0,5.5,24.6)
const_rhs_so<-c(1,1,1,1,25000,3,3,3,2,2,2,2,2,2,2)
x = lp ("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin=TRUE, all.int=TRUE
)
> x
Success: the objective function is 68.1
> x$solution
[1] 1 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0
虽然上面产生了一个解决方案,但这不是我想要的解决方案,因为我实际上希望最后七个约束为 >=2 或 0。我不知道如何编写代码或是否可能。任何帮助,将不胜感激。我不是线性规划高手,所以请原谅对这种方法的任何误解。
我的理解是,最后 7 个约束中的每一个都大于 2 或等于 0,即不是 1。
1) 只有 7 个这样的约束,所以有 2^7 = 128 种可能性,这足够小,我们可以 运行 每个使用公式的可能性他在没有过多 运行 时间的情况下提出问题,然后找到其中的最大值。
dec2bin
以 10 为基数(即十进制)并将其转换为 0 和 1 的二进制向量。 运行 它在 0 到 127 之间的每个数字上给出二进制数,使得 1 对应于 >= 2 的约束(其余等于 0)。
dec2bin <- function(dec, digits = 7) {
# see
tail(rev(as.integer(intToBits(dec))), digits)
}
runLP <- function(i) {
bin <- dec2bin(i)
n <- length(const_rhs_so) # 15
ix <- seq(to = n, length = length(bin)) # indexes of last 7 constraints, i.e. 9:15
const_dir_so[ix] <- ifelse(bin, ">=", "=")
const_rhs_so[ix] <- 2*bin
lp("max", objective_so, const_mat_so, const_dir_so, const_rhs_so, all.bin = TRUE)
}
lpout <- lapply(0:127, runLP)
ixmax <- which.max(sapply(lpout, "[[", "objval"))
ans <- lpout[[ixmax]]
ans
ans$solution
tail(c(const_mat_so %*% ans$solution), 7)
给予:
> ans
Success: the objective function is 62
> ans$solution
[1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
> tail(c(const_mat_so %*% ans$solution), 7) # last 7 constraint values
[1] 0 0 0 0 0 0 0
2) 在@Erwin Kalvelagen 的第二个替代方案中,它指的是约束变量,但我认为他的回答中的 x
是 LHS 的值最后 7 个约束之一。也就是说,如果 C
是最初的最后 7 个约束的矩阵,则用这 14 个约束替换原来的 7 个约束:
Cx + D1 y <= 0
Cx + D2 y >= 0
其中D1是对角元素为任意足够大的负数的对角矩阵,D2是对角元素全为-2的对角矩阵。在这里,我们正在优化二进制变量的 x
和 y
向量。 x
变量与问题中的一样,并且有 7 个新的 y
二进制变量,使得 y[i] 为 0 以将最后 7 个原始约束中的第 i 个约束为 0 或 1 以将其约束为2个或更多。 y
变量在 (1) 中称为 bin
。 objective中y
个变量的系数全为零
lpSolve R 代码方面:
objective_so2 <- c(objective_so, numeric(7))
const_mat_so2 <- cbind(rbind(const_mat_so, const_mat_so[9:15, ]),
rbind(matrix(0, 8, 7), diag(-100, 7), diag(-2, 7)))
const_dir_so2 <- c(const_dir_so, rep(">=", 7))
const_rhs_so2 <- c(const_rhs_so[1:8], numeric(14))
x2 = lp ("max", objective_so2, const_mat_so2, const_dir_so2, const_rhs_so2, all.bin = TRUE)
给出与 (1) 相同的值 62。 y
变量(最后 7 个)均为 0,也对应于 (1)。这也提供了双重检查,因为两种方法现在给出了一致的答案。
> x2
Success: the objective function is 62
> x2$solution
[1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
我相信 LpSolve 支持半连续变量。具有下限 L 和上限 U 的半连续变量可以假定值为 0 或介于 L 和 U 之间。我不确定 R 包 lpSolve 是否支持此变量类型。
但是我们可以用一个额外的二进制变量 y
和额外的约束来模拟它。所以你需要让你的 x
变量连续(或者整数,如果你只想要整数值)并添加约束:
2*y <= x <= U*y
其中 U
是 x
的上限。
lpSolveAPI
软件包为 "lp_solve" 提供了更高级的接口。正如@Erwin Kalvelagen 提到的,"lp_solve" 和 lpSolveAPI
支持半连续变量(半连续决策变量可以在其上限和下限之间取允许值,也可以取零)。约束矩阵使您能够将第 9-15 个约束公式的输出转换为第 18-24 个变量。例如(关于第9个约束),当x6 + x11 + x14 + x16 - x18 = 0
,x6 + x11 + x14 + x16 = x18
。所以我认为你可以通过半连续变量 x18
.
x6 + x11 + x14 + x16
library(lpSolveAPI)
## add 18-24th cols to define the 18-24th variables
const_mat_so2 <- cbind(const_mat_so, rbind(matrix(0, nrow = 8, ncol = 7), diag(-1, 7)))
## [EDITED] make a model and set a constraint matrix and objective coefs
model <- make.lp(nrow(const_mat_so2), 0)
for(i in 1:ncol(const_mat_so2)) add.column(model, const_mat_so2[,i])
set.constr.type(model, c(const_dir_so[-c(9:15)], rep("=", 7)))
set.rhs(model, c(const_rhs_so[-c(9:15)], rep(0, 7))) # each original output - 18-24th = 0
set.objfn(model, c(objective_so, rep(0, 7))) # 18-24th are 0
## define semi-continuous and bounds.
set.semicont(model, col = 18:24)
set.bounds(model, lower = rep(1.9, 7), col = 18:24) # default upper is Inf.
## define other things
set.type(model, col = 1:17, type = "binary") # original variable
set.type(model, col = 18:24, type = "integer") # outputs of original constraint formulas
lp.control(model, sense = "max") # do maximize
# write.lp(model, "filename.lp", "lp") # if you want to watch the whole model
solve(model)
get.variables(model)
# [1] 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 [18] 0 0 0 0 0 0 0
get.objective(model)
# [1] 62
t(const_mat_so %*% res[1:17])
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
# [1,] 1 1 1 1 22300 1 0 0 0 0 0 0 0 0 0