Packing Circles into the Smallest Square

This example formulates a simple nonlinear optimization problem that seeks to pack a number of circles into the smallest possible square box.  The problem is solved with the KNITRO for Mathematica application package.  The KNITRO solver, developed by Ziena Optimization, Inc. (www.ziena.com), is the premier algorithm for large-scale nonlinear optimization.
This note shows how to formulate the problem in Mathematica, and solve it numerically using the KNITRO package. The formulation is easily adapted to similar packing problems.
Suppose we have 20 circles of known size.  The Mathematica listing below shows how 20 radii were generated for this example, ranging in size from 0.1 to 1.2.

n = 20;
    maxRanRad = 1.1;
    r = Table[0.1 + maxRanRad * Random[], {n}];

We wish to find a non-overlapping arrangement of the circles inside a square of minimum area.  The circle centers are unknowns, plus the size of the bounding square.  Define w to be the half-width of a square centered at the origin.  The problem requires inequalities that prevent the circles from overlapping, and inequalities that confine the circles within the square.  Formally, we write

\(TraditionalForm\`\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(min\ \ \ \ \ \ \ w

    	subject\ \ to\ \ \ \ \ \ \ \((x\_i\  - \ x\_j)\)\^2 + \ \((y\_i\  - \ y\_j)\)\^2 ≥ \ \((r\_i\  + \ r\_j)\)\^2\ , \ \ \ i = 1.  . N, \ 
    j = i + 1.  . N

    			x\_i + \ r\_i\  - \ w ≤ \ 
      0, \ \ \ x\_i - \ r\_i\  + \ w ≥ \ 0\ , \ \ \ \ i = 
      1.  . N

    			y\_i + \ r\_i\  - \ w ≤ \ 
      0, \ \ y\_i - \ r\_i\  + \ w ≥ \ 0\ , \ \ \ i = 
      1.  . N

    			\(-B\) ≤ x\_i ≤ B\ , \ \ \(-B\) ≤ 
      y\_i ≤ B\ , \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ i = 
      1.  . N

    			0 ≤ w ≤ B\)\)\)

This optimization problem has 2N+1 unknowns, ((N^2)-N)/2 quadratic inequalities, and 4N  linear inequality constraints.  The constant B is a simple "worst case" estimate of the bounding box, specified as N (0.1+maxRanRad).  We assemble the constraints and call KNITRO to compute a solution.

Here is Mathematica code that constructs the inequality and bound constraints:
    (* ---- CONSTRAINTS THAT PREVENT THE PACKED CIRCLES FROM OVERLAPPING. *)
    conOverlap = Last[Last[Reap[
          Do[ Sow[ (x[i] - x[j])^2 + (y[i] - y[j])^2 >= (r[[i]] + r[[j]])^2 ],
              {i,1,n}, {j,i+1,n}] ] ] ];
    (* ---- CONSTRAINTS THAT KEEP THE CIRCLES PACKED IN THE SQUARE. *)
    conPack = Last[Last[Reap[
          Do[ (Sow[ x[i] + r[[i]] - w <= 0 ];
               Sow[ x[i] - r[[i]] + w >= 0 ];
               Sow[ y[i] + r[[i]] - w <= 0 ];
               Sow[ y[i] - r[[i]] + w >= 0 ]; ), {i,1,n} ] ] ] ];
    (* ---- SIMPLE BOUNDS ON THE VARIABLES. *)
    worstWidth = n * (0.1 + maxRanRad);
    lobnds = Append[ Table[-worstWidth, {2 n}], 0 ];
    upbnds = Append[ Table[ worstWidth, {2 n}], worstWidth ];
    bnds = {lobnds, upbnds};
Here is Mathematica code that calls KNITRO:
    knOptions = {{"multistart",1},{"ms_maxsolves",5}};
    result = KnitroMinimize[{w, c}, bnds, vars, KnitroOptionList->knOptions];

KNITRO returns the minimum square size and the locations of the packed circles.  The plot below visualizes the solution.

The solution must have degeneracies because some small circles have room to "float" within the square without disturbing other circles.  Nevertheless, KNITRO found the solution without difficulty.
All variables in the problem were bounded by the constant B.  KNITRO still solves the problem if the bounds are removed (i.e., if the variables are left free); however, it takes longer to converge.  This is because the default KNITRO algorithm is an interior-point method, and interior-point methods generally perform better if reasonably tight bounds can be specified for the variables. 
The optimization problem formulated in this example may have many local minima; that is, points which satisfy the constraints and minimize the objective function in a local neighborhood of (R^(2  N + 1)).  Of course we want the best local minimum (the global minimum), but global optimization of nonconvex problems can be extremely time consuming when the number of unknowns is large.  KNITRO provides a basic "multi-start" procedure that chooses random initial values for the variables from which KNITRO computes a nearby local minimum.  The options {{"multistart",1},{"ms_maxsolves",5}} enable multi-start and tell it to try 5 different start points.  The best of the 5 local minima is returned.  Note that multi-start performance is also improved if tight bounds can be specified for the variables.

A full listing of the Mathematica code used for this example is available on request from info@ziena.com.  For more information about KNITRO, go to http://www.ziena.com/

Copyright (c) 2006 by Ziena Optimization, Inc.  All rights reserved.