The calculations are done with FindMinimum, using EvaluationMonitor to get the coordinates of the circle centers at every step. The objective function is radius of the circumscribing circle, r0. { x
, y } is the center of circle i. A simplified version of the Mathematica code is n = 10;
Do[r[i] = i^(-1/2), {i, n}];
vars = Join[{r0}, Flatten[Table[{x[i], y[i]}, {i, n}]]];
minimaxcons = Thread[r0 >= Table[r[i] + Sqrt[x[i]^2 + y[i]^2], {i, n}]];
overlapcons = Flatten[Table[(r[i] + r[j] )^2 <= (x[i] - x[j])^2 + (y[i] - y[j])^2, {i, n - 1}, {j, i + 1, n}]];
FindMinimum[{r0, Sequence @@ Join[minimaxcons, overlapcons]}, Thread[{vars, start}]]
A "minimax" approach is used. start is the the starting values of the variables, which I chose to illustrate how a small difference could cause a much bigger difference in the result.
I think this illustrates the difficulties of doing global optimization in a problem with a large number of local minimum. I believe it occurs due to the line search nature of the constrained optimization. However, any other method would be much slower.
