Я хочу минимизировать следующее уравнение:
F=SUM{u 1:20}sum{w 1:10} Quw(ruw-yuw)
со следующими ограничениями:
yuw >= yu,w+1
yuw >= yu-1,w
y20,0 >= 100
y0,10 >= 0
У меня есть матрица 20*10 ruw и 20*10 quw, теперь мне нужно создать матрицу yuw, которая придерживается ограничений. Я пишу код на R и знаком с пакетами lpsolve и optimx, но не знаю, как их использовать для этого конкретного вопроса.
Поскольку Quw и ruw являются данными, все ограничения, а также цель являются линейными в yuw переменных решения. В результате это задача линейного программирования, к которой можно подойти с помощью пакета lpSolve.
Чтобы немного абстрагироваться, предположим, что R=20 и C=10 описывают размерность входных матриц. Затем есть R*C переменных решения, и мы можем присвоить им порядок y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC, считывая столбцы матрицы переменных.
Коэффициент каждой переменной yuw в задаче равен -Quw; обратите внимание, что Quw*ruw членов в суммировании являются константами (то есть не зависят от значений, которые мы выбираем для переменных решения) и, следовательно, не вводятся в решатель линейного программирования. Интересно, что это означает, что ruw фактически не влияет на решение модели оптимизации.
Первые R*(C-1) ограничений соответствуют yuw >= yu,w+1 ограничениям, а следующие (R-1)*C ограничениям соответствуют yuw >= yu-1,w ограничениям. Последние два ограничения соответствуют ограничениям y20,1 >= 100 и y1,10 >= 0.
Мы можем ввести эту модель в пакет lpsolve с помощью следующей команды R, взяв в качестве входных данных простую матрицу Q, где каждая запись равна -1 (результирующее решение должно иметь все переменные решения, установленные на 0, кроме нижнего левого угла, который должен быть равен 100). ):
# Sample data
Quw <- matrix(-1, nrow=20, ncol=10)
R <- nrow(Quw)
C <- ncol(Quw)
# Build constraint matrix
part1 <- matrix(0, nrow=R*(C-1), ncol=R*C)
part1[cbind(1:(R*C-R), 1:(R*C-R))] <- 1
part1[cbind(1:(R*C-R), (R+1):(R*C))] <- -1
part2 <- matrix(0, nrow=(R-1)*C, ncol=R*C)
pos2 <- as.vector(sapply(2:R, function(r) r+seq(0, R*(C-1), by=R)))
part2[cbind(1:nrow(part2), pos2)] <- 1
part2[cbind(1:nrow(part2), pos2-1)] <- -1
part3 <- rep(0, R*C)
part3[R] <- 1
part4 <- rep(0, R*C)
part4[(C-1)*R + 1] <- 1
const.mat <- rbind(part1, part2, part3, part4)
library(lpSolve)
mod <- lp(direction = "min",
objective.in = as.vector(-Quw),
const.mat = const.mat,
const.dir = rep(">=", nrow(const.mat)),
const.rhs = c(rep(0, nrow(const.mat)-2), 100, 0))
Теперь мы можем получить доступ к решению модели:
mod
# Success: the objective function is 100
matrix(mod$solution, nrow=R)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] 0 0 0 0 0 0 0 0 0 0
# [2,] 0 0 0 0 0 0 0 0 0 0
# [3,] 0 0 0 0 0 0 0 0 0 0
# [4,] 0 0 0 0 0 0 0 0 0 0
# [5,] 0 0 0 0 0 0 0 0 0 0
# [6,] 0 0 0 0 0 0 0 0 0 0
# [7,] 0 0 0 0 0 0 0 0 0 0
# [8,] 0 0 0 0 0 0 0 0 0 0
# [9,] 0 0 0 0 0 0 0 0 0 0
# [10,] 0 0 0 0 0 0 0 0 0 0
# [11,] 0 0 0 0 0 0 0 0 0 0
# [12,] 0 0 0 0 0 0 0 0 0 0
# [13,] 0 0 0 0 0 0 0 0 0 0
# [14,] 0 0 0 0 0 0 0 0 0 0
# [15,] 0 0 0 0 0 0 0 0 0 0
# [16,] 0 0 0 0 0 0 0 0 0 0
# [17,] 0 0 0 0 0 0 0 0 0 0
# [18,] 0 0 0 0 0 0 0 0 0 0
# [19,] 0 0 0 0 0 0 0 0 0 0
# [20,] 100 0 0 0 0 0 0 0 0 0
Обратите внимание, что ваша модель может легко стать неосуществимой, если Quw изменится (например, если мы заполним ее 1 вместо -1). В этих случаях модель выйдет со статусом 3 (это можно увидеть, запустив mod или mod$status).
ruw. Вы что-то упустили в своей формулировке, например. абсолютное значение в цели? - person josliber♦ schedule 20.06.2015