如何在R中求解线性规划模型 [英] How to solve linear programming model in R
问题描述
我需要解决以下微观经济问题:
I need to solve the following microeconomic problem:
- 我有五年(2011年至2015年)可以生产的六种资产(资产1-6).
- 每个资产只能在一年内生产.
- 每个资产都必须在我的五年期内生产.
- 生产不是相互排斥的；我一年可以生产一种以上的商品，而不会影响任何一种的生产.
- 每项资产的固定生产成本等于30. 我每年必须有非负利润；收入必须至少为30.
- I have six assets I can produce (asset 1 - 6) across five years (2011 - 2015).
- Each asset can only be produced during one year.
- Each asset must be produced in my five year period.
- Production is not mutually exclusive; I can produce more than one good in a year without affecting the production of either.
- Each asset has a fixed cost of production equal to 30.
- I must have non-negative profit in each year; revenues must be at least 30.
下面是一个矩阵，表示我在给定年份(j)中生产每种资产(i)的潜在收入.
Below is a matrix representing my potential revenue for producing each asset (i) in a given year (j).
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*
星号( * )代表最佳解决方案集.
The asterisks (*) represent what should be the optimal solution set.
在概述的约束下，我如何使用R来解决使我的收入(以及利润)最大化的生产计划.我的输出应该是类似的6x5矩阵，分别是 0 和 1 ，其中 1 表示选择在商品中生产商品.给定年份.
How can I use R to solve for the production plan that maximizes my revenue (and therefore profit) subject to the constraints outlined. My output should be a similar 6x5 matrix of 0's and 1's, where 1's represent choosing to produce a good in a given year.
推荐答案
这是一个经典问题，需要重新制定.
This is a classic problem, and one that needs to be reformulated.
从重新制定问题开始
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] 的向量开始.为方便起见，让我们为其命名(尽管未使用)，只是为了便于跟踪.
In the lpSolve package the maximization problem is given in a linear representation, eg. in non-matrix format. Lets start by making a vector representing our x_[i,t]. For ease let's name it (although this is not used), just so we can keep track.
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
现在，我们需要创建条件.从第一个条件开始
Now now we need to create our conditions. Starting with the first condition
sum_t x_[i,t] = 1 [ for all i ]
我们可以相当简单地创建它.需要注意的是，尺寸必须正确.我们有一个长度为30的向量，因此我们需要我们的条件矩阵具有30列.此外，我们有6个资产，因此对于这种情况我们将需要6行.再次让我们命名行和列以跟踪自己.
we can create this rather simply. The thing to watch out for here, is that the dimension has to be right. We have a vector of length 30, so we'll need our conditions matrix to have 30 columns. In addition we have 6 assets, so we'll need 6 rows for this condition. Again lets name the rows and columns to keep track ourself.
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循环是解决此问题的最简单方法
Next we fill our the correct fields. x_[1,1] + x[1, 2] + ... = 1 and x_[2,1] + x_[2,2] + ... = 1 and so forth. Using a for loop is the simplest for this problem
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并指定方向，但我现在将等待它.
所以接下来让我们为第二个条件创建矩阵
We still have to create the RHS and specify direction but I'll wait with this for now.
So next lets create our matrix for the second condition
sum_i x_[i,t] >= 30 [ for all t ]
此过程非常相似，但是现在每个周期都需要一行，因此矩阵的尺寸为5x30.这里的主要区别是，我们需要插入 x_ [i，t]
The process for this one is very similar, but now we need a row for each period, so the dimension of the matrix is 5x30. The main difference here, is we need to insert the values of 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 ，并阅读该内容，指出
Note that I'm printing the result for x_[1, t] to illustrate we've got it right.
Last we have the final condition. For this we note the ?lpSolve::lp has an argument all.bin, and reading this, it states
逻辑:所有变量都应为二进制吗?默认值:假.
Logical: should all variables be binary? Default: FALSE.
因此，由于所有变量均为1或0，我们只需将此值设置为 TRUE .在继续之前，让我们将条件合并到一个矩阵中
So since all variables are either 1 or 0, we simply set this value to TRUE. Before continuing lets combine our conditions into one matrix
cond <- rbind(cond1, cond2)
现在，RHS和方向都简单地从这两个条件中得出.从 const.dir 参数
Now both the RHS and the direction are simply taken from the 2 conditions. From the documentation on the const.dir argument
给出约束方向的字符串向量:每个值应该是"<"，-"或-"中的一个.<=，""=""=="，"，或"=".(在每对中，两个值是相同的.)
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)次> = .
In our conditions we have 6 rows representing the first condition, and rows represeting condition 2. Thus we need n (6) times == and t (5) times >=.
cond_dir <- c(rep('==', n), rep('>=', t))
RHS是用类似的方式创建的
The RHS is created in a similar fashion
RHS <- c(rep(1, n), rep(30, t))
就是这样！现在，我们准备使用 lpSolve :: lp 函数解决我们的问题.
And that's it! Now we're ready to solve our problem using the lpSolve::lp function.
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
The weights for the solution are stored in 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
我们很快就会发现这是正确的解决方案.:-)
Which we quickly see is the correct solution. :-)
一个人可能已经注意到费用到底去了哪里?".在这种特定情况下，成本是固定的，不是很有趣.这意味着我们可以在计算期间忽略这些，因为我们知道总成本将为 30 * 6 = 180 (必须减去目标值).但是，成本取决于各种因素并可能影响最佳解决方案的情况并不少见.为了说明，我将在此示例中包括如何合并成本.
首先，我们必须扩展目标向量，以合并每个时期每个产品的成本
One may have noticed "Where the hell did the costs go?". In this specific case, costs are fixed and not very interesting. This means we can ignore these during the calculations because we know the total cost is going to be 30 * 6 = 180 (which has to be substracted from the objective value). However it is not uncommon that costs depend on various factors, and might affect the optimal solution. For illustration, I'll include how we could incorporate costs in this example here.
First we'll have to extend our objective vector to incorporate the costs for each product at each period
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 行，每个成本和期间对应一个.另外，我们可以使用第一个约束，而实际上只与一个常数生活在一起
This constraint ensures that if x_[i,t] = 1 then the relevant cost is added to the problem. There's 2 ways to create this constraint. The first is to have a matrix with n * t rows, one for each cost and period. Alternatively we can use our first constraint and actually live with only a single constrant
sum_[i,t] x_[i,t] - C_[i,t] = 0
因为我们的第一个约束条件确保 x [1，1]！= x [1，2] .所以我们的第三个约束变成
because our first constraint makes sure x[1, 1] != x[1, 2]. So our third constraint becomes
cond3 <- c(rep(1, n * t), rep(-1, n * t))
最后，我们必须扩展RHS以及条件1和2矩阵.只需在条件矩阵上加0即可使尺寸合适.
Lastly we have to extend our RHS and condition 1 and 2 matrices. Simply add 0's to the condition matrices to make the dimensions fit.
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 .
which is equal to our previous value 275 minus our fixed costs Fixed_C * n = 180.
这篇关于如何在R中求解线性规划模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
