# Maximal circle packing inside a given square

How to pack n equal disks(radius a/2) with maximal area inside a given rectangle ?

I'm trying to solve this problem using Mathematicas Region functions.

Here my approach:

Maximize the area of the union of all disks regKand require all disks to lie inside the square RegionWithin[rect,regK].

sys[a_?NumericQ, pts_ ?(MatrixQ[#, NumericQ] & )]:=Block[{
rect = Rectangle[{-1/2, -1/2}, {1/2, 1/2}],
regK = RegionUnion[Map[Disk[#, a/2] &, pts]]}
,
{ Area[regK],
Join[{RegionWithin[rect, regK] , 0 < a < 1/2},
Map[-1/2 < # < 1/2 &, Flatten[ pts]]]}]


Unfortunately

n=3;  (*number of circles*)
NMaximize[ sys[a, Table[{x[i], y[i]}, {i, 1, n}]]
, Join[{a}, Flatten[Table[{x[i], y[i]}, {i, 1,n}]]]]
(*NMaximize::nnum: The function value {-0.835066,{-False,-False,-False,-True,-False,-True,-False,-False}}
is not a number at
{a,x[1],x[2],x[3],y[1],y[2],y[3]} = {0.677613,0.614816,0.51572,0.66946,-0.0824906,0.251149,0.6107}.*)


NMaximize stops evaluation with an error message.

What's wrong with my code? Thanks!

• The formulation is not clear to me: is a fixed in "radius a/2 maximal"? Commented Apr 19, 2021 at 10:18
• I'm looking for maximal radius a/2. Commented Apr 19, 2021 at 10:20
• See that Maple application as an example (Preview of its as PDF file is available.). Good luck! Commented Apr 19, 2021 at 10:40
• @user64494 Thanks, very interesting! Commented Apr 19, 2021 at 10:45
• Commented Apr 19, 2021 at 14:43

Edit

It seems that RegionWithin is better than the original method.

reg = Rectangle[];
n = 5;
sol = NMaximize[{r, r > 0,
Table[RegionWithin[reg, Disk[{x[i], y[i]}, r]], {i, n}],
Table[SignedRegionDistance[RegionBoundary@reg]@{{x[i], y[i]}} >=
r, {i, n}], Table[{x[i], y[i]} \[Element] reg, {i, 1, n}],
Table[EuclideanDistance[{x[i], y[i]}, {x[j], y[j]}] >= 2 r, {i,
n}, {j, i - 1}]}, {r, Table[{x[i], y[i]}, {i, n}]} // Flatten]
Graphics[{EdgeForm[Cyan], FaceForm[],
Rectangle[], {FaceForm[Red], Table[Disk[{x[i], y[i]}, r], {i, n}]} /.
sol[[2]]}]


{0.207107, {r -> 0.207107, x[1] -> 0.207107, y[1] -> 0.792893, x[2] -> 0.207107, y[2] -> 0.207107, x[3] -> 0.792893, y[3] -> 0.792893, x[4] -> 0.792893, y[4] -> 0.207107, x[5] -> 0.5, y[5] -> 0.5}}

Original

For the $$n$$ points $$p_i,i=1\cdots n$$, we set $$d(p_i,p_j)\geq 2r,i\not=j$$ and all the distance to the boundary of region $$d(p_i,Boundary)\geq r$$

reg = Rectangle[];
n = 5;
sol = NMaximize[{r,
SignedRegionDistance[RegionBoundary@reg] /@
Table[{x[i], y[i]}, {i, n}] >= r,
Table[{x[i], y[i]} ∈ reg, {i, 1, n}],
Table[EuclideanDistance[{x[i], y[i]}, {x[j], y[j]}] >= 2 r, {i,
n}, {j, i - 1}]}, {r, Table[x[i], {i, n}], Table[y[i], {i, n}]} //
Flatten]
Graphics[{EdgeForm[Cyan], FaceForm[],
Rectangle[], {FaceForm[Red],
Table[Disk[{x[i], y[i]}, r], {i, n}]} /. sol[[2]]}]


• Thank you for your helpful answer. I tried your solution for n=3 and got a "L-shaped" result which isn't optimal. Commented Apr 19, 2021 at 11:50
• @UlrichNeumann I am using version 12.2 in Win and Linux. Commented Apr 19, 2021 at 11:58
• The conditions Flatten[{0 < r < 1, Table[{r <= x[i] <= 1 - r, r <= y[i] <= 1 - r}, {i, 1, n}], Table[Total[({x[i], y[i]} - {x[j], y[j]})^2] >= (2 r)^2, {i, n}, {j, i - 1}]}] runs faster and gives a larger radius for n == 5 but a smaller radius for n == 6. Commented Apr 19, 2021 at 12:02
• @cvgmt My version 12.2/Win gives a different "optimal" solution, don't know why. Your modified answer shows what I would expect as an optimum! Commented Apr 19, 2021 at 12:03
• @ChipHurst Thanks, I test in my Linux , both of the two codes get the same radius. For n=6, r -> 0.187681, for n=5, r -> 0.196438 Commented Apr 19, 2021 at 12:14

Not an answer, but to long for a comment:

Here is "my" result evaluated with Mathematica v12.2 and the approach from @cvgmt

reg = Rectangle[];
n = 3;

sol = NMaximize[{r,
SignedRegionDistance[RegionBoundary@reg] /@
Table[{x[i], y[i]}, {i, n}] >= r,
Table[{x[i], y[i]} \[Element] reg, {i, 1, n}],
Table[EuclideanDistance[{x[i], y[i]}, {x[j], y[j]}] >= 2 r, {i,
n}, {j, i - 1}]}, {r, Table[x[i], {i, n}], Table[y[i], {i, n}]} //
Flatten]
(*{0.25, {r -> 0.25, x[1] -> 0.75, x[2] -> 0.749409, x[3] -> 0.25,
y[1] -> 0.25, y[2] -> 0.75, y[3] -> 0.250385}}*)

Graphics[{EdgeForm[Cyan], FaceForm[],
Rectangle[], {FaceForm[Red], Table[Disk[{x[i], y[i]}, r], {i, n}]} /.
sol[[2]]}]


Obviously the result r=0.25 isn't optimal and deviates from @cvgmt answer. Don't know why!

The following optimum is only evaluated, if an additional constraint x[1] == r,y[1] == 1 - r  is added

reg = Rectangle[];
n = 8;
sol = NMaximize[{r,
Map[SignedRegionDistance[RegionBoundary[reg]],
Table[{x[i], y[i]}, {i, n}]] >= r,
Table[{x[i], y[i]} \[Element] reg, {i, 1, n}],
Table[EuclideanDistance[{x[i], y[i]}, {x[j], y[j]}] >= 2 r, {i, n}, {j, i -1}]
, x[1] == r,y[1] == 1 - r }, {r, Table[x[i], {i, n}], Table[y[i], {i, n}]} //Flatten,
Method -> {Automatic, "DifferentialEvolution", "RandomSearch","SimulatedAnnealing"; "NelderMead"}[[2]]]
(*{0.17054, {r -> 0.17054, x[1] -> 0.17054, x[2] -> 0.5, x[3] -> 0.5,x[4] -> 0.829459, x[5] -> 0.258819, x[6] -> 0.17054,x[7] -> 0.829459, x[8] -> 0.741182,y[1] -> 0.82946,y[2] -> 0.741181, y[3] -> 0.258817, y[4] -> 0.170541,y[5] -> 0.499999, y[6] -> 0.17054, y[7] -> 0.829459, y[8] -> 0.5}}*)

Graphics[{EdgeForm[Cyan], FaceForm[],Rectangle[], {FaceForm[Red], Table[Disk[{x[i], y[i]}, r], {i, n}]} /.sol[[2]]}]


• Try Method->"DifferentialEvolution" option. Commented Apr 19, 2021 at 12:50
• Thanks, "DifferentialEvolution" and "SimulatedAnnealing" evaluate the "right" optimum! Commented Apr 19, 2021 at 12:57