Maximum area of projection of cuboid onto plane

Question

For didactic purposes I solve with Mathematica the following problem (see Vasilyev's box on p. 17): when the area of the projection of a rotating cuboid in the three-dimensional space with coordinates x,y,zonto the xOy plane is maximum? (The problem has the exact symbolic answer, but I don't know an accurate proof of that.) I consider a concrete cuboid 1 x 3 x 2 since Mathematica (and other CASes) does not operate with abstract geometric objects. I rotate its vertices by EulerMatrix. Then I find their projections, cutting the third coordinate, and find the area of the convex hull of the projected vertices. Here is my code.

ClearAll["Global`*"];
f[\[Alpha]_,\[Beta]_,\[Gamma]_]:= Area[ConvexHullMesh[Table[(EulerMatrix[{\[Alpha],\[Beta],\[Gamma]}].x)[[1;;2]], 
{x, {{0, 0, 0}, {1, 0, 0}, {0, 3, 0}, {1, 3, 0}, {0, 0, 2}, {1, 0, 2}, {0, 3, 2}, {1, 3, 2}}}]]]
f[-Pi/4, Pi/3, Pi/6]
(*6.86603*)
NMaximize[{f[a, b, c],a >= -Pi && a <= Pi && b >= 0 && b <= Pi && c >= -Pi && 
c <= Pi}, {a, b, c}, Method -> "DifferentialEvolution"]
(*{7.,{a->-2.55692,b->1.12789,c->-0.321751}}*)

and a lot of error communications. To be sure,

NMaximize[{f[a, b, c],a >= -Pi && a <= Pi && b >= 0 && b <= Pi && c >= -Pi && 
c <= Pi}, {a, b, c}, Method -> "RandomSearch"]
(*{7., {a -> -1.73782, b -> 2.01371, c -> 0.321751}}*)

(a lot of error communications as a bonus) and

NMinimize[{f[a, b, c],a >= -Pi && a <= Pi && b >= 0 && b <= Pi && c >= -Pi && 
c <= Pi}, {a, b, c}, Method -> "RandomSearch"]
(*{2., {a -> -2.01037, b -> 1.5708, c -> 1.5708}}*)

The question arises: is it possible to obtain the exact solution with help of Mathematica?

Answer 1

There seems to be some bug in "Area". If I write my own area function, I get at least a result without crashing the kernel:

Note, the first rotation about the z axis has no effect on the area in the x/y plane and need not be maximized.

area[a1_, a2_ /; NumericQ[a2], a3_] := Module[{},
  cub = {{0, 0, 0}, {1, 0, 0}, {0, 3, 0}, {1, 3, 0}, {0, 0, 2}, {1, 0,
      2}, {0, 3, 2}, {1, 3, 2}};
  heron[p1_, p2_, p3_] := 
   Module[{n1 = Norm[p1 - p2], n2 = Norm[p2 - p3], n3 = Norm[p3 - p1],
      s},
    s = Total[n1 + n2 + n3]/2;
    Sqrt[s (s - n1) (s - n2) (s - n3)]
    ];
  pts = ConvexHullMesh[
      Transpose[EulerMatrix[{a1, a2, a3}].Transpose@cub][[All, 
        1 ;; 2]]]["BoundaryPolygons"][[1, 1]];
  heron[First[pts], ##] & @@@ Partition[Rest@pts, 2, 1] // Total
  ]
NMaximize[area[0, a2, a3], {a2, a3}]