如何解决 R 中的线性规划模型
How to solve linear programming model in R
我需要解决以下微观经济问题：
- 我可以在五年内 (2011 - 2015) 生产六项资产（资产 1 - 6）。
- 每个资产只能在一年内生产。
- 每项资产必须在我的五年期限内生产。
- 生产不互斥；我可以在一年内生产不止一种商品，而不会影响任何一种商品的生产。
- 每项资产的固定生产成本等于 30。
- 我每年必须有非负利润；收入必须至少为 30。
下面是一个矩阵，表示我在给定年份 (j) 生产每项资产 (i) 的潜在收入。
2011 2012 2013 2014 2015
Asset1 35* 37 39 42 45
Asset2 16 17 18 19 20*
Asset3 125 130 136*139 144
Asset4 15 27 29 30* 33
Asset5 14 43* 46 50 52
Asset6 5 7 8 10 11*
星号 (*) 代表最佳解决方案集。
我如何使用 R 解决生产计划，使我的收入（以及利润）在所列约束条件下最大化。我的输出应该是类似的 0 和 1 的 6x5 矩阵，其中 1 代表选择在给定年份生产商品。
这是一个经典问题，需要重新表述。
首先重新表述您的问题
Max( sum_[i,t] (pi_[i,t] - C_[i,t]) * x_[i,t])
Sd.
sum_t x_[i,t] = 1 [ for all i ]
sum_i x_[i,t] >= 30 [ for all t ]
x_[i,t] >= 0 [for all i, t]
在 lpSolve 包中，最大化问题以线性表示形式给出，例如。非矩阵格式。让我们从制作一个代表我们的 x_[i,t] 的向量开始。为了方便起见，我们给它起个名字（虽然没有用到），这样我们就可以跟踪。
n <- 6
t <- 5
#x ordered by column.
x <- c(35, 16, 125, 15, 14, 5, 37, 17, 130, 27, 43, 7, 39, 18, 136, 29, 46, 8, 42, 19, 139, 30, 50, 10, 45, 20, 144, 33, 52, 11)
# if x is matrix use:
# x <- as.vector(x)
names(x) <- paste0('x_[', seq(n), ',', rep(seq(t), each = n), ']')
head(x, n * 2)
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2] x_[2,2] x_[3,2] x_[4,2] x_[5,2] x_[6,2]
35 16 125 15 14 5 37 17 130 27 43 7
length(x)
[1] 30
现在我们需要创造条件。从第一个条件开始
sum_t x_[i,t] = 1 [ for all i ]
我们可以很简单地创建它。这里要注意的是尺寸必须正确。我们有一个长度为 30 的向量，因此我们需要条件矩阵有 30 列。此外，我们有 6 个资产，因此我们需要 6 行来满足此条件。再次让我们命名行和列以跟踪我们自己。
cond1 <- matrix(0, ncol = t * n,
nrow = n,
dimnames = list(paste0('x_[', seq(n), ',t]'),
names(x)))
cond1[, seq(n + 1)]
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2]
x_[1,t] 0 0 0 0 0 0 0
x_[2,t] 0 0 0 0 0 0 0
x_[3,t] 0 0 0 0 0 0 0
x_[4,t] 0 0 0 0 0 0 0
x_[5,t] 0 0 0 0 0 0 0
x_[6,t] 0 0 0 0 0 0 0
接下来我们填写正确的字段。 x_[1,1] + x[1, 2] + ... = 1 和 x_[2,1] + x_[2,2] + ... = 1 等等。对于这个问题，使用 for 循环是最简单的
for(i in seq(n)){
cond1[i, seq(i, 30, n)] <- 1
}
cond1[, seq(n + 1)]
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2]
x_[1,t] 1 0 0 0 0 0 1
x_[2,t] 0 1 0 0 0 0 0
x_[3,t] 0 0 1 0 0 0 0
x_[4,t] 0 0 0 1 0 0 0
x_[5,t] 0 0 0 0 1 0 0
x_[6,t] 0 0 0 0 0 1 0
我们仍然需要创建 RHS 并指定方向，但我现在等一下。
所以接下来让我们为第二个条件创建我们的矩阵
sum_i x_[i,t] >= 30 [ for all t ]
这个过程非常相似，但现在每个周期需要一行，所以矩阵的维度是 5x30。这里的主要区别是我们需要插入 x_[i, t]
的值
cond2 <- matrix(0, ncol = t * n,
nrow = t,
dimnames = list(paste0('t=', seq(t)),
names(x)))
for(i in seq(t)){
cond2[i, seq(n) + n * (i - 1)] <- x[seq(n) + n * (i - 1)]
}
cond2[, seq(1, n * t, n)]
x_[1,1] x_[1,2] x_[1,3] x_[1,4] x_[1,5]
t=1 35 0 0 0 0
t=2 0 37 0 0 0
t=3 0 0 39 0 0
t=4 0 0 0 42 0
t=5 0 0 0 0 45
请注意，我正在打印 x_[1, t] 的结果以说明我们做对了。
最后我们有最终条件。为此，我们注意到 ?lpSolve::lp 有一个参数 all.bin，读到这里，它表示
Logical: should all variables be binary? Default: FALSE.
因此，由于所有变量不是 1 就是 0，我们只需将此值设置为 TRUE。在继续之前，让我们将我们的条件组合成一个矩阵
cond <- rbind(cond1, cond2)
现在 RHS 和方向都简单地取自 2 个条件。来自 const.dir 参数
的文档
Vector of character strings giving the direction of the constraint: each value should be one of "<," "<=," "=," "==," ">," or ">=". (In each pair the two values are identical.)
在我们的条件中，我们有 6 行代表第一个条件，行代表条件 2。因此我们需要 n (6) 次 == 和 t (5) 次>=.
cond_dir <- c(rep('==', n), rep('>=', t))
RHS 以类似的方式创建
RHS <- c(rep(1, n), rep(30, t))
就是这样！现在我们准备好使用 lpSolve::lp 函数来解决我们的问题。
sol = lpSolve::lp(direction = 'max',
objective.in = x,
const.mat = cond,
const.dir = cond_dir,
const.rhs = RHS,
all.bin = TRUE)
sol$objval
[1] 275
解决方案的权重存储在 sol$solution
names(sol$solution) <- names(x)
sol$solution
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2] x_[2,2] x_[3,2] x_[4,2] x_[5,2] x_[6,2] x_[1,3] x_[2,3] x_[3,3]
1 0 0 0 0 0 0 0 0 0 1 0 0 0 1
x_[4,3] x_[5,3] x_[6,3] x_[1,4] x_[2,4] x_[3,4] x_[4,4] x_[5,4] x_[6,4] x_[1,5] x_[2,5] x_[3,5] x_[4,5] x_[5,5] x_[6,5]
0 0 0 0 0 0 1 0 0 0 1 0 0 0 1
matrix(sol$solution,
ncol = t,
dimnames = list(rownames(cond1),
rownames(cond2)))
t=1 t=2 t=3 t=4 t=5
x_[1,t] 1 0 0 0 0
x_[2,t] 0 0 0 0 1
x_[3,t] 0 0 1 0 0
x_[4,t] 0 0 0 1 0
x_[5,t] 0 1 0 0 0
x_[6,t] 0 0 0 0 1
我们很快就会发现这是正确的解决方案。 :-)
关于费用的附注
人们可能已经注意到“成本到底去哪儿了？”。在这种特定情况下，成本是固定的并且不是很有趣。这意味着我们可以在计算期间忽略这些，因为我们知道总成本将为 30 * 6 = 180（必须从 objective 值中减去）。然而，成本取决于各种因素并可能影响最佳解决方案的情况并不少见。为了便于说明，我将在此处包括我们如何将成本合并到此示例中。
首先，我们必须扩展我们的 objective 向量以包含每个产品在每个时期的成本
Fixed_C <- -30
x <- c(x, rep(Fixed_C, n * t))
接下来我们将添加一个伪约束
x_[i,t] - C_[i,t] = 0 [for all i, t]
此约束确保如果 x_[i,t] = 1 则将相关成本添加到问题中。有两种方法可以创建此约束。第一个是有一个包含 n * t 行的矩阵，每个行对应一个成本和期间。或者，我们可以使用我们的第一个约束条件，实际上只使用一个约束条件
sum_[i,t] x_[i,t] - C_[i,t] = 0
因为我们的第一个约束确保x[1, 1] != x[1, 2]。所以我们的第三个约束变成
cond3 <- c(rep(1, n * t), rep(-1, n * t))
最后我们必须扩展我们的 RHS 和条件 1 和 2 矩阵。只需将 0 添加到条件矩阵即可使尺寸合适。
cond1 <- cbind(cond1, matrix(0, nrow = n, ncol = n * t))
cond2 <- cbind(cond2, matrix(0, nrow = n, ncol = n * t))
cond <- rbind(cond1, cond2, cond3)
cond_dir <- c(cond_dir, '==')
RHS <- c(RHS, 0)
现在我们可以再次使用lpSolve::lp
找到最优解
solC = lpSolve::lp(direction = 'max',
objective.in = x,
const.mat = cond,
const.dir = cond_dir,
const.rhs = RHS,
all.bin = TRUE)
solC$objval
[1] 95
这等于我们之前的价值 275 减去我们的固定成本 Fixed_C * n = 180。
我需要解决以下微观经济问题：
- 我可以在五年内 (2011 - 2015) 生产六项资产（资产 1 - 6）。
- 每个资产只能在一年内生产。
- 每项资产必须在我的五年期限内生产。
- 生产不互斥；我可以在一年内生产不止一种商品，而不会影响任何一种商品的生产。
- 每项资产的固定生产成本等于 30。
- 我每年必须有非负利润；收入必须至少为 30。
下面是一个矩阵，表示我在给定年份 (j) 生产每项资产 (i) 的潜在收入。
2011 2012 2013 2014 2015
Asset1 35* 37 39 42 45
Asset2 16 17 18 19 20*
Asset3 125 130 136*139 144
Asset4 15 27 29 30* 33
Asset5 14 43* 46 50 52
Asset6 5 7 8 10 11*
星号 (*) 代表最佳解决方案集。
我如何使用 R 解决生产计划，使我的收入（以及利润）在所列约束条件下最大化。我的输出应该是类似的 0 和 1 的 6x5 矩阵，其中 1 代表选择在给定年份生产商品。
这是一个经典问题，需要重新表述。
首先重新表述您的问题
Max( sum_[i,t] (pi_[i,t] - C_[i,t]) * x_[i,t])
Sd.
sum_t x_[i,t] = 1 [ for all i ]
sum_i x_[i,t] >= 30 [ for all t ]
x_[i,t] >= 0 [for all i, t]
在 lpSolve 包中，最大化问题以线性表示形式给出，例如。非矩阵格式。让我们从制作一个代表我们的 x_[i,t] 的向量开始。为了方便起见，我们给它起个名字（虽然没有用到），这样我们就可以跟踪。
n <- 6
t <- 5
#x ordered by column.
x <- c(35, 16, 125, 15, 14, 5, 37, 17, 130, 27, 43, 7, 39, 18, 136, 29, 46, 8, 42, 19, 139, 30, 50, 10, 45, 20, 144, 33, 52, 11)
# if x is matrix use:
# x <- as.vector(x)
names(x) <- paste0('x_[', seq(n), ',', rep(seq(t), each = n), ']')
head(x, n * 2)
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2] x_[2,2] x_[3,2] x_[4,2] x_[5,2] x_[6,2]
35 16 125 15 14 5 37 17 130 27 43 7
length(x)
[1] 30
现在我们需要创造条件。从第一个条件开始
sum_t x_[i,t] = 1 [ for all i ]
我们可以很简单地创建它。这里要注意的是尺寸必须正确。我们有一个长度为 30 的向量，因此我们需要条件矩阵有 30 列。此外，我们有 6 个资产，因此我们需要 6 行来满足此条件。再次让我们命名行和列以跟踪我们自己。
cond1 <- matrix(0, ncol = t * n,
nrow = n,
dimnames = list(paste0('x_[', seq(n), ',t]'),
names(x)))
cond1[, seq(n + 1)]
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2]
x_[1,t] 0 0 0 0 0 0 0
x_[2,t] 0 0 0 0 0 0 0
x_[3,t] 0 0 0 0 0 0 0
x_[4,t] 0 0 0 0 0 0 0
x_[5,t] 0 0 0 0 0 0 0
x_[6,t] 0 0 0 0 0 0 0
接下来我们填写正确的字段。 x_[1,1] + x[1, 2] + ... = 1 和 x_[2,1] + x_[2,2] + ... = 1 等等。对于这个问题，使用 for 循环是最简单的
for(i in seq(n)){
cond1[i, seq(i, 30, n)] <- 1
}
cond1[, seq(n + 1)]
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2]
x_[1,t] 1 0 0 0 0 0 1
x_[2,t] 0 1 0 0 0 0 0
x_[3,t] 0 0 1 0 0 0 0
x_[4,t] 0 0 0 1 0 0 0
x_[5,t] 0 0 0 0 1 0 0
x_[6,t] 0 0 0 0 0 1 0
我们仍然需要创建 RHS 并指定方向，但我现在等一下。
所以接下来让我们为第二个条件创建我们的矩阵
sum_i x_[i,t] >= 30 [ for all t ]
这个过程非常相似，但现在每个周期需要一行，所以矩阵的维度是 5x30。这里的主要区别是我们需要插入 x_[i, t]
cond2 <- matrix(0, ncol = t * n,
nrow = t,
dimnames = list(paste0('t=', seq(t)),
names(x)))
for(i in seq(t)){
cond2[i, seq(n) + n * (i - 1)] <- x[seq(n) + n * (i - 1)]
}
cond2[, seq(1, n * t, n)]
x_[1,1] x_[1,2] x_[1,3] x_[1,4] x_[1,5]
t=1 35 0 0 0 0
t=2 0 37 0 0 0
t=3 0 0 39 0 0
t=4 0 0 0 42 0
t=5 0 0 0 0 45
请注意，我正在打印 x_[1, t] 的结果以说明我们做对了。
最后我们有最终条件。为此，我们注意到 ?lpSolve::lp 有一个参数 all.bin，读到这里，它表示
Logical: should all variables be binary? Default: FALSE.
因此，由于所有变量不是 1 就是 0，我们只需将此值设置为 TRUE。在继续之前，让我们将我们的条件组合成一个矩阵
cond <- rbind(cond1, cond2)
现在 RHS 和方向都简单地取自 2 个条件。来自 const.dir 参数
Vector of character strings giving the direction of the constraint: each value should be one of "<," "<=," "=," "==," ">," or ">=". (In each pair the two values are identical.)
在我们的条件中，我们有 6 行代表第一个条件，行代表条件 2。因此我们需要 n (6) 次 == 和 t (5) 次>=.
cond_dir <- c(rep('==', n), rep('>=', t))
RHS 以类似的方式创建
RHS <- c(rep(1, n), rep(30, t))
就是这样！现在我们准备好使用 lpSolve::lp 函数来解决我们的问题。
sol = lpSolve::lp(direction = 'max',
objective.in = x,
const.mat = cond,
const.dir = cond_dir,
const.rhs = RHS,
all.bin = TRUE)
sol$objval
[1] 275
解决方案的权重存储在 sol$solution
names(sol$solution) <- names(x)
sol$solution
x_[1,1] x_[2,1] x_[3,1] x_[4,1] x_[5,1] x_[6,1] x_[1,2] x_[2,2] x_[3,2] x_[4,2] x_[5,2] x_[6,2] x_[1,3] x_[2,3] x_[3,3]
1 0 0 0 0 0 0 0 0 0 1 0 0 0 1
x_[4,3] x_[5,3] x_[6,3] x_[1,4] x_[2,4] x_[3,4] x_[4,4] x_[5,4] x_[6,4] x_[1,5] x_[2,5] x_[3,5] x_[4,5] x_[5,5] x_[6,5]
0 0 0 0 0 0 1 0 0 0 1 0 0 0 1
matrix(sol$solution,
ncol = t,
dimnames = list(rownames(cond1),
rownames(cond2)))
t=1 t=2 t=3 t=4 t=5
x_[1,t] 1 0 0 0 0
x_[2,t] 0 0 0 0 1
x_[3,t] 0 0 1 0 0
x_[4,t] 0 0 0 1 0
x_[5,t] 0 1 0 0 0
x_[6,t] 0 0 0 0 1
我们很快就会发现这是正确的解决方案。 :-)
关于费用的附注
人们可能已经注意到“成本到底去哪儿了？”。在这种特定情况下，成本是固定的并且不是很有趣。这意味着我们可以在计算期间忽略这些，因为我们知道总成本将为 30 * 6 = 180（必须从 objective 值中减去）。然而，成本取决于各种因素并可能影响最佳解决方案的情况并不少见。为了便于说明，我将在此处包括我们如何将成本合并到此示例中。
首先，我们必须扩展我们的 objective 向量以包含每个产品在每个时期的成本
Fixed_C <- -30
x <- c(x, rep(Fixed_C, n * t))
接下来我们将添加一个伪约束
x_[i,t] - C_[i,t] = 0 [for all i, t]
此约束确保如果 x_[i,t] = 1 则将相关成本添加到问题中。有两种方法可以创建此约束。第一个是有一个包含 n * t 行的矩阵，每个行对应一个成本和期间。或者，我们可以使用我们的第一个约束条件，实际上只使用一个约束条件
sum_[i,t] x_[i,t] - C_[i,t] = 0
因为我们的第一个约束确保x[1, 1] != x[1, 2]。所以我们的第三个约束变成
cond3 <- c(rep(1, n * t), rep(-1, n * t))
最后我们必须扩展我们的 RHS 和条件 1 和 2 矩阵。只需将 0 添加到条件矩阵即可使尺寸合适。
cond1 <- cbind(cond1, matrix(0, nrow = n, ncol = n * t))
cond2 <- cbind(cond2, matrix(0, nrow = n, ncol = n * t))
cond <- rbind(cond1, cond2, cond3)
cond_dir <- c(cond_dir, '==')
RHS <- c(RHS, 0)
现在我们可以再次使用lpSolve::lp
solC = lpSolve::lp(direction = 'max',
objective.in = x,
const.mat = cond,
const.dir = cond_dir,
const.rhs = RHS,
all.bin = TRUE)
solC$objval
[1] 95
这等于我们之前的价值 275 减去我们的固定成本 Fixed_C * n = 180。