Shape Optimization of Oscillating Buoy Wave Energy Converter Based on the Mean Annual Power Prediction Model

: In order to improve the energy capture efﬁciency of an oscillating buoy wave energy converter (WEC), a buoy-shape optimization design method based on the mean annual power prediction model is proposed. According to the statistical data of long-term wave characteristics in the Chinese sea area, the optimal design space is determined. Sixty-three sample points were randomly selected in the optimized space. Based on simulation, the mean annual power corresponding to each sample point is calculated to quantitatively describe the energy capture ability. The response surface method (RSM), radial basis function neural network (RBFNN), and elliptical basis functions neural network (EBFNN) are used to establish the mean annual power prediction models, respectively. By combining the prediction model with the multi-island genetic algorithm (MIGA), the optimal solution in the design space is easily obtained. The reliability of the optimal solution is further proved by quantitative analysis about the inﬂuence of optimization parameters on the mean annual captured power. Compared with the common RSM and RBFNN methods, the prediction model established by the EBFNN method has a higher prediction accuracy. In the optimization process, the simulation calculation is replaced by a prediction model, which can effectively solve the problem of high simulation calculation cost.


Introduction
In recent years, the energy crisis and environmental problems have become increasingly serious. Vigorously developing and utilizing renewable energy has become a universal consensus of the international community [1,2]. Wave energy is a kind of clean and renewable energy. The theoretical reserve of the global wave energy resource is 3 TW, and further development can effectively alleviate the pressure of the global energy shortage [3]. Therefore, various wave energy converters (WECs) have been developed after decades of research and exploration. The WEC has gradually developed from the original concept of various and different types into five of the most typical device forms, namely, oscillating, raft, wave crossing, pendulum and oscillating water column [4][5][6]. The oscillating WEC has the characteristics of a simple structure, high power generation efficiency, and is not limited by water depth. It has received extensive attention, and many prototypes of oscillating WECs have been developed, such as the Archimedes Wave Swing, Power Buoy, Aqua Buoy, and Wavebob [7][8][9][10].
However, these prototypes have not been implemented into commercial operation, because of the low power generation efficiency and poor stability of WECs. In order to achieve the goal of high efficiency and stable output, researchers have carried out a lot of work, such as buoy optimization, power take-off (PTO) improvement, control strategy, and so on. Optimizing the buoy shape is one of the important ways to improve the wave power generation efficiency [11]. Therefore, many scholars have carried out related research. Some scholars have studied the hydrodynamic characteristics of the oscillating buoy, and analyzed the influence of buoy-shape parameters on energy capture efficiency to optimize the buoy shape. In order to efficiently extract the energy, Goggins et al. [12] established the optimal structural geometry and analyzed the optimal configuration for devices with a different radius. Shi et al. [13] numerically simulated the heave motion of conical bottom buoys with different masses and cross-sectional surfaces, and analyzed the effects of mass and cross-sectional surfaces on the energy capture. In order to analyze the effect of shape change on wave energy capture efficiency, Berenjkoob et al. [14] constructed 10 different geometric models under the same conditions. The results show that the proper shape was more efficient in capturing wave energy. Shadman et al. [15] took the cylindrical buoy as the research object and used the design of experiments (DOE) method to optimize the design parameters. Under the condition of mooring, Wu [16] obtained the optimal design parameters of the buoy by limiting the maximum displacement of the buoy and changing the draft ratio of the buoy radius. Some scholars also have used the heuristic algorithm to optimize buoy shape, such as genetic algorithm (GA), particle swarm optimization (PSO) et al. McCabe et al. [17] used a B-spline surface to parameterize the buoy geometry, and automatically optimized shapes based on GA. After 50,436 simulation calculations, a better design scheme was obtained. Garcia-Teruel et al. [18] further analyzed the adaptation of different algorithms such as GA and PSO, and used the non-dominated sorting genetic algorithm (NSGA-II) algorithm for multi-objective optimization design. Neshat et al. [19] optimized buoy shape and PTO based on bio-inspired swarm-evolutionary optimization algorithms. In general, the optimal buoy is mainly attained by analyzing the response relationship between the buoy parameters and the captured power, or by searching the design space through the heuristic optimization algorithm. Whether analyzing the response relationship or searching in the design space, enough sampling points are required. The hydrodynamic parameters of each sampling point are obtained by simulation to calculate the buoy captured power. Obviously, the hydrodynamic calculations take the majority of time during the buoy optimization. The significant time cost is considered a challenge for engineering design.
Therefore, reducing the simulation calculation times can effectively improve the optimization efficiency and reduce the optimization time cost. The prediction model is one of the best choices because it requires only a small number of simulation calculations. Prediction models have been used to solve the WEC optimization problems in engineering design and have been proven to be useful. Ji et al. [20] applied the response surface method (RSM) to build a hydrodynamic coefficients prediction model, which was used to calculate the captured power. Koh et al. [21] used the RSM prediction model to replace the simulation calculation, which determined the WEC geometric parameters. Colby et al. [22] developed a neural network function approximator and combined it with an evolutionary algorithm to optimize the WEC ballast configuration. The power output was maximized under a specific wave climate. Avila et al. [23] have predicted the energy generated from different WECs in the Canary Islands by using the Monte Carlo method. Avila et al. [24] proposed two mathematical models to forecast the wave energy based on fuzzy inference systems (FIS) and artificial neural networks (ANN), respectively. Wang et al. [25] proposed a prediction model of WEC buoy energy absorption power based on a back-propagation neural network (BPNN). Liu et al. [26] trained a radial basis function neural network (RBFNN) prediction model, which took the combination of design parameters as input and the capture factor as output. The prediction model was used with the genetic algorithm to optimize the converters, and the high-performance converter was obtained. Commonly, the prediction model indicates the mapping relations between the input (e.g., design parameter combinations of the buoy) and the output (e.g., capture factors, capture power, etc.), which is established by means of mathematical models. Combining the prediction model with the optimization algorithm to find the optimal solution in the design space becomes an optimization method that can quickly obtain the optimization results. So the accuracy of the prediction model is directly related to the reliability of the optimization results.
In recent years, methods such as RSM, BPNN, and RBFNN have been applied to WECs optimization design, but there are some shortcomings. The RSM cannot guarantee that the response surface passes through all the sample points, and the prediction effect of the highly complex functional relationship is not as good as that of the RBFNN. The BPNN has a slow convergence speed and might fall into local optimal solutions. The RBFNN has strong fault tolerance, but the modeling time is much longer than the RSM. Compared with the above methods, the elliptical basis functions neural network (EBFNN) has the advantages of a fast training speed and high prediction accuracy. Therefore, in order to improve the optimization efficiency, this paper proposes a mean annual power prediction model based on EBFNN. Compared with RSM and RBFNN, this method can establish a more accurate prediction model. Combining the mean annual power prediction model with MIGA can effectively reduce the simulation calculation time and improve the optimization design efficiency that an oscillating buoy with high performance can be easily designed. In this study, the optimization method not only optimizes the oscillating buoy, but also has the potential to solve other types of WECs optimization problems.
The remainder of this paper is organized as follows. Section 2 introduces the optimization method, including the geometric model, mathematical models, and optimization procedures. Section 3 is based on the wave statistics of the ocean near Chengshantou, Weihai City, Shandong Province, China; the mean annual power of the sample points is calculated. The prediction models with optimized parameters as input and mean annual power as output are established, respectively, by using RSM, RBFNN and EBFNN. Section 4 introduces the optimization results and discussion. Section 5 presents the conclusions.

Geometric Model
The schematic diagram of the oscillating buoy WEC is shown in Figure 1. In this study, a generalized WEC is considered, which is a single buoy with an ideal linear damper as the PTO system. The buoy only does heave motion under wave excitation. In the case of slack mooring, the mooring force has little effect on the wave energy capture by the buoy [20,21]. Therefore, this paper ignores the mooring force acting on the buoy. The origin of the Cartesian coordinate system is at the calm waterline. The x-axis points to the direction of wave propagation. The z-axis passes vertically upward through the center of gravity of the buoy. The y-axis is perpendicular to the former according to the right-hand rule. power, etc.), which is established by means of mathematical models. Combining the prediction model with the optimization algorithm to find the optimal solution in the design space becomes an optimization method that can quickly obtain the optimization results. So the accuracy of the prediction model is directly related to the reliability of the optimization results.
In recent years, methods such as RSM, BPNN, and RBFNN have been applied to WECs optimization design, but there are some shortcomings. The RSM cannot guarantee that the response surface passes through all the sample points, and the prediction effect of the highly complex functional relationship is not as good as that of the RBFNN. The BPNN has a slow convergence speed and might fall into local optimal solutions. The RBFNN has strong fault tolerance, but the modeling time is much longer than the RSM. Compared with the above methods, the elliptical basis functions neural network (EBFNN) has the advantages of a fast training speed and high prediction accuracy. Therefore, in order to improve the optimization efficiency, this paper proposes a mean annual power prediction model based on EBFNN. Compared with RSM and RBFNN, this method can establish a more accurate prediction model. Combining the mean annual power prediction model with MIGA can effectively reduce the simulation calculation time and improve the optimization design efficiency that an oscillating buoy with high performance can be easily designed. In this study, the optimization method not only optimizes the oscillating buoy, but also has the potential to solve other types of WECs optimization problems.
The remainder of this paper is organized as follows. Section 2 introduces the optimization method, including the geometric model, mathematical models, and optimization procedures. Section 3 is based on the wave statistics of the ocean near Chengshantou, Weihai City, Shandong Province, China; the mean annual power of the sample points is calculated. The prediction models with optimized parameters as input and mean annual power as output are established, respectively, by using RSM, RBFNN and EBFNN. Section 4 introduces the optimization results and discussion. Section 5 presents the conclusions.

Geometric Model
The schematic diagram of the oscillating buoy WEC is shown in Figure 1. In this study, a generalized WEC is considered, which is a single buoy with an ideal linear damper as the PTO system. The buoy only does heave motion under wave excitation. In the case of slack mooring, the mooring force has little effect on the wave energy capture by the buoy [20,21]. Therefore, this paper ignores the mooring force acting on the buoy. The origin of the Cartesian coordinate system is at the calm waterline. The x-axis points to the direction of wave propagation. The z-axis passes vertically upward through the center of gravity of the buoy. The y-axis is perpendicular to the former according to the right-hand rule.

Mean Annual Power
The modified Pierson-Moskowitz (P-M) spectrum is used to describe the distribution of wave energy relative to each component wave [27]. The spectrum can be expressed as: where H s is the significant wave height, T av is the average wave period, ω is the wave angular frequency. Divide the wave spectrum into 200 parts, the frequency ω i is from 0.02 to 4 rad/s, then the frequency interval is 0.02 rad/s [28]. The energy per unit area of an irregular wave is given by: where ρ is the density of the fluid, g is the acceleration due to gravity, ∆ω i is the frequency interval, N is the number of component angular frequencies (200), and α(ω i ) are the component wave amplitudes. Hence, the component wave amplitudes can be expressed as: Each sea state is assigned a set of proprietary random number generator seeds, from which repeatable normal distribution phase sets can be generated. The wave surface function of the i-th regular wave can be expressed as: where ε i,j is the initial phase uniformly random distributed between 0 and 2 π, and is associated with angular frequency ω i . j is the number of sea states. The corresponding incident wave is then synthesized by the linear superposition of all frequency components, as: Assuming an irrotational, non-viscous, and incompressible flow in the fluid domain, the waves are of small amplitude. Potential flow theory and three-dimensional radiation/diffraction theory are applied in the analysis. The velocity potential of the flow motion can be expressed as: where t is the time, ϕ → X is the spatial velocity potential function, ω is the wave frequency.
where ϕ I is the incident wave potential, ϕ D is the diffracted wave potential, ϕ rj is the radiation wave potential. ϕ → X satisfies the Laplace equation that can be expressed as: The radiation force acting on the buoy can be expressed as follows: where n is the direction vector, A m is the added mass, B C is the radiation damping. The complex amplitude of the exciting force acting on the buoy can be expressed as: Based on the linear potential flow theory and Newton's second law, the heave response in the frequency domain is: where M is the mass of the buoy, A m is the added mass, B C is the radiation damping, C PTO is damping coefficient of the PTO system, K is the restoring force coefficient, Z(ω) is the buoy displacement, F ex (ω) is the excitation force.
In irregular waves, the buoy heave displacement can be expressed as: The instantaneous power can be expressed as: With a spectral frequency interval of 0.02 rad/s, the maximum duration of any nonrepeating time series derived from this spectrum is t = 2π/0.02 ≈ 314 s. Then, the mean power over time series can be expressed as: Considering the random nature of waves, the mean power value for each sea state is the average of ten values, which can be expressed as: The mean annual power is obtained by taking the weighted average of the mean power values for all the sea states.
where H(H s , T av ) is the probability that each sea state occurs within a year.

Response Surface Method
RSM has a wide range of applications in the field of optimization design, and also has successful experience in the optimization of the WEC buoy shape [20]. RSM uses the knowledge of mathematics and statistics to analyze the response of a series of sample points in the design variable space, and usually uses the least square method to establish a prediction function, which is generally a low-order polynomial. Commonly, the relationship between the response value y(x) and the design variable x can be expressed as: where y(x) is the response function of the input parameter x, ε is the observation error. The second-order RSM is the most commonly used and is expressed as: where k is the total number of design variables, and x i is an input design variable, α 0 , α i , α ii , and α ij are the unknown regression coefficients of the constant, linear, quadratic and interaction terms, respectively. The least square method is employed to estimate the unknown coefficients α, which is expressed as:

Radial Basis Functions Neural Network
RBFNN is a feedforward network model which consists of input layer, hidden layer and output layer, and its basic structure is shown in Figure 2. From the input layer to the hidden layer is a fixed non-linear transformation that maps the input vector directly into a new space. The mapping from the hidden layer space to the output layer space is linear, and the output layer implements a linear weighted combination in the new linear space, where the weight is the network adjustable parameter. The RBFNN does not need any mathematical assumptions and has the characteristics of a black box. It has a fast learning speed, a strong ability to predict complex nonlinear functions and excellent generalization ability.
Energies 2022, 15, x FOR PEER REVIEW 6 of 21 a prediction function, which is generally a low-order polynomial. Commonly, the relationship between the response value ( ) and the design variable can be expressed as: where �( ) is the response function of the input parameter , is the observation error. The second-order RSM is the most commonly used and is expressed as: where is the total number of design variables, and is an input design variable, 0 , , , and are the unknown regression coefficients of the constant, linear, quadratic and interaction terms, respectively. The least square method is employed to estimate the unknown coefficients α, which is expressed as:

Radial Basis Functions Neural Network
RBFNN is a feedforward network model which consists of input layer, hidden layer and output layer, and its basic structure is shown in Figure 2. From the input layer to the hidden layer is a fixed non-linear transformation that maps the input vector directly into a new space. The mapping from the hidden layer space to the output layer space is linear, and the output layer implements a linear weighted combination in the new linear space, where the weight is the network adjustable parameter. The RBFNN does not need any mathematical assumptions and has the characteristics of a black box. It has a fast learning speed, a strong ability to predict complex nonlinear functions and excellent generalization ability. For n dimensional input variables 1 , 2 , ⋯ , ∈ , the RBFNN can be expressed as: where is the number of neurons in the hidden layer, is the weight factor, ( ) is the p-th basis function of the hidden layer, and � − � is the euclidean distance between the test point and the sample point.
The Gaussian basis function is used as the radial basis function of the neural network, which is expressed as: For n dimensional input variables x 1 , x 2 , · · · , x n ∈ R, the RBFNN can be expressed as: where p is the number of neurons in the hidden layer, λ p is the weight factor, ϕ p (x) is the p-th basis function of the hidden layer, and x − x j is the euclidean distance between the test point and the sample point. The Gaussian basis function is used as the radial basis function of the neural network, which is expressed as:

Elliptical Basis Functions Neural Network
Kavuri et al. [29] have constructed the EBFNN based on the study of pattern recognition. The EBFNN uses hyperellipse to divide the design variable space, and the hidden layer neurons use the ellipse function as the basis function. The EBFNN is the same as the structure of the RBFNN. The EBFNN takes the Mahalanobis distance between the point to be measured and the sample point as the independent variable. For n dimensional input variables x 1 , x 2 , · · · , x n ∈ R, the basis function can be expressed as: where S is the covariance matrix. The EBFNN prediction model is constructed by linear superposition, and the input and output layers can be expressed as: where c ji , b ji , respectively, represent the center value and semi-axis length of the ellipse element j on the i-th dimension. β is the adjustment parameter for the slope of the sigmoid curve. The larger the value of β, the faster the output value of the elliptic basis function tends to zero.

Error Analysis
The complex correlation coefficient R 2 is used to evaluate the consistency between the estimated value of the prediction model and the real value. It can be expressed as: where n is the number of sample points, y i is the actual value of the response,ŷ i is the model estimated value, y is the average value of the actual value. When R 2 is closer to 1, the model accuracy is higher. It is generally considered that when the R 2 is above 0.9, the accuracy of the prediction model can meet the engineering needs [30].

Optimization Algorithms and Processes
The origins of GA can be traced back to the early 1960s. The algorithm is designed and proposed according to the evolution law of organisms in nature. The GA imitates the genetic reproduction mechanism in the process of biological evolution, encodes the individuals in the solution space of the optimization problem, and then carries out genetic operations on the encoded individual population (such as selection, crossover, variation, etc.), and finds the combination containing the optimal solution or better solution from the new species population through iteration. GA is explained in detail in the works of Eiben and Smith [31].
The MIGA is an improvement of the parallel distributed genetic algorithm by Hiroyasu [32]. It divides a large population into several subpopulations, which are called "islands" vividly, and uses traditional GA algorithms to perform subpopulation evolution on each island. The schematic diagram is shown in Figure 3. The genetic evolution among the subpopulations does not interfere with each other. Every certain number of generations, some individuals are randomly selected to transfer to other islands, which helps complete the exchange of individuals between populations, and increase the diversity of individuals. Compared with GA, it can effectively suppress the occurrence of precocious phenomenon, which is beneficial to obtain the global optimal result [33]. the subpopulations does not interfere with each other. Every certain number of generations, some individuals are randomly selected to transfer to other islands, which helps complete the exchange of individuals between populations, and increase the diversity of individuals. Compared with GA, it can effectively suppress the occurrence of precocious phenomenon, which is beneficial to obtain the global optimal result [33]. The optimization process is illustrated in Figure 4. Some sample points are randomly selected in the design space by the optimal Latin hypercube design (OLHD). Calculating the hydrodynamic parameters (i.e., the wave excitation force, the additional mass, the radiation damping), the mean annual power for each sample point is calculated. The prediction model is established by analyzing the input and output variables of the sample library. When the prediction model meets the accuracy requirements, it is loaded into the optimization process. Then, the simulation calculation is replaced by the prediction model in the next optimization process. After initialization, according to the results of the previous generation, the parameters of the next-generation population individuals are obtained by genetic evolution. The individual parameters are input into the prediction model to obtain the energy capture power of the individual. The above process is repeated until the maximum generations and the population will contain the best individual after the maximum generations. The optimization process is illustrated in Figure 4. Some sample points are randomly selected in the design space by the optimal Latin hypercube design (OLHD). Calculating the hydrodynamic parameters (i.e., the wave excitation force, the additional mass, the radiation damping), the mean annual power for each sample point is calculated. The prediction model is established by analyzing the input and output variables of the sample library. When the prediction model meets the accuracy requirements, it is loaded into the optimization process. Then, the simulation calculation is replaced by the prediction model in the next optimization process. After initialization, according to the results of the previous generation, the parameters of the next-generation population individuals are obtained by genetic evolution. The individual parameters are input into the prediction model to obtain the energy capture power of the individual. The above process is repeated until the maximum generations and the population will contain the best individual after the maximum generations.

Determining the Design Space
In this paper, the WEC will be deployed in the sea near Chengshantou, Weihai, Shandong, China (123° E, 37.5° N). The wave resource statistics data were obtained from 2011 to 2020 by ERA5 reanalysis data [34]. The probability of occurrences of significant wave height and average wave period is shown in Table 1. According to the statistics data, the

Determining the Design Space
In this paper, the WEC will be deployed in the sea near Chengshantou, Weihai, Shandong, China (123 • E, 37.5 • N). The wave resource statistics data were obtained from 2011 to 2020 by ERA5 reanalysis data [34]. The probability of occurrences of significant wave height and average wave period is shown in Table 1. According to the statistics data, the average wave period concentrates in 3-6 s. the significant wave height mostly is less than 2 m. The WEC captures wave energy through the movement of the buoy with the wave. Energy capture efficiency mainly depends on the shape of the buoy and the response of the PTO system [35,36]. Therefore, the radius R and the ratio of draft radius D and the damping coefficient C PTO are selected as the optimization parameters.
(1) The diameter of buoy is preferably 5-10% of the main wavelength [37]. Therefore, the value range of the buoy diameter can be expressed as follows: (2) In order to ensure the rationality of the buoy shape, the draft of the buoy is normalized, which is expressed as D = H/R. The value range of D can be expressed as follows: (3) The traditional design method for the PTO system damping is based on the spectral peak frequency resonance, which leads to low energy capturing efficiency [20]. Therefore, the global search method is used to find the optimal PTO system damping which matches the wave resources. The value range of C PTO can be expressed as follows: 50 KNs/m < C PTO < 300 KNs/m (30) According to the wave resources in the target sea area and considering the cost, the goal of this paper is to maximize the mean annual power per unit submerged surface area of the buoy [38]. Constraints and objective functions are as follows: 0.6 m ≤ 2R ≤ 6 m 0.5 ≤ D ≤ 1 50 KNs/m ≤ C PTO ≤ 300 KNs/m (31)

Hydrodynamic Calculation Verification
In order to verify the accuracy of simulation calculations and avoid the influence of the number of meshes on the simulation results, the numerical calculation results are compared with the results published in this paper by Mei [39]. The radius of the cylinder buoy is 4.05 m, the draft is 4.05 m. The wave frequency varies from 0.1 to 1.9 rad/s with an interval of 0.1 rad/s. k is the corresponding wave number. The water depth d is 15 m. The calculation meshes are 1310, 3945, 5825 and 8000, respectively. The calculated results are shown in Figure 5. When the total number of meshes are 1310, the response amplitude operator (RAO) peak calculation error is larger and the root mean square error (RMSE) is 0.021. When the number of meshes are 3945, 5825 and 8000, the RMSE is 0.017, 0.016 and 0.016, respectively. The calculation results tend to be consistent, and continuous refinement of the mesh has little effect on the simulation results. Therefore, in order to ensure high numerical simulation accuracy, the same mesh density is used. max � ⁄ . . 0.6 m ≤ 2 ≤ 6 m 0.5 ≤ ≤ 1 50 KNs m ⁄ ≤ ≤ 300 KNs m ⁄

Hydrodynamic Calculation Verification
In order to verify the accuracy of simulation calculations and avo the number of meshes on the simulation results, the numerical calculati pared with the results published in this paper by Mei [39]. The radius of is 4.05 m, the draft is 4.05 m. The wave frequency varies from 0.1 to interval of 0.1 rad/s. k is the corresponding wave number. The water de calculation meshes are 1310, 3945, 5825 and 8000, respectively. The cal shown in Figure 5. When the total number of meshes are 1310, the re operator (RAO) peak calculation error is larger and the root mean squa 0.021. When the number of meshes are 3945, 5825 and 8000, the RMSE 0.016, respectively. The calculation results tend to be consistent, and ment of the mesh has little effect on the simulation results. Therefore, high numerical simulation accuracy, the same mesh density is used.

Calculating the Sample Points
The OLHD makes all test points as evenly distributed in the design space as possible, and has very good space filling and balance. The sample library is randomly selected in the design space by using OLHD, which covers 63 sample points [40]. The sample points are shown in Table 2, and the scatter diagrams of these samples are shown in Figure 6.  According to the sample points parameters, the simulation models are establish wave excitation force, additional mass and radiation damping are calculated for ea ple point based on ANSYS-AQWA. According to Equation (17), the mean annual p each sample point was calculated, and the scatter diagrams are shown in Figure 7.

Training and Testing of the Prediction Model
The sample library is divided into two groups, respectively named and the validation set. The first 52 groups of sample points are used and the last 11 groups of sample points are used as the verification s models are trained based on RSM, RBFNN and EBFNN, respectively. W dius is constant, the mapping relationship between the PTO system d radius ratio and mean annual power is shown in Figure 8.

Training and Testing of the Prediction Model
The sample library is divided into two groups, respectively named as the training set and the validation set. The first 52 groups of sample points are used as the training set, and the last 11 groups of sample points are used as the verification set. The prediction models are trained based on RSM, RBFNN and EBFNN, respectively. When the buoy radius is constant, the mapping relationship between the PTO system damping, draft-to-radius ratio and mean annual power is shown in Figure 8. The verification set is used to test the prediction performance of the three prediction models. The predicted value is compared with the numerical simulation calculation, and the results are shown in Figure 9. It can be seen intuitively from Figure 9 that the predicted values of the three prediction models are in good agreement with the actual values, and each point is distributed around the line = . Among the three prediction models, the EBF prediction model has the smallest error, and the RSM prediction model has the largest error. In addition to the use of test set verification, other evaluation methods are used. By Equation (27), the complex correlation coefficient R 2 of the training set is calculated. The R 2 of RSM prediction model is 0.99269. The R 2 of BFNN prediction model is 0.99913. The R 2 of EBFNN prediction model is 0.99963. The test results show that the obtained model could be used to predict the mean annual power. The verification set is used to test the prediction performance of the three prediction models. The predicted value is compared with the numerical simulation calculation, and the results are shown in Figure 9. It can be seen intuitively from Figure 9 that the predicted values of the three prediction models are in good agreement with the actual values, and each point is distributed around the line y = x. Among the three prediction models, the EBF prediction model has the smallest error, and the RSM prediction model has the largest error. In addition to the use of test set verification, other evaluation methods are used. By Equation (27), the complex correlation coefficient R 2 of the training set is calculated. The R 2 of RSM prediction model is 0.99269. The R 2 of BFNN prediction model is 0.99913. The R 2 of EBFNN prediction model is 0.99963. The test results show that the obtained model could be used to predict the mean annual power.

Optimization Results and Discussion
Based on three prediction models, the MIGA is used to optimize the shape of the buoy. The population size is 60, the subpopulation size is 20, the number of subpopulations are 3, the number of generations are 100, the rate of crossover is 0.85, the rate of mutation is 0.01, the rate of migration is 0.01, and the interval generation of migration is 5. The optimization results are shown in Figure 10.

Optimization Results and Discussion
Based on three prediction models, the MIGA is used to optimize the shape of the buoy. The population size is 60, the subpopulation size is 20, the number of subpopulations are 3, the number of generations are 100, the rate of crossover is 0.85, the rate of mutation is 0.01, the rate of migration is 0.01, and the interval generation of migration is 5. The optimization results are shown in Figure 10.
During the optimization process, the mean annual power is calculated by the prediction model. According to the calculation results of the previous generation, the optimization variables are continuously updated to form new generations, and the iterations are repeated until the maximum number of generations is reached. After 100 generations of evolution, the optimization results are shown in Table 3.  During the optimization process, the mean annual power is calculated by the prediction model. According to the calculation results of the previous generation, the optimization variables are continuously updated to form new generations, and the iterations are repeated until the maximum number of generations is reached. After 100 generations of evolution, the optimization results are shown in Table 3. As can be seen from Table 3, the comparison between optimization results and simulation results shows that the RSM prediction model has the lowest accuracy between the predicted value and the actual value, with an error of 9.45%. The RBFNN prediction model has the better accuracy by comparing the predicted value with the actual value, with an error of 5.54%. The EBFNN prediction model has the highest accuracy by comparing the predicted value with the actual value, with an error of 0.13%. The accuracy As can be seen from Table 3, the comparison between optimization results and simulation results shows that the RSM prediction model has the lowest accuracy between the predicted value and the actual value, with an error of 9.45%. The RBFNN prediction model has the better accuracy by comparing the predicted value with the actual value, with an error of 5.54%. The EBFNN prediction model has the highest accuracy by comparing the predicted value with the actual value, with an error of 0.13%. The accuracy comparison of the three prediction models is the same as that of the corresponding complex correlation coefficient R 2 . The accuracy of the prediction model established by the EBFNN method is the best, while the accuracy of the prediction model established by the RSM method is the worst.
In order to further verify the accuracy of the optimization results, the three design parameters are qualitatively analyzed respectively. The influence of these three parameters on P/S is analyzed by some simulation calculations. Figure 11 illustrates the response of P/S with buoy-shape parameters under optimal C PTO . In Figure 11a, when the D is same, the P/S first increases and then decreases with the increase of the buoy radius. The P/S reaches peak when the radius is around 1.33 m. In Figure 11b, when the R is the same, the P/S decreases with the increase of the D. According to the change trend in Figure 11, it can be judged that the optimal parameters of the buoy are about when the R is 1.33 m and the draft is 0.665 m. ofP/S with buoy-shape parameters under optimal CPTO. In Figure 11a, when the D is same, theP/S first increases and then decreases with the increase of the buoy radius. TheP/S reaches peak when the radius is around 1.33 m. In Figure 11b, when the R is the same, the P/S decreases with the increase of the D. According to the change trend in Figure 11, it can be judged that the optimal parameters of the buoy are about when the R is 1.33 m and the draft is 0.665 m. When the buoy is determined, the effect of CPTO on the capture power of WEC is analyzed, as shown in Figure 12. In Figure 12a, when D is the same, the optimal CPTO increases with the increase of R, and the smaller D, the faster the optimal CPTO increases with the radius. In Figure 12b, when the R is the same, the optimal CPTO decreases with the increase of the D, and the larger the radius, the faster the optimal CPTO decreases with the increase of D. Since the wave energy is mainly distributed on the sea surface, wave force decays rapidly with increasing water depth. In Figure 12c, theP/S first increases and then decreases with the increase of CPTO. In the design space, each set of shape parameters has an optimal CPTO. According to the change trend in Figure 12, it can be judged that the optimal CPTO is 50 KNs/m. When the buoy is determined, the effect of C PTO on the capture power of WEC is analyzed, as shown in Figure 12. In Figure 12a, when D is the same, the optimal C PTO increases with the increase of R, and the smaller D, the faster the optimal C PTO increases with the radius. In Figure 12b, when the R is the same, the optimal C PTO decreases with the increase of the D, and the larger the radius, the faster the optimal C PTO decreases with the increase of D. Since the wave energy is mainly distributed on the sea surface, wave force decays rapidly with increasing water depth. In Figure 12c, the P/S first increases and then decreases with the increase of C PTO . In the design space, each set of shape parameters has an optimal C PTO . According to the change trend in Figure 12, it can be judged that the optimal C PTO is 50 KNs/m.
According to the comparison between the quantitative analysis results and the optimization results, it is found that the error of the optimization results of the RSM prediction model and the RBFNN prediction model is relatively large, and the error of the optimization results of the EBFNN prediction model is small, which further proves the reliability of the optimization results of the EBFNN prediction model.
The comparison of optimization results, numerical simulation results and qualitative analysis results shows that the EBFNN prediction model can reflect the function characteristics of input parameters and mean annual power most effectively. Combining the EBFNN prediction model with the MIGA to perform global optimization in the design space, it is easy to obtain high-performance buoy design parameters. When the prediction model method is adopted, only the sample points need to be calculated. After the prediction model meets the accuracy, the prediction model is used to replace the simulation calculation in the optimization process, and a total of 63 simulation calculations are required. When the direct optimization design method is adopted, it is necessary to perform simulation calculations on the newly generated population individuals, and then the genetic variation generates the next-generation population. This iteration is repeated until the optimization generations are reached, which requires a total of 6300 simulation calculations. Obviously, the number of repeated modeling and simulation calculations required by the prediction model method is only 1% of the direct optimization design, and correspondingly, 99% of the original workload can be saved. The comparison between the prediction model optimization design method and the direct optimization design shows that the prediction model replacing the simulation calculation method can effectively reduce the number of repeated simulation calculations in the optimization process, effectively shortening the optimization design cycle, and improving the optimization design efficiency. According to the comparison between the quantitative analysis results and the optimization results, it is found that the error of the optimization results of the RSM prediction model and the RBFNN prediction model is relatively large, and the error of the optimization results of the EBFNN prediction model is small, which further proves the reliability of the optimization results of the EBFNN prediction model.
The comparison of optimization results, numerical simulation results and qualitative analysis results shows that the EBFNN prediction model can reflect the function characteristics of input parameters and mean annual power most effectively. Combining the EB-FNN prediction model with the MIGA to perform global optimization in the design space, it is easy to obtain high-performance buoy design parameters. When the prediction model method is adopted, only the sample points need to be calculated. After the prediction model meets the accuracy, the prediction model is used to replace the simulation calculation in the optimization process, and a total of 63 simulation calculations are required. When the direct optimization design method is adopted, it is necessary to perform simulation calculations on the newly generated population individuals, and then the genetic variation generates the next-generation population. This iteration is repeated until the optimization generations are reached, which requires a total of 6300 simulation calculations. Obviously, the number of repeated modeling and simulation calculations required by the prediction model method is only 1% of the direct optimization design, and correspondingly, 99% of the original workload can be saved. The comparison between the prediction model optimization design method and the direct optimization design shows that the prediction model replacing the simulation calculation method can effectively reduce the

Conclusions
In order to further improve the capture efficiency of WEC, this paper focuses on the optimization design of the buoy shape. In view of the shortcomings of traditional optimization methods such as long-time consumption and low optimization efficiency, this paper discusses the feasibility of using prediction models to replace simulation calculations. The method uses a small number of experimental samples to build a prediction model. In the optimization process, it replaces the actual model, which greatly reduces the evaluation and computation of the objective function and constraint function. Aiming at the shortcomings of RSM and RBFNN to establish the prediction model, the EBFNN is proposed to train the prediction model. The prediction models obtained by the three methods are combined with MIGA to conduct global optimization in the design space. The obtained optimization results are compared with the simulation calculation results and qualitative analysis results, and the following conclusions can be obtained: (1) The comparison shows that the prediction model established by the RSM method has the worst accuracy, the prediction model trained by the RBFNN method has better accuracy, and the prediction model trained by the EBFNN method has the best accuracy. The mean annual power prediction model trained by the EBFNN method can more accurately reflect the mapping relationship between the input and output. According to the shape parameters of the buoy, the mean annual power can be accurately predicted. (2) Taking the wave statistics data of the Chengshantou sea area near Weihai City, Shandong Province, China, as an example, the method of combining MIGA and the mean annual prediction model is adopted to obtain a high-performance design scheme, which provides a reference for engineering design. In the optimization process, the mean annual power prediction model replaces the simulation calculation, which can reduce a lot of workload (i.e., repeated modeling, simulation, and calculation). Compared with optimization design based on simulation results, this method can save considerable time and cost, effectively shorten the optimization design cycle, and improve the optimization efficiency. This optimization method can also be extended to the optimal design of other sea areas or types of WECs. In the future, we will continue to explore WEC array optimization methods. In order to promote the development of WEC commercialization, WEC array optimization methods will be investigated in the future. (3) The three optimization parameters have a significant impact on the energy capture performance. When the buoy shape is determined, the P/S first increases and then decreases with the increase of damping, and there is an optimal C PTO . The optimal C PTO is significantly affected by buoy radius and draft, which is positively correlated with the buoy radius and negatively correlated with the buoy draft. When C PTO is the optimal damping, the P/S increases first and then decreases with the increase of radius and the P/S decreases with the increase of draft.

Conflicts of Interest:
The authors declare no conflict of interest.