Это классическая c проблема, и ее необходимо переформулировать.
Начните с переформулирования вашей проблемы
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 и так далее. Использование для l oop является самым простым для этой проблемы.
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, и, читая это, он сообщает
Логический: все ли переменные должны быть двоичными? По умолчанию: FALSE.
Итак, поскольку все переменные равны 1 или 0, мы просто устанавливаем это значение на TRUE. Прежде чем продолжить, давайте объединим наши условия в одну матрицу
cond <- rbind(cond1, cond2)
Теперь и правая часть, и направление просто взяты из двух условий. Из документации по аргументу const.dir
Вектор символьных строк, задающих направление ограничения: каждое значение должно быть одним из "<," "<=," "=," "= =, "">, или "> =". (В каждой паре два значения идентичны.)
В наших условиях у нас есть 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
Что мы быстро видим, правильное решение. : -)
Замечание по затратам
Можно было заметить «где, черт возьми, затраты go?». В этом конкретном случае c затраты фиксированы и не очень интересны. Это означает, что мы можем игнорировать их во время расчетов, потому что мы знаем, что общая стоимость будет 30 * 6 = 180 (которую нужно вычесть из объективного значения). Однако нередко затраты зависят от различных факторов и могут повлиять на оптимальное решение. Для иллюстрации я расскажу, как мы могли бы учесть затраты в этом примере. Сначала нам нужно расширить наш целевой вектор, чтобы включить в него затраты для каждого продукта за каждый период
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, то соответствующая стоимость добавляется к проблеме. Есть 2 способа создать это ограничение. Первый - это матрица с 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.