I have a list of locations and their weights (calculated distances apart) in a matrix. I would like the optimal solution for each location having 3 connections, minimizing total distance.
costs6 <- matrix(c(0,399671,1525211,990914,1689886,1536081,399671,0,1802419,1128519,1964930,1603803,1525211,1802419,0,814942,164677,943489,990914,1128519.4,814942.7,0,953202,565712,1689886,1964930,164677,953202,0, 1004916,1536081,1603803,943489,565712,1004916,0),ncol=6,byrow=TRUE)
plantcap <- rep(3,6)
citydemand <- rep(3,6)
plant.signs <- rep("=",6)
city.signs <- rep("=",6)
lptrans <- lp.transport(costs6,"min",plant.signs,plantcap,city.signs,citydemand)
lptrans$solution
lptrans
This LP solver returns
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3 0 0 0 0 0
[2,] 0 3 0 0 0 0
[3,] 0 0 3 0 0 0
[4,] 0 0 0 3 0 0
[5,] 0 0 0 0 3 0
[6,] 0 0 0 0 0 3
I am wondering if there is a way to max out any Xij at 1, so that the solver will give me three ones in each column/row, rather than one 3 in each column/row? If not, is there another solver I can use to find the solution?
Something like this, setting it up as an LP problem (assuming a symmetric solution matrix)?
library(lpSolve)
costs6 <- matrix(c(0,399671,1525211,990914,1689886,1536081,
399671,0,1802419,1128519,1964930,1603803,
1525211,1802419,0,814942,164677,943489,
990914,1128519.4,814942.7,0,953202,565712,
1689886,1964930,164677,953202,0, 1004916,
1536081,1603803,943489,565712,1004916,0),ncol=6,byrow=TRUE)
nLoc <- nrow(costs6)
nParams <- sum(1:(nLoc - 1L))
# set up the constraint matrix
# columns are parameters corresponding to the lower triangular of costs6 (by column)
# the first six constraints are for the row/column sums
# the last 15 constraints are for the maximum number of times each path can be used (1)
nConst <- sum(1:nLoc)
mConst <- matrix(0L, nConst, nParams)
mConst[matrix(c(c(combn(1:nLoc, 2)), rep(1:nParams, each = 2)), ncol = 2)] <- 1L
mConst[(nLoc + 1L):nConst,] <- diag(nParams)
lpSol <- lp(
direction = "min",
objective.in = unlist(costs6[lower.tri(costs6)]),
const.mat = mConst,
const.dir = c(rep("=", nLoc), rep("<=", nParams)),
const.rhs = c(rep(3L, nLoc), rep(1L, nParams)),
all.int = TRUE
)
lpSol
#> Success: the objective function is 8688039
# convert the solution to a transport matrix
mSol <- matrix(0, nLoc, nLoc)
mSol[lower.tri(mSol)] <- lpSol$solution
mSol[upper.tri(mSol)] <- t(mSol)[upper.tri(mSol)]
mSol
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0 1 1 1 0 0
#> [2,] 1 0 0 1 1 0
#> [3,] 1 0 0 0 1 1
#> [4,] 1 1 0 0 0 1
#> [5,] 0 1 1 0 0 1
#> [6,] 0 0 1 1 1 0