A novel model-order reduction technique for the solution of the fine-scale equilibrium problem appearing in computational homogenization is presented. The reduced set of empirical shape functions is obtained using a partitioned version that accounts for the elastic/inelastic character of the solution - of the Proper Orthogonal Decomposition (POD). On the other hand, it is shown that the standard approach of replacing the nonaffine term by an interpolant constructed using only POD modes leads to ill-posed formulations. We demonstrate that this ill-posedness can be avoided by enriching the approximation space with the span of the gradient of the empirical shape functions. Furthermore, interpolation points are chosen guided, not only by accuracy requirements, but also by stability considerations. The approach is assessed in the homogenization of a highly complex porous metal material. Computed results show that computational complexity is independent of the size and geometrical complexity of the Representative Volume Element. The speedup factor is over three orders of magnitude - as compared with finite element analysis - whereas the maximum error in stresses is less than 10%.