r - p 色散(maxmin)问题的最佳线性化？

Question

与我的另一个问题部分相关here .

在我的例子中，“最初”的目标是从 N=292 个对象中选择 n=50 个对象，这样所选对象之间所有成对距离的总和就会最大化(最大和或 p 分散和)。

感谢提供建议的用户，我做了一些进一步的阅读，现在我明白这个问题确实是最简单形式的二次问题，并且像 CPLEX 这样的求解器可能能够解决它。

然而，this article by Kuby 指出，maxsum 结果并不能保证不会有彼此非常接近的对象；事实上，从我通过蛮力对模拟的较小情况进行的一些测试中，我发现具有非常高最大值和的解决方案有时包含非常接近的对象。

所以现在我认为 p 色散 (maxmin) 方法可能更适合我想要实现的目标。这本来也是一个二次问题。

由于我还没有 CPLEX，我无法尝试二次公式，所以我研究了线性化方法。这两篇文章对我来说似乎很有趣:
Franco, Uchoa
Sayah, 2015

后者指向另一篇文章，我也觉得很有趣:
Pisinger, 2006

我的下一步是尝试以下操作:

1. 根据 Kuby/Erkut 线性化 p 色散，对象有 N 个二进制变量，最大最小距离有 1 个连续变量，在距离矩阵中的最小距离和最大距离之间有界
2. 蛮力，从 N 中枚举 n 个对象的所有组合，并找到最小距离最大的一个
3. 类似于 1，但使用 Sayah/Pisinger 的方法为连续变量设置了更严格的上限
4. 根据 Sayah 的线性化 p 色散，对象有 N 个二元变量，成对距离有多达 N*(N-1)/2 个额外的二元变量

我没有尝试收紧下限或添加更多不等式，因为文章中建议的方法超出了我的数学水平。

令我困惑的是，本应是“紧凑”的方法 4 实际上具有大量二进制变量和随之而来的约束，并且在我运行的测试中，它的性能比方法 1 和 2 差得多。拧紧另一方面，上限产生了巨大的影响，事实上，目前方法 2 似乎是唯一能够在合理时间内解决大型问题的方法。
但我确实没有完全实现 Sayah 论文中的方法，所以我的观察可能是不正确的。

问题:您如何看待这些文章中描述的各种线性化方法？你能推荐更好的吗？您是否认为像 Kuby 的公式那样将最大最小距离保持为连续变量比像 Sayah 的公式那样将其“量化”更好？

事实上，在此期间出现了进一步的并发症和发展，例如“强制”对象的存在以及对每个对象使用分数的需要，但我想先解决上述问题。

我在下面粘贴了我用来测试它的 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))

}

Answer 1

您没有提及是否可以容忍非最佳解决方案。但是您应该能够做到，因为您不能指望能够普遍找到该问题的最佳解决方案。在这种情况下，有一个因子为 2 的近似值。

Let V be the set of nodes/objects
Let i and j be two nodes at maximum distance
Let p be the number of objects to choose
p = set([i,j])
while size(P)<p:
  Find a node v in V-P such that min_{v' in P} dist(v,v') is maximum
  \That is: find the node with the greatest minimum distance to the set P
  P = P.union(v)
Output P

此近似算法保证找到一个值不超过最优值两倍的解，并且除非 P=NP，否则没有多项式时间启发式算法可以提供更好的性能保证。

最优界在 White (1991) 中得到证明和 Ravi et al. (1994) .后者证明启发式是最好的。

作为引用，我针对 p=50、n=400 运行了完整的 MIP。 6000 秒后，最优性差距仍为 568%。近似算法用了 0.47s 获得了 100%(或更小)的最优性差距。

近似算法的 Python(抱歉，我不在 R 中建模)表示如下:

#!/usr/bin/env python3

import numpy as np

p = 50
N = 400

print("Building distance matrix...")
d = np.random.rand(N,N) #Random matrix
d = (d + d.T)/2             #Make the matrix symmetric

print("Finding initial edge...")
maxdist  = 0
bestpair = ()
for i in range(N):
  for j in range(i+1,N):
    if d[i,j]>maxdist:
      maxdist = d[i,j]
      bestpair = (i,j)

P = set()
P.add(bestpair[0])
P.add(bestpair[1])

print("Finding optimal set...")
while len(P)<p:
  print("P size = {0}".format(len(P)))
  maxdist = 0
  vbest = None
  for v in range(N):
    if v in P:
      continue
    for vprime in P:
      if d[v,vprime]>maxdist:
        maxdist = d[v,vprime]
        vbest   = v
  P.add(vbest)

print(P)

而 Gurobi Python 表示可能如下所示:

#!/usr/bin/env python
import numpy as np
import gurobipy as grb

p = 50
N = 400

print("Building distance matrix...")
d = np.random.rand(N,N) #Random matrix
d = (d + d.T)/2             #Make the matrix symmetric

m = grb.Model(name="MIP Model")

used  = [m.addVar(vtype=grb.GRB.BINARY) for i in range(N)]

objective = grb.quicksum( d[i,j]*used[i]*used[j] for i in range(0,N) for j in range(i+1,N) )

m.addConstr(
  lhs=grb.quicksum(used),
  sense=grb.GRB.EQUAL,
  rhs=p
)

# for maximization
m.ModelSense = grb.GRB.MAXIMIZE
m.setObjective(objective)

# m.Params.TimeLimit = 3*60

# solving with Glpk
ret = m.optimize()