A Comparative Study of Metaheuristic Algorithms for Wave Energy Converter Power Take-Off Optimisation: A Case Study for Eastern Australia

: One of the most encouraging sorts of renewable energy is ocean wave energy. In spite of a large number of investigations in this ﬁeld during the last decade, wave energy technologies are recognised as neither mature nor broadly commercialised compared to other renewable energy technologies. In this paper, we develop and optimise Power Take-off (PTO) conﬁgurations of a well-known wave energy converter (WEC) called a point absorber. This WEC is a fully submerged buoy with three tethers, which was proposed and developed by Carnegie Clean Energy Company in Australia. Optimising the WEC’s PTO parameters is a challenging engineering problem due to the high dimensionality and complexity of the search space. This research compares the performance of ﬁve state-of-the-art metaheuristics (including Covariance Matrix Adaptation Evolution Strategy, Gray Wolf optimiser, Harris Hawks optimisation, and Grasshopper Optimisation Algorithm) based on the real wave scenario in Sydney sea state. The experimental achievements show that the Multiverse optimisation (MVO) algorithm performs better than the other metaheuristics applied in this work.


Introduction
The recent uprising demand for more energy resources is imposing huge pressure on energy supply industries [1]. The most sustainable solution to offset the current depletion of fossil energy resources is using renewable energy technologies. Marine wave energy has been considered as a reliable source for coastal demands since 1799 [2]. Thus, its development has occurred vastly-compared to the other resources-due to the following benefits: (i) It is another source of sustainable energy contributing to the mix of energy resources that leads to more diversity and attraction for coastal cities and suppliers [3]. (ii) Wave energy could be exploited offshore and does not need any land, resulting in a lower cost of in-land site selection and reducing unfavourable visual impact [4]. (iii) By considering the best layout and offshore site location, the permanent generation of energy will be accessible (compared to solar energy, for instance, which is time-dependent) [5]. By taking all the mentioned reasons into account, recent studies have focused on enhancing the performance of wave energy technologies by developing optimisation-based solutions in the section of geometric design, PTO parameters, and layout. The second division proved to have a crucial role in energy transmission and exploitation in electricity networks [6]. While the energy sector is developing different technologies in this area, some challenges still need to be dealt with, such as design limitations [7], operation and maintenance difficulties in offshore area, power take-off efficiency under heavy currents, etc. [8]. Hence, it is not far-fetched to assume that optimal solutions are always required to achieve better performance. A summary of recent studies has been classified in Table 1.  [30] In the following, this article is classified into five sections. We describe a summary of the hydrodynamic interaction model and wave resource, and present the computation of the produced power in Section 2. Section 3 shows the optimisation strategies and manifests the details of the Multiverse optimisation (MVO) algorithm. In Section 4, we discuss the power take-off optimisation results in terms of efficiency and produced power output. Finally, in Section 5, we summarise the major findings of this comparative study, and give our suggestion for future studies.

Model Setup
In this section, we introduce the wave scenario in the selected site on the eastern coasts of Australia. Then, we describe the geometry of the converter and analyse the governing equation of the motion of the converter's body. Finally, by considering all implied forces into account, the power output of the WEC is numerically simulated.

Wave Resource
In this study, a site on the eastern coast of Australia (Sydney) is investigated according to the 2016 real wave dataset [31]. One of the main ways to represent wave climate is the wave rose, which is a graphical demonstration of significant wave height over its various directions [16].
For instance, the dominant sea state specifies Sydney with a peak period (T p ) around 9 s and a Significant Wave Height (H s ) of 2 m; the Pierson-Moskowitz wave spectrum is helpful for identifying the probability distribution of the wave frequency ( f w ). In fact, the zone of plausible incident wave angle (β) at Sydney is 180. Due to the diversity in directional distributions and probable cumulative energy, wave regimes may have differences with one another. Therefore, Sydney has been chosen to assess the performance of optimisation methods. In this research, a modified spectrum called Bretschnider is used, which needs two parameters (T p and H s ) to model irregular waves. It is worth considering the fetch length and the storm's location together along with wind speed and wind direction as they can change the spectrum on a different level.
In Equation (1), f p and H m0 are peak wave period and significant wave height, respectively.

Wave Energy Converter Modelling
A fully submerged, symmetric and spherical buoy is used as the WEC in this study. To keep the buoy in location with specific tolerances against both severe storm load and normal load conditions, the buoy has equipped with three tethers, each one is connected to a mooring point on the seabed with an angle of 55 degrees (α). Details of the buoy parameters can be seen in Table 2 and Figure 1 demonstrates the buoy perfectly [16].  The Frequency-domain simulation provides a faster approximation, which is applied in this research. This simulation enjoys a significant advantage when it comes to running many model evaluations in designing different scenarios. In the frequency-domain, a linear PTO model is applied. We require modelling of only translational degrees of freedom sway, surge, and heave. The buoy can harness energy from all modes of tethers [16,32].
The buoy's three tethers have power take-off mechanisms, which can be modelled as an electric generator or a hydraulic circuit. A PTO system is connected to each tether and the loading force from each of them is assumed to have a linear spring and a damping effect relating to the tether lengthening and rate of change of its lengthening [33].
The fully submerged buoy bears forces that play important roles in numerical modelling. It includes exerted force from the PTO system and excitation force. The excitation forces are felt when the body stands in front of incoming waves, while radiation forces are sensed when the body moves in a steady-state. To calculate the excitation force, incident and scattered waves in hydrodynamic pressure need to be summed up. On the other side, Radiation forces arise on account of body motions, and the radiation resistance, which may be presumed as wave damping force, also has a relation with the average energy exchanging among the WEC's body and sea. This force informs us about the amount of harnessed energy from incoming waves. Another force is added-mass, which is the inertia of water entrained with body motion. It should be considered that the optimal device parameters contingent upon the wave frequency would change relentlessly [4]. The mooring system is designed according to the linear wave theory, considering small motion amplitude compared with tethers' length. The following S i indicates the spatial arrangement of tethers.
where in the reference coordinate frame O xyz , r is the vector of buoy's position and R is the radiation matrix. The position vector of the anchor point of tether i is related to G, which is the centre of mass, and d i is the position vector of the anchor point of tether i on the seafloor in O xyz . The length of tether can be calculated as and the change of the length is modelled as ∆l i = l i − l 0 [34]. The below equation is used for modelling the PTO force.
where D pto and K pto are the damping and stiffness matrices, respectively. C pto is a counteracting force of hydrostatic force, which is obtained from the following equation: where m w is the mass of displaced water, m b is the mass of the buoy, g is the gravitational acceleration constant, and α is assumed to be 55 degrees as proposed in [34]. The buoy velocities can be mapped in a Cartesian coordinate frame considering the change rate of the tether length, which is provided by the inverse kinematic Jacobian in [34]. Other forces including viscous forces and the end stop forces are neglected and kinematics are linearised in this research. Among the discussed forces, the most important one is the excitation force, which is calculated by Equation (6) [16].
where F exc is the frequency-dependent vector of excitation forces, M is the mass matrix, and I is the identity matrix. We use a constant 3 value because there are three degrees of freedom.Ẍ is a vector of body acceleration in the surge, heave, and sway directions.
Added-mass and radiation forces are defined through A and B matrices, respectively. There are three tethers attached to the buoy modelled as an oscillating spring. There are 50 individual K pto and D pto parameters; therefore, in each tether, there are 100 individuals, and in total it will be 300 for a buoy. Hydrodynamic reactions can be calculated according to the semianalytical solution described in [35].
The power output of the buoy is calculated in an irregular wave frequency domain according to Equation (7) [16]. In this study, the applied wave energy converter simulator was developed and published by Sergiienko on Matlab-R2019b [36].

Optimisation
In this section, we describe the general formulation of the problem regarding the optimisation viewpoint, then, we briefly introduce optimisation algorithms. Finally, the MVO optimiser has been reviewed as a novel, population-based optimiser in the application of WEC's PTO assessment.

Optimisation Formulation
The applied formulation of the optimisation problem in order to maximise the produced power output of the WEC can be seen in the following: where Power represents the annual average power generated for given PTO settings of a WEC where the location is fixed: Power Take

Optimisation Algorithms
In this study, we use five metaheuristic swarm optimisation algorithms to achieve the optimal values of power take-off parameters in terms of either damping or stiffness variables. These algorithms have been chosen based on recent reviews of their application in WECs' power take-off performance assessments [37,38]. The collected algorithms include the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [39], Gray Wolf optimiser [40], Harris Hawks optimisation [41], and Grasshopper Optimisation Algorithm [42]. More specifically, we go through the application details of another recently released metaheuristic optimisation algorithm called Multiverse optimiser. All algorithms have been used for 10,000 total evaluations, carried out by 25 search agents at 400 iterations. We replicate each experiment ten times to find the minimum, maximum, mean, and standard deviation for each optimiser in this problem. Table 3 presents the settings of the five employed algorithms. Table 3. The settings of optimisation approaches. The maximum evaluation number is 10 5 .

Algorithms Settings
CMA-ES [39] with the default settings and nPop = 25; MVO [43] with nPop = 25, (decreased with a damping ratio w f = 0.99 exponentially); GWO [40] with nPop = 25 and the de , α = 2 (linearly decreased to zero) GOA [42] with nPop = 25, and the default settings HHO [41] with nPop = 25, and the default settings As the application of CMA-ES, GOA, HHO, and GWO algorithms have been studied individually in the previous works [17,44], we will go through details of the state-of-the-art Multiverse optimiser for this problem as follows.

Multiverse Optimiser (MVO)
In a recent released study [43], a promising and novel population-based optimisation algorithm is introduced entitled "Multiverse optimiser (MVO)" that divides the search mode toward the principal two steps: exploration and exploitation. In the MVO method, the main inspiration is related to the theories of black holes and, inversely, white holes to perform the exploitation and exploration processes, respectively. Furthermore, wormholes are effective in order to develop the search abilities as a local search strategy. In spite of what preceded, some distinct concepts are implemented in the MVO including a universe that represents a solution, an individual in the universe corresponds to each variable of the solution, and the inflation rate accords with the value of the fitness function. The main strategies of the MVO algorithm can be seen in the following:

1.
The larger values of the inflation rate are regarded as performing a white hole with a higher probability rate than that of a black hole appearing.

2.
A universe with a large inflation degree leads to transfer candidates. However, a low inflation rate results in receiving candidates within black and white spirits.

3.
Total elements in the universe will stochastically turn around the best candidate within the wormhole, notwithstanding the statistical values of the inflation rate.
In each generation of the MVO approach, a roulette-wheel technique is applied to choose a white hole from the total universes based on the inflation rate. The main objective of this mechanism is to improve candidates transfer in various universes and also to develop exploration abilities. Thus, the MVO scheme assumes that, as follows: where n is the number of solutions in the population and d is the length of decision variables. Therefore, each solution can be defined as follows: In Equation (10), r 1 shows a random number in the interval of [0, 1], x j i symbolises the jth variable of ith universe, where x j k shows the jth variable of kth universe chosen by a roulette-wheel selection algorithm. Meanwhile, U i explains the ith universe and N I(Ui) denotes the rate of normalised inflation of the ith universe, so updating the new location of the candidates is as follows: where r 2 , r 3 , and r 4 are generated random numbers within [0, 1], x j i mentions the jth variable of ith universe, and x j intimates the jth variable of the best universe developed so far. Both TDR and WEP are applied as the coefficients, lb j determines the lower bound of jth variable and ub j denotes the upper bound of jth variable. In the MVO algorithm, there are two significant coefficients including wormhole existence probability (WEP) and travelling distance rate (TDR). Both coefficients are developed as adaptive formulas as follows: where min and max are constant values at 0.2 and 1, respectively. In the original MVO, iter is the current iteration and Max iter is the maximum number of iterations.
Finally, ρ implies the exploitation accuracy beyond the iterations, which is chosen as 6 in the original MVO paper. Where the ρ value progresses higher, it is likely to receive more advanced and precise exploitation. Moreover, the Quick-sort algorithm is proposed to sort the universe after each iteration. This process facilitates finding optimal answers in this optimisation problem in terms of power take-off variables.

Experimental Optimisation Results
In this section, the experimental results of optimising parameters are presented. The output records of the model revealed the convergence curve of the converter's power output to the maximum value. Then, the minimum, maximum, mean value, and standard deviation of the experiments have been calculated. Furthermore, the overall trend of PTO coefficient changes has been recorded, followed by statistical analysis of achieved results based on the PTO's coefficient configuration. As mentioned in Section 3.2, each experiment comprises 10 5 evaluation numbers. The numerical analysis has been carried out using five different metaheuristic promising optimisation algorithms, including MVO, CMA-ES, GWO, GOA, and HHO. To begin with, Figure 2 tracks the convergence of WEC's power output to its maximum amount using five above mentioned optimisation algorithms over 10 5 iterations.
By looking at Figure 2, we can perceive the competence of each algorithm to optimise this parametric problem. As it can be easily seen, the Multiverse optimiser succeeds other algorithms to extract higher power. According to the search process attributes of the HHO algorithm, it seems to be stuck in local optima, thereby losing its chance to perform an exhaustive search throughout the whole search space. Second to the MVO, the CMA-ES approach achieves the best results, followed by the other three algorithms by a non-negligible distance. In brief, both MVO and CMA-ES were observed to be the most productive methods of optimising in this problem. In order to examine the performance of studied methods to find the optimal number of PTO coefficients, we analyse the trajectory of damping and spring coefficients for each of the three PTO systems in three different frequency samples-including 1, 25, and 50-in Figure 3. By looking at Figure 3, we can perceive the competence of each algorithm to optimise this parametric problem. As it can be easily seen, the Multi-verse optimizer succeeds other algorithms to extract higher power. According to the search process attributes of the HHO algorithm, it seems to be stuck in local optima, thereby losing its chance to perform an exhaustive search throughout the whole search space. Second to the MVO, the CMA-ES approach achieves the best results, followed by the other three algorithms by a non-negligible distance. In brief, both MVO and CMA-ES were observed to be the most productive methods of optimising in this problem.
In order to examine the performance of studied methods to find the optimal amount of PTO coefficients, we analyse the trajectory of damping and spring coefficient for each three PTO systems in three different frequency samples, including 1, 25, and 50 in Figure 4.   Looking at Figure 3, we can observe that the parameter coefficients for each PTO stiffness of the damping system change nonmonotonically over the course of iterations as its frequency changes. Moreover, high chaotic fluctuations are observed at the beginning of the search process in all frequencies, specifying the exploration phase of the searching process.
Having a closer look at variables, most damping coefficients tend to reach higher values than stiffness coefficients in lower frequencies, like 1 and 25. On the other hand, this trend is vice versa in the 50th frequency. Having said that, the maximum changes from the initial amount are experienced by the stiffness coefficient. This leads us to the fact that the objective function is affected by this variable more than those related to the damping part of the PTO system. Figure 4 presents the box-and-whisker plot for the power output of the best solution per run for all search metaheuristic algorithms for the Sydney wave scenario on the eastern Australian coasts. Considering Figure 4, a comparative-statistical analysis can be found. The best mean power output of the WEC is observed using the MVO method. It is followed by the CMA-ES, GWO, GOA, and HHO methods, respectively. Furthermore, the dispersion of the achieved data over the mean value is high in the GWO and HHO methods compared to CMA-ES and GOA. Overall, as the power output is just under 300 kW considering MVO, this method is found to be the most compatible algorithm for this optimisation problem.
Finally, Table 4 summarises the mean, minimum, maximum, and standard deviation of each algorithm's results over 10 runs.  Table 4 reveals that the maximum power output of the converter can be produced by letting the MVO algorithm use optimum PTO coefficients. This maximum result is followed by the results of CMA-ES, GWO, HHO, and GOA, respectively. Regarding the mean and minimum amount of power outputs, these trends are found to be repeated. As for the standard deviation of the results, we can sort the algorithms in order as HHO, MVO, GWO, GOA, and CMA-ES. Additionally, the convergence of the results is more precise in HHO and MVO while, on the other hand, results from GWO, GOA, and CMA-ES are much more spread out.

Conclusions
In this paper, we developed a comparative optimisation framework in order to evaluate and optimise the Power Take-off (PTO) settings of a fully-submerged three-tether wave energy converter (WEC). This WEC was introduced and developed by Carnegie Clean Energy Company in Australia.
Due to a large-scale dimension of PTO parameters and complex interactions among three tethers, optimising the WEC's PTO parameters is a challenging engineering problem. In this way, we developed and compared the effectiveness of five state-of-the-art metaheuristic algorithms on the Sydney wave scenario, Australia. The preliminary results confirm that the Multiverse optimisation (MVO) algorithm can considerably outperform the other four optimisation algorithms applied in this paper. Second to the MVO, the highest mean value and maximum value of the WEC's power output were achieved using CMA-ES, GWO, GOA, and HHO, respectively. As another concluding remark of this study, it worth mentioning that the stiffness coefficients proved to have more effect on objective functions than damping coefficient in the studied converter's power take-off system.
In future work, we will develop a more comprehensive comparative framework to analyse the performance of a wide range of Genetic, Evolutionary, and Swarm optimisation algorithms in order to optimise the PTO settings. Furthermore, developing and assessing the impact of hybridisation, alternating, and cooperative optimisation techniques can be considered.  Institutional Review Board Statement: Ethical review and approval were waived for this study, due to studies not involving humans or animals.