Частично связан с другим моим вопросом здесь .
В моем случае «первоначальной» целью было выбрать n = 50 объектов из N = 292, так чтобы сумма всехпарные расстояния между выбранными объектами максимальны (максимальная сумма или сумма p-дисперсии).
Благодаря пользователям, предоставившим рекомендации, я немного углубился в чтение, и теперь я понимаю, что проблема действительно квадратичная в самом простомформа, и решатель, такой как CPLEX, может решить ее.
Однако эта статья Кьюби указывает, что результаты maxsum не гарантируют, что не будет объектов, очень близких кдруг с другом;и действительно, из некоторых тестов, которые я проводил методом грубой силы в смоделированных небольших случаях, я обнаружил, что решения с очень высоким максимальным значением иногда содержат очень близкие объекты.
Так что теперь я думаю, что p-дисперсия (maxmin)подход может быть более подходящим для того, чего я хочу достичь.Это также изначально квадратичная проблема.
Поскольку у меня еще нет CPLEX, я не могу попробовать квадратичную формулировку, поэтому я посмотрел на подходы линеаризации.Эти 2 статьи кажутся мне довольно интересными:
Франко, Учоа
Сая, 2015
Последняя указывает на другую статью, которую я нахожу оченьтоже интересно:
Пизингер, 2006
Следующим моим шагом было опробовать следующее:
- линеаризованная p-дисперсия по Кьюби / Эркуту, с N двоичными переменными для объектов и 1 непрерывной переменной для максимального минимального расстояния, ограниченного между наименьшим и наибольшим расстоянием в матрице расстояний
- грубой силой, перечисляя все комбинации из n объектов из N и находятот, который имеет наибольшее минимальное расстояние
- как 1, но устанавливает более жесткую верхнюю границу для непрерывной переменной, используя метод Сая / Пизингера
- линеаризованной p-дисперсии по Сая, с N-двоичнымпеременные для объектов и до N * (N-1) / 2 дополнительных двоичных переменных для парных расстояний
Я не пытался ужесточить нижнюю границуили добавить еще неравенство, потому что методы, предложенные в статьях, выходят за рамки моего уровня математики.
Меня удивляет то, что метод 4, который должен быть «компактным», на самом деле имеет огромное количествобинарные переменные и последующие ограничения, и в тестах, которые я запускал, он работал намного хуже, чем методы 1 и 2. Затягивание верхней границы, с другой стороны, имело огромный эффект, и на самом деле метод 2 на данный момент является единственным, который кажетсяуметь справляться с большими проблемами в разумные сроки.
Но это правда, что я не реализовал точно метод в статье Сая, поэтому, возможно, мои наблюдения не верны.
Вопросы : что вы думаете о различных методах линеаризации, описанных в этих статьях?Можете ли вы предложить лучшие?Считаете ли вы, что сохранение максимального минимального расстояния в качестве непрерывной переменной, как в формулировке Кьюби, лучше, чем ее «квантование», как в формулировке Сая?
Фактически за это время возникли дополнительные сложности и события, например, наличие«принудительных» объектов и необходимости использовать оценки для каждого объекта, но я хотел бы сначала рассмотреть вышеупомянутые.
Я вставил ниже код R, который я использовал для тестирования этого.
Спасибо!
#Test of linearized methods for the solution of p-dispersion (maxmin) problems
#-----------------------------------------------------------------------------
#Definitions
#Given N objects, whose distance matrix 'distmat' is available:
#p-dispersion (maxmin): select n (n >= 2, n < N) objects such that the minimal distance between any two objects is maximised
#p-dispersion sum (maxsum): select n (n >= 2, n < N) objects such that the sum of all the pairwise distances between them is maximised
#Literature
#Kuby, 1987: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1538-4632.1987.tb00133.x
#Pisinger, 1999: https://pdfs.semanticscholar.org/1eb3/810077c0af9d46ed5ff2b0819d954c97dcae.pdf
#Pisinger, 2006: http://yalma.fime.uanl.mx/~roger/work/teaching/clase_tso/docs_project/problems/PDP/cor-2006-Pisinger.pdf
#Franco, Uchoa: https://pdfs.semanticscholar.org/4092/d2c98cdb46d5d625a580bac08fcddc4c1e60.pdf
#Sayah, 2015: https://download.uni-mainz.de/RePEc/pdf/Discussion_Paper_1517.pdf
#Initialization
require(Matrix)
if (length(find.package(package="Rsymphony",quiet=TRUE))==0) install.packages("Rsymphony")
require(Rsymphony)
par(mfrow = c(2,2))
#0. Choose N, n and which methods to run
N = 20
n = ceiling(0.17*N)
run_PD_Erkut = TRUE
run_PD_brute_force = TRUE
run_PD_Erkut_UB_Sayah = TRUE
run_PD_Sayah = TRUE
#1. Make random distance matrix for testing
set.seed(1)
coords <- cbind(runif(N,-5,5),runif(N,-5,5))
distmat <- t(as.matrix(dist(coords,diag=T)))
distmat[lower.tri(distmat)] <- 0
distmat <- Matrix(distmat,sparse=T)
N.i <- NROW(distmat)
colnames(distmat) <- paste("j",1:N.i,sep="_")
rownames(distmat) <- paste("i",1:N.i,sep="_")
#2. Make a 2D representation of the points using classic multidimensional scaling
cmds <- cmdscale(as.dist(t(distmat)))
#3. Link the pairwise distances to the rows and columns of the distmat
distmat_summary <- summary(distmat)
N.ij <- NROW(distmat_summary)
distmat_summary["ID"] <- 1:(N.ij)
i.mat <- xtabs(~ID+i,distmat_summary,sparse=T)
j.mat <- xtabs(~ID+j,distmat_summary,sparse=T)
ij.mat <- cbind(i.mat,0)+cbind(0,j.mat)
colnames(ij.mat)[[N.i]] <- as.character(N.i)
zij.mat <- .sparseDiagonal(n=N.ij,x=1)
#4. MaxMin task by Kuby/Erkut (N binary variables + 1 continuous variable for max Dmin)
if (run_PD_Erkut == TRUE) {
#4a. Building the constraint matrix (mat), direction (dir), right-hand-side (rhs) and objective (obj) for the LP task
dij <- distmat_summary$x
M <- max(dij)
m <- min(dij)
#Erkut's condition: for each i,j i<j, D (min distance to maximise) + M*xi + M*xj <= 2*M + dij
constr.dij <- cbind("D"=1,ij.mat*M)
dir.dij <- rep("<=",N.ij)
rhs.dij <- 2*M+dij
constr.D <- c(1,rep(0,N.i))
dir.DM <- "<="
rhs.DM <- M
dir.Dm <- ">="
rhs.Dm <- m
#constraining the total number of objects to be n
constr.n <- c(0,rep(1,N.i))
dir.n <- "=="
rhs.n <- n
#assembling the constraints
mat <- rbind(constr.n,constr.dij,constr.D,constr.D)
dir <- c(dir.n,dir.dij,dir.DM,dir.Dm)
rhs <- c(rhs.n,rhs.dij,rhs.DM,rhs.Dm)
#objective
obj <- setNames(c(1,rep(0,N.i)), c("D",colnames(ij.mat)))
#4.b. Solution
st <- system.time(LP.sol <- Rsymphony_solve_LP(obj,mat,dir,rhs,types=c("C",rep("B",N.i)),max=TRUE,verbosity = -2, time_limit = 5*60))
ij.sol <- names(obj[-1])[as.logical(LP.sol$solution[-1])]
items.sol <- rownames(distmat)[as.numeric(ij.sol)]
Dmin <- LP.sol$solution[1]
#4.c. Plotting the results
plot(cmds,main=paste(c("p-dispersion (Erkut), N =",N,", n =",n,"\nUB =",round(M,2),", time =",round(st[3],2),"s, Dmin =",round(Dmin,2)),collapse=" ") )
points(cmds[as.numeric(ij.sol),],pch=16,col="red")
text(cmds[as.numeric(ij.sol),],ij.sol,cex=0.9,col="red",adj=c(0,1))
}
#5. MaxMin task by brute force
if (run_PD_brute_force == TRUE) {
if (choose(N,n) <= 200000) {
st <- system.time({combs <- as.data.frame(t(combn(N,n)))
combs["maxmin"] <- apply(combs, 1, function(x) {min(distmat_summary[(distmat_summary$j %in% x) & (distmat_summary$i %in% x),"x"])})
combs["maxsum"] <- apply(combs, 1, function(x) {sum(distmat_summary[(distmat_summary$j %in% x) & (distmat_summary$i %in% x),"x"])})
combs_maxmin_max <- combs[combs$maxmin == max(combs$maxmin),][1,]})
ij.sol <- as.character(combs_maxmin_max[,1:n])
items.sol <- rownames(distmat)[as.numeric(ij.sol)]
Dmin <- combs_maxmin_max[1,"maxmin"]
plot(cmds,main=paste(c("p-dispersion (brute force), N =",N,", n =",n,"\ntime =",round(st[3],2),"s, Dmin =",round(Dmin,2)),collapse=" ") )
points(cmds[as.numeric(ij.sol),],pch=16,col="red")
text(cmds[as.numeric(ij.sol),],ij.sol,cex=0.9,col="red",adj=c(0,1))
}
}
#6. MaxMin task by Erkut with Sayah's upper bound
if (run_PD_Erkut_UB_Sayah == TRUE) {
#6a. Building the constraint matrix (mat), direction (dir), right-hand-side (rhs) and objective (obj) for the LP task
m <- min(distmat_summary$x)
M <- sort(sapply(1:(N.i), function(it) {min((sort(distmat_summary[(distmat_summary$i == it) | (distmat_summary$j == it),"x"],decreasing = TRUE)[1:(n-1)]))}),decreasing=TRUE)[n]
#Erkut's condition: for each i,j i<j, D (min distance to maximise) + M*xi + M*xj <= 2*M + dij
constr.dij <- cbind("D"=1,ij.mat*M)
dir.dij <- rep("<=",N.ij)
rhs.dij <- 2*M+dij
constr.D <- c(1,rep(0,N.i))
dir.DM <- "<="
rhs.DM <- M
dir.Dm <- ">="
rhs.Dm <- m
#constraining the total number of objects to be n
constr.n <- c(0,rep(1,N.i))
dir.n <- "=="
rhs.n <- n
#assembling the constraints
mat <- rbind(constr.n,constr.dij,constr.D,constr.D)
dir <- c(dir.n,dir.dij,dir.DM,dir.Dm)
rhs <- c(rhs.n,rhs.dij,rhs.DM,rhs.Dm)
#objective
obj <- setNames(c(1,rep(0,N.i)), c("D",colnames(ij.mat)))
#6.b. Solution
st <- system.time(LP.sol <- Rsymphony_solve_LP(obj,mat,dir,rhs,types=c("C",rep("B",N.i)),max=TRUE,verbosity = -2, time_limit = 5*60))
ij.sol <- names(obj[-1])[as.logical(LP.sol$solution[-1])]
items.sol <- rownames(distmat)[as.numeric(ij.sol)]
Dmin <- LP.sol$solution[1]
#6.c. Plotting the results
plot(cmds,main=paste(c("p-dispersion (Erkut, UB by Sayah), N =",N,", n =",n,"\nUB =",round(M,2),", time =",round(st[3],2),"s, Dmin =",round(Dmin,2)),collapse=" ") )
points(cmds[as.numeric(ij.sol),],pch=16,col="red")
text(cmds[as.numeric(ij.sol),],ij.sol,cex=0.9,col="red",adj=c(0,1))
}
#7. MaxMin task by Sayah (N binary variables + binary variables from unique values of dij)
if (run_PD_Sayah == TRUE) {
#7a. Building the constraint matrix (mat), direction (dir), right-hand-side (rhs) and objective (obj) for the LP task
#7a.1. Finding the upper (M) and lower (m) bound for the minimal distance
m <- min(distmat_summary$x)
M <- sort(sapply(1:(N.i), function(it) {min((sort(distmat_summary[(distmat_summary$i == it) | (distmat_summary$j == it),"x"],decreasing = TRUE)[1:(n-1)]))}),decreasing=TRUE)[n]
dijs <- unique(sort(distmat_summary$x))
dijs <- dijs[dijs <= M]
N.dijs <- length(dijs)
z.mat <- .sparseDiagonal(N.dijs,1)
#Sayah's formulation:
#applying z[k] <= z[k-1]
constr.z <- cbind(rep(0,N.i*(N.dijs-1)),cbind(0,z.mat[-1,-1])-z.mat[-NROW(z.mat),])
dir.z <- rep("<=",N.dijs-1)
rhs.z <- rep(0,N.dijs-1)
#applying x[i]+x[j]+z[k] <= 2
constr.ijk <- NULL
for (k in 2:N.dijs) {
IDs <- distmat_summary[distmat_summary$x < dijs[k],"ID"]
constr.ijk <- rbind(constr.ijk,cbind(ij.mat[IDs,,drop=F],z.mat[rep(k,length(IDs)),,drop=F]))
}
dir.ijk <- rep("<=",NROW(constr.ijk))
rhs.ijk <- rep(2,NROW(constr.ijk))
#constraining the total number of objects to be n
constr.n <- c(rep(1,N.i),rep(0,N.dijs))
dir.n <- "=="
rhs.n <- n
#assembling the constraints
mat <- rbind(constr.n,constr.z,constr.ijk)
dir <- c(dir.n,dir.z,dir.ijk)
rhs <- c(rhs.n,rhs.z,rhs.ijk)
#objective
obj <- setNames(c(rep(0,N.i),1,diff(dijs)), c(colnames(ij.mat),paste("z",1:N.dijs,sep="_")))
#7.b. Solution
st <- system.time(LP.sol <- Rsymphony_solve_LP(obj,mat,dir,rhs,types="B",max=TRUE,verbosity = -2, time_limit = 5*60))
ij.sol <- names(obj[1:N.i])[as.logical(LP.sol$solution[1:N.i])]
items.sol <- rownames(distmat)[as.numeric(ij.sol)]
Dmin <- sum(LP.sol$solution[(1+N.i):(N.dijs+N.i)]*obj[(1+N.i):(N.dijs+N.i)])
#7.c. Plotting the results
plot(cmds,main=paste(c("p-dispersion (Sayah), N =",N,", n =",n,"\nUB =",round(M,2),", time =",round(st[3],2),"s, Dmin =",round(Dmin,2)),collapse=" ") )
points(cmds[as.numeric(ij.sol),],pch=16,col="red")
text(cmds[as.numeric(ij.sol),],ij.sol,cex=0.9,col="red",adj=c(0,1))
}