Abstract. Wind farm layout optimization is usually subjected to boundary constraints of irregular shapes. The analytical expressions of these shapes are rarely available, and consequently, it can be challenging to include them in the mathematical formulation of the problem. This paper presents a new methodology to integrate multiple disconnected and irregular domain boundaries in wind farm layout optimization problems. The method relies on the analytical gradients of the distances between wind turbine locations and boundaries, which are represented by polygons. This parameterized representation of boundary locations allows for a continuous optimization formulation. A limitation of the method, if combined with gradient-based solvers, is that wind turbines are placed within the nearest polygons when the optimization is started in order to satisfy the boundary constraints, thus the allocation of wind turbines per polygon is highly dependent on the initial guess. To overcome this and improve the quality of the solutions, two independent strategies are proposed. A study case is presented to demonstrate the applicability of the method and the proposed strategies. In this study, a wind farm layout is optimized in order to maximize the AEP in a non-uniform wind resource site. The problem is constrained by the minimum distance between wind turbines and five irregular polygon boundaries, defined as inclusion zones. Initial guesses are used to instantiate the optimization problem, which is solved following three independent approaches: (1) a baseline approach that uses a gradient-based solver, (2) approach 1 combined with the relaxation of the boundaries, which allows for a better design space exploration, and (3) the application of a heuristic algorithm, smart-start, prior to the gradient-based optimization, improving the allocation of wind turbines within the inclusion polygons based on the potential wind resource and the available area. The results show that the relaxation of boundaries combined with a gradient-based solver achieves on average +10.2 % of AEP over the baseline, whilst the smart-start algorithm, combined with a gradient-based solver, finds on average +20.5 % of AEP with respect to the baseline and +9.4 % of AEP with respect to the relaxation strategy.