For a honeycomb structure used for absorbing crash energy and protecting the safety of human or instruments, the bigger the specific energy absorption (SEA) is, the more popular it would be when the peak crushing stress (σp) was retained small enough. In order to improve the energy absorption capacity, crashworthiness optimization for honeycomb structures with various cell specifications are studied in this paper. Detailed numerical models are established for those honeycomb structures by using an explicit finite element method code LS-DYNA. The numerical simulation results are then used as the design samples for constructing metamodels. The optimal Latin hypercube design (OLHD) method is employed for the selection of sampling design points in the design space, and the polynomial functions, radial basis functions (RBF), Kriging, multivariate adaptive regression splines (MARS), and support vector regression (SVR) are utilized to formulate the two optimal objectives SEA and σp. It is found that the polynomial function is the most efficient in constructing the crashworthiness metamodels of honeycombs among the above-mentioned methods. Then, the polynomial function models of SEA and σp are chosen as the surrogate models in the crashworthiness optimization. In order to further validate the polynomial function models, the polynomial function models of SEA and σp are compared with the analytical solutions based on Wierzbicki's theory and Kunimoto and Yamada's theory, respectively. An excellent correlation has been established. As such, the multi-objective particle swarm optimization algorithm (MOPSOA) is applied to obtain the Pareto front of SEA with σp of the honeycomb structures with various cell specifications, which has resulted in a range of optimal designs of honeycomb structures by the multi-objective optimization.