The performance of metasurfaces measured experimentally often discords with expected values from numerical optimization. These discrepancies are attributed to the poor tolerance of metasurface building blocks with respect to fabrication uncertainties and nanoscale imperfections. Quantifying their efficiency drop according to geometry variation are crucial to improve the range of application of this technology. Here, we present a novel optimization methodology to account for the manufacturing errors related to metasurface designs. In this approach, accurate results using probabilistic surrogate models are used to reduce the number of costly numerical simulations. We employ our procedure to optimize the classical beam steering metasurface made of cylindrical nanopillars. Our numerical results yield a design that is twice more robust compared to the deterministic case.