This paper investigates the use of constrained surrogate models to solve the multi-design optimization problem of a flexible hydrofoil. The surrogate-based optimization (EGO) substitutes the complex objective function of the problem by an easily evaluable model, constructed from a limited number of computations at carefully selected design points. Associated with ad-hoc statistical strategies to propose optimum candidates within the estimated feasible domain, EGO enables the resolution of complex optimization problems. In this work, we rely on Gaussian processes (GP) to model the objective function and adopt a probabilistic classification method to treat non-explicit inequality constraints and non-explicit representation of the feasible domain. This procedure is applied to the design of the shape and the elastic characteristics of a hydrofoil equipped with deformable elements providing flexibility to the trailing edge. The optimization concerns the minimization of the hydrofoil drag while ensuring a non-cavitating flow, at selected sailing conditions (boat speed and lifting force). The drag value and cavitation criterion are determined by solving a two-dimensional nonlinear fluid-structure interaction problem, based on a static vortex lattice method with viscous boundary layer equations, for the flow, and a nonlinear elasticity solver for the deformations of the elastic components of the foil. We compare the optimized flexible hydrofoil with a rigid foil geometrically optimized for the same sailing conditions. This comparison highlights the hydrodynamical advantages brought by the flexibility: a reduction of the drag over a large range of boat speeds, less susceptibility to cavitation and a smaller angle of attack tuning range.