Maximal circle packing inside a given square

Question

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!

Answer 1

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]]}]

Answer 2

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!

addendum case n=8

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]]}]