Advances and Challenges in Wave Energy Park Optimization—A Review

A commercial wave energy system will typically consist of many interacting wave energy converters installed in a park. The performance of the park depends on many parameters such as array layout and number of devices, and may be evaluated based on different measures such as energy absorption, electricity quality, or cost of the produced electricity. As wave energy is currently at the stage where several large-scale installations are being planned, optimizing the park performance is an active research area, with many important contributions in the past few years. Here, this research is reviewed, with a focus on identifying the current state of the art, analyzing how realistic, reliable, and relevant the methods and the results are, and outlining directions for future research.


INTRODUCTION
Wave energy has the potential to contribute significantly to the world's electricity consumption. To produce electricity in the range above a few MW, most wave energy concepts require that wave energy converters (WECs) are deployed together in arrays, or parks.
Interaction between the devices will affect the full performance, reliability, cost, and life-time of the park. The interaction can be hydrodynamical (scattered and radiated waves), mechanical (shared mooring and foundations), electrical (sea cables, substations, grid connection), and economic (shared capital and operational costs). As a result, many parameters will affect the interaction and the park performance, reliability, and costs. The number of devices and their separation distance, the park layout, mooring configurations, electrical and power take-off (PTO) systems, rated power of individual devices, constraints of subsystems, and so on-all parameters should be tuned to obtain the optimal design of the wave energy park before installation. In addition, environmental parameters have a large impact on the park. Wave climate, wave direction and variability, water depth, currents, and distance from shore are some factors that all affect the wave energy system, and different sites and environmental conditions will require different optimal solutions. Optimization is the procedure of identifying the best solution from some set of available alternatives, under given constraints. In the simplest case, an optimization problem consists of maximizing or minimizing a real function by systematically choosing input values from an allowed set and computing the value of the function. In a more complex situation, optimization includes finding the best available values of some objective functions given a defined parameter space.
Many papers on wave energy parks have claimed to carry out optimization, whereas, in reality, they have only compared a few distinct configurations or tuned one parameter to obtain a minimum or maximum point. In recent years, there has been a vast increase in the research field of wave energy park optimization, and several global optimization algorithms have been developed and applied. Ideally, optimization of a wave energy park should find the best available solution to an objective function that considers all aspects of the park, including all costs and total revenue over the lifetime, reliability, constraints regarding available ocean area, deployment, and maintenance, allowed power fluctuations, water depth, etc. In short, the goal of the research should be to provide wave energy developers with clear, unbiased and reliable answers on how they should best design their park, given the constraints they are facing.
The aim of the current paper is to ask how far we are from providing such answers. The paper poses questions on how realistic the current models and results for wave energy park optimization are, how reliable they are, and if they are relevant. Further, the paper aims to outline possible routes on how to reach that goal. The focus of the paper is on the state of the art in this research field. Most of the referenced papers have been published during the last five years, although older key references have also been included for consistency. All figures have been reproduced with permission from the journal copyright holders.

Structure of the Paper
The objective of this review paper is to review the state of the art of wave energy park optimization and to analyze how realistic, reliable, and relevant the current methods and results of the research area are. To arrive at this, first, the modeling methods and the optimization algorithms are reviewed. Modeling methods used in wave energy park optimization are reviewed in section 2, with focus on hydrodynamic modeling and WEC dynamics. Optimization methods are reviewed in section 3, and both simple "optimization" procedures such as comparing distinctive configurations or tuning single parameters are included, as well as global optimization algorithms such as evolutionary strategies. In section 4, we then return to the technical aspects reviewed in sections 2-3 and analyze how realistic (section 4.1), reliable (section 4.2), and relevant (section 4.3) the different methods and output are. Different approaches and systems are also compared to find general trends. Challenges and constraints for the research field are discussed in section 5, and future routes are outlined. Finally, some general conclusions are presented in section 6.

MODELING METHODS
As mentioned in the introduction, a large wave energy system consists of hydrodynamical, mechanical, and electrical subsystems, all with their own complexities and requiring different modeling techniques. In the current section, different methods used to model hydrodynamics, wave-structure interaction, and WEC dynamics, including power take-off models, are reviewed.

Hydrodynamic Modeling
The vast majority of works on hydrodynamical interactions in wave energy parks have been carried out based on linear potential flow theory. In other words, the fluid is assumed to be non-viscous, non-rotational, and incompressible, such that the governing equation of the fluid velocity reduces to the Laplace equation = 0, where (x, t) is the fluid velocity potential. By further assuming non-steep waves, the boundary constraints at the free surface can be linearized and the first-order linear approximation taken. Although there is some progress in studying the hydrodynamical interactions in wave energy arrays using non-linear potential theory and even computational fluid dynamics methods (Devolder et al., 2017(Devolder et al., , 2018, the optimization of arrays typically involves a large number of evaluations. This demands fast methods for the hydrodynamic computations, requiring that linear potential flow theory, and sometimes additional approximations, are used. For further reviews of the hydrodynamics modeling methods of wave energy parks than the one presented here; see Li and Yu (2012), Babarit (2013) and Folley et al. (2013).

Analytical Methods
Hydrodynamic modeling of wave energy parks by means of analytical methods was initiated already in the early works on wave energy. Budal (1977) studied a system of buoys under the assumption that the floats were so small that the interaction due to scattered waves could be neglected. The same point-absorber approximation was used by Evans (1980) and Falnes (1980) in early works on wave energy arrays and is sometimes used even today to facilitate modeling of large parks or modeling with low computational costs (Göteman et al., 2015b;McGuinness and Thomas, 2019). Similarly, in the plane-wave approximation, the structures are assumed to be so far apart that they influence each other only by plane waves, and the evanescent modes are neglected (Simon, 1982;McIver and Evans, 1984). The pointabsorber approximation has been found to perform well at low frequencies and for large separating distances, whereas the planewave approximation works better at high frequencies. However, the error is expected to grow with an increase in the number of interacting units (McIver, 1994;Mavrakos and McIver, 1997).
By applying an acoustic multi-body diffraction theory to water waves, the iterative multiple scattering method was developed by Ohkusu (1974) to study offshore platforms and further developed by Mavrakos and Koumoutsakos (1987) and Mavrakos (1991) for wave energy applications. In the iterative multiple scattering method, the diffraction and radiation properties of isolated bodies are used, and the reflected waves within an array are added iteratively until convergence is achieved. This method was combined with the direct matrix method of Simon (1982) and Kagemoto and Yue (1986). The resulting method, sometimes referred to as the direct matrix method and sometimes as the multiple scattering method, solves the wave amplitude around each body simultaneously without iteration and is exact within the assumptions of linear potential flow theory. For floating bodies, the infinite matrices in the solution must be truncated, and the method is semi-exact. The method was later combined with the single-body diffraction solution of Garrett (1971) and Yılmaz and Incecik (1998) and extended to independent radiation by Siddorn and Eatock Taylor (2008). Recently, the analytical multiple-scattering method was coupled by McNatt et al. (2015) to a numerical method to allow for arbitrary geometry (see Figure 1) and used to study the park interaction and existence of trapped wave modes by Flavià et al. (2018). Semi-analytical methods were also developed for oscillating wave surge converters (OWSCs) and used to study wave energy parks by Dias (2012, 2013), and Noad and Porter (2015). To improve the computational speed and allow for the optimization of large arrays, methods using resonant FIGURE 1 | (Left) The analytical multiple-scattering method was extended by McNatt et al. (2015) to allow for arbitrary geometry and used to study array interactions of attenuator WECs. The figure shows wave fields produced with the BEM software WAMIT and the analytical method (denoted IT in the figure). (Right) Wave propagation in a nine-WEC array in regular (top) and irregular (bottom) waves computed using NEMOH (Verao Fernández et al., 2018). modes (De Chowdhury and Manasseh, 2017;Wolgamot et al., 2017), interaction distance cut-off (Göteman et al., 2015a), a nearest neighbor approach (Sarkar et al., 2016), and Haskind's relation (Flavià and Clément, 2017) have been introduced. Both the iterative and non-iterative versions of the multiple scattering theory have been further developed and used to model the hydrodynamics of wave energy parks (Ji et al., 2015;Konispoliatis and Mavrakos, 2016;Göteman, 2017;Ruiz et al., 2017;Fang et al., 2018;Giassi and Göteman, 2018;Zheng et al., 2018Zheng et al., , 2019Liu et al., 2019). Several other analytical methods have been developed and used for wave energy park applications, including matched asymptotic expansions (McIver and Evans, 1988), multipole expansions (Linton and Evans, 1993), and Bragg scattering (Li and Mei, 2007).

Numerical Methods
Hydrodynamic modeling in offshore engineering is most often performed using numerical methods. The most commonly used approach is the boundary element method (BEM). At the EWTEC conference in 2015, over 37% of the papers in the wave energy tracks explicitly stated that a BEM software had been used (Penalba et al., 2017). In BEM, the boundaries of the fluid domain are discretized and the integral representation of the fluid velocity potential is used. The boundary condition on the body and free surface are applied and the fluid potential can then be determined using Green's functions anywhere in the fluid domain. Several commercial (Lee, 1995;AQWA, 2013) and open-source BEM software packages (Babarit and Delhommeau, 2015) exist. BEM software has also been coupled to wave propagation models such as MILDwave or OceanWave3D to study the perturbed wave field in a larger ocean area with varying bathymetry (Verbrugghe et al., 2017;Verao Fernández et al., 2018, 2019 and are further developed and used extensively to model the hydrodynamical interaction in wave energy parks (Sinha et al., 2016;Ruiz et al., 2017;Tay and Venugopal, 2017;Bosma et al., 2019;Lyu et al., 2019).
In the above-mentioned methods, the underlying theory is based on potential flow, to linear or second non-linear order. When steep waves are considered, the assumptions of linearization are no longer valid, and in addition, viscosity may have a large impact on the dynamics of floating bodies. Computational fluid dynamics (CFD) methods solve the full Navier-Stokes equations using numerical methods. Different commercial and open-source CFD software packages are used increasingly to study the hydrodynamic properties of WECs, but the computational cost is still too high to allow for wave park optimization studies, even if Devolder et al. (2017Devolder et al. ( , 2018 have recently studied interaction between a few WECs with CFD.

Motion and PTO Modeling
The dynamics of WECs differ fundamentally from those of traditional offshore structures, mainly due to the small scale of the WEC structures, the existence of a PTO, and the aim to absorb energy from the waves, which often requires large-amplitude motion at resonance with the waves. Many different PTO systems and methods of modeling them exist in the literature (see Folley et al., 2013;Day et al., 2015;Forehand et al., 2015;Penalba and Ringwood, 2016). For instance, the PTO can be a hydraulic system driving a rotating generator, a hydro-or air-turbine, or a linear direct-driven generator. The PTO extracts energy from the system, which can then be converted into electricity and connected to the grid. This extraction of energy affects the motion of the device, which again affects the hydrodynamical damping and the dynamics. Thus, the modeling of WECs requires that both the hydrodynamics and the PTO and the coupling between them are modeled in a realistic way. The PTO is often modeled as a spring-damper system, but it is only the damping part of the force that is responsible for power take-off, F PTO (t) = γ (t)ẋ(t), where γ (t) can be the pressure differential of a hydraulic system or the generator damping of a direct-driven generator, and x(t) is the dynamical part of the PTO. In simple cases, the damping γ can be approximated as constant. For oscillating water columns (OWCs), the PTO is typically an enclosed air chamber with air flow driving an air turbine, which can be modeled as a thermodynamic pressure model but has also been modeled successfully as a point-absorber by Sharp et al. (2017).
The equations of motion for a system of floating bodies can be described by Cummins' equation (Cummins, 1962), which includes a convolution term between the radiation impulse response function and the velocity of the float. A standard approach when optimizing wave energy parks is to solve the equations of motion in the frequency domain, which dramatically reduces the computational cost but excludes the possibility of modeling non-linear and time-dependent PTO systems. With the growing field of control algorithms to steer the WEC dynamics, the complexity of PTO modeling is growing with advanced time-domain models and wave-to-wire models. Realtime simulations of WECs are computationally demanding due to the radiation convolution integral, and real-time array simulations have therefore been restricted to small arrays or have been realized by approximating the radiation convolution term with a state-space model . WEC-Sim is an open-source MATLAB/Simulink software package that solves the dynamics of WECs in the time domain. Hydrodynamical force parameters are required as input and can be computed with numerical BEM software such as the one described above. WEC-Sim has been used recently by Bosma et al. (2019) to study the performance of an array of five OWC devices and by Mankle et al. (2019) to study an array of five pitching WaveStar WECs in a staggered layout. The SESAM package developed by DNV-GL was used to model the hydrodynamic and structural response of arrays of 10 WECs by Yang et al. (2019). Real-time simulation was carried out by Balitsky et al. (2014) to model the optimal configuration for a six-body array using measured wave data from the west coast of Ireland. Control models have been introduced both in experiments and numerical modeling of arrays Li and Belmont, 2014;Mercadé Ruiz et al., 2017;Nader et al., 2017;Thomas et al., 2018a), but it is beyond the scope of the current paper to go into the details of advanced generator simulations or wave-to-wire models. Instead, we refer to reviews such as Penalba and Ringwood (2016) and Wang et al. (2018) and focus on the most commonly used PTO models used in wave energy park optimization.
The WEC dynamics can also be affected by the mooring system. For floating WECs, the main purpose of the mooring system is station-keeping. For several wave energy concepts, however, the mooring system is integrated into the PTO. This is the case, for example, when the WEC consists of a buoy connected to a direct-driven linear generator or when the buoy is tethered to pumps driving the PTO. The mooring system of WEC arrays is most commonly modeled as linear springs, as in the works by Vicente et al. (2009), Konispoliatis and Mavrakos (2014), and Yang et al. (2019). Some examples of non-linear mooring models for wave energy parks exist. A coupled timedomain method was used by Gao and Moan (2009) to solve the WEC motion and mooring line tension simultaneously for a system of nine WECs. Vicente et al. (2010) modeled the mooring cables as catenary lines in a quasi-static analysis, and a finite element cable model was used by Bailey et al. (2015) to study arrays of different WECs with a stochastic method, focusing on power variations. For a comprehensive review of modeling methods for WEC mooring systems, see Davidson and Ringwood (2017). In section 4.1, we will return to the numerical models of WEC arrays and discuss how well they represent the physical systems studied.
By solving the equations of motion and obtaining the dynamics of the system, the absorbed power can be computed. Once the absorbed power has been computed, the hydrodynamical interaction can be evaluated using the interaction factor, or q-factor. This is usually defined as the ratio between the time-averaged absorbed power of the full farm, P tot , and the power that would have been absorbed by the same number N of isolated WECs, P i isolated , An interaction factor q > 1 implies that the hydrodynamical interaction in the park is constructive, i.e., the gain is larger than the destructive interactions. The interaction factor is often used to evaluate the total performance of the park, although this is a rather restrictive measure, as will be discussed in section 4.3.

Physical Experiments
Both analytical and numerical modeling are always connected with some approximations and uncertainty, and for an accurate understanding of the systems and for reliable results, physical experiments are required. However, experiments with wave energy arrays are both expensive and complex to carry out, and it is only in recent years that several large-scale physical experiments have been conducted and their results published.
Numerical predictions of the response of an array of five heaving floats in regular waves were compared to experimental measurements by Thomas et al. (2008) and in irregular waves by Weller et al. (2010). Magagna et al. (2011) measured the interaction factor in an array of OWCs experimentally. Trapped modes within an array of eight fixed truncated cylinders were demonstrated experimentally by Wolgamot et al. (2016), and near-trapping was also observed in an array of fixed cylinders in short-crested waves by Ji et al. (2015). A triangular array of three spar-buoy WECs at a scale of 1:32 was studied experimentally by Correia da Fonseca et al. (2016), and it was found that, in wave climates with large energy periods, the array performed better than three isolated devices, i.e., an interaction factor larger than one was observed. Two layouts of five fixed OWCs were modeled experimentally and numerically in regular and irregular waves by Bosma et al. (2019). The layouts were chosen based on earlier layout optimization, and a maximal power increase of 12% from the non-optimal to optimal layout was found. However, when taking the average, only a modest increase could be established.
Recently, experiments with more units and increased complexity in terms of multiple degrees of freedom dynamics or with advanced control algorithms have been conducted. An experimental campaign of large arrays of up to 25 heaving pointabsorber WECs was carried out in the shallow water wave basin of DHI in Denmark in the WECwakes project and has been presented by Troch et al. (2013) and Stratigaki et al. (2014); see Figure 2. The PTO was modeled as friction damping, and both long-crested and short-crested waves were studied. It was found that the wave height was reduced after the park by up to 18% due to park interactions (Stratigaki et al., 2014) and that it could be increased by up to 11% up-wave of the array due to radiated waves (Troch et al., 2013). Similarly, experiments with arrays of up to 24 point-absorber WECs were conducted at Queen's University Belfast, and hydrodynamical array interactions were found experimentally to cause losses of up to 26% in energy yield compared to WECs in isolation (Child and Laporte Weywada, 2013). Arrays of five WaveStar 1:20 scale models equipped with linear control PTO systems were studied experimentally by Mercadé Ruiz et al. (2017) (Figure 2), and the data were used to validate the hydrodynamics tool implemented in the DTOcean software.
Arrays of up to six point-absorber WECs, each moving in six degrees of freedom were carried out both at the Australian Maritime College by Nader et al. (2017) and at the COAST laboratory at the University of Plymouth, UK, by Thomas et al. (2018a) and Giassi et al. (2019b). In the former experiment, the surface elevation was tracked by videogrammic measurements, and the experimentally measured interaction factor was presented for 1-2 floats moving in heave and surge (Nader et al., 2017). In the latter, both a linear damping and an advanced control algorithm based on machine learning and artificial neural networks were used as PTO systems (Thomas et al., 2018b). The collaborative control algorithm required no previous knowledge of the incident waves but still increased the energy absorption through communication between the WECs (Thomas et al., 2018a), and the performance of the WECs with linear damping was compared for three array layouts by Giassi et al. (2019b).
The experimental works discussed above were all carried out in wave tanks. Experimental results from real sea tests of arrays of WECs are even more rare. Three fullscale WECs were deployed by Uppsala University at the Lysekil offshore site in 2009, and array interactions were studied (Rahm et al., 2012). Several commercial developers have claimed to have deployed and grid-connected arrays of WECs, but no data have been made publicly available.

Comparing Distinctive Configurations
Most studies on designing optimal wave energy parks have been based on comparisons of a few different configurations (Babarit, 2013), often with little or no arguments as to why these particular configurations were chosen for the study. Whereas these works were not optimizations per se, they lay the grounds for the current optimization strategies.
Several layouts and parameter settings of arrays of heaving cylinders and surging barges were studied by Borgarino et al. (2012), and it was concluded that grouping of WECs into Frontiers in Energy Research | www.frontiersin.org arrays had a constructive effect in general. Vicente et al. (2013) compared different park layouts consisting of 12 WECs each with three different numerical methods, and separation distances between the WECs were identified that led to constructive interactions between the units. Similarly, Engström et al. (2013) compared circular, rectangular, and randomized array layouts of 32 point-absorbers for power production and power fluctuations. Three different array layouts and three different separation distances were analyzed for life-cycle performance by López-Ruiz et al. (2018a) (see Figure 3). Five different array layouts consisting of 12 WECs were compared by Sinha et al. (2016) (see Figure 4), and their performance in sea states with different incident wave directions was studied. Several array layouts were compared for two-body heaving WECs by De Andrés et al. (2014). Triangular arrays were found to be optimal for multidirectional wave regimes, whereas square arrays were better for unidirectional waves, and tuning the separation distance was found to be important to reach higher interaction factors. Different array layouts such as linear, rectangular, and staggered arrays were studied for Oyster and Pelamis WECs in 10year wave spectra by Gorr-Pozzi et al. (2019) and for pointabsorber WECs at four Italian sites by Bozzi et al. (2017). Power production, costs, and mooring fatigue were compared for four array layouts of 10 WECs by Yang et al. (2019) (see Figure 4).

Parameter Sweep
The next step after comparing distinctive configurations is to vary one parameter to find the optimal value, keeping all other parameters constant.  number of devices in the array. Layout optimization by sweeping over parameters such as separation distance and angles between devices was performed by Snyder and Moarefdoost (2015) and Göteman et al. (2015a). The power production of arrays of 12 OWSCs was studied in both regular, long-crested, and shortcrested waves by Tay and Venugopal (2017), and the interaction function was evaluated as a function of the separation distance between devices. Also, parameter sweeps to determine the impact of wave parameters such as wave number, wave direction, and significant wave height has also been conducted in a number of works, including Konispoliatis and Mavrakos (2016)

Global Optimization Algorithms
As has been discussed above, many parameters may affect the performance of a wave energy park: the layout of the park, the separation distance between devices, the number and size of the devices in the park, the individual PTO settings of each unit, mooring configurations, control strategies, wave climate and wave direction, etc. Simple sweeps over single parameters are not sufficient to find optima in the broad and complex solution space of wave energy parks, and more advanced optimization algorithms are needed.

Non-linear Programming Optimization
Non-linear programming optimization is useful for handling optimization problems that cannot be handled by simple parameter sweeps but which require optimization across a moderately large parameter space. These methods can also be adapted to handle non-linear constraints on the WEC motion, which was applied already in early works by Evans (1981) and Pizer (1993) on wave energy optimization. The efficiency and accuracy of the solution to the constrained optimization problem depend on the number of constraints and design variables and on the characteristics of the objective function and constraints. When constraints are linear and the objective function is quadratic or linear, reliable solutions are readily available. If instead the objective function and the constraints are both non-linear, the problem can be transformed into a subproblem that can be solved and used as the core of an iterative process to find the full solution to the optimization problem.
Optimization of the array layout for point-absorber, spherical WECs operating in heave only was performed by McGuinness and Thomas (2015Thomas ( , 2016. To simplify the optimization, the WECs were constrained to be positioned on a straight line or along a circle, and no other constraints on the motion were imposed. More generally, the optimization problem for five WECs without restrictions on the array layout was studied by McGuinness and Thomas (2019). In these papers, the layout optimization algorithm identified the minimal value of an objective function using a sequential quadratic programming method, which resulted in the optimal layouts shown in Figure 5.
Optimization in terms of maximizing the energy absorption under WEC motion constraints was considered by  and  for three arrays of 2-3 WECs. Similarly, an iterative quasi-Newton optimization algorithm was applied by Wu et al. (2017) to identify the maximum captured power of an array of Duck WECs under motion constraints. Optimization of the absorbed power for point-absorbers operating in heave in a fixed array layout but with individually chosen constant damping of the WECs was considered by Wang et al. (2016a,b). The impact of motion constraints of the individual WECs was also taken into account for different scenarios. An improvement in power absorption of up to 18% was found with individual control of the damping compared to a situation with optimal but same damping for all WECs. Here, the non-linear problem was transformed into quadratic programming using a penalty function for constraints near the constraint boundary, implemented in MATLAB in the "fmincon" function. Optimization of arrays with individually chosen time-varying damping was attempted by Wang et al. (2016c) but for arrays consisting of two units only due to the heavy computational requirements. Moarefdoost et al. (2017) used an iterative quadratic programming to optimize farm layout for up to 15 point-absorbers; see Figure 5. However, their method also used a metaheuristic algorithm (see section 3.3.2) to choose the starting points of the optimization process.
More examples of array optimizations coupled to WEC control algorithms exist (Garcia-Rosa et al., 2015;Sharp et al., 2017). Advanced control methods fall outside the scope of the current paper, and the reader is referred to review papers such as Penalba and Ringwood (2016) and Faedo et al. (2017) for further details.

Metaheuristic Algorithms
Metaheuristic optimization methods search the solution space for sufficiently good solutions, given specified constraints and convergence criteria. They are useful when the optimization problem is too large for all solutions to be evaluated and when the solution space is multi-peaked, and several different algorithms have been developed and applied to wave energy park optimization.
Genetic algorithms. Genetic algorithms (GAs) originate from the theory of evolution studies but have since been applied in a variety of areas, ranging from string theory to pharmaceutical research (Holland, 1992). The algorithm mimics the process of natural selection and incorporates inheritance, mutation, and selection among states in the solution space to find the global minima of stated problems. Several versions of GAs exist, but the idea is the same: a certain percentage of the best performing solutions is kept in each generation based on a specified fitness value. In the next generation, the surviving solutions (parents) are combined two and two to generate new solutions (children), using different crossover mechanisms. As the population evolves, individuals with better and better fitness values are generated. In addition, mutation introduces random FIGURE 5 | (Left) Optimal layout for arrays of 3-7 WECs, obtained through sequential quadratic programming where the starting point for the layout has been chosen based on a heuristic optimization algorithm (Moarefdoost et al., 2017). (Right) Optimal layouts for arrays with five WECs with constrained motion amplitude δ (McGuinness and Thomas, 2019). In the upper figure, the incident waves propagate in directions close to the y-axis, and in the lower figure, the incident waves propagate in a wider range of directions. All layouts have a WEC at the origin. solutions in the population and helps avoid local minima in the solution space. A flowchart of the algorithm is shown in Figure 6. Covariance matrix adaptation. The covariance matrix gives a geometric interpretation of the solution distribution in a generation. Covariance matrix adaptation (CMA) is a stochastic metaheuristic evolution strategy that adaptively increases or decreases the search space for the next generation based on the shape of the covariance matrix (Hansen, 2006). A large search space is used when the solutions are far from the global minima but is automatically reduced when we are close to a minimum and the solution space just needs fine-tuning. Differential evolution. Differential evolution algorithms also belong to the family of evolutionary optimization algorithms. Like the GAs, they allow solutions to evolve by inheriting strengths from previous generations. The main difference from GAs is that real-valued vectors are used to represent the individuals in the population. Mutation/crossover is performed by adding the weighted difference between two solutions to a third solution. This difference vector can be used with a scaling factor to traverse the search space (Storn and Price, 1997). Particle swarm algorithms. A particle swarm algorithm mimics swarm behavior in birds flocking and fish schooling to guide the particles to search for global optimal solutions (Eberhart and Kennedy, 1995). Based on the fitness values of all particles in the swarm, the velocities and positions of the particles are updated, and the swarm propagates toward optimal solutions. In the related glowworm optimization algorithm (Krishnanand and Ghose, 2009), the particles are equipped with a luciferin level, depending on their fitness values. Each glowworm is attracted by the brighter glow of other neighboring glowworms, bringing the swarm toward the optimal solutions. Other. Several other global and metaheuristic algorithms have been used for wave energy park optimization, for instance, the gray wolf algorithm, another optimization strategy inspired by animals living in groups (Mirjalili et al., 2014). The parabolic intersection method was designed to optimize layouts by placing WECs iteratively on the parabolic scattered wave fields produced by existing WECs in the array (Child and Venugopal, 2010).
GAs were first applied for the layout of wave energy arrays in Child and Venugopal (2010) and further developed by Child et al. (2011) to obtain optimal layouts and PTO settings for point-absorber arrays. The GA toolbox from MATLAB was used and was compared to a parabolic intersection method.  coupled the analytical multiple scattering method of Göteman et al. (2015a) to a GA optimization tool, and  used it to find optimal layouts of arrays with WECs of different sizes. This work was further continued by Giassi and Göteman (2018) to obtain optimal arrays with 4 to 14 WECs where the interaction factor was larger than one. The layout of five point-absorber WECs was optimized with a binary GA by Sharp and DuPont (2018), and constraints in terms of minimal distance between the WECs implied different optimal layouts, although the general trend was to align perpendicular to the wave direction. Sharp et al. (2017) coupled the array layout optimization to active device control. An adaptive mutation factor was introduced by Fang et al. (2018) to improve the convergence rate of the array optimization. Faraggiana et al. (2019) used a GA to optimize the separation distance between devices, the PTO tuning, and the rated power of a system of three WaveSub WECs.
To evaluate the reliability and efficiency of different metaheuristic optimization strategies for wave energy applications, several recent works have carried out array optimizations using different algorithms and compared their results. A particle swarm optimization was compared to a GA by Faraggiana et al. (2019) and was found to outperform the latter for two WECs, but not for larger arrays. In the research project DTOcean (optimal design tools for ocean energy arrays), a set of optimization tools for ocean energy systems have been developed (DTOcean, 2017). Ruiz et al. (2017) compared the CMA evolution optimization algorithm of the DTOcean tool to a GA and a particle swarm optimization algorithm. The yearly absorbed power was maximized as a function of the number of WECs and the park layout, under the constraints of minimal interaction factor q ≥ 0.9, minimal separation distance 65 m, and maximal ocean area 500 × 500 m. It was found that the GA and the particle swarm algorithms performed slightly better than the CMA algorithm but that the latter had lower computational cost. Similarly, the accuracy and computational cost of the GA in the DTOcean modeling tool were compared to a meta-model optimization method by Ferri (2017), with the conclusion that the former was more accurate, but with a higher computational cost. A GA was compared to a Monte-Carlo method to obtain optimal arrangements of 40 OWSCs by Sarkar et al. (2016) (see Figure 7). The hydrodynamical interaction was restricted between nearest neighbors of WECs to allow for fast evaluations, and a machine learning approach was used to identify optimal three-WEC clusters, after which the GA was applied to optimize the full park layout. With the same number of 16 million evaluations, the GA was able to find more optimal layouts than the Monte-Carlo simulations. Neshat et al. (2018Neshat et al. ( , 2019a compared many different metaheuristic optimization algorithms for wave energy park applications; see Figure 6. Neshat et al. (2018Neshat et al. ( , 2019a optimized the coordinates of the WECs to achieve maximal absorbed power under the constraints of maximal ocean area and a minimal separation distance between the units. The parks consisted of 4 and 16 submerged CETO WECs, and the wave climates at several Australian sites were considered. Neshat et al. (2019b,d) extended the optimization problem to include the individual PTO parameters as well. Several optimization strategies were compared: either all parameters were optimized simultaneously, using several algorithms such as the CMA strategy, differential evolution, gray wolf optimization, or particle swarm optimization, or the WEC coordinates and PTO parameters were optimized in an alternated fashion, or the WECs were positioned one-by-one and a local search was applied to find the best coordinates for the next WEC. In terms of convergence speed, the last approach outperformed the former, but possible limitations were that no backtracking was allowed to further optimize WEC positions once they had been placed, which was, however, not believed to provide any large gain in the power absorption (Neshat et al., 2019c). Metaheuristic evolutionary algorithms were also applied by Arbonès et al. (2016Arbonès et al. ( , 2018 for multi-objective optimization and were seen to outperform single-objective algorithms; see further details in section 4.3.

CURRENT STATE OF THE ART
The research area of wave energy park optimization is still young, with many important developments in the past few years. But as the wave energy industry is currently moving toward ocean deployments and fullscale systems, it is important to identify the state of the art and where to focus future work. In this section, we return to the technical aspects of wave energy park modeling and optimization reviewed in sections 2, 3 and analyze how realistic, reliable, and relevant the methods and results are. These aspects are discussed in sections 4.1-4.3, respectively.

Realistic Representation
In this section, we analyze how well the WECs and environmental systems are modeled, i.e., whether the models are realistic representations of the studied physical systems.

Realistic Waves
As discussed before, most results on wave energy parks have been obtained under the assumption of potential flow theory, where the viscosity and the rotational behavior of the fluid have been neglected. Whereas some advances have recently been presented to model the hydrodynamical interaction between WECs with CFD methods (see section 2.1), it is still unthinkable to carry out optimization of large wave energy systems using high-fidelity viscous and rotational fluid models. Instead, physical experiments will have to be used to validate the analytical/numerical methods used for wavestructure modeling, which will be further discussed in section 4.2. But even within the assumptions of potential flow theory, further simplifications on the waves have often been assumed, and most arrays have been studied in a few sea states with long-crested waves. To study the performance in more realistic wave conditions, wave directionality and long-term wave spectra have been incorporated in recent works on park optimization.

Wave directionality and short-crested waves
The impact of the wave directionality on wave energy park performance has been studied in a number of works. As could be expected, the park performance of array layouts with rotational symmetry is less affected by wave directions than that of corresponding parks with rectangular or linear layouts (De Andrés et al., 2014;Göteman et al., 2015a). Yang et al. (2019) considered two sea states and included wind and water current velocity aligned in the same directions as the waves. McGuinness and Thomas (2019) performed layout optimization over an interval of incident wave directions for unidirectional waves. Not surprisingly, it was found that better performance could be obtained for narrow-banded wave directionality, as the layout could then be tuned for those particular wave directions. The directional wave spectra at two Australian sites were represented by seven wave directions by Neshat et al. (2019a,b,c). The identified optimal layouts were different for different sites: at one of the sites, where the directional wave spread was small, the WECs were placed roughly perpendicular to the predominant wave direction (see Figures 11, 12), whereas no obvious geometric arrangement was found at another site with a larger interval of incident wave directions (Neshat et al., 2019b) (see Figure 7, right).
In reality, ocean waves are usually not long-crested but short-crested, i.e., they consist of many waves traveling in different directions simultaneously. Recently, some works have considered the performance and optimization of wave energy arrays in short-crested, or omni-directional, waves. Wave runup and trapped wave modes were studied both experimentally and analytically in an array of bottom-mounted cylinders by Ji et al. (2015). The computed and experimentally obtained results agreed well, and it was seen that unidirectional wave modeling would overpredict the normal force and underpredict the transversal force on the cylinders. An array of 25 WEC models was studied experimentally in short-crested waves by Stratigaki et al. (2014), and it was found that the wave height decrease occurred earlier than in corresponding long-crested waves. Tay and Venugopal (2017) studied the interaction factor of a park of 12 OWSCs in short-crested wave spectra, and it was found to be slightly lower than in long-crested waves. The performance of an array of 16 WECs in short-crested waves was studied by Göteman et al. (2018), and it was concluded that the energy absorption was comparable to the situation in long-crested waves, whereas the power fluctuations were considerably lower.

Long-term and large domain wave propagation
Most array optimization studies have focused on optimization in one or a few wave conditions. The relevant objective should rather be the full life-cycle performance of the arrays, which requires a long-term assessment. However, the wave propagation models that are typically used to model waves over long time period or large domains cannot be used directly to model the dynamics and instant power of the WECs, and simplified estimations of the absorbed energy have previously been used to assess the potential of arrays over large domains (Defne et al., 2009;Arinaga and Cheung, 2012;Akpınar and Kömürcü, 2013;Vicinanza et al., 2013;Soomere and Eelsalu, 2014).
The wave climate over 25 years, comparable to the life-time of the devices, was used by López- Ruiz et al. (2018a,b) to compare layouts and separation distances for arrays with nine floating overtopping WECs; see model as the basis for waves to calculate the instant and average absorbed energy from a park of 200 generic point-absorbers in a 1 km 2 grid over the Swedish exclusive zone. A ten-year hindcast model was used by Gorr-Pozzi et al. (2019) to study three wave parks with 5-25 WECs and their impact on the near-shore environment in long-term wave spectra.
Efforts to couple accurate hydrodynamic modeling close to the WEC with wave propagation models a distance away from the WEC have been made both for potential flow (McNatt et al., 2015;Verbrugghe et al., 2017;Verao Fernández et al., 2018, 2019 and CFD models (Paulsen et al., 2014;Verbrugghe et al., 2018). Comparison of the predicted wave height from a coupled MILDwave-NEMOH model with experimental data was presented by Verao Fernández et al. (2019) (see Figure 8). Decomposition of the fluid domain into connected volumes that require different modeling techniques can be a useful approach both when non-linear effects should be studied in an array, when studying the far-field effects of parks, and when large-scale parks over several kilometers are designed.

Realistic WEC Systems
A wave energy system consists of many parts, often including multi-body floats, moorings, power take-off systems, foundations, and electrical systems in generators, substations, sea cables, and grid connection points. Each of the subsystems is complex, and it is impossible to include all of the complexity of each subsystem in the optimization of a full wave energy park. But as the obtained results can be misleading if important aspects of the full system are neglected, proper attention must be given to the approximations made and how well they reflect the realistic wave energy system.

Power take-off and electrical subsystems
As discussed in section 2.2, the standard approach when optimizing wave energy parks is to model the PTO as a linear spring-damper system, where the power absorption is proportional to the Coulomb damping. In addition, the dynamical degrees of freedom are usually restricted; for example, point-absorbers are considered to move in heave only. Comparison of numerical PTO models to experimental data of scaled models (with simplified PTO systems) has been carried out in a number of works. Linear and non-linear PTO models for the Wavepiston WEC were compared in frequency and time domain by Read and Bingham (2018), both with and without viscosity. The overall agreement was good, although the frequency domain model predicted slightly higher energy absorption when the PTO was non-linear. Similarly, experimental data for a two-body floating WEC with eddy current brake PTO was compared to a frequency and time domain model in Kurniawan et al. (2019). Although the general agreement was good, the time domain model was able to predict non-linear instability behavior that the frequency model was not able to capture. In the WECwakes project (Troch et al., 2013;Stratigaki et al., 2014), the predicted total power from a linear time-domain model overestimated the measured total power for single WECs by around 20% (Troch et al., 2013). Mankle et al. (2019) compared WEC-Sim with experimental data from five pitching WaveStar WECs in a staggered layout.
Comparisons of numerical PTO models to fullscale, realistic WECs in offshore operations are rarer. Almost all tests with fullscale WECs have been conducted by companies and are kept within the company, but a few exceptions exist. A linear timedomain model for a single point-absorber WEC was validated against experimental data from a fullscale WEC by Eriksson et al. (2007). The PTO in the numerical model consisted of a realistic circuit with generator, sea cable, and resistive load. The comparison showed good agreement for time-averaged energy absorption. However, the PTO model had a time-dependent damping coefficient, in contrast to the constant damping usually used in wave energy park modeling. A different approach was used by Ulvgård et al. (2016), where onshore tests were conducted on a fullscale WEC of point-absorber type by simply pulling the translator of the direct-driven generator with a crane. It was seen that the generator damping coefficient was indeed approximately FIGURE 8 | Comparison of the wave height obtained experimentally and numerically with a coupled MILDwave-NEMOH model (Verao Fernández et al., 2019). The K d coefficient is defined as the ratio between the numerically calculated and the target significant wave height, and the wave height has been recorded by wave gauges (WG) in an array of nine heaving damped WECs. constant when the translator was fully within the stator and reduced with decreasing active area for the translator.

Mooring
Mooring lines of floating wave energy systems may exhibit nonlinear dynamics and snap loads, which will affect both the absorbed power and the life-time of the energy system. Whereas numerical models exist to study this dynamical behavior for isolated WECs (Bhinder et al., 2015;Palm et al., 2016), the common assumption in wave park optimization is to model the mooring as linear springs, as discussed in section 2.2. In many cases of wave energy park optimization, the mooring dynamics are simply neglected. The performance of a system of two WECs was found to be significantly affected by the presence of the mooring lines by Vicente et al. (2009). A system of three floating OWCs was studied by Konispoliatis and Mavrakos (2014), and it was concluded that the power absorption was higher in the case when mooring was considered than for freely floating OWCs. The time-domain method used by Gao and Moan (2009) to model the mooring of nine WECs showed that in certain situations, large mooring line tension was observed in the array. Vicente et al. (2010) found that the performance of the array was significantly affected by the mooring system, as the occurrence of low-frequency horizontal oscillations of large amplitude was observed. The discussed results show that mooring dynamics may play an important role in the dynamics of wave energy parks and that careful consideration must be taken before assumptions on the mooring system are made.

Reliability of the Results
As for all sound science, the methods and the results must be reliable; the methods should produce stable and consistent results. Methods should be validated by, for example, experimental data, and methods applied on the same object a number of times should produce the same results, i.e., the results should be repeatable. Also, the methods and studied systems should be chosen without bias to make sure that the research questions and methods are thoughtfully chosen. In addition, the modeling errors and how they manifest in the results should be quantified with adequate methods.

Validation With Experiments
With the physical experiments of wave energy parks carried out in the past few years, validation of analytical or numerical modeling methods has become possible. An analytical multiple scattering method was compared to experimental data for an array of bottom-mounted cylinders in short-crested waves by Ji et al. (2015), and the agreement on wave run-up was generally good, but the computations underpredicted the wave forces. Simulations using BEM software and analytical methods were compared to experimental data for five WaveStar models at a scale of 1:20 by Mercadé Ruiz et al. (2017). The power take-off was modeled as constant damping, and the predicted results from the simulations were in good agreement with the experimental results (see Figure 9). The largest deviations occurred at resonance with 40% deviation, or 20% for panchromatic waves. The agreement between prediction and measurement was shown to be dependent on the floats' position in the array, which was identified to be caused by the difficulty of reproducing twodimensional plane parallel waves in the physical models due to undesired three-dimensional wave effects. Mankle et al. (2019) validated that the float response computed by WEC-Sim for five floats moving in pitch with experimental data had an acceptable error of 11% or less, although there was a larger variability in the experimental results, and only results from three regular waves were presented. When comparing the frequency domain solver and the spectral domain solver with experimental data of 24 point-absorber WECs, Child and Laporte Weywada (2013) found that the frequency domain solver was better at predicting the hydrodynamic interactions alone (ignoring absolute WEC performance) than the spectral solver, while the spectral solver was better at predicting long-term average power due to its ability to use non-linear WEC characteristics, and the spectraldomain method was expected to improve when more realistic array definitions would be considered.
Physical tank tests and numerical modeling using WAMIT for OWCs were compared by Sharp et al. (2017) to validate the numerical model used in wave park optimizations. In the physical models, a butterfly was used as PTO instead of the turbine used in fullscale devices, and in the numerical model, the mass of water moving inside the OWC was modeled as a point-absorber. The response amplitude operators obtained by the experimental and numerical methods agreed rather well; see Figure 10. The work was extended by Bosma et al. (2019), who compared a numerical model based on WAMIT and WEC-Sim to experimental results. The numerical predictions agreed well with the experimental, although the phase and amplitude were off in a few of the tested cases; see Figure 10.
Most of the experiments on arrays have been carried out using heaving buoys (Nader et al., 2017), despite the fact that the fullscale WECs may be designed for operation in more degrees of freedom. Several experimental works with waveactivated WECs free to move in multiple degrees of freedom have reported motion instabilities (Tarrant and Meskell, 2016;Gomes et al., 2017), which can often be attributed to the PTO settings and large amplitude motions and can be understood in terms of Mathieu equations (Orszaghova et al., 2019). Whereas experimental results in recent years have enabled validation of numerical models used to model and optimize wave energy parks, more work remains to be done in terms of quantifying the uncertainties in the numerical models, and with experiments using more realistic and advanced PTO and WEC dynamics. Sinha et al. (2016) compared several different park layouts, and it was found that the circular layout gave the most uniform power output and better power quality compared to linear and grid arrays. Similar studies with similar results had been presented before by Engström et al. (2013) and Göteman et al. (2014), which shows consistency in the results. A circular array is less dependent on wave direction, although it should be kept in mind that these simulations were carried out with long-crested waves, and the benefits of circular layouts as compared to linear and grid layouts should be less in omnidirectional waves.  Arrays with WECs of different dimensions were studied by Göteman (2017), , and Lyu et al. (2019), where  and Lyu et al. (2019) optimized both array layout and buoy dimensions using GAs. In all three works, it was concluded that the full park performance can be improved if the park contains WECs of different individual dimensions.

Repeatability
Many wave energy park optimization algorithms have started from prescribed layouts and optimized parameters such as separation distance between devices. A more unbiased approach was taken by McGuinness and Thomas (2019), where the layout of arrays with five point-absorber devices was optimized with the only constraint that the devices should be located within a half-sphere of a given area. Similar optimization routes with few a priori constraints on the park configuration have been enabled by the global optimization methods described in section 3.3, and some comparison of their results will be presented here.
Optimal configurations through single-or multi-objective optimization routines such as evolutionary algorithms were studied by Child and Venugopal (2010), Sharp and DuPont (2018), Giassi and Göteman (2018), Child et al. (2011, Fang et al. (2018), Arbonès et al. (2018), andNeshat et al. (2019b). All studies, although based on different devices, assumptions, parameters, models, and optimization algorithms, show similar layouts trends, with WECs aligned perpendicular to the incident wave field; see Figures 11-13. It also appears that, if physically possible, the WECs align in a single row perpendicular to the wave field direction, while they are staggered in two rows when the area or the separation distance constraints do not permit full alignment, as shown in Figure 12. The results are comparable, as all the authors model waves that are long-crested or, in the case of real wave climates, incoming mostly from one predominant direction. In some publications, the WECs were not aligned perpendicular to the wave direction but along a line with a small offset angle; see Figure 13. In wave climates where waves have a larger directional spread, the optimal layouts have been found to take less intuitive form (Neshat et al., 2019b); see Figure 7 (right).

Relevance of the Results
It is of little value that the methods and results are realistic and reliable if the results are not relevant for the wave energy industry. In other words-are we optimizing the correct quantities? If the optimization algorithm maximizes the power output of a park, but this produces power peaks that are too large for any realistic electrical system or grid connection, the results are of little use for the development of large-scale wave energy systems. Or, if a relevant objective function is being optimized but at the expense of immensely increased costs, this again will not be a  . (C) Multi-objective optimization using an evolutionary algorithm (Arbonès et al., 2018). (D) Optimal layout of four WECs obtained with a symmetric local search combined with a Nelder-Mead algorithm (Neshat et al., 2019b). (G) Optimal layout obtained with a GA in Sharp and DuPont (2018) for the case of a 4-m minimum separation distance. Figure 11 but for the case when the allowed area for the park does not allow for the WECs to all align perpendicular to the predominant wave direction. The waves are propagating along the x-direction in all figure parts except for (C), where the wave direction is a narrow interval around 240 • . (A,B) GA optimization in long-crested waves with a gridded and continuous layout, respectively . (C) Optimal layout of 16 WECs (Neshat et al., 2019b). (D) GA optimization for the case of a 3-m minimum separation distance (Sharp and DuPont, 2018). (E) Optimal eight-WEC array layout obtained in Fang et al. (2018) with an evolutionary algorithm and regular waves. (F) Multi-objective optimization of array with nine WECs (Arbonès et al., 2018).

FIGURE 13 | (A)
Optimal five-WEC array layout obtained by Fang et al. (2018) with an evolutionary algorithm and regular waves. (B) Optimal layout of five WECs in regular waves by a GA optimization (Child and Venugopal, 2010). (C) Optimal layout of 10 WECs in irregular waves propagating along the x-direction, obtained with the array design tool in WaveFarmer (Child et al., 2011). relevant option for the industry. In order to compete with other energy technologies, the objective function should maximize the total revenue over the lifetime of the devices within given constraints, such as maximal ocean area that can be used or maximal allowed power fluctuations. Many of the constraints can also be connected to economic values, for example, that increased separation distance between WECs requires more use of sea cable in the park, which increases the costs. Although most array optimization works have mainly focused on maximizing the absorbed power or the interaction factor, recent works have started to include more comprehensive objective functions.

Delivering Grid-Compatible Electricity
Rapidly varying voltage magnitudes, known as flicker, is a major problem when integrating the produced electricity into the grid. Different voltage levels define different flicker severity over the short term (10 min) and long term (over 2 h) (Penalba and Ringwood, 2016). Flicker evaluation, or electricity quality, is therefore a crucial aspect to consider when designing wave energy parks. Smoothing of the power from wave energy arrays was studied in Tissandier et al. (2008). As expected, it was concluded that the power peaks reduce with increased number of devices, which was later also found in other wave park simulations by Engström et al. (2013), Vicente et al. (2013), Göteman et al. (2014), and Bailey et al. (2015) and in offshore experiments by Rahm et al. (2012). By arranging the park layout such that the WECs are not excited simultaneously by the incident waves, the power output can be further smoothened (Engström et al., 2013;Sjolte et al., 2013;Göteman et al., 2014). Armstrong et al. (2015), Kovaltchouk et al. (2016), and Blavette et al. (2016) used approximate methods for hydrodynamics and power production in the array to study flicker severity with a focus on transmission and grid aspects. The impact of a park consisting of ten OWCs and ten heaving buoys on the transmission system was studied by Armstrong et al. (2015), and it was concluded that the transformers and transmission lines would remain below 45% capacity. Kovaltchouk et al. (2016) studied limits to voltage fluctuations for arrays of different WEC types, and compensatory actions were proposed to allow grid integration. Blavette et al. (2016) estimated the flicker level as a function of grid impedance angle. Hydrodynamics and electro-mechanic aspects were included in the works by Tedeschi and Santos-Mugica (2013) and Parwal et al. (2018). A 20-MW park at the Bimet facility was chosen as a test case by Tedeschi and Santos-Mugica (2013), and a control strategy and simplified shortterm energy storage system were proposed to reduce the power variability. Parwal et al. (2018) proposed an energy management system in terms of a combined battery and supercapacitor to minimize the power fluctuations from a wave energy park to the grid. Other dedicated power electronics and filtering equipment can also be used to reduce flicker and obtain acceptably smooth electricity output, but these methods fall outside the scope of this paper.

Economical Objective Functions
Simplified economic cost functions have been implemented in wave energy park optimizations in several recent works. Sharp and DuPont (2018) defined the objective function as the ratio of the "costs" and the produced power, where the costs were approximated by a simple formula obtained from Sandia National Lab's reference model project (Previsic, 2012), where N is the number of WECs in the park.  defined and evaluated different objective functions in a GA to optimize park layouts with WECs of different sizes. The functions were defined as ratios of the total absorbed power and the total mass of the units, as this is often used as a crude estimate of the cost of the technology, and MWh/tons is commonly used as a measure to compare different wave energy concepts . Sea cables and electrical subsystems contribute significantly to the costs of the park, both in terms of installed capital costs and costs due to energy losses (Henfridsson et al., 2007;Sharkey et al., 2011Sharkey et al., , 2013Engström et al., 2019). Cable dimensions and lengths, distance from grid connection point, number of substations and related infrastructure all affect the total revenue. Sharkey et al. (2011Sharkey et al. ( , 2013 studied how the costs related to the electrical network can be reduced by optimizing the spacing and capacity factors of the individual WECs. Optimization of parks with respect to cable length and costs was done by Arbonès et al. (2016) and Giassi et al. (2019a) (see Figure 14). The cable length was minimized by Arbonès et al. (2016) by means of multiobjective optimization algorithms. As can be seen in Figure 14, the cable costs will increase with the number of WECs, but the increase is not linear due to the increased number of electrical substations that must be installed in the park at discrete intervals (Giassi et al., 2019a).
The levelized cost of electricity (LCOE) computes the cost to generate electricity over the life-time of an energy system and is a useful measurement to compare different energy sources (IEA and NEA, 2015). It is calculated by taking the ratio between the total capital (CAPEX) and operational (OPEX) costs of the device, discounted to present day value, and the amount of electricity E out delivered to the grid throughout the device's lifetime (Short et al., 1995), where t = 1, . . . n is the years of the lifespan, and r is the discount rate. However, identifying economic parameters and values for a technology that has not yet reached full maturity is difficult, both because very little data is available and because it is expected that the costs and uncertainties will be reduced as the technology matures (Astariz and Iglesias, 2015;Farrell et al., 2015;Ocean Energy Systems, 2015;De Andres et al., 2016). As discussed by Piscopo et al. (2018), the EU Strategic Energy Technology Plan expects the LCOE of wave energy to decrease to 0.15 EUR/kWh by 2030 and 0.10 EUR/kWh by 2035.
Despite these challenges, some authors have attempted to include economic parameters in their optimization algorithms. Teillant et al. (2012) modeled a wave farm of 100 units at the AMETS test site on the west coast of Ireland, and an economic assessment based on production, location, and maintenance strategies was presented. Rinaldi et al. (2016b) presented a computational tool focused on maintenance operations that computes key performance indicators such as annual electricity generation and total gross revenue, with inputs such as the failure rates, wave climate, and power matrix of the device. A reliability and financial model was developed by Macadré et al. (2015) FIGURE 14 | Wave energy park optimization from the perspective of minimizing cable costs. (Left) Arbonès et al. (2016) performed both the minimization of cable length and park area and the maximization of power output in a multi-objective optimization. The figure shows the minimization of cable length, computed as a Euclidean minimum spanning tree, with three different optimization algorithms. (Right) Cable length as a function of the number of WECs and associated costs for a park of point-absorbers up to 100 WECs. The sharp increase in cable lengths and costs corresponds with the additional number of electrical substations that must be installed in the park with an increasing number of WECs (Giassi et al., 2019a). and used to study a combined wind and wave power plant with respect to availability and maintenance. An LCOE optimization in terms of the device dimensions was carried out by Piscopo et al. (2018). Faraggiana et al. (2019) minimized the LCOE for a system of three WaveSub WECs. The LCOE was computed based on a formula from Marine Power Systems Ltd, but the input values were not stated explicitly in the paper. The LCOE of four parks of 10 WaveEL WECs each were compared for two wave climates by Yang et al. (2019). It was concluded that the LCOE values could be reduced or increased by drastic amounts if the hydrodynamical and mechanical coupling between the WECs was taken into account. The LCOE of wave energy parks of up to 50 WECs was optimized using a GA by Giassi et al. (2019a). Results from Faraggiana et al. (2019), Giassi et al. (2019a), andTeillant et al. (2012) are shown in Figure 15. As can be seen from the figure, the LCOE decreases with the number of generations in the optimization algorithm (Faraggiana et al., 2019). The net present value (NPV) increases in a park with an increase in the number of units (Teillant et al., 2012) unless the shadowing effect in the park is too large, in which case the NPV decreases (Giassi et al., 2019a).

Multi-objective Functions
In the works above, the economic optimization of the full system with the LCOE as the objective function is still performed as a single-objective optimization. To the authors' knowledge, the only attempt to study the park optimization problem in a pure multi-objective way was presented by Arbonès et al. (2016Arbonès et al. ( , 2018 for a submerged three-tether heaving buoy. Arbonès et al. (2016) optimized three competing objectives (the total energy production, the cable length, and the marine area needed to place the buoys) in long-crested waves with two different evolutionary multi-objective optimization algorithms (SMS-EMOA, MO-CMA-ES). The results showed that the multiobjective method improved the power output of a singleobjective optimization by 3.8%. Arbonès et al. (2018) extended the work by considering realistic wave conditions and incoming waves from multiple directions. Improvements in terms of power output, minimized cable length, and marine area were achieved through the methodology for parks of up to 36 converters.

CHALLENGES AND FUTURE WORK
The state of the art of wave energy array modeling was reviewed by Babarit (2013). He concluded that most of the work concerned small arrays and that more research should be carried out for large arrays to increase the reliability of the results. Since then, many important steps have been taken to achieve reliable, realistic, and relevant optimizations of large wave energy parks, as has been discussed in the previous sections. In this section, we will discuss some challenges and constraints to this work and list some ideas for future directions in the research area.

Computational Costs
As the wave energy systems increase in size, so do the computational costs. Both with an increased number of WECs and an increased number of parameters to optimize, the time required for the optimization algorithms to converge quickly grows out of hand. Evaluation of each configuration is expensive, and the search space is non-convex and multi-modal (Neshat et al., 2019c). Run times of 1-3 days on 20 parallel processors were reported by Faraggiana et al. (2019) to optimize several parameters of an array of three WECs. The approach of Neshat et al. (2018Neshat et al. ( , 2019a was to compare different modeling and optimization methods and constrain the computational budget to three days on a moderately high-performance sharedmemory platform, using 12 processors in parallel; see Figure 16.  (Teillant et al., 2012). Lower right: Comparison of the NPV for different park sizes up to 25 WECs when assuming no vs. full hydrodynamic interactions among the devices. The incoming wave directions are either perpendicular to the array layout ("y-waves") or along the direction of the array layout ("x-waves") (Giassi et al., 2019a).
A method able to find optimal solutions within this time frame would enable developers to analyze and evaluate two scenarios per week.
Much work has been carried out to speed up the optimization process. Ruiz et al. (2017) reported that the optimization time could be reduced if the hydrodynamical interactions were evaluated with an analytical method instead of BEM software. Also, a comparison was made on the computational cost for three global optimization algorithms, where the CMA method achieved convergence much faster than GAs or particle swarm optimization. Several developments to enable modeling of large arrays have been reported where the hydrodynamical interaction is restricted to nearest neighboring WECs (Sarkar et al., 2016) or within a specified interaction distance (Göteman et al., 2015a). Reducing the computational cost of optimizing large wave energy parks was one of the main objectives of Wu et al. (2016), where an approximation in the buoy interactions model resulted in a 350-fold computational speed-up, enabling optimization of parks with 100 WECs with acceptable accuracy. Two global optimization algorithms were compared by Ferri (2017) with respect to computational time and accuracy. The separation distance and skew-angle in the array layout were studied, and it was found that the meta-model-based optimization was 3-5 times faster than the studied GA but not as accurate. Neshat et al. (2019a) equipped the optimization method with an adaptive neural-surrogate model based on a machine learning approach, which estimated the absorbed power of the park instead of carrying out time-consuming exact computation. This resulted in faster convergence compared to the off-the-shelf evolutionary algorithms. A different strategy was introduced by Neshat et al. (2019d) in terms of a combined symmetric local search, a Nelder-Mead and a cooperative co-evolution algorithm, which improved the convergence rate and resulted in significantly increased absorbed power by the parks. FIGURE 16 | (Left) Convergence rates of three optimization algorithms (CMA, GA and particle swarm, GSO) for the optimization of the number of WECs and their positions in a park, showing that the CMA requires fewer iterations to converge but results in slightly lower park performance . (Right) Convergence rate for many algorithms to optimize the layout and PTO parameters for a system of four WECs (Neshat et al., 2019d).
It can be expected that future works will continue to investigate and develop methods based on machine learning or other advanced algorithms to enable faster optimization of large parks. To ensure that the models are still realistic representations of the wave energy systems, the same systems should be modeled with several methods and the impact of approximations or constraints thoroughly evaluated.

Reliable Data
Lack of available realistic data is a problem for the full wave energy sector, not only in the optimization of wave energy parks (Michelez et al., 2010;Sharp and DuPont, 2018). Data are required regarding performance, failure rates, power ratings, lifetime expectancy, capital and operational costs, and so on. In the traditional offshore oil and gas industry, the OREDA (Offshore REliability DAta) database was developed as a collaboration between global oil and gas companies with the purpose of collecting and exchanging reliability data among the participating companies (OREDA, 1984). The entries in the database were anonymous and could not be tracked back to the companies, and the process developed was then formalized in the international standard ISO 14224. Similar approaches have also been taken for offshore wind in the System Performance, Availability and Reliability Trend Analysis database (SPARTA, 2014) and also recently for tidal current turbines (ORE Catapult, 2015). For wave energy, similar initiatives have not yet been taken, mainly due to a lack of convergence on the technological concept, lack of experimental data (in particular long-term offshore data), and the competitive nature of the wave energy industry in its early stage of technology development.
For completeness, in addition to the papers discussed above, economic values regarding wave energy systems can be found in Teillant et al. (2012) Piscopo et al. (2018), although these papers and reports are not primarily focused on wave energy park optimizations but rather on risk assessments and reliability.

Repeatable and Unbiased Studies
To reduce the uncertainty in numerical simulations of wave energy systems using CFD, several blind test workshops and paper series have recently been presented under the name CCP-WSI Blind Test Series (CCP-WSI, 2019). The aim of the series has been to invite the research community to use their numerical codes to simulate a series of specific wave structure interaction problems, using only relevant input data, and where the key outputs can be compared to experimental results only after the community have submitted its numerical results. Similar blind test programs with the aim of reducing biases and achieving repeatability in the results produced by different research groups would be of high value for wave energy park optimization research.

CONCLUSIONS
In this paper, we have reviewed the state of the art of wave energy park optimization, with a focus on developments in the last five years. After discussing the different modeling and optimization methods that are used, we have analyzed whether the methods and results are realistic, reliable, and relevant.
It is not feasible to include all aspects of the waves when modeling wave energy parks. Whereas viscosity and rotation are important fluid properties when considering the detailed response and loads of single WECs, potential flow theory is the dominant assumption for large arrays of WECs. However, wave directionality, waves propagating in several directions simultaneously, and waves propagating over a long period of time should be considered when designing parks, and important steps have been taken in the last five years on realistic representation of the waves. Also, coupling high-fidelity fluid models close to the WECs to low-fidelity wave models further away from the bodies is a promising direction for a more realistic representation of the fluid-structure interactions. Optimization of large wave energy parks is still often carried out in the frequency domain with linear spring-damper systems as PTO. Some validation of this assumption has been carried out, but several works have also shown that linear frequency domain analysis cannot capture the non-linearities that are inherent in the system. Also, mooring dynamics have been shown to significantly affect the performance of parks and should not be neglected, and several important contributions have been published lately.
Straightforward comparison of different research groups' results is challenging, since the wave energy parks studied by the groups usually differ in terms of significant characteristics. However, some tendencies can be seen when comparing similar works.
• Optimal layouts will position the WECs along lines perpendicular to wave direction at sites with a narrow interval of incident wave directions, but less obvious optimal layouts will emerge at sites where waves are coming from multiple directions. If there is not sufficient ocean area for perpendicular lines, the WECs will align along several lines, where the separation distance will depend on the predominant wave conditions. • However, these "perpendicular" layouts were in general obtained using single-objective functions, optimizing only the power output or interaction factor. Other optimal solutions might emerge when economic or multi-objective functions are considered, for example, when also minimizing the cable costs and losses. • Several works have shown that the total performance might be improved if devices of different individual dimensions are installed in the park. • The performance of large parks might be slightly lower in short-crested waves as compared to long-crested but with desired lower power fluctuations. • Some results have been published in recent years that have provided experimental data on large arrays and/or arrays with complex dynamics or PTO systems, but there are still very few publications with validation of park modeling and optimization. Although reasonable agreement has been shown, some works have also reported non-linear dynamics or deviations between numerical and experimental results.
As discussed above, the large, non-convex, and multimodal optimization problem of wave energy parks might be too large for conventional population-based optimization algorithms, and new approximate or machine learning methods might be useful complements for the evaluation of park performance, in particular for large parks. Simultaneously, careful consideration must be taken when making simplifying assumptions, as it has been shown that non-linear dynamics, mooring systems, and other complexities cannot be neglected. Validation of wave energy park modeling and optimization is still only in its initial phases. To increase reliability, more systems with increasing complexity should be validated, and the modeling uncertainties should be adequately quantified. In addition, biases are reduced and the reliability improved when different methods or systems are compared and shown to give similar results. The relevance of the modeling results increases when the correct quantities are optimized. Multi-objective and/or economic cost functions should be used to provide useful guidelines on how wave energy parks should be designed, given a site and other constraints. If the research area is to move from developing methods for wave park optimization to actually making useful predictions that can guide developers, it is of crucial importance that data and economic figures are made available and are used in the optimization. In recent years, several important contributions have been made that not only optimize the power output or interaction factor but also consider grid connection, cost and length of installed cables, and life-time costs and revenue. This is a direction in which we will hopefully see more work in the future.

AUTHOR CONTRIBUTIONS
The work in this review paper was initiated and lead by MGö. MGö wrote the manuscript with contributions from MGi, JE, and JI. All authors have proofread the final paper and are accountable for the content of the work.

FUNDING
The research in this paper was supported by the Swedish Research Council (VR, grant number 2015-04657), STandUP for Energy, and the Swedish Energy Authority (project number 40421-1).