In this work, surrogate models based on support vector regression (SVR) of a multi-resonant unit cell in a geometrical 4-D parallelotope domain are trained and used in a reflectarray antenna design. The multiple sharp resonances of the unit cell prevent a suitable training process in the whole orthotope defined by the available degrees of freedom (DoF). Thus, a strategy to improve the training process and obtain highly accurate models is devised. It consists in defining a parallelotope around a rectangle of stability, which is in turn defined at a lower dimensionality. The SVR models with four geometrical DoF obtained in this parallelotope are shown to provide highly accurate results for the design of a large contoured-beam reflectarray for space applications. The direct optimization with the surrogate models allows to improve the cross-polarization performance several dB while considerably increasing computational performance. Furthermore, compared to lower dimensionality models, the 4-D models offer better results when applied to wideband and dual-band reflectarray direct optimization.