Multi-Mode Wave Energy Converter Design Optimisation Using an Improved Moth Flame Optimisation Algorithm

Abstract

Ocean renewable wave power is one of the more encouraging inexhaustible energy sources, with the potential to be exploited for nearly 337 GW worldwide. However, compared with other sources of renewables, wave energy technologies have not been fully developed, and the produced energy price is not as competitive as that of wind or solar renewable technologies. In order to commercialise ocean wave technologies, a wide range of optimisation methodologies have been proposed in the last decade. However, evaluations and comparisons of the performance of state-of-the-art bio-inspired optimisation algorithms have not been contemplated for wave energy converters’ optimisation. In this work, we conduct a comprehensive investigation, evaluation and comparison of the optimisation of the geometry, tether angles and power take-off (PTO) settings of a wave energy converter (WEC) using bio-inspired swarm-evolutionary optimisation algorithms based on a sample wave regime at a site in the Mediterranean Sea, in the west of Sicily, Italy. An improved version of a recent optimisation algorithm, called the Moth–Flame Optimiser (MFO), is also proposed for this application area. The results demonstrated that the proposed MFO can outperform other optimisation methods in maximising the total power harnessed from a WEC.

1. Introduction

Global energy demand has increased tremendously during the last decade and is expected to surge by $4.6\%$ in 2021 [1] and rise to 45% by 2040 [2]. Renewable energy sources are considered alternatives and are one of the fastest-growing technologies globally. Renewables are set to produce around $30\%$ of generated electricity globally in 2021. The most significant contribution is from solar and wind renewable energy. The electricity produced from wind power and solar PV is predicted to increase by 275 terawatt-hours and 145 terawatt-hours in 2021, respectively. More than $50\%$ of the world’s green power has been generated by China, followed by the United States, the European Union (EU) and India [1].

Ocean wave energy is currently considered the cleanest, safest and most reliable source of renewable energy [3]. Furthermore, the density of wave energy is significantly higher than that of solar and wind energy [4]. Nevertheless, this new kind of renewable energy needs more development and has not fully matured; its commercial perception is still shallow. The main reason for this is that the cost associated with generating energy using sea waves is much higher than that of other kinds of renewable and clean power, such as wind energy [5].

To handle these challenges, numerous studies have been conducted in order to optimise different aspects of wave energy converters (WECs), such as geometry parameters [6,7,8,9,10,11], power take-off settings (PTO) [12,13] and the arrangement of WECs [14,15,16,17,18]. The optimisation of the shape of WECs can be significant for wave energy conversion models, such as oscillating water columns, where the optimum parameters achieved an immersion depth of $0.45$ of the water depth and a chamber diameter of $0.92$ of the water depth [19]), and over-topping designs [20], where the optimum length to opening ratio should be 2.5–3 m.

Most of the studies on optimising the geometry parameters of WECs have been developed to examine a few predefined shapes, because the computational cost of exploring the multi-parameter search space and assessing all feasible designs is high.

This paper develops an improved Moth–Flame Optimisation (IMFO) method in order to maximise the harnessed power output of a WEC. To evaluate the quality of the best-found configurations using the IMFO, we developed a comparison framework including five robust and popular optimisation algorithms including Grey Wolf Optimiser [21] (GWO), Whale Optimisation Algorithm (WOA) [22], standard Moth–Flame Optimisation Algorithm (MFO) [23], Particle Swarm Optimisation (PSO) [24] and Covariance matrix adaptation evolution strategy (CMA-ES) [25]. The experimental results confirm that the proposed MFO method performed the best compared with the other meta-heuristics applied in this work, in terms of higher produced power.

The structure of the article is as follows: Section 2 presents a review of the previous developments of the shape and PTO parameters. Next, Section 3 describes the details of the WEC design and the formulas that are used in order to simulate the incorporated power. Then, the optimisation problem is explained, and Section 5 outlines the proposed and applied meta-heuristic approaches. The optimisation performances are introduced and studied in Section 6. Subsequently, Section 7 presents the conclusions of this research.

2. Related Works

One of the initial efforts to optimise the shape of WECs was proposed by Vantorre et al. [26]; an assortment of geometries for a heaving point absorber were evaluated and compared based on the Belgian coastal area. These shapes incorporated a hemisphere and two different conical geometries. In [26], the most reliable power efficiency was achieved with a cylindrical cone at 90${}^{\circ }$.

Goggins and Finnegan [27] investigated different heights and radii for a point absorber with a vertical-shape cylinder under wave conditions off the west coast of Ireland. Goggins and Finnegan [27] proposed that a trimmed cylinder with a hemisphere, with a full draft at an aspect ratio of 2.5, can produce the maximum power. However, Hager et al. [28] suggested that a concave-shape buoy can extract more power from ocean waves than other shapes, such as a flat or convex model, after evaluating several asymmetrical converter designs.

Another recent published study recommends using a bi-cubic B-spline wave converter shape [29], which can produce more energy than other conventional WEC models. Recently, Neshat et al. [30] proposed a bi-level optimisation algorithm that consists of Grey Wolf Optimiser [21] (GWO) and a self-adaptive differential evolution with ensemble sinusoidal parameter adaptation called LSHADE-EpSin [31]). The authors suggested a large vertical cylinder converter with a radius and height of 14.51 m and 30 m, respectively.

3. Modelling

3.1. Wave Energy Converter

This study focuses on optimising a generic multi-mode wave energy converter that absorbs power from multiple degrees of freedom due to the arrangement of mooring lines. The buoy is positively buoyant, and three inclined tethers should always be under tension to keep it fully submerged. The power is generated using a rotary generator connected to each tether to convert the orbital motion of the WEC to electricity. A simplified sketch of the system is shown in Figure 1. The buoy sits 2 m below the sea water level, has a cylindrical shape, and its geometry is defined by a radius of a and a height of H. The angles ${\alpha }_t$ and ${\alpha }_{ap}$ define how the tethers are inclined and attached to the buoy hull.

3.2. Site Location and Wave Resource

The wave characteristics dataset was collected using a point absorber located at 37.96${}^{\circ }$ N, 12.04${}^{\circ }$ E, Aegadian sea (Italy), near Marettimo island, over the 32 years between 1979 and 2010. The climate of Marettimo island is Mediterranean; the number of cold months is less than four and the average temperature is below 10 ${}^{\circ }$C. Generally, summer is warm (temperature >22 ${}^{\circ }$C), dry, and calm in terms of wind speed. Most of the storms happen in the fall and winter seasons. The maximum wind speed is approximately 28 Km/h, blowing from North to North West.

The potential location for this WEC deployment is demonstrated in Figure 2 and the detailed parameters are specified in

Table 1. Parameters of the Marettimo test site.

Location	37.96${}^{\circ }$ N, 12.04${}^{\circ }$ E
Type of data	Real sea measurement [32]
Water depth	10 m
Mean wave power density	6.38 kW/m
WXSD resource class [33]	Class 1

.

Figure 3a shows the wave scatter diagram of the deployment site. The wave climate is not energetic, with a wave power density of 6.38 kW/m. The dominant wave period is approximately 6 s, and the dominant wave height is less than 1 m.

As optimisation requires a large number of objective function evaluations, it is possible to replace the wave climate in Figure 3a with a subset of sea states in order to evaluate the performance of the WEC. Following the procedure explained in [34], a cluster analysis was applied to extract 10 representative sea states using the k-means clustering method, as shown in Figure 3b. The probability of each centroid was then recalculated to make sure that the weighted wave power flux from each centroid corresponded to the sum of all of the probabilities within each cluster.

3.3. Equations of Motion

For this work, it was assumed that waves are mono-directional and propagate along the x axis. The wave structure interaction was modelled using linear wave theory assuming inviscid, irrotational and incompressible fluid. The only nonlinear effects included in the WEC model were viscous forces. It was assumed that the PTO machinery behaves as a linear spring-damper system, where spring and damper parameters that are identical for all three generators are optimised for each sea state of interest.

The dynamics of the WEC were modelled using a spectral-domain approach, where all nonlinear forces were replaced by equivalent linear matrices that were calculated iteratively, applying the statistical linearisation technique [35]. The WEC motion is described by the vector $\widehat{\mathbf{u}}(\omega )$, which represents the complex amplitude of the buoy velocity in all six degrees of freedom:

$$\left[i\omega \left(\mathbf{M}+{\mathbf{A}}_{rad}(\omega )\right)+\left({\mathbf{B}}_{rad}(\omega )+{\mathbf{D}}_{pto}+{\mathbf{B}}_{visc}\right)-i\frac{{\mathbf{K}}_{pto}}{\omega }\right]\widehat{\mathbf{u}}(\omega )={\widehat{\mathbf{F}}}_{exc}(\omega ),$$

(1)

where $\mathbf{M}$ is the diagonal mass matrix, ${\mathbf{A}}_{rad}(\omega )$ is the frequency-dependent matrix of added mass coefficients, ${\mathbf{B}}_{rad}(\omega )$ is the frequency-dependent radiation damping matrix, ${\mathbf{K}}_{pto}$ and ${\mathbf{D}}_{pto}$ are linearised stiffness and damping matrices that represent the load from the power take-off machinery (refer to [36]), ${\mathbf{B}}_{visc}$ is the equivalent damping term that corresponds to the viscous drag force, and ${\widehat{\mathbf{F}}}_{exc}(\omega )=A_w(\omega ){\widehat{\mathbf{f}}}_{exc}(\omega )$ is the excitation force vector ($A_w$ is the wave amplitude, and ${\widehat{\mathbf{f}}}_{exc}(\omega )$ is the vector of excitation force coefficients).

In this work, the hydrodynamic parameters ${\mathbf{A}}_{rad}(\omega )$, ${\mathbf{B}}_{rad}(\omega )$ and ${\widehat{\mathbf{f}}}_{exc}(\omega )$ were calculated using a semi-analytical model [37,38]. The mass matrix $\mathbf{M}$ is defined by the buoy geometry, assuming that the mass of the buoy is equal to half of the displaced water mass. The PTO matrices ${\mathbf{K}}_{pto}$ and ${\mathbf{D}}_{pto}$ are defined by the arrangement of tethers and PTO control parameters. The value of the equivalent damping term ${\mathbf{B}}_{visc}$ depends on the WEC motion and should be calculated iteratively according to [35]:

$${\mathbf{B}}_{visc}=-〈\frac{\partial {\mathbf{F}}_{visc}(t)}{\partial \mathbf{u}(t)}〉,$$

(2)

where $〈·〉$ denotes the mathematical expectation, and the viscous force in time domain has the form:

$${\mathbf{F}}_{visc}(t)=-\frac{1}{2}{\rho }_w{\mathbf{C}}_d{\mathbf{A}}_d(|\left|\mathbf{u}(t)\right||\odot \mathbf{u}(t)).$$

(3)

${\rho }_w$ is the water density, ${\mathbf{C}}_d$ is the diagonal matrix of drag coefficients, ${\mathbf{A}}_d$ is the diagonal matrix of the cross-section areas of the buoy perpendicular to the direction of motion, and ⊙ denotes the Hadamard product (element-wise multiplication).

The following iterative procedure should be followed to estimate ${\mathbf{B}}_{visc}$ and the WEC response in any irregular sea state of interest:

Step 1.

Select the sea state of interest and calculate the incident wave spectrum $S_{\eta }(\omega )$ for the given $H_s$ and $T_p$. A Pierson–Moskowitz wave spectrum was used in this study [39].

Step 2.

Calculate the power spectral density (PSD) matrix of the excitation force:

$${\mathbf{S}}_{\mathbf{F}}(\omega )=S_{\eta }(\omega ){\widehat{\mathbf{f}}}_{exc}(\omega ){\widehat{\mathbf{f}}}_{exc}^{\ast }(\omega ),$$

(4)

where ${()}^{\ast }$ denotes the conjugate transpose of a vector/matrix.

Step 3.

Evaluate the WEC transfer function:

$$\mathbf{H}(\omega )={\left[i\omega \left(\mathbf{M}+{\mathbf{A}}_{rad}(\omega )\right)+\left({\mathbf{B}}_{rad}(\omega )+{\mathbf{D}}_{pto}+{\mathbf{B}}_{visc}\right)-i\frac{{\mathbf{K}}_{pto}}{\omega }\right]}^{-1},$$

(5)

where ${\mathbf{B}}_{eq}={\mathbf{0}}_{6\times 6}$ is the first iteration.

Step 4.

Calculate the PSD matrix of the buoy velocity:

$${\mathbf{S}}_{\mathbf{u}}(\omega )=\mathbf{H}(\omega ){\mathbf{S}}_{\mathbf{F}}(\omega ){\mathbf{H}}^{\ast }(\omega ).$$

(6)

Step 5.

Compute the covariance matrix of the buoy velocity:

$${\sigma }_{\mathbf{u}}^2=cov[\mathbf{u},\mathbf{u}]={\int }_0^{\infty }{\mathbf{S}}_{\mathbf{u}}(\omega )d\omega .$$

(7)

Step 6.

Approximate the viscous damping matrix ${\mathbf{B}}_{visc}$ as [35]:

$$\begin{array}{c}\hfill {\mathbf{B}}_{visc}=-〈\frac{\partial {\mathbf{F}}_{visc}}{\partial \mathbf{u}}〉=\frac{1}{2}\sqrt{\frac{8}{\pi }}{\rho }_w{\mathbf{C}}_d{\mathbf{A}}_d{\sigma }_{\mathbf{u}}^2.\end{array}$$

(8)

Step 7.

Check the convergence:

$$|{\mathbf{B}}_{visc}\left[n\right]-{\mathbf{B}}_{visc}[n-1]|<\delta ,$$

(9)

where n is the iteration number and $\delta$ is set to 0.01. Return to Step 3 if Step 7 is not satisfied.

Equation (1) describes the motion of the WEC in the Cartesian coordinate system. However, the power is generated based on the tether elongations and velocities that can be calculated using the inverse kinematic Jacobian as $\dot{\mathbf{q}}={\mathbf{J}}^{-1}\mathbf{u}$ ($\dot{\mathbf{q}}$ is the vector containing the rate of change of each tether length).

If the damping coefficient of each PTO unit is equal to $D_{pto}$, the average wave power converted by the WEC in the irregular wave is defined by $H_s$ and $T_p$ is [35]:

$$\overline{P}(H_s,T_p)=D_{pto}\sum _{k=1}^3{\sigma }_{{\dot{q}}_k}^2(H_s,T_p),$$

(10)

where k corresponds to the PTO number, and the variance matrix of the tether velocity $\dot{\mathbf{q}}$ is calculated as:

$${\sigma }_{\dot{\mathbf{q}}}^2={\int }_0^{\infty }{\mathbf{J}}^{-1}{\mathbf{S}}_{\mathbf{u}}(\omega ){\mathbf{J}}^{-T}d\omega .$$

(11)

For a location with a wave climate described by the matrix $\mathbf{O}(H_s,T_p)$, the annual average power production is:

$$P_{AAP}=\sum _{H_s}\sum _{T_p}O(H_s,T_p)·\overline{P}(H_s,T_p).$$

(12)

4. Optimisation Setup

With regard to defining the fitness function, we consider that the most significant optimisation variables of the applied cylinder are as follows: (1) The cylinder radius (a); (2) the aspect ratio, which is the dimension of the cylinder height over the radius ($H/a$); (3) tether attachment angle (${\alpha }_{ap}$); (4) tether inclination angle (${\alpha }_t$); (5) a vector of parameters of damping PTO; and (6) a vector of stiffness coefficients’ PTO values, which are shown as follows:

$${\mathbf{d}}_{pto}={[D_{pto}^{(1)},D_{pto}^{(2)},\dots ,D_{pto}^{(N_d)}]}^T,\mathrm{and}{\mathbf{k}}_{pto}={[K_{pto}^{(1)},K_{pto}^{(2)},\dots ,K_{pto}^{(N_k)}]}^T.$$

(13)

The variable numbers of both $k_{pto}$ and $d_{pto}$ vectors are the same at $N_k=N_d=10$, and also the total length of optimisation designs is 24. These decision variables must be optimised in what follows:

$$\mathbf{z}=[a,H,{\alpha }_t,{\alpha }_{ap},{\mathbf{k}}_{pto}\in {\mathbb{R}}^{N_k\times 1},{\mathbf{d}}_{pto}\in {\mathbb{R}}^{N_d\times 1}].$$

(14)

The objective function we considered in this study is:

$$f_O=\underset{\mathbf{z}}{argmax}{P}_{AAP}(\mathbf{z}),\mathrm{subject}\mathrm{to}:\mathbf{z}\in [{\mathbf{z}}_{min},{\mathbf{z}}_{max}].$$

(15)

Equation (15) shows the generated power output, on average, of the WEC in one year, and this power should be maximised.

Table 2. The feasible range of the cylinder hyper-parameters settings.

Parameter	Unit	Min	Max	Length
a	m	1	10	1
H	m	1	10	1
$(H/a)$		0.4	2	1
${\alpha }_t$	deg	10	80	1
${\alpha }_{ap}$	deg	10	80	1
$K_{pto}$	N/m	${10}^3$	${10}^8$	10
$D_{pto}$	N/(m/s)	${10}^3$	${10}^8$	10

presents the feasible ranges of the applied optimisation variables.

5. Bio-Inspired Optimisation Algorithms

To optimise the performance of a three-tether converter that is fully-submerged, in this paper we mainly considered a broad and well-known optimisation strategy called ‘all-at-once’. In this strategy, we performed the optimisers on all decision variables simultaneously. The mentioned variables are as follows: the radius a, aspect ratio ($H/a$), both tether inclinations (${\alpha }_t$), attachment angle (${\alpha }_{ap}$), both series of the spring stiffness (${\mathbf{k}}_{pto}$) and damping coefficients (${\mathbf{b}}_{\mathbf{pto}}$ ) of the power take-off features.

According to the previous study [30], by analysing the characteristics of the search space of this optimisation problem, we can see that these combined parameters formed a multi-modal, constrained and large-scale search space that is sensitive to PTO parameters alteration. These challenges provide the main motivation for applying the meta-heuristic optimisation methods, which are similar to genetic, evolutionary and swarm optimisation methods.

The meta-heuristics performed and compared in this work incorporate Covariance matrix adaptation evolution strategy (CMA-ES) [40], Grey Wolf Optimiser (GWO) [21], Particle Swarm Optimisation (PSO) [41], Whale Optimisation Algorithm (WOA) [22], standard Moth–Flame Optimisation Algorithm (MFO) [23] and the developed MFO.

Table 3. Details of the control parameters of the applied optimisation algorithms. All algorithms are limited to the same number as the fitness function evaluation.

Methods	Settings
CMA-ES [40]	$\lambda =25$ with the default settings recommended in reference [40];
PSO [41]	with $\lambda =25$, $c_1=1.5$, $c_2=2$, $\omega =1$ (reduced by a damping ratio $w_f=0.99$ exponentially);
GWO [21]	$\lambda$= 25, $\alpha =2$ (linearly decreased to zero)
WOA [22]	$\lambda =25$, ${\alpha }_1=2$ (declines linearly from 2 to 0), ${\alpha }_2=-1$ (linearly reduced from −1 to −2), $\beta =1$
MFO [23]	$\lambda =25$, $\alpha =-1$ (linearly drops from −1 to −2), $\beta =1$, $t=(\alpha -1)\times rand+1$;
IMFO	$\lambda =25$, the same MFO’s control parameters

lists the applied optimisation methods and their control parameters, where $\lambda$ shows the initial population size; $c1$ and $c_2$ are cognitive and social acceleration coefficients, respectively; and $\omega$ signifies the inertia weight. In GWO, $\alpha$ shows the relative coefficient for balancing the exploration and exploitation of the search process. Moreover, both ${\alpha }_1$ and ${\alpha }_2$ explain the convergence coefficients of WOA. $\beta$ denotes a constant in order to determine the shape of the logarithmic spiral.

5.1. Improved Moth–Flame Optimisation (IMFO)

The Moth–Flame Optimisation (MFO) algorithm [23] is a state-of-the-art, nature-inspired meta-heuristic method that demonstrates a considerable performance in optimising various numerical [42,43] and real engineering [44,45] optimisation problems. The principal motivation of the Moth–Flame Optimisation algorithm is the navigating strategy of the moth, called ’transverse orientation’. Indeed, some groups of insects, like the moth family, explore their environment based on a navigation model that involves flying a route at a specific angle based on the distance between them and the light source, such as natural or man-made lights. If the natural light sources are followed by the moths, due to the long distance this navigation system assists with keeping to a straight flight path. However, for the man-made lights, the situation is different and leads to the moths flying in a spiral pattern. Equation (16) shows the mathematical model of this logarithmic spiral movement in one dimension, as follows:

$$S=D{e}^{bt}cos(2\pi t)+F,$$

(16)

where F denotes the situation of each flame in the search space; M signifies the location of the moth, $D=|M-F|$ is the absolute interval between each moth and the existing flame. Moreover, $t\in (-1,1)$ equals a time step. Based on Equation (16), a route from the moth to a flame is formed while t linearly reduces between $-1.5$ and 1; b is a constant and restricts the speed of approach to a flame.

In the MFO, particles cross in the d dimensional search space and their situations are depicted as feasible or infeasible solutions. Hence, a population of moths and their situations can be identified as follows:

$$\begin{array}{c}\hfill M=\left[\begin{array}{ccccc}{M}_{11}& M_{12}& \cdots & \cdots & M_{1d}\\ M_{21}& M_{22}& \cdots & \cdots & M_{2d}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ M_{n1}& M_{n2}& \cdots & \cdots & M_{nd},\end{array}\right]\end{array}$$

(17)

where n denotes the size of the population, and the fitness value of the solutions are recorded as a vector in the following:

$$\begin{array}{c}\hfill f(M)=\left[\begin{array}{c}f(M_1)\\ f(M_2)\\ ⋮\\ f(M_n).\end{array}\right]\end{array}$$

(18)

To show the flames, the same matrix that can be seen in Equation (19) is used, and each flame shows the best-found solution so far.

$$\begin{array}{c}\hfill F=\left[\begin{array}{ccccc}{F}_{11}& F_{12}& \cdots & \cdots & F_{1d}\\ F_{21}& F_{22}& \cdots & \cdots & F_{2d}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ F_{n1}& F_{n2}& \cdots & \cdots & F_{nd}\end{array}\right]\begin{array}{c}\hfill ⟶f(F)=\left[\begin{array}{c}f(F_1)\\ f(F_2)\\ ⋮\\ f(F_n).\end{array}\right]\end{array}\end{array}$$

(19)

Thus, the area around the flames in the search space is explored by a formal moth in order to search for better solutions. The flames’ situations will be updated when a better configuration is observed. Furthermore, Function I initialises a population of random solutions and then determines the fitness values of the solutions. The following approach is used in this way:

$$\begin{array}{c}{M}_{ij}=(u{b}_i-l{b}_i)\ast rand()+l{b}_i⟶i\in [1,n],j\in [1,d],\hfill \end{array}$$

(20)

where both $ub$ and $lb$ are vectors of decision variables’ upper and lower boundary constraints as follows:

$$\begin{array}{cc}\hfill ub=& [u{b}_1,u{b}_2,\dots ,u{b}_d],\hfill \end{array}\begin{array}{cc}\hfill lb=& [l{b}_1,l{b}_2,\dots ,l{b}_d].\hfill \end{array}$$

(21)

The second step is the P function, which signifies the principal component of the MFO algorithm. This function transfers solutions throughout the search space and also renews the status of the recorded flames during the optimisation process. The state of the $i$-th solution in the population with regard to the $j$-th flame will be updated as follows:

$$\begin{array}{cc}\hfill M_i=& S(M_i,F_j)=D_{ij}{e}^{bt}cos(2\pi t)+F_j⟶\begin{array}{cc}\hfill D_{ij}=& |M_i-F_j|.\hfill \end{array}\hfill \end{array}$$

(22)

In the MFO, the exploration procedure strengthens the candidate solutions to improve the flying direction in general. As a result, the exploration strategy supports the optimisation process to evade trapping in local optima. Furthermore, the exploitation procedure runs while the new location extends inside the search space between the moth and the flame. The next step of the MFO is to sort the candidate solutions with regard to their quality from the best performers to the worst performers. Later, the moth’s location will be updated again, depending on the corresponding flame. Furthermore, in order to obtain a good balance between the exploitation and exploration processes, an adaptive strategy was recommended. This formula reduces the number of flames through iteration, as follows:

$$\begin{array}{c}\hfill F_n=round\left(N-iter\ast \frac{N-1}{Ma{x}_{iter}}\right),\end{array}$$

(23)

where $F_n$ and N are the number of flames and the size of the population, respectively.

5.2. Diversification Strategy

In the MFO, in order to avoid local optima and premature convergence, and to improve the exploration ability, we proposed a diversification strategy for the MFO. In this mechanism, we alllow the MFO to generate a small percentage of new random solutions in each iteration, based on a probability threshold ($\mu =0.1$). These newly generated solutions keep the diversity between the solutions in the population and lead to strengthening the exploration ability until the last iterations. This is mainly because, in the last iterations of the standard MFO, the diversity rate of the population has been reduced and the search algorithm is not able to explore properly. Equation (24) shows the formulation of the diversification strategy.

$$M_{ij}=(u{b}_i-l{b}_i)\ast rand()+l{b}_i⟶ifrand<\mu .$$

(24)

As in the standard MFO, the number of both moths and flames dwindles, and the exploration ability gradually weakens. The proposed improved MFO was developed using a simple and effective diversification strategy. Therefore, keeping the population diversity extends the algorithm’s global search capability and prevents premature convergence. However, it is expected that the convergence rate of IMFO decreases due to the generation and addition of a few random design solutions in the population. This improved algorithm will be applied to case studies in the next section.

6. Results and Discussions

In this paper, in order to evaluate the searchability of the proposed IMFO algorithm, we developed a comparative framework consisting of a well-known evolutionary algorithm—called CMA-ES [40]—and four popular swarm intelligence methods—Particle Swarm Optimisation (PSO) [41], Grey Wolf Optimiser (GWO) [21], Whale Optimisation Algorithm (WOA) [22] and standard Moth–Flame Optimisation Algorithm (MFO) [23].

The applied control parameters for all optimisation methods are listed in Table 3, as suggested by their developers. In this way, we developed a fair comparison between the applied optimisation methods.

In the first step of this comparison for maximising the absorbed power output of the cylinder, each optimisation method was run ten times and the minimum, maximum and median of the best-found solutions are reported as a box-and-whisker diagram (Figure 4). In Figure 4, we can see that the IMFO could propose the best settings for the cylinder in terms of harnessing sea wave energy. Furthermore, the performance of the standard MFO is competitive and, in the following, the GWO performed better than the WOA, PSO and CMA-ES.

The second comparison metric for the optimisation methods’ performance is the convergence rate. In this study, we set the maximum evaluation number at $5\times {10}^3$, which was suggested by a prior study [30]. According to the results represented in Figure 5, it is obvious that the standard WOA could converge substantially faster than the other meta-heuristics. However, WOA could not be released from the local optima and was faced with a premature convergence issue. Interestingly, CMA-ES shows the same convergence behaviour as WOA, but the performance of the WOA method is better than CMA-ES in terms of best-found configurations. In the following convergence trajectories, the original MFO algorithm outperforms both WOA and CMA-ES after consuming only $25\%$ of the whole computational budget, and maintains this superiority over the others until near the last iterations. However, surprisingly, the improved MFO surpasses all optimisation algorithms and can propose better solutions on average. This achievement comes from the high exploration ability, using the proposed diversification strategy.

Table 4. Optimisation performance of IMFO and the five other optimisation methods based on the best-found configuration per each experiment.

	CMA-ES	PSO	GWO	WOA	MFO	IMFO
Mean	8.0621 × 104	7.9192 × 104	8.1757 × 104	8.1134 × 104	8.3430 × 104	8.4448 × 104
Min	8.0003 × 104	7.1264 × 104	7.9016 × 104	7.9956 × 104	8.1406 × 104	8.1534 × 104
Max	8.1032 × 104	8.5783 × 104	8.5985 × 104	8.4410 × 104	8.5933 × 104	8.8274 × 104
STD	3.6621 × 102	5.1481 × 103	2.3205 × 103	1.3638 × 103	1.7451 × 103	2.5527 × 103

presents the average, minimum, maximum and standard deviation of the best-found solutions proposed by the improved MFO and the five other optimisation algorithms. It is noted that IMFO found the best optimisation settings for the cylinder in order to harness the maximum power from the sea wave. Meanwhile, both the GWO and MFO heuristics worked well as the second and third best performing wave power optimisers at finding solutions with the highest absorbed power at 8.5985 × 10⁴ and 8.5933 × 10⁴, respectively.

Table 5. The best-found geometric and PTO settings of the cylinder.

	CMA-ES	PSO	GWO	WOA	MFO	IMFO
a [m]	10	10	9.97	9.99	10	10
H [m]	10	10	10	10	10	10
${\alpha }_t$ [deg]	80	36	10	79.99	80	79.85
${\alpha }_{ap}$ [deg]	79	33	80	76.98	29.85	11
${\overline{)\sum }}_{i=1}^{N_K}$$K_{pto}$$(\times {10}^7)$	8.85	4.9415	4.4482	7.5518	5.8939	5.0907
${\overline{)\sum }}_{i=1}^{N_B}$$D_{pto}$$(\times {10}^7)$	9.575	6.3568	5.5574	9.4748	7.6913	6.7046
Power (Watt)	8.10 × 104	8.58 × 104	8.60 × 104	8.44 × 104	8.59 × 104	8.83 × 104

shows the list of the best-found design parameters, including the cylinder radius (a), height (H), both tether attachment angles (${\alpha }_{ap}$), the tether inclination angle, (${\alpha }_t$) and, finally, the damping-spring PTO values. The first and foremost observation from Table 5 is that both the radius and height values of the best-found solutions are converged to the upper bound, and it is expected that a bigger cylinder can produce more power. However, most of the time the costs of developing, installing and maintaining the large wave converters restricted the developers. Furthermore, it can be seen from Table 5 that the optimised values of the tether angles are within a wide range between 10${}^{\circ }$ and 80${}^{\circ }$. This evidence demonstrates that there is a straightforward method for setting these design parameters. This is mainly because there are strong dependencies between design parameters (tether angles and PTO parameters) and the hydrodynamic model that dominates the power extraction (surge, pitch or heave).

In order to provide deeper insight into the optimisation achievements, we compared the performance of the proposed optimisation method with other meta-heuristics, which can be seen in

Table 6. The percentage improvement of the WEC absorbed power using IMFO compared with the other optimisers applied in this study.

	CMA-ES	PSO	GWO	WOA	MFO
Mean	4.75%	6.64%	3.29%	4.08%	1.22%
Worst	1.91%	14.41%	3.19%	1.97%	0.16%
Best	8.94%	2.90%	2.66%	4.58%	2.72 %

. In terms of the average best-found WEC’s settings, IMFO can perform substantially better than PSO by $6.64\%$, and the percentage improvements of IMFO can be considerable compared with CMA-ES, WOA and GWO at $4.75\%$, $4.08\%$ and $3.29\%$, respectively. Furthermore, according to the best-found optimal WEC’s configurations, IMFO is able to outperform the standard MFO by $2.72\%$.

The sizes of the cylinder’s radius and height have a substantial impact on the power conversion from the converter. Figure 6 and Figure 7 show the process of the cylinder size optimisation using (a) the standard MFO and (b) the proposed MFO. Each line represents the search process of the optimisation algorithm, which is repeated ten times. As there is a semi-direct relationship between the radius and height size and the produced power output, both methods could converge rapidly and have approximately the same performance.

The inclination angle of the tethers ($\alpha$) of the developed cylinder can modify the absorbed power of the ocean wave. As the tether angle cannot be adjusted after installing the converter, optimising this parameter plays an important role in the average performance of the cylinder during its lifetime. As a result, for more powerful sea states, it is suggested that the tether angles be larger because more power will be absorbed from the surge than the heave of the incident wave. The comparison of both optimisation methods—(a) MFO and (b) IMFO—in order to explore the best configurations of the inclination and attachment angles can be seen in Figure 8 and Figure 9, respectively.

Furthermore, from Figure 8, we can observe that for both MFO and IMFO methods, after exploring the search space of the tether’s inclination angle ($\alpha$), most of the experiments are converged to the upper bound in order to absorb the maximum level of wave energy. It can be seen in Figure 8a that the standard MFO rapidly converged to a specific range (perhaps a local optimum) and was not able to keep the exploration ability. However, the proposed MFO kept its exploration searchability even among the last iterations (Figure 8b).

Figure 9 demonstrates the trajectory of the (a) MFO and (b) IMFO methods to optimise the attachment angles of the three tethers. According to the optimisation results, MFO proposed a specific range of angles from 50${}^{\circ }$ to 80${}^{\circ }$ for ten independent runs. However, we can see that the IMFO algorithm suggested a different range of attachment angles during the optimisation process, with a higher exploration ability compared to the MFO method.

In Figure 10, we can see clearly that the best-found Kpto parameters using the MFO method are converged to less than $5\times {10}^7$. Nevertheless, IMFO proposed the values of Kpto in a wide range of search domains. As a result, these found Kpto configurations can produce more power compared with the MFO performance. Additionally, Figure 10b shows that the search space of the PTO parameters is multi-modal with highly dynamic behavior.

In order to clarify and provide an observable comparative landscape for both MFO and IMFO, Figure 11 is presented. To study the exploration and exploitation abilities of IMFO carefully and in great detail, the last ten iterations of the proposed method are plotted. We can notice clearly that the best-found configurations for both Kpto and Dpto using IMFO are distributed throughout the search space; however, the proposed settings of MFO are concentrated in a narrow range.

Moreover, we can see that it is more challenging to find the optimal values for the tether angles than for the radius and height values. Finally, it seems that the most challenging part of the optimisation is related to finding the appropriate values for the PTO parameters.

7. Conclusions

In this study, an improved Moth–Flame Optimisation Algorithm is proposed for the purpose of enhancing the power output extracted from ocean waves by a fully-submerged, cylindrical wave energy converter with three tethers. The applied wave climate in this optimisation work was a Mediterranean sea site that is located in Sicily, Italy. As optimising the shape features, PTO parameters and tether angles simultaneously is a complex and time-consuming optimisation problem, we proposed a robust and reliable development of the Moth–Flame Optimisation (MFO) Algorithm using a diversification strategy. With regard to evaluating the developed method’s effectiveness, we take into account five state-of-the-art swarm and evolutionary optimisation algorithms for comparison. According to the findings, the improved MFO method is able to outperform recent meta-heuristic methods in terms of the absorbed power output of the best-found WEC configurations. In order to develop the current method in terms of dealing with multi-directional wave conditions, in the future we will focus on the adaptive and self-adaptive meta-heuristics to achieve a good search balance that will result in a fast and effective wave energy converter optimiser. Future plans also include considering whether the speed of the convergence rate can be beneficial in terms of saving on computational costs.