# Packing unequal spheres into minimal cuboid

There are several non-overlapping spheres in 3D. How to find a cuboid Cuboid[{a,b,c}] containing these spheres with minimal a+b+c? Here are my trials with two and three concrete spheres.

(i) The spheres Sphere[{x,y,z},13] ang Sphere[{r,s,t},12] are done. Here is my code

res2 = NMinimize[{a + b + c,  x + 13 <= c && x - 13 >= 0 && y + 13 <= c && y - 13 >= 0 &&
z + 13 <= c && z - 13 >= 0 && r + 12 <= c && r - 12 >= 0 &&
s + 12 <= c && s - 12 >= 0 && t + 12 <= c &&
t - 12 >= 0 && (x - r)^2 + (y - s)^2 + (z - t)^2 >= (13 + 12)^2 &&
a >= b && b >= c}, {a, b, c, x, y, z, r, s, t},
WorkingPrecision -> 20, MaxIterations -> 400]


118.30134496888302115, {a -> 39.433781777141207591, b -> 39.433781662216129186, c -> 39.433781529525684370, x -> 26.433781527577711100, y -> 26.433781527543169115, z -> 26.433781527524494303, r -> 12.000024797809696713, s -> 12.000024797806752675, t -> 12.000024797806996020}}

The verification is as follows.

x + 13 <= a && x - 13 >= 0 && y + 13 <= b && y - 13 >= 0 &&
z + 13 <= c && z - 13 >= 0 && r + 12 <= a && r - 12 >= 0 &&
s + 12 <= b && s - 12 >= 0 && t + 12 <= c && t - 12 >= 0 && a >= b &&
b >= c && (x - r)^2 + (y - s)^2 + (z - t)^2 >= (13 + 12)^2 /. res[[2]]


True

However, the above solution is not optimal Cuboid[{50,26,26}]. How to find the optimal solution with Mathematica?

(ii) The spheres Sphere[{x,y,z},13] ang Sphere[{r,s,t},12] and Sphere[{k,n,h},9/2] are done. The same question for

res3 = NMinimize[{a + b + c, x + 13 <= c && x - 13 >= 0 && y + 13 <= c && y - 13 >= 0 &&
z + 13 <= c && z - 13 >= 0 && r + 12 <= c && r - 12 >= 0 &&
s + 12 <= c && s - 12 >= 0 && t + 12 <= c &&
t - 12 >= 0 && (x - r)^2 + (y - s)^2 + (z - t)^2 >= (13 + 12)^2 &&
k - 9/2 >= 0 && k + 9/2 <= c && h - 9/2 >= 0 && h + 9/2 <= c &&
n - 9/2 >= 0 &&
n + 9/2 <=
c && (x - k)^2 + (y - n)^2 + (z - h)^2 >= (13 +
9/2)^2 && (r - k)^2 + (s - n)^2 + (t - h)^2 >= (12 + 9/2)^2 &&
a >= b && b >= c}, {a, b, c, x, y, z, r, s, t, k, n, h},
WorkingPrecision -> 20, MaxIterations -> 400]


{118.30129311042972531, {a -> 39.433764379149049236, b -> 39.433764371039300303, c -> 39.433764360241375774, x -> 26.433763996833655112, y -> 13.000010247112034705, z -> 13.000010247112188711, r -> 12.000000035837539518, s -> 27.433764329963458967, t -> 27.433764358870261060, k -> 9.1229513952118878652, n -> 10.988012916666727197, h -> 15.679144217401481948}}

I may present two references

Stoyan,Y.G.,Scheithauer,G.& Yaskov,G.N.Packing Unequal Spheres into Various Containers.Cybern Syst Anal 52,419-426 (2016). https://doi.org/10.1007/s10559-016-9842-1

Labib Yousef,Contribution à la résolution des problèmes de placement en trois dimensions,Doctorate University de Picardie Jules Verne (France),2017,172 pages.Downloadable from:https://www.theses.fr/2017AMIE0020.pdf

Unfortunately, the linked article is not accessible to me and I don't read French.

Edit. x+13<=c instead of x+13<=a in res2. Sorry for the typo.

• I'm afraid to even think about 9 spheres. Commented Oct 6, 2023 at 18:30

A generalization of @Ulrich Neumann's code:

findMinimalBox[rList_] :=
(
nSpheres = Length@rList;
coordsList = Transpose[Array[#, nSpheres] & /@ {x, y, z}];
dM = Simplify@
DistanceMatrix[coordsList,
DistanceFunction -> SquaredEuclideanDistance];
rSquareDists = Outer[Plus[#1, #2]^2 &, rList, rList, 1];

eqs = MapThread[GreaterEqual, {dM, rSquareDists}, 2];
p0 = Position[dM, 0, {2}];
notOverlapping = Flatten[Delete[eqs, p0]] // DeleteDuplicates;

spheres = MapThread[Sphere[#1, #2] &, {coordsList, rList}];

rb = RegionBounds[RegionUnion @@ spheres];

abc = Map[#[[2]] - #[[1]] &, rb];

mini = NMinimize[{Total[abc], notOverlapping}, Flatten@coordsList];
abcSoln = abc /. mini[[2]];
graph = Graphics3D[
Thread[{Table[ColorData[3, i - 1], {i, nSpheres}],
spheres /. mini[[2]]}], Axes -> True, Axes -> {x, y, z}];

{abcSoln, graph}
)

findMinimalBox[{12, 13,9/2}]

(*{{26., 49.96, 26.} ,and the graph below}*)


I rotated the graph by hand a little to see the three spheres better:

Trying with 9 spheres (takes about 15 mins on my laptop)

SeedRandom[1234];
rList =  RandomInteger[{5, 15}, 9];
findMinimalBox@rList

(*{34.5274, 33.7461, 30.9396}, and graph below}*)


• Thank's for your support and assistance. Of course I know the generalization but because of the fruitless discussion there was no reason for me to show it. Commented Oct 6, 2023 at 19:56
• Thanks, my original implementation was going to be a lot uglier, but then I saw how you used RegionBounds and that my code a lot better looking. Also it doesn't really matter, but I could also fix one of the spheres at the origin to reduce the number of variables by 3 I guess.
– ydd
Commented Oct 6, 2023 at 20:06
• Just caught an error.  eqs = MapThread[Equal, {dM, rSquareDists}, 2];  was supposed to be  eqs = MapThread[GreaterEqual, {dM, rSquareDists}, 2];  This makes the code much slower, but not unbearable for 9 spheres (15 mins)
– ydd
Commented Oct 6, 2023 at 22:26
• What is your Mathematica version? Commented Oct 12, 2023 at 8:32
• 13.3.1 for Mac OS X ARM
– ydd
Commented Oct 12, 2023 at 15:19

Try RegionBounds to get the enclosing cuboid a,b,c of the two spheres

rb = RegionBounds[RegionUnion[Sphere[{x, y, z}, 13], Sphere[{r, s, t}, 12]]];
abc = Map[#[[2]] - #[[1]] &, rb];
mini = NMinimize[{Total[abc], (x - r)^2 + (y - s)^2 + (z - t)^2 >= (13 +12)^2}, {x, y,z, r, s, t}] // Quiet
abc /. mini[[2]] (* sidelengths of the cuboid*)
(*{26., 49.96, 26.}*)

Graphics3D[{Sphere[{x, y, z}, 13], Sphere[{r, s, t}, 12]} /.mini[[2]], Axes -> True, Axes -> {x, y, z}]


• @user64494 How should your code work with undefined a,b,c? By the way the solution in my answer seems to be close to the optimum Cuboid[{0,0,0},{50,26,26}] which you expects. Commented Oct 6, 2023 at 19:24
• The exact result mini = Minimize[{Total[ abc], (x - r)^2 + (y - s)^2 + (z - t)^2 >= (13 + 12)^2}, {x, y, z, r, s, t}] // Quiet abc /. mini[[2]] performs {77 + Sqrt[623], {x -> -Sqrt[623], y -> 0, z -> 0, r -> 0, s -> 1, t -> 1}} and{25 + Sqrt[623], 26, 26} which is too small. Also how about more than two spheres? Commented Oct 6, 2023 at 19:25
• @user64494 I don't know the exact result, though your reply "result is to small" isn't helpful. My solution idea seems to be not so bad after all. Generalization for more spheres is straight forward (and left up to you) Commented Oct 6, 2023 at 19:38
• I think about it and draw the conclusion that you do find the optimal solution in the case of two spheres, better than my hypothesis. +1. I don't see a possible generalization on more than two spheres. Commented Oct 6, 2023 at 19:52
• @user64494 The connecting line between the centers of the spheres is not parallel to the cuboid sides! Therefore, the optimum must be less than 50! That's not emotional, that's basic geometry in 3D space. I'm out of this fruitless discussion! Commented Oct 6, 2023 at 19:53