Visualizing optimization process

One of the approaches to graph drawing is application of so called force-directed algorithms. In its simplest form the idea is to layout the nodes on plane so that all edges in the graph have approximately equal length. This problem has very intuitive visualization so it is a nice case for showing how different optimization algorithms behave in high dimensions.
I want to position several balls (8 in the example below) on a plane in such a way that distance between all balls is approximately 0.5. Precisely my objective is to minimize the highest deviance from 0.5 of distance between any pair of balls.

The interesting feature of this problem is that visualization of 8 balls on a plane is simple and intuitive - it is important because it allows me to visually assess the solution quality - not having to rely on objective function value only.
On the other hand the optimization involves 16 decision variables for 8 balls and the objective function is not smooth and has many local minima - a perfect case to compare different optimization algorithms.

The code that generates the plots for BFGS optimization algorithm is given below:

library(animation)
b.init <- function(balls) {
    best.fit <<- +Inf
    pos <- runif(2 * balls)
    names(pos) <- paste(c("x", "y"), rep(1:balls, each = 2),
                  sep = "")
    return(pos)
}

b.plot <- function(pos) {
    par(mar = c(2, 2, 1.5, 0.5))
    # extract x and y coordinates of balls
    x <- pos[c(T, F)]
    y <- pos[c(F, T)]
    # plot balls
    plot(x, y, cex = 5, pch = 20, main = b.fit(pos),
         xlim = c(min(x) - 0.1, max(x) + 0.1),
         ylim = c(min(y) - 0.1, max(y) + 0.1))
    # and lines connecting them
    for (i in 1:(length(x) - 1)) {
        for (j in (i + 1):length(x)) {
            lines(x[c(i, j)], y[c(i, j)])
        }
    }
}

b.fit <- function(pos) {
    # our goal is to make distance between all balls equal to 0.5
    max(abs(dist(t(matrix(pos, 2))) - 0.5))
}

f <- function(pos) {
    fit <- b.fit(pos)
    # update plot if improvement is found
    if (fit < best.fit) {
        best.fit <<- fit
        b.plot(pos)
    }
    return(fit)
}

run <- function() {
    set.seed(1)
    pos <- b.init(8)
    for (i in 1:20) {
        b.plot(pos)
    }
    pos <- optim(pos, f, method = "BFGS", control = list(maxit = 1000))$par
    for (i in 1:20) {
        b.plot(pos)
    }
}

ani.options(interval = 0.05)
saveGIF(run())

It generates the following animation:

It is interesting to run this simulation for different initial positions of balls as BFGS might converge to different local minima or to try different optimization algorithms either from optim function (like SANN) or any other optiization routine (function f keeps track of best solution found independently from optimization procedure).