obj <- function(x, g, np){
  mean((tapply(g, x, sum) - sum(g) / np)^2)
}

gg <- sample(1:30, 2e3, TRUE)
numPack <- 20

idx <- matrix(NA, numPack, length(gg) / numPack)
for (i in 1:ncol(idx)) {
  if (i %% 2 == 0) {
    idx[ , i] <- numPack:1
  } else {
    idx[ , i] <- 1:numPack
  }
}

obj(as.vector(idx), sort(gg), numPack) # 10432.6

tapply(sort(gg), as.vector(idx), sum)
#    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19   20
# 1532 1532 1533 1533 1532 1535 1536 1533 1533 1534 1536 1535 1535 1532 1532 1532 1532 1531 1531 1531
sum(gg) / numPack # 1533