First of all you can compute the convex hull of the meshes so to drastically reduce the number of computed points. Indeed, the minimal bounding box of the convex hull is guaranteed to be the same than the target mesh.
There are ways to compute the optimal volume for 1 mesh by computing the diameter of the point set of the convex hull but AFAIK this cannot be easily transposed to N meshes. The rotating calipers can be used to compute the diameter of a point set (the supporting parallel lines of the antipodal pair can be used to form the bounding box).
One simple way to address this problem is by using basic optimizations methods (AFAIK there is no exact algorithm to solve this in linear time nor any simple exact algorithm).
Here is an example of metric function to minimize (total volume):
function metric(α, β, θ)
R_α <- RotationMatrixX(α)
R_β <- RotationMatrixY(β)
R_θ <- RotationMatrixZ(θ)
R <- R_α * R_β * R_θ
sum <- 0
for each mesh M:
K <- Points(PrecomputedConvexHull(M)) * R
B <- BoundingBox(K)
sum <- sum + volume(B)
return sum
α, β and θ are the rotation angles of the X, Y and Z axis to optimize so to minimize the total volume. There are many optimization methods that can be used to minimize this metric function. For example, you can use the L-BFGS algorithm (it should be relatively good since the metric function should be "smooth").