One geometric remark (which is reminiscent of a symmetrization
principle used in isoperimetric problems): a minimizing polycube
may be chosen to lie in the first octant, with three sides tangent
to the coordinate planes, and so that the number of cubes in each column is
monotonic in each coordinate, like a lozenge tiling:
The point is that you can choose some direction to be the
surface of the earth, and let gravity pull all of the cubes
down. This does not increase the surface area: the result
ends up with at most one bottom and top face in each column.
In adjacent columns of height $m$ and $n$, the number of
exposed faces is at least $|m-n|$, which is minimized when
they are sitting flush. Repeat this in each direction, until
the configuration is monotonic (like a 3-dimensional Young diagram,
but I don't know if there's a mathematical label).
Added: This argument gives an interpretation for the
desired minimization problem. One has a sequence of Young diagrams
in three directions such that the total number of squares
adds up to $n$, and which are monotonic in the number of
squares in each diagram (with some compatibility conditions).
Let the number of squares in the young diagram adjacent to
each coordinate plane be $x_1, y_1, z_1$ respectively. Then
the surface area is $2(x_1+y_1+z_1)$. This is counting twice
the number of blue, red, and green squares (excluding the
corners) in the picture (pictured are the top squares, to
which there is a corresponding bottom square in each column).
We have $x_1\geq x_2 \geq \cdots \geq x_k$, and $x_1+x_2+\cdots+x_k=n$,
and corresponding formulae for $y_i$ and $z_i$ counting
the number of squares in each tableau parallel to a given
coordinate plane. It shouldn't be too hard to solve this
combinatorial optimization problem, but I don't have time
to work it out now.
Addendum 2: I think I can describe the optimal solution, but
I don't have time to write down the complete argument now.
It is helpful to consider the 2-dimensional case of Young diagrams.
A length $n$ Young diagram is a sequence $x_1\geq x_2 \geq \cdots \geq x_{y_1}$, such that $\sum x_j=n$.
The dual diagram is $y_1 \geq y_2 \geq \cdots y_{x_1}$, where $y_i=\max\{ j | x_j\geq i\}$. The length of the boundary is $2(x_1+y_1)$. Then the symmetrization principle says that the minimal "polysquare" is a Young diagram. Let $P(2k)$ be the maximal size of a diagram with perimeter $2k$ (the "isoperimetric function"). Since a Young diagram fits inside a rectangle of sides $x_1, y_1$ with the same perimeter, clearly $P(2k) = \max_{a+b=k} ab = \lfloor \frac{k}{2} \rfloor \lceil \frac{k}{2}\rceil$. Also, the Young diagram which minimizes perimeter for a given $n$ clearly has $P(2(x_1+y_1-1))< n \leq P(2(x_1+y_1))$. Let $k=x_1+y_1$. Then the minimal perimeter will be realized by a Young diagram which, if $k$ is even, is a $k/2(k/2-1)$ rectangle with a row of $n-k(k-2)/4$ squares added on (so fitting in a $(k/2)^2$ square and containing a $k/2(k/2-1)$ rectangle). If $k$ is odd, then the Young diagram will be obtained from an $((k-1)/2)^2$ square with a row of $n-(k-1)^2/4$ squares added (so containing a $(k-1)/2\times (k-1)/2$ square, and contained in a rectangle of size $(k-1)/2\times (k+1)/2$). Note: The perimeter minimizer is far from unique. For example, if $n=8$, then the minimizer will be realized by the Young sequences $3,3,2$ and $4,4$. We are only describing a preferred choice of minimizers.
Now we may use this description of perimeter minimizing Young diagrams to describe minimal 3-D Young diagrams. As described above, we have a monotonic sequence of Young diagrams, each fitting inside the other, of sizes $x_1\geq x_2 \geq \cdots \geq x_k$. Then one checks that the perimeter is $2x_1$ plus twice the sums of the perimeters of the Young diagrams. Replace each Young diagram with a perimeter-minimizing one described above (one may observe that the choices are nested for $x_i \geq x_{i+1}$). Repeat this in the $y-$ and $z-$ directions to obtain a 3-D Young diagram in which the cross sections are perimeter minimizing 2-D Young diagrams. Then the 3-D Young diagram fits into an $a\times b \times c$ box, with $a\leq b\leq c$, and $|c-a|=1$ (so $b= a$ or $c$). A similar isoperimetric argument to the above implies that it contains an $a^3$ cube and is contained in a $c^3$ cube, where $a=\lfloor n^{1/3} \rfloor$. If $ n < a^2(a+1)$, then it is obtained from an $a^3$ cube with a $n-a^3$ optimal Young diagram attached to one face. If $a^2(a+1)\leq n < a(a+1)^2$, then it is obtained from an $a^2\times (a+1)$ box by attaching an optimal $n-a^2(a+1)$ Young diagram to an $a\times (a+1)$ face. And if $a(a+1)^2 \leq n$, then it is obtained from an $a\times (a+1)^2$ box by attaching an $n-a(a+1)^2$ Young diagram to an $(a+1)^2$ face. Again, these are not the unique perimeter minimizers, but a description of optimal ones. It's possible that I've made a mistake, but at least the symmetrization principles underlying the argument seem solid.
