Meta‐heuristic optimisation approach for wave energy converter design by means of a stochastic hydrodynamic model

Correspondence Marcos Blanco, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avda. Complutense 40, 28040 Madrid, Spain. Email: marcos.blanco@ciemat.es Abstract A methodology to automate the sizing of a wave energy converter (WEC) is presented. The preliminary design of WEC geometry is addressed as a mathematical program by means of a meta-heuristic optimisation, which includes as key aspects the use of a multi-objective differential evolutionary (DE) algorithm and the implementation of a stochastic hydrodynamic model of the WEC. In this methodology, the WEC concept, the control strategy, the power take-off characteristics and the location have been taken into consideration due to their strong influence on the WEC design. As a case example, a real project has been selected, and the methodology has been applied to two test locations. Stochastic models are used to handle irregular sea states and to obtain statistical information while reducing the calculation effort when compared to time-domain models. The output derived from the proposed methodology is a set of feasible solutions so-called Pareto frontier. These solutions can be used as inputs for a detailed design process. Data from real projects could be used to further develop the sizing methodology by improving the algorithm convergence, refining the requirements, or including new restrictions.


INTRODUCTION
Electricity generation from ocean-wave energy has been developed for more than two centuries, experiencing a growing interest since the 1970s [1][2][3]. One of its main attractions is the considerable energy resource available in the ocean waves. In absolute terms, the global wave power incident on the world's coastline is estimated at ∼2.1 TW [4]. Development of wave energy is still considered to be at a pre-commercial stage, but showing a significant market potential [5]. Indeed, industrial companies have started to show interest in the wave energy sector. In recent years, some examples of developments and collaborations between industrial partners and technology developers are the following: FCC has participated in the UNDIGEN project with Wedge Global; Iberdola participates as an investor in the start-ups Oceantec and CorPower; and ENEL has collaborated with waves for energy in the development of the ISWEC prototype. Other signs of maturity in this sector can be found in the development of project guidelines, (EMEC [6], Carbon Trust Flow chart of marine power project development according to EMEC guidelines, adapted from [6] order to maximise the energy capture. In this context, automatic WEC-sizing methodologies, as the meta-heuristic methodology presented in this paper, might become a particularly helpful assist-tool at different WEC design project stages. In fact, the automation of the prime-mover sizing becomes specially useful in WEC design projects due to the aforementioned dependency on the location that often requires device re-design in each new location. According to the design and development flowchart for marine energy projects proposed by EMEC [6] (see Figure 1) the applicability of the proposed methodology extends along several stages of project development. First, a basic version of the sizing methodology, as outlined in [16], could be applied to tasks such as 'site screening' at Stage 1 or 'outline design' at Stage 2 in order to assist in the conceptual design process. Second, the automated WEC-sizing methodology could be applied at project Stages 3-4 to provide a set of preliminary dimensions, according to 'project design' or 'preliminary design' requirements, which are inputs to the 'detailed design'.
Optimisation algorithms have been proposed as automatised approaches to tackle different WEC design tasks [17]. For the sake of example, they can be applied to the design of specific components such as moorings [18] or PTOs [19][20][21]; or to WEC design steps such as wave farm layout design [22,23], or WEC structural reliability studies [24]. In the particular design step relative to WEC sizing, optimisation algorithms have been applied to different WEC technologies. Examples of point-absorber [25] shape optimisation studies are: [26], where a heave oscillating water column (OWC) point-absorber geometry is optimised by two multi-objective algorithms; [27], where a heave IPS point absorber is optimised by a multi-objective genetic algorithm; or [28] where a one-body surge-pitch point absorber is optimised by a single-objective genetic algorithm. Additional examples of the application of optimisation algo-rithms to alternative WEC technologies can be found in the literature, such as the shape optimisation of a one-body inertial rotating-mass device by means of a multi-objective algorithm based on a hybrid genetic/gradient-descent method [29]; or the shape optimisation of a submerged pitch-oscillating device [30] by means of a multi-objective evolutionary algorithm. Finally, examples of WEC sizing or WEC designs based on brute-force search, parametric analysis or other approaches can also be found in the literature [31][32][33][34][35][36].
The developed methodology addresses WEC sizing as a mathematical multi-objective program (or mathematical optimisation problem). The optimisation problem is solved by means of a differential evolutionary (DE) algorithm, where location and the PTO technical boundaries act as constraints. The complex mathematical models used to represent the WEC dynamics and involved in the WEC-sizing process, prompt to use of meta-heuristic optimisation algorithms [17,22,28], and DE has been selected from all the potentially suitable options-for example, genetic algorithms-because of its robustness and rapid convergence [26,37]. A stochastic hydrodynamic model is implemented in order to evaluate the WEC performance under ocean irregular-wave scenarios. This model improves the accuracy of the estimated energy extraction when compared to a previously used frequency-domain model using a single regular wave to represent the sea conditions [16]. It incorporates more detailed information-to be used at the 'detailed design' stage in Figure 1-from various WEC sub-systems by means of the obtained statistical information, such as the likelihood of exceedance of the PTO stroke constraints. This stochastic model uses a linearised hydrodynamic model and considers the irregular sea state spectrum [26,38,39].
The basic sizing method proposed in this paper is based on several previous works which are summarised in [16]. The methodology employed in [16] is modified in this work by incorporating a stochastic hydrodynamic model of the WEC and integrating the statistical analysis of its results. For example, percentage of occurrence of maximum amplitudes in velocity, displacement or force as calculated by means of the Longuet-Higgins equations [40]. In terms of computational cost, stochastic models-evaluated with irregular wav-despite being more demanding than frequency-domain models-evaluated with single-frequency regular waves [38]-are still a better option than computational fluid dynamics (CFD) programs or timedomain models [41]. Moreover, stochastic models provide statistical information and, when evaluated under irregular waves, offer more realistic results than frequency-domain models that consider regular waves. These more realistic results provide an additional level of detail that enables the utilisation of the methodology at more advanced stages of a WEC development project, such as at 'project design' or 'detailed design' stages than the proposed ones for the previous methodology; for example, at 'project design' or 'detailed design' stages. In addition, to facilitate a proper utilisation of the derived statistical information, all the constraints have been redefined in order to incorporate the available information on the different sub-systems: for example, information about the over-stroke probability can be linked with the impact absorption capacity of the end-stops.

FIGURE 2
Simplified prime-mover geometry and its dimensional variables In summary, the paper presents a WEC-sizing methodology based on a meta-heuristic algorithm using a stochastic hydrodynamic model. This methodology is applied for a particular study case described in Section 2; the stochastic model of the specific WEC technology of the study case is described in Section 3; the analysis of the wave climate at the location considered in the study case is presented in Section 4; the sizing methodology, including the definition of the optimisation problem and the optimisation algorithm description, is presented in Section 5; at last, results and conclusions are summarised in Sections 6 and 7, respectively.

CASE STUDY CHARACTERISTICS
This preliminary sizing method is framed within the 'project design' task (Stage 3), whereas the WEC design characteristics would have been pre-set in previous project stages: location is defined at the 'site screening' task (Stage 1) and prime mover (i.e. the WEC part that directly interacts with the ocean waves) and PTO concept are defined in 'technology selection' and 'outline design' tasks (Stage 2). For the sake of example, a case study based on the project UNDIGEN [42] is described below, although the methodology could be directly applied to other WEC development projects: 1. The prime-mover concept is a two-body point absorber [1] shown in Figure 2. 2. Two sites are considered to evaluate the impact of the location on the geometry design, both of them being wave energy testing facilities: PLOCAN test site in the Canary Islands (Spain), and BiMEP test site in Biscay (Spain). 3. The PTO is a linear switched reluctance generator (LSRG) [43,44]. Table 1 shows PTO rated characteristics. 4. An R&D pilot plant with a single WEC connected to the grid is assumed. According to the classification established by MaRINET [9], the case presented in this paper would cor- During the 'project design', the geometry of the prime mover, the equipment and the rest of the WEC parts should be defined and specified. In the proposed preliminary sizing method, the definition of the prime-mover geometry consists in the extraction of the values of the basic dimensional parameters shown in Figure 2 (see Section 5.1). It provides the necessary information to continue with the specification of the WEC elements, such as PTO, moorings, ballast tanks, etc. Additionally, it provides the basic dimensional parameters to be used as the starting point in the 'detailed design' (Stage 4). At this stage, the complete fabrication drawings are ready to be generated from these basic dimensions. Apart from the previously described use, this basic sizing method could be used in parametric pre-analysis where the WEC concept would be evaluated at different locations, taking into account alternative control strategies, or in order to obtain the WEC dimensions for a range of target rated powers. The WEC concept and its associated dynamic equations and models are presented in Section 3. Models and equations allow for the evaluation of each WEC geometry solution. The WEC location and its ocean-wave climate is described in Section 4.

MODEL OF THE WAVE ENERGY CONVERTER CONCEPT
The preliminary sizing method calculates the dimensions of a simplified WEC prime-mover geometry. This simplified geometry has been generated from cylindrical bodies. The prime mover is the physical structure that interacts with the waves. In this case study, as seen in Figure 2, the two-body point-absorber prime mover is composed of a toroidal float body (body 1) and a semi-submerged body (body 2). Body 2 is composed of two coaxial rigidly connected cylindrical bodies: a semi-submerged one and a fully submerged one, which would be moored to the sea bottom. The only permitted relative motion between 'body 1' and 'body 2' is the heave movement, since 'body 1' sliding motion along 'body 2' is mechanically guided to block the rest of the degrees of freedom. The PTO linear generator is composed of two parts, with one of them being solidly attached to 'body 1' of the prime mover, and the other part solidly attached to 'body 2'. This PTO is able to generate electric energy from only this heave relative movement between the two bodies.

WEC prime-mover frequency-domain model for single regular-wave analysis
WEC-wave interaction is a complex, non-linear physical process. However, assuming that wave amplitudes are sufficiently small compared to the wavelength, the linear water wave theory can be applied. This assumption is reasonably valid for most of the time during which a WEC is operating [45]. The assumed linear model adopted has as its main limitation the exclusion of non-linear effects such as viscous dissipation or end-stop behaviour.
According to Newton's second law, WEC frequency-domain equations are shown in (1) and (2) in terms of phasors for a motion in a single degree of freedom in heave. This WEC frequency-domain model is evaluated for regular waves. The mooring effects are not considered since no significant impact on the power performance is expected on the supposition that moorings are properly designed [46,47]: for example, in [47], body 2 is slack moored with a sufficient length to accommodate displacements. Since the relative motion between the two bodies is restricted to heave, the surge and pitch degrees of freedom can be assumed to be decoupled from the former [21]. Therefore, surge and pitch motions can be neglected in the energy production assessment [48] of the considered WEC. In consequence, power performance has been evaluated only in heave for the WEC concept of the case study, which is the only motion the PTO can take advantage of in order to extract energy.
where subscripts 1 and 2 stand for body 1 and 2, respectively; represents a phasor variable (complex number); j is the complex number √ −1; F e,i is the excitation force of body i; D rii is the radiation resistance of body i (real part of the hydrodynamic radiation impedance); m adii is the added mass of body i (imaginary part of the hydrodynamic radiation impedance); D ril and m adil are the corresponding hydrodynamic impedance terms of body i with respect to body l ; D vi is the additional friction resistance of body i (associated to mechanical losses); D vPTO is the PTO mechanical friction resistance; m i , mass of body i; S i is the hydrostatic restoring force of body i; and F PTO is the mechanical force exerted by the PTO [45][46][47][48][49]. The corresponding hydrodynamic impedance terms are reciprocal, hence D ril =D li and m adil =m adli .
Equations (1) and (2) can be expressed in the matrix form (3) and (4): [ where Z il is the total impedance of the body i dynamics with respect to the body l velocity; Z PTO is the PTO impedance (F PTO = Z PTO ⋅ ( 1 − 2 )); and Z ′ il is the total impedance including the PTO impedance. The excitation force F e,i is defined as the hydrodynamic parameter f e,i multiplied by the wave amplitudeÂ. The total impedance matrix is symmetrical, The evaluation of the above expressions requires calculating the values of all parameters. The WEC hydrodynamic parameters f e1 , f e2 , D r11 , D r22 , D r12 , m ad 11 , m ad 22 , m ad 12 and S 1 are frequency dependent, uniquely defined by the WEC geometry. In the developed methodology, they are calculated by applying a Matching Eigenfunction Expansion method [49], but other methods such as BEM codes (e.g. NEMOH) are equally valid.
Based on this frequency-domain model-presented here to be evaluated with regular waves-a stochastic model has been set out. The stochastic model can handle irregular sea states like a frequency-domain model, and both models shorten the calculation time with respect to time-domain modelling [50].
However, the stochastic model shows the extra benefit of being capable of evaluating statistical information about the mechanical variables.

WEC prime-mover stochastic model for irregular-wave analysis
The stochastic model allows to take into account irregular sea states by means of their energy spectrum information. This allows to evaluate the power performance without needing time-domain models, which are usually more computationally demanding and, therefore, more suitable to be used in a detailed engineering framework where the non-linear effects should be analysed. The stochastic model is based on the fact that ocean waves can be modelled as zero-mean random Gaussian variables [39,40]. The ocean free surface movement ( ) can be characterised by its frequency harmonics (5). The contribution of each frequency to a wave amplitude in a certain sea state is obtained from the wave energy spectrum, S ( ) [51], as shown by (6).
where ℝ(Y ) is the real part of a generic phasor Y ;Â k is the amplitude amplitude phasor of the discrete frequency k , with amplitude A k and angle k ; Δ , discrete frequency step used in S ( ). In addition, WEC stochastic models are based on the ocean waves' statistical properties and on the frequency-domain transfer functions between wave amplitude and certain dynamic variables such as: the velocity of each WEC body (7,8); the displacement of each body (9, 10); the distance from floating body mass center to sea surface (11); and the exerted PTO force (12).
where H X stands for the frequency-domain transfer function, that is the relation between the WEC dynamic variable X , output, and the wave amplitudeÂ k , input. Consequently, mechanical variables can be characterised as spectral distributions according to their transfer functions, in the same way as the incident wave oscillations. Equation (13) shows a generic spectral function S X of a dynamic variable X with H X as the transfer function, and (14) shows its generic (nth order) spectral momentum (m Xn ).
Ocean surface oscillations are considered Gaussian variables and the WEC is modelled linearly. According to these premises, WEC oscillating mechanical variables can be statistically characterised similarly to ocean waves [38]. The joint probability density function (PDF) of the incident wave oscillations (specific values of period and amplitude) for each sea state can be defined for a generic dynamic variable X as per (15) [40]. Equation (17) evaluates its expected value E (X ).
where A X is the amplitude of X ; p X is the density probability of A X ; L X , normalisation factor; x , spectral width parameter. The expected value of the WEC generated electric power (P elec ) is obtained from the WEC mechanical power (P mec ) and the PTO power losses (P loss ). P mec is evaluated from the relative velocity and PTO force variance [38], as per (18). The damping PTO value (R PTO ) is obtained according to the control strategy detailed in section 3.3.2, which relates PTO force to relative velocity as per (19). P loss is evaluated as detailed in Section 3.3.1.
where D PTO and B PTO are real and imaginary parts of Z PTO , respectively; r is the relative velocity of body 1 with respect to body 2 ( r = 1 − 2 ); and m X,0 is the zeroth order spectral momentum (equal to variance) of X .
The consideration of irregular waves-used in a stochastic model in this case-provides more accurate results than the consideration of regular waves--used in a frequency-domain model in this case. Therefore, the former is preferable in a design process. Such accuracy could be further improved through representative sea state selection [50]. In addition, PDFs of the mechanical variables are useful to design auxiliary devices such as end-stops, brakes or a back-up energy systems [52]; for example, statistical information about the overstroke probability gives information about the expected number of impacts per unit of time.

WEC electric power train model
The PTO considered in the analysis is a LSRG developed by WedgeGlobalSL [44], one of the partners of the case study project UNDIGEN. Its rated characteristics or target main parameters, which have an important impact on the WEC design process and energy-extraction control strategy, are presented in Table 1 [44].

Power loss model of the PTO
The LSRG can establish its nominal force in tens of milliseconds, whereas the WEC oscillates with periods in the order of seconds (few seconds to dozens of seconds) [53]. Consequently, the LSRG force command, defined by the energy-extraction strategy, can be considered to be exerted instantaneously, neglecting the LSRG electro-mechanic dynamics. As a result, the PTO losses model in frequency domain given by (20) is used to evaluate the electric to mechanical power conversion. This assumes that the mechanical losses and the iron losses can be neglected due to the relatively low LSRG velocity. Copper losses are assumed predominant due to the high force over velocity ratio, and thus, due to the high current value [54,55].
where F PTO denotes the PTO force; R cu is the LSRG electric resistance; I PTO is the LSRG current; R ′ PTO is the losses coefficient of the LSRG; K I −F is the relationship between current and force in the LSRG; and R ′ PTO is the losses coefficient of the PTO.

Energy extraction control strategy
The control strategy of the energy extraction determines the instantaneous F PTO value and it is based on optimum control [45] with two modifications: (I) a force constraint is included as explained in [56] and (II) electric power is optimised instead of mechanical power. In order to define the control strategy, the two-body WEC model is simplified to a one-body WEC equivalent model. This equivalent model is represented by an equivalent excitation force and an equivalent intrinsic mechanical impedance, which are obtained as per [57]. The equivalence is obtained by means of the application of the Thévenin theorem to the WEC electric analogue or equivalent circuit [45]. Equations (21) and (22) present equivalent force (F ′ TH ) and impedance (Z ′ TH ), respectively, including PTO losses.
Second, the expression of |F PTO | is obtained in order to maximise the electric power and to derive the maximum F PTO constraint (rated LSRM force). Electric power in terms of F PTO magnitude and angle is shown in (23). Optimum F PTO magnitude (24) and angle (25) are obtained by applying Lagrange multiplier theorem.
where |F PTO-OPT | and PTO-OPT are magnitude and angle of optimum F PTO phasor. The previous mathematical expressions allow the calculation of the PTO control strategy parameters. Therefore, the values of the control parameters are evaluated from the values of the hydrodynamic parameters, which in turn are obtained from the WEC dimensions. This implies that, for every WEC geometry solution, the optimal values of F PTO are obtained for each frequency. Subsequently, Z PTO optimal values can be derived from (18). Given Z PTO , the rest of variables can be easily obtained, such as velocities (7,8), or electrical power generation (19,20). The computational cost for a single-solution evaluation with the presented stochastic model results in, approximately, 0.24 s of computing time.

Scatter diagram and wave spectrum analysis
The ocean-wave scatter diagrams represent the occurrence of each sea state energy spectrum S • PLOCAN seems like an adequate site for testing the energyextraction capacity of scaled prototypes, while BiMEP seems to be adequate for testing higher-TRL WEC prototypes and verifying survivability requirements, due to their different wave energy levels. From the 1-year data of the hourly measured S ( f ), it is possible to generate the occurrence scatter diagrams of Figure 3. In addition, from this data, it is possible to select the most suitable parametric wave energy spectrum. Pierson-Moskowitz (PM), Jonswap (JP) and Scott (SC) [51] have been considered as candidate parametric spectra. The minimum mean-square error of the difference between each candidate parametric spectrum and each hourly measured spectrum have been evaluated. PM and JP spectra were found to be the best options (i.e. the option

WEC power operation analysis
Due to the uncertainty of ocean-wave energy resource, WEC operation should adapt to wave energy level. Four operational states can be envisaged [7].
1. Stand-by: low wave energy level fails to compensate for PTO losses and therefore the PTO remains inactive. 2. Resonance: wave energy level generates electric energy with PTO, below rated power, using reactive mechanical power to maximise power extraction (reactive control [45]). 3. Saturated: wave energy level leads the PTO to generate its rated power, managing just active mechanical power; that is the PTO force is always proportional to the velocity (damping control [45,59]). Here 'active power' refers to the component of the instantaneous power that contributes to the time-averaged power. 4. Survival: under extreme wave conditions, the relative movement of the device is blocked.
Sea states, described in terms of H s − T p , have been associated with operational ranges in the locations as follows (Figure 4): • Sea states are sorted according to their wave power level (J ) [45], which is proportional to T p ⋅ H 2 s . • From the highest energetic sea state, and continuing in descending order, the sea states are assigned to state 3 (saturated) until the cumulative sum of their annual occurrence percentage reaches 30% (marked in Figure 4 with a green line). • From the last sea-state assigned to Stage 3, and continuing in descending order, the sea states are assigned to state 2 (resonance) until the cumulative sum of the annual occurrence percentage (of the sea-states assigned to the state 2) reaches 90% (marked in Figure 4 with a red line). • The rest of the sea states are associated with ranges 1 (standby) and 4 (survival). The distinction between these ranges implies a further analysis of the survival mode, beyond the scope of the project-stage context considered in this work.

PRELIMINARY SIZING METHOD
The preliminary sizing methodology obtains suitable basic WEC geometries according to certain basic given characteristics, such as location or PTO rated values. The method, implemented in programming code, automates the preliminary WECsizing process, which is set up as a mathematical optimisation problem. To formulate this problem, the following concepts need to be defined: • Search space variables: free variables, varied by the algorithm, that define each candidate solution. • Objective functions (cost functions): output variables to be maximised or minimised. • Constraints: rules that certain variables must comply with.
They define the zone of feasible candidate solutions. • Model: equations used to evaluate the objective function and constraints from specified values of search-space variables. The WEC mathematical model was previously described in Section 3. • Optimisation algorithm: mathematical strategy used to solve the optimisation problem, optimising objective functions and managing constraints.

Search space variables
Each candidate solution is composed of simplified cylindrical prime-mover geometries, defined by their radii and heights (see Figure 2). The spar radius (R 3 ) is fixed to host the PTO and therefore should not be varied. Thus, the search-space variables that define each solution are: R 1 , d 1 , R 2 , d 2 , and d 3 . In this case, the search space is a subset of a five-dimensional Euclidean space delimited by the constraints described in Section 5.3.

Objective function
In mature technologies, the function to be optimised (through the so-called 'a priori' methods [60]) will be purely economic: for example single-objective optimisation problems such as the minimisation of the levelised cost of energy or the maximisation of the profits. The use of this type of variables requires detailed techno-economic assessments [50] in order to maximise energy extraction while minimising total costs, such as cost of materials, production, operation, etc. However, for a technology under development such as wave energy extraction (and considering that the methodology is framed under a preliminary design stage of a WEC prototype), setting up a multi-objective optimisation mathematical problem and using an 'a posteriori' method could be desirable. It permits to obtain as much information as possible from the method in a first instance, and to make design decisions about the solutions found at later design stages: for example detailed design stage according to Figure 1.
The following simplified objective functions-representative of the total cost and profits-are used in this work: • The objective function associated with costs is the WEC volume, V, (26), since it will have an effect on the amount of material needed, the size of moorings, etc.
• The objective function associated with income is the annual generation of electric energy (E elec ) (27).
where E elec is the annual energy production; E[P elec (T p , H s )] is the expected WEC electric power matrix (20) evaluated for each sea state; and h waves (T p , H s ) is the annual sea-state occurrence matrix for the location. It is worth mentioning that choosing the volume as proxy for the WEC cost is relatively common at preliminary-stage developments [28,61]. However, other authors prefer to use surface as the proxy since it is directly related to the amount of material required to construct the WEC external surface [30,62]. In addition, in complex WEC geometries (such as complete hull geometries) volume-driven optimisation might lead to complex, less manufacturable and more numerous wet surfaces [62]. On the contrary, simple shapes do not show significant discrepancies between the use of the two proxies [16]. In any case, the proxy for the WEC cost can be easily updated and, hence, in future works, further studies would be performed in order to ascertain what the most suitable option is.

Constraints functions
The constraints are equality or inequality expressions that define the feasibility of candidate solutions. Non-feasible solutions are discarded. In the present case study, these constraints are formulated in terms of variables from the WEC model described in Section 3, which is fed with the values of the search-space variables, in this case, the WEC dimensions. In this work, four inequality constraints were defined with a strong focus on WEC suitability for the given location/PTO and on the minimisation of the amount of critical situations such as impacts with the end-stop devices due to large over-excursions. The implemented constraints have adjustment parameters that modify the limit values of the constraints. These parameters can be used to relax the constraints and then enlarge the feasible search space.
In addition, the parameters can be used to adapt the limit values in the constraints to the rated values of the equipment. In the analysis of the case study presented in Section 6, the adjusting parameters values have been chosen as one, and hence have no influence. Moreover, ranges of the search-space variables present additional direct restrictions, with the aim of gaining control of the size of the search space. Nevertheless, it must be remarked that this kind of constraints have no influence on the optimal solutions since bounds are broad enough so that they are never reached. Variables' upper and lower bounds are shown in Table 2. R 3 is fixed to 0.75m in order to host the PTO.

Minimum expected electric power
This constraint (28) ensures that the PTO works in a regime near to its rated operating conditions, increasing the equivalent (full-load) hours. It is imposed only in the range of 30% of the most frequent sea states (defined as range 3: Saturated State; as per Section 4.2). The expected extracted power (19) should exceed a minimum value P min (a percentage of P rtd ,see Table 1, defined by K p ; P min = K p ⋅ P rtd ). The parameter K p is the adjustment parameter of the power constraint, which is set to one as previously mentioned.

Maximum expected velocity
The aim of the velocity constraint is to avoid excessive heating and associated damage in the linear bearings that guide the LSRM parts: stator, attached to body 2, and translator, attached to body 1. The heating process is slow enough to be represented by the stochastic expected value. This constraint is imposed in operating ranges 2 and 3 (see Section 4.2), which represents 90% of the most commonly occurring sea states. The WEC dynamic of these sea states should comply with (29) where the expected relative velocity (see (7), (8) and (15)) should not exceed a maximum value max . max is defined as a percentage of rtd , see Table 1, defined by the parameter K ; max = K ⋅ rtd . The parameter K v is the adjustment parameter of the velocity constraint.

Probability of maximum displacement amplitude
Firstly, the maximum displacement constraint establishes a limit value (n smax ) to the probability of amplitude oscillations exceed a maximum value. This probability value is calculated from the stochastic WEC model (15). The maximum value of the amplitude-oscillation is the rated PTO value s maxSTO (Table 1). This constraint (30) limits the over-stroke probability evaluated in sea-states belonging to ranges 2-3 (see Section 4.2).
where n smax is the maximum probability of over-stroke and s max PTO (equal to s rtd in Table 1) is the limiting distance in the relative movement. Second, the distance between the position of body 1 (floating body) and the position of the free water surface is limited, considering that the position of body 1 is measured from its equilibrium position. The frequency-domain transfer function of this variable is shown in (11). This constraint (31) is set in terms of probability of occurrence, similarly to the previous constraint (30), evaluated in ranges 2-3 (see Section 4.2).
where n SL max is the maximum probability of slamming; and d max STO (floating body draught d 1 ) is the limit value of the previously defined distance. This constraint reduces the probability of slamming [63] which would damage the WEC structure. The stochastic model, as linear model, could overestimate the displacement amplitudes since non-linear effects (e.g. viscous effects) are neglected. However, this only happens under large ocean-wave amplitudes, which do not usually occur during the majority of the WEC generation time. In consequence, it can be concluded that these constraints impose upper limit values in the displacements, which implies working with a high safety factor in the project design process (conservative approach).

Multi-objective optimisation algorithm
The implemented optimisation algorithm is based on a differential evolution (DE) algorithm upgraded and modified according to NSGA-II to manage several objective functions [64] and to manage constraints according to Deb's rules [65]. The initial data of each iteration (generation t ) is a set, called parent population P t , of WEC design solutions previously evaluated (see Figure 5). Each solution is defined by specific values of the search-space variables.
In the first step, the DE algorithm generates the new candidate population (Q t ) from P t by means of a mutation operation, followed by a crossover operation [66]. The main difference of the DE algorithm versus other bio-inspired algorithms (such as genetic algorithms) is the way the mutation mechanism works. To generate each candidate solution of Q t , first, three solutions belonging to P t are randomly selected. The mutant solution is the result of the vectorial sum of the first solution and the subtraction of the two remaining solutions multiplied by a certain factor (which has been set at 0.8 for robustness [66]). Second, In a second step, Q t is evaluated to obtain objective functions' and constraint' values. These values are calculated by means of the WEC mathematical model for the new combinations of search-space variable values. Third, the combined population R t (P t and Q t solutions) is sorted according to dominance rules into groups or frontiers. In addition, each frontier of R t is sorted according to the crowding distance [64]. The dominance rules allow for obtaining solutions closer to the real Pareto frontier (PF), while the crowding distance sorting increases the diversity of the solutions.
In summary, the multi-objective DE algorithm has to solve an optimisation problem with two objective functions (see Section 5.2), fulfilling four inequality constraints (see Section 5.3), in a search space of five dimensions (see Section 5.1). The optimisation algorithm begins to generate a randomly distributed initial population P 1 over the search space. As in the rest of the iterations, in this initial iteration the objective functions and constraints of each solution are evaluated by means of the WEC dynamic model (described in Section 3). From this point, the algorithm follows the previously explained steps: DE algorithm generates the candidates solutions of Q 1 ; the best solutions of P 1 and Q 1 are selected to conform P 2 by applying dominance and crowding distance rules; and the process continues iteratively until the maximum number of iterations (1000 iterations in this case) is reached. The final information obtained from the optimisation algorithm is a set of feasible solutions and its socalled Pareto frontier. These solutions are Pareto optimal, that is no alternative solutions has been found that would surpass them in all objective functions at the same time.

RESULTS OF THE CASE STUDY
The results of the optimisation process in both locations are presented in Figures 6 and 7. An analysis of a suitable WEC geometry is shown in Figures 9 and 10. In Figure 6, the feasible solutions are projected over the plane defined by each pair of combinations of the five search-space variables. These projections are shown as a scatter plot matrix of the search-space variables for each location-PLOCAN and BiMEP. In addition, the case of the common feasible solutions, that is WEC geometries that comply with constraints in both locations, is also shown. The whole BiMEP solution's zone and its PF solutions are contained in the common feasible zone while a large number of PLOCAN solutions, even PF solutions, fall outside. Therefore, BiMEP Pareto frontier solutions could be located at the PLOCAN site without violating the constraints. The generated energy maximisation tends to increase R 1 in order to capture as much energy as possible. Hence, large WEC solutions would be interesting in locations with less energetic sea states such as PLOCAN, and these solutions, located in a more energetic location such as BiMEP, can go through larger excursions than those allowed by the design constraints.
The variable d 2 remains at a very narrow range around its minimum value for the two locations. This low value reduces as much as possible the volume since it does not impact on the energy generated significantly. Additionally, in the two locations the draught of Body 2 (d 3 ) increases with R 1 to reduce the excitation force of Body 2 and, therefore, to reduce excursions. Figure 7 shows the obtained Pareto frontier (dark grey curve) and feasible solutions (light grey points) with respect to the two objective functions for each location: (a) for PLOCAN and (b) for BiMEP. The solutions optimised at one location have also been evaluated at the other one. In Figure 7, this has been highlighted in terms of Pareto frontier (dark blue line) and feasible solutions (light blue points), that is the solutions of the common feasible area and its PF are shown in blue colour. As can be observed in the figure, the area under the grey curve in PLO-CAN is quite separated from the area under the blue curve, whereas both areas almost overlap in BiMEP (similarly to Figure 6). The left-hand side part of the PF curve (less energetic and lower-volume solutions) are bounded by the capacity of the WEC to extract energy, while the right-hand side part (more energetic and higher-volume solutions) is bounded for the fulfilment of the constraints, especially velocity and displacement constraints. For this reason, the WEC geometries on the left part of the PF are prone to be employed in alternative additional locations, since they are not affected by the site-related constraints. This is a direct consequence of the conclusion highlighted in the analysis of Figure 6: the solutions obtained for only BiMEP are contained in the common feasible area for both locations. The design of a suitable WEC for more than one location implies a reduction in the extracted energy, with a more severe reduction in the least energetic locations. The common Pareto frontier overlaps the specific Pareto frontiers in the lowvolume and low-energy solution range; therefore, these solutions are suitable for both locations. For the sake of example, a suitable versatile solution in this range is highlighted in Figure 7, with approximately 350m 3 of volume. Nevertheless, WEC solutions outside the common area could be selected, ignoring the 'versatility' criterion. An acceptable alternative criterion could be to select the WEC solutions located at the 'corner point' of the PF, where the slope of the energy versus volume plot presents an abrupt change (at WEC volumes of aroSd 2000m 3 for PLOCAN and 1000m 3 for BiMEP). Compared to the 'versatile' solution, the 'corner point' solution in PLOCAN is capable of extracting 50.15% of extra power, but with 488.2% of extra volume. In the case of BiMEP, the 'corner point' solution is able to extract 23.8% of extra power, but thanks to an increase of 208.2% in extra volume (compared to the 'versatile' solution).
The constraints taken into account in the methodology, in particular the constraints related with displacements or velocities, lead to small WECs in highly energetic locations. Therefore, the higher the energy content of the location, the greater the impact of the constraints on the feasible space. In the present analysis, this issue implies that the PLOCAN PF solutions are outside the common feasible solutions zone (see Figure 7a), and vice versa, the BiMEP solutions are practically totally included in the common feasible zone (see Figure 7b). However, it is worth highlighting that the predominant periods of the two locations are similar (8 vs. 9 s), which has an effect on the size of the common feasible zone. The common feasible zone could be considerably reduced if studying two more dissimilar locations; that is with clearly differentiated ocean-wave climates.
These solutions can be used as an input in a detailed design process, so it is important that the aforementioned set of feasible solutions is as wide as possible. The detailed design process would include extra analyses and constraints that may change the feasible solutions zone and therefore modify the Pareto frontier, such as reinforcement structure design, pitch stability analysis, mooring design, etc. A high slope in the PF curve implies that, for an increment of volume (investment), a worthy increment in generated energy (income) is obtained. In Figure 7, the ratio energy/volume decreases as the volume and energy increase, presenting an abrupt change (or 'corner point') at volume values of around 2000m 3 for PLOCAN and 1000m 3 for BiMEP. At higher volume values, the ratio energy/volume decreases significantly, and the PF enters an area where the increase in extracted energy is almost negligible. One factor influencing this behaviour is the problem definition itself. In the present study, the WEC size optimisation is taking into account pre-defined PTO characteristics: the maximum force condition is directly applied through the control strategy, so that the PTO never exerts forces over its rated value. In addition, there is one additional constraint related to the PTO rated velocity. In conclusion, the rated power of the PTO is restraining the WEC's maximum achievable power extraction. For the sake of example, Figure 8 shows the constraint values of the WEC solutions obtained for BiMEP location. This figure shows that the velocity and maximum displacement constraint values are activated in the zone of low ratio energy/volume (high-volume and highenergy solution range).
Therefore, the most worthwhile zone of the Pareto frontier corresponds to volumes under the aforementioned 'corner point' values. However, as previously mentioned, considering lower volumes leads to an over-sizing of the PTO because the operating time at rated power will be limited. In other words, constraints, whose aim is to adapt the WEC geometry to the nominal characteristics of the PTO, will have a reduced effect in the zone of smaller volumes. In consequence, solutions with volumes around the corner point zone could be suitable tradeoff solutions in order to maximise the energy/volume ratio and to limit the PTO over-sizing.
For the sake of illustration, and following the previous comments on the Parent Frontier shape, a WEC geometry has been selected (in the absence of any additional criteria) and marked with a yellow point in Figure 7a and b (see Figure 7c for the 3D corresponding design). This is a trade-off solution following the aforementioned criteria: feasibility for the two locations, maximisation of the volume-energy slope, and maximisation of the utilisation of PTO at rated power. The corresponding WEC dynamic behaviour is shown in Figures 9 and 10 in PLOCAN and BiMEP, respectively. These two figures show the matrices of the expected values of: mechanical extracted power (upper left), electric extracted power (upper -right), relative velocity (lower left), and probability of not satisfying the relative velocity constraint (lower right).
Stochastic models are considered accurate enough to evaluate WEC performance [38,50], and they provide useful statistic information for the detailed design process: for example the probability of over-stroke probability could help in end-stop specification. In addition, data from real projects could be used to improve the sizing method in aspects such as design parameter matching or the inclusion of new constraints. Results can be used as baseline information for a detailed WEC-sizing process (supported by CFD programs), thus helping to reduce the detailed design workload.

CONCLUSIONS
The proposed WEC-sizing methodology automates the calculation of a WEC prime-mover dimensions according to its given characteristics. This method considers WEC concept, control strategy, PTO characteristics, and location. Correlation between the WEC design and the location often requires device redesign at every new location. Hence, automation of the aforementioned prime-mover sizing is extremely useful in WEC design projects. Such correlation, consequence of the energy resource nature, suggests that future extraction of energy from waves could be probably linked to bespoke WEC prime-mover designs, potentially combined with certain off-the-shelf (and/or modular) components such as the PTO or the end-stops. The method presented in this paper automates the production of suitable baseline solutions, which become, in turn, inputs to the detailed sizing. The stochastic hydrodynamic model improves the energy capture estimations-due to the consideration of irregular-wave conditions-and provides statistical information about the WEC mechanical behaviour with a low computational cost. The statistical information allows to integrate extra information on different components and sub-systems, such as the end-stops, that would otherwise be lost. For these reasons, the method is suitable to be used framed at the project design stage of a WEC. Although in this case the preliminary sizing method is presented at the project design stage, as a future work option, this method could be integrated into the detailed engineering stage, along with a CFD program, to undertake a more detailed sizing process. Additional future work options for development include the consideration of alternative WEC concepts, other energy-extraction strategies, the extension of the method to wave farm planning, or the joint optimisation of the PTO and the prime mover. In addition, further developments could include: the increase in degrees of freedom of the model in order to include the pitch and surge motions for stability analysis or to evaluate their influence on WEC performance; structural analysis by means of CFD programs; integration of the mooring effects in the dynamic model and the analysis of its characteristics as constraints; or survival mode analysis with non-linear numerical models. These improvements expand the applicability of the method to other design project stages, such as: site screening, outline design or the technology selection, where the method would permit to automate the process of parametric analysis in terms of WEC concept, PTO control or location evaluation.