In this paper, we apply a semi-analytical method and numerical simulation to study the hydrodynamic behavior of large arrays of point-absorber wave energy converters (WECs) and the multi-objective optimization of the arrays layout using particle swam optimization (PSO) algorithm based on surrogate models. To investigate wave interactions among multiple WECs among an array, a semi-analytical model is developed based on the linear potential flow theory and the matched eigen-function expansions method. The fluid domain is divided into two kinds of regions, interior regions underneath the cylinders and an exterior region surrounding all the cylinders. The matched eigen-function expansions method is used to solve the radiation potential problem in each domain, and the hydrodynamic coefficients and motion response of the cylinders in the array are evaluated. To validate the accuracy of our semi-analytical method, the hydrodynamic software WAMIT is adopted to simulate the wave energy park numerically. Results obtained by the semi-analytical method are compared with the numerical simulation results. The hydrodynamic characteristics and power absorption performance of the WECs within the wave energy park are discussed. The power performance of a wave energy park is studied as functions of global geometry, incident wave direction and separating distance between WECs, respectively. Finally, a multi-objective optimization algorithm based on a surrogate model is used to optimize the layout of wave energy array. A surrogate model is established based on radial basis function (RBF) to improve the computing efficiency. With the power output and power smoothing of wave energy array as two optimization objectives, the safety of cylinders and the occupied sea area of wave energy park as constraint conditions, the multi-objective particle swarm optimization (MOPSO) was used to search for the Pareto-optimal solution, and the feasible optimization scheme of wave energy array was given.