Packing Equal Size Circles in a Triangle

This – not quite simple – example below illustrates both the modeling power of Mathematica and the solver quality of LGO. For a given triangle, we wish to maximize the radius of identical size, non-overlapping circles that all need to be packed in the triangle. Such models will pose a challenge to most nonlinear solvers (in any modeling environment).

Area of the triangle:

area[tripts_List] := (1/2) 
        Det[{tripts[[2]] - 
              tripts[[1]], 
            tripts[[3]] - 
              tripts[[1]]}]

Distance squared from a point {xc,yc} to a line defined by points {x1,y1} and {x2,y2}:

dist2[{x1_, y1_}, {x2_, y2_}, {xc_, 
            yc_}] = ((((xc (-y1 + y2) + x2 (y1 - yc) + x1 (-y2 + yc)))^2)/((((x1 - x2))^2) + (((y1 - y2))^2)));

Inequalities satisified by any point {x,y} inside a triangle made up of 3 points {x1,y1}, {x2,y2} and {x3,y3}

triIneqs[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}, {x_, 
          y_}] :=  
{(x2 (y - y1) + 
                x (y1 - y2) + x1 (-y + y2)) *
            Sign[x3 (-y1 + y2) + x2 (y1 - y3) + 
                x1 (-y2 + y3)] ≤ 
          0, 
(x3 (y - y1) + x (y1 - y3) + 
                x1 (-y + y3)) *
            Sign[x3 (y1 - y2) + x1 (y2 - y3) + 
                x2 (-y1 + y3)] ≤ 
          0, 
 (x3 (-y + y2) + 
                x2 (y - y3) + x (-y2 + y3))*
            Sign[x3 (y1 - y2) + x1 (y2 - y3) + 
                x2 (-y1 + y3)] ≤ 0}

Define now tripts as the vertices of the triangle, and n as the number of circles

packTriangle[tripts_List, n_Integer]:=
Module[{cons1, parpts, cons2, cons3, allcons, xmin, xmax, ymin, 
            ymax},

  
cons1 = 
              Flatten [
                Table[ 2  
                      r  ≤   (Sqrt[(((x[i] - x[j]))^2) + (((y[i] - y[j]))^2)]), {i, n - 1}, {j, i + 1, n}]];

  
partpts = 
              Partition[
                tripts[[{1, 2, 3, 1}]], 
                2, 1];

  
cons2 = 
              Thread[Flatten[(dist2[#[[1]], #[[2]], {x, 
                                y}] &  /@  partpts) /.  
                      Table[{x  ->  x[i], y  ->  y[i]}, {i, 
                          n}]]  ≥  (r^2)];

  
cons3 = 
              Flatten [Table[triIneqs[tripts, {x[i], y[i]}], {i, n}]];

allcons = Join[cons1, cons2, cons3];

xmin = 
              Min[tripts[[All, 1]]];xmax = 
              Max[tripts[[All, 1]]];
ymin = 
              Min[tripts[[All, 2]]];ymax = 
              Max[tripts[[All, 2]]];

vars = 
              Join[Table[{x[i], xmin, xmax}, {i, n}], 
                Table[{y[i], ymin, ymax}, {i, n}], {{r, 
                    0, (Sqrt[area[tripts]/(π n)])}}];

{-r, allcons, vars}
]

<<Graphics`ArgColors`

showTri[tripts_List, n_Integer, 
        soln_] := 
Module[{p1, p2, xmin, xmax, ymin, 
          ymax}, 
xmin = 
          Min[tripts[[All, 1]]]; 
        xmax = Max[
            tripts[[All, 
              1]]]; 
ymin = 
          Min[tripts[[All, 2]]]; 
        ymax = Max[
            tripts[[All, 
              2]]]; 
p1 = Show[#, 
                Axes  ->  True, 
                AspectRatio  ->  ((ymax - ymin)/(xmax - xmin)), 
                PlotRange  ->  {{xmin, xmax}, {ymin, ymax}}, 
                DisplayFunction  ->  Identity] & @ Graphics @ 
              Line[{tripts[[1]], 
                  tripts[[2]], 
                  tripts[[3]], 
                  tripts[[1]]}]; 
p2 = 
          Show[Graphics[
              Table[{ColorCircle[2 π (i - 1)/n], 
                    Disk[{x[i], y[i]}, r]}, {i, n}] /.  
                soln[[2]]], 
            DisplayFunction  ->  Identity]; 
Show[
          p1, p2, DisplayFunction  ->  $DisplayFunction];
]

triPts = {{ .3,  .3}, { .4, 1.1}, {- .3, 1.3}};

The 7-circle model instance is set up by

t7 = packTriangle[triPts, 7];

Next, it is solved by MathOptimizer Professional:

AbsoluteTiming[trisoln = callLGO[Sequence  @@  t7]]

Finally, the solution is visualized:

showTri[triPts, 7, trisoln]

In this example, we found the optimized radius r->0.0921645 in less than 1 second, on an AMD 64 computer running Windows XP Professional, and using Mathematica 5.2. Let us remark that the 7-circle model generated above has 15 decision variables and 21 non-convex (no overlaps) constraints, in addition to 21 + 7 simpler (convex) constraints.