Как использовать R для решения программы с кусочно-постоянной целевой функцией?
Я хотел бы решить проблему минимизации, используя R, которая имеет кусочно-постоянную целевую функцию. Идея состоит в том, что при более низких значениях моей (целочисленной) переменной решения х возникают более высокие штрафные затраты, чем при более высоких значениях. Я хочу минимизировать общие штрафные расходы, учитывая некоторые ограничения.
Итак, моя программа выглядит следующим образом:
min P(x)
s.t. A x <= b
x >= 0
У меня проблема с программированием кусочно-постоянной целевой функции P(x), где P(x) - сумма всех элементов вектора x. Я знаю, что это не может быть использовано в сочетании с функцией lp() от linprog библиотека. Однако я не могу найти способ сделать это без указания огромного количества дополнительных переменных. Кроме того, обширный поиск в Интернете не дал никаких полезных идей.
Позвольте мне привести пример того, как эта функция P выглядит
[,1] [,2] [,3] [,4] [,5] [,6]
21 11 9 9 0 0
45 28 17 17 6 0
14 0 0 0 0 0
17 11 11 5 5 0
26 11 0 0 0 0
38 18 18 13 0 0
Это следует читать следующим образом: если x1=2, штраф в размере 11 понесется. Если x6 = 4, штраф составляет 13. То есть для x=c(2, ..., 4)у нас есть это P=c(11, ..., 13) и общие штрафные расходы (объективная стоимость) sum(11, ..., 13),
Моя матрица A (она абсолютно унимодулярна) и вектор b выглядят следующим образом.
A <- matrix(c(1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1),nrow=6)
b <- c(4,5,1,5,2,4),
Итак, мой вопрос:
Как я могу найти минимум кусочно-постоянной целевой функции, используя R?
1 ответ
У меня есть два подхода к вашему вопросу. Кажется, ваша проблема в целочисленном программировании, где x1 ... x6 может принимать значение в {1, 2, 3, 4, 5, 6}.
Pmat <- matrix(c(25, 11, 9, 9, 0, 0,
45, 28, 17, 17, 6, 0,
14, 0, 0, 0, 0, 0,
17, 11, 11, 5, 5, 0,
26, 11, 0, 0, 0, 0,
38, 18, 18, 13, 0, 0),
byrow = TRUE, nrow = 6, ncol = 6)
A <- matrix(c(1,0,1,0,0,1,1,0,0,1,0,0,0,1,0,1,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1),nrow=6)
b <- c(4,5,1,5,2,4)
Грубая сила
Поскольку у вас есть 6 таких переменных, есть только 6^6 = 46656 кандидатов. Следовательно, это не займет много времени, чтобы решить это грубой силой; для каждой комбинации проверьте, удовлетворяет ли она ограничению и значению.
d <- expand.grid(x1 = 1:6, x2 = 1:6, x3 = 1:6, x4 = 1:6, x5 =1:6, x6 = 1:6)
condition <- logical(nrow(d))
value <- numeric(nrow(d))
for (i in seq(nrow(d))) {
x <- as.matrix(as.numeric(d[i, 1:6]), ncol = 1)
d$condition[i] <- all(A %*% x <= b)
d$value[i] <- sum(Pmat[cbind(1:6, x)])
}
d <- d[order(d$value, decreasing = FALSE), ]
d <- d[order(d$condition, decreasing = TRUE), ]
print(head(d))
Это получает
# x1 x2 x3 x4 x5 x6 condition value
#7789 1 3 1 1 1 2 TRUE 117
#7783 1 2 1 1 1 2 TRUE 128
#229 1 3 1 2 1 1 TRUE 131
#13 1 3 1 1 1 1 TRUE 137
#223 1 2 1 2 1 1 TRUE 142
#7777 1 1 1 1 1 2 TRUE 145
Таким образом, решение [1, 3, 1, 1, 1, 2] со значением 117.
Формулировка LP
Грубая сила неэффективна, и вы, очевидно, ищете способ решить проблему с помощью решателя LP. Чтобы сформулировать LP для вашей задачи, я бы переопределил набор двоичных переменных, как показано ниже (хотя это может быть то, что вы называете "очень большим количеством дополнительных переменных"):
y_{ij} = 1(x_i = j) for each i = 1,...,6 and j = 1,...,6
Здесь 1() - это функция, которая принимает 1 в том и только в том случае, если утверждение внутри истинно.
Затем векторизовать эти переменные как
y = [y_{11}, y_{12}, ..., y_{16}, y_{21}, ..., y_{66}]
Это новая управляющая переменная для LP. Обратите внимание, что размер y 36
Затем вам нужно векторизовать P матрица в вектор весов для y, Вам также необходимо конвертировать матрицу A, Например, если исходное ограничение
x_1 + x_2 <= b
Тогда это эквивалентно
y_{11} + 2*y_{12} + 3*y_{13} + 4*y_{14} + 5*y_{15} + 6*y_{16}
+ y_{21} + 2*y_{22} + 3*y_{23} + 4*y_{24} + 5*y_{25} + 6*y_{26} <= b
Продукт Kronecker удобен для этого преобразования.
Наконец, для каждого iдолжен быть только один j такой, что y_{ij} = 1, Следовательно, вам нужно
y_{i1} + y_{i2} + ... + y_{i6} = 1 for each i.
Продукт Kronecker также полезен для этого.
library(lpSolve)
Pvec <- as.numeric(t(Pmat))
A1 <- kronecker(A, matrix(1:6, nrow = 1))
b1 <- b
A2 <- kronecker(diag(6), matrix(rep(1, 6), nrow = 1))
b2 <- rep(1, 6)
o <- lp(direction = "min",
objective.in = Pvec, const.mat = rbind(A1, A2),
const.dir = c(rep("<=", 6), rep("=", 6)), const.rhs = c(b1, b2),
all.bin = TRUE)
# show the variable names
tmp <- expand.grid(1:6, 1:6)
vname <- paste("x", tmp[,2], tmp[,1], sep = "")
names(o$solution) <- vname
print(o)
print(o$solution)
Это дает то же решение, что и грубая сила.
Success: the objective function is 117
x11 x12 x13 x14 x15 x16 x21 x22 x23 x24 x25 x26 x31 x32 x33 x34 x35 x36
1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0
x41 x42 x43 x44 x45 x46 x51 x52 x53 x54 x55 x56 x61 x62 x63 x64 x65 x66
1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0