Проблема, которую я пытаюсь решить, в основном та же, что и в этом посте:
https://stats.stackexchange.com/questions/339935/python-library-for-combinatorial-optimization
И моя текущая реализация действительно использует оптимизатор на основе генетического алгоритма.
Тем не менее, я хотел бы решить ее как бинарную задачу линейного программирования (по крайней мере, попытаться, хотя, по-видимому, это 'NP-hard').
Мой вопрос заключается в том, как сформулировать LPнаилучшим образом, потому что я не уверен, что делаю это правильно.
Ниже приведена упрощенная версия того, с чем я имею дело, которая, однако, точно показывает, в чем проблема.
Мы создаем m * n (в данном случае 6) объектов комбинаторным процессом, беря m (3) объектов типа 'R1' (скажем, {A, B, C}) и n (2) объектов типа 'R2' (скажем {X, Y}).
6 объектов {AX, AY, BX, BY, CX, CY} оцениваются, и каждый получает оценку D, в данном случае {0.8,0.7,0.5,0.9,0.4,0.0}, в этом порядке.
CL <- cbind(expand.grid(R2=LETTERS[24:25],R1=LETTERS[1:3],stringsAsFactors = FALSE),D=c(0.8,0.7,0.5,0.9,0.4,0.0))
Теперь мы хотим выбрать 2 различных R1 и 1 R2 так, чтобы сумма D была максимальной.
В этом примере ответом является R1 = {A, B}, R2 = {Y}.
Однако нельзя прийти к такому выводу, принимая, например, 2 R1 и R2 ссамый высокий средний D.
Это будет работать для R1, но не для R2:
aggregate(D~R1,CL,mean)
# R1 D
#1 A 0.75
#2 B 0.70
#3 C 0.20
aggregate(D~R2,CL,mean)
# R2 D
#1 X 0.5666667
#2 Y 0.5333333
Я знаю, как сформулировать это как задачу линейного программирования;только я не уверен, что моя формулировка эффективна, потому что в основном это приводит к проблеме с m n + m + n переменными и 2 (m + n) +2 ограничениями.
Основная трудность заключается в том, чтоМне нужно как-то подсчитать количество выбранных R1 и R2, и я не знаю никакого способа сделать это, кроме того, что я покажу ниже (и также описано в моем другом посте здесь ).
Вот что я бы сделал:
CL["Entry"] <- seq_len(dim(CL)[[1]])
R1.mat <- table(CL$R1,CL$Entry)
R2.mat <- table(CL$R2,CL$Entry)
N_R1 <- dim(R1.mat)[[1]]
N_R2 <- dim(R2.mat)[[1]]
N_Entry <- dim(CL)[[1]]
constr.mat <- NULL
dir <- NULL
rhs <- NULL
constr.mat <- rbind(constr.mat,cbind(R1.mat,-diag(table(CL$R1)),matrix(0,N_R1,N_R2)))
dir <- c(dir,rep("<=",N_R1))
rhs <- c(rhs,rep(0,N_R1))
constr.mat <- rbind(constr.mat,cbind(R2.mat,matrix(0,N_R2,N_R1),-diag(table(CL$R2))))
dir <- c(dir,rep("<=",N_R2))
rhs <- c(rhs,rep(0,N_R2))
constr.mat <- rbind(constr.mat,constr.mat)
dir <- c(dir,rep(">=",N_R1+N_R2))
rhs <- c(rhs,1-table(CL$R1),1-table(CL$R2))
constr.mat <- rbind(constr.mat,c(rep(0,N_Entry),rep(1,N_R1),rep(0,N_R2)))
dir <- c(dir,"==")
rhs <- c(rhs,2)
constr.mat <- rbind(constr.mat,c(rep(0,N_Entry),rep(0,N_R1),rep(1,N_R2)))
dir <- c(dir,"==")
rhs <- c(rhs,1)
obj <- c(aggregate(D~Entry,CL,c)[["D"]],rep(0,N_R1+N_R2))
Что можно решить, например, с помощью lpSolve:
sol <- lp("max", obj, constr.mat, dir, rhs, all.bin = TRUE,num.bin.solns = 1, use.rw=FALSE, transpose.constr=TRUE)
sol$solution
#[1] 0 1 0 1 0 0 1 1 0 0 1
, показывающего, что продукты {AY, BY}были выбраны, соответствующие R1 = {A, B} и R2 = {Y}:
CL[as.logical(sol$solution[1:N_Entry]),]
# R2 R1 D Entry
#2 Y A 0.7 2
#4 Y B 0.9 4
Я обнаружил, что при больших проблемах lpSolve застревает на века;Rsymphony, кажется, работал лучше.
Но опять же, мой главный вопрос: эффективен ли этот способ формулировки LP?Должен ли я сделать это по-другому?
Спасибо!
РЕДАКТИРОВАТЬ
А пока, работая над несколько связанной проблемой, я обнаружил , что может быть достаточно только одного набора ограничений, если добавить «затраты» (в данном примере, отрицательные) к целевому вектору для «различных переменных R1 и R2».
Здесь вместо:
obj <- c(aggregate(D~Entry,CL,c)[["D"]],rep(0,N_R1+N_R2))
Я бы сделал:
obj <- c(aggregate(D~Entry,CL,c)[["D"]],rep(-1,N_R1+N_R2))
Это сделало бы ненужными ограничения m + n .
Это все еще остается огромной проблемой для решения,даже для относительно небольших m, n , поэтому, если кто-нибудь может посоветовать, как это сделать лучше ...
Я посмотрел на lp.transport, но это будет ограничено двумя измерениями (т.е. толькоR1 и R2, а не R1, R2, R3, например), и я не думаю, что вы можете ограничить количество различных объектов на категорию в этом типе решателя.