A New Bi-Level Optimisation Framework for Optimising a Multi-Mode Wave Energy Converter Design: A Case Study for the Marettimo Island, Mediterranean Sea

Abstract

To advance commercialisation of ocean wave energy and for the technology to become competitive with other sources of renewable energy, the cost of wave energy harvesting should be significantly reduced. The Mediterranean Sea is a region with a relatively low wave energy potential, but due to the absence of extreme waves, can be considered at the initial stage of the prototype development as a proof of concept. In this study, we focus on the optimisation of a multi-mode wave energy converter inspired by the CETO system to be tested in the west of Sicily, Italy. We develop a computationally efficient spectral-domain model that fully captures the nonlinear dynamics of a wave energy converter (WEC). We consider two different objective functions for the purpose of optimising a WEC: (1) maximise the annual average power output (with no concern for WEC cost), and (2) minimise the levelised cost of energy (LCoE). We develop a new bi-level optimisation framework to simultaneously optimise the WEC geometry, tether angles and power take-off (PTO) parameters. In the upper-level of this bi-level process, all WEC parameters are optimised using a state-of-the-art self-adaptive differential evolution method as a global optimisation technique. At the lower-level, we apply a local downhill search method to optimise the geometry and tether angles settings in two independent steps. We evaluate and compare the performance of the new bi-level optimisation framework with seven well-known evolutionary and swarm optimisation methods using the same computational budget. The simulation results demonstrate that the bi-level method converges faster than other methods to a better configuration in terms of both absorbed power and the levelised cost of energy. The optimisation results confirm that if we focus on minimising the produced energy cost at the given location, the best-found WEC dimension is that of a small WEC with a radius of 5 m and height of 2 m.

1. Introduction

Renewable energy is the fastest-growing new energy source globally. As an example, in the United States, the growth rate of this technology increased by 100% between 2000 and 2018 [1]. On a global scale, renewable energy technologies produced 26.2% of the global electricity demand in 2018, and this is expected to climb to 45% by 2040 [1]. A large number of investigations have been applied in order to optimise various characteristics of renewable energy systems such as dealing with the uncertainty in renewable energy accessibility, support decision-making in the built environment [2] and the appropriation of energy storage operations for dampening the chaotic problems [3]. Among the different renewable energy sources, ocean wave energy is the cleanest, safest, most reliable and predictable source of renewable energy [4] with a power density significantly higher than that of solar and wind [5]. However, wave energy technology is not fully developed, and their commercial penetration is still shallow. This is because the costs involved in producing energy using ocean waves are currently much higher than those for other renewables [6]. Therefore, in the last decade, a large number of investigations have been carried out to optimise wave energy converter (WEC) design and dimensions [7, 8, 9, 10, 11, 12] power generation settings (PTO) [13,14], and the position of WECs in a wave farm [15, 16, 17, 18, 19].

The wave energy resource around the globe has been divided into six major classes depending on the wave energy potential, directional and spectral characteristics, and extreme waves [20]. However, it has been noted [20] that while wave energy developers mainly target wave climates with the highest energy content (class 5 and 6), other resource classes can provide additional benefits to the technology development. For example, the Mediterranean Sea due to its enclosed nature has low wave power availability [21, 22, 23] and belongs to the resource class 1 but the absence of extreme wave heights makes this region attractive for the initial prototype testing.

Shape optimisation is important for all types of wave energy conversion systems, including oscillating water columns [24]), and over-topping designs [25]. The majority of efforts, to date, have been restricted to analysing a few specific shapes. The main reason for this is that the computational demands of searching and evaluating all feasible designs are high. Vantorre et al. [26] evaluated and compared the performance of a set of geometries for a heaving point absorber in a Belgian coastal area. These included a hemisphere and some conical geometries. The authors proposed that the best power efficiency was related to a cylindrical extension with a 90° cone. Later work by Goggins and Finnegan [27] contemplated a vertical cylinder of various heights and radii under wave conditions off the west coast of Ireland. They found that the most substantial significant heave velocity response was that of a trimmed cylinder with a hemisphere joined to its foundation, with a whole draft to the aspect ratio of 2.5. In other recent publications, a wide range of asymmetrical buoy designs has been proposed, including a concave buoy face which is better able to absorb power than a flat or convex model [28]. Another recommendation of a surface described by bi-cubic B-spline [29] outperforms conventional WEC models. However, in these studies, the main objective was to maximise the harnessed power of the WEC, and the authors did not consider the design, installation and maintenance costs of these asymmetric converters.

Other work has taken into account the trade-offs between absorbed power and the cost of building and deploying the WECs. These analyses have considered the cost-efficiency or levelised cost of energy (LCoE) [30]. This metric is one of the most reliable indices for the evaluation of energy investments. Recently, Piscopo et al. [31] combined an LCoE minimisation with a power take-off (PTO) control optimisation based on point-absorber dimensions in five Mediterranean Sea sites. This refined earlier work, optimising LCoE through optimisation of both WEC geometry and PTO settings [32,33].

In this work, we consider a single fully submerged, three-tether, cylindrical wave energy converter. This WEC is under development by Carnegie Clean Energy Limited, Australia. Two initial attempts [12,34] were performed to investigate the impact of different geometries and PTO parameters on power efficiency and the LCoE. However, in these prior works, only some predefined geometries were studied, and the results showed that in the cylinder-shaped WEC, an optimal tethers angle depends on the ratio between the buoy height and radius. However, optimisation procedures were not adequately outlined [34]. In another study [12], the performance of a few conventional optimisation methods was investigated in order to maximise the absorbed power and minimise the LCoE.

This paper improves upon previous research by expanding the findings of [12] to include another two state-of-the-art meta-heuristics including the Grey Wolf Optimiser [35] (GWO) and a self-adaptive version of differential evolution (LSHADE-EpSin [36]). Moreover, we propose two novel bi-level optimisation methods consisting of a global search method that works in the upper-level combined with a local search method in the lower-level. In total, nine optimisation methods are applied and compared in order to maximise the absorbed power and minimise the LCoE in a real wave regime from the southern coast of Marettimo (an island in the Mediterranean Sea). We also improve previous research by modelling waves regimes with a higher granularity of wave-directions.

The experimental outcomes show that a bi-level optimisation technique consisting of a self-adaptive differential evolution search (LSHADE-EpSin) interleaved with Nelder–Mead (NM) simplex direct search outperforms previous heuristic methods used in prior works in terms of convergence rate, higher absorbed power output, and lower levelised cost of energy.

The paper is structured as follows. Section 2 outlines the design of the WEC and the model that is applied to simulate both the absorbed power and LCoE. In the next section, the optimisation problem is described, and Section 4 represents the proposed meta-heuristic methods. The optimisation achievements are presented and considered in Section 5. Finally, Section 6 presents the conclusions of this work and canvasses future work.

2. Modelling

2.1. Wave Energy Converter

A wave energy converter chosen for this case study is a fully submerged cylindrical buoy connected to three tethers to absorb wave power from its motion in multiple degrees-of-freedom (or multiple modes). namely surge, heave and pitch. As shown in Figure 1, the geometry of this WEC is determined by the radius a and height H of the cylinder, tethers inclination angle ${\alpha }_t$, and the angle ${\alpha }_{ap}$ that defines the tether attachment point (from the centre of mass of the buoy). The submergence depth (distance from the undisturbed water level to the top of the buoy) is considered fixed and equal to 2 m regardless of the buoy size. The mass of the buoy is taken as half the displaced mass of water $m_b=0.5{\rho }_{w}V$ (the density of water is ${\rho }_w=1025$ kg/m${}^3$, and the buoy volume is $V=\pi a^{2}H$). The hollow buoy houses three direct mechanical drive power take-off units (each connected to the tether). Each PTO acts as a spring-damper system where stiffness and damping coefficients can be adjusted for each sea state.

2.2. Wave Climate

A potential wave energy development site located near the west coast of Marretimo Island (Italy) in the Mediterranean Sea is chosen for this analysis. According to the WXSD classification [20], this wave climate belongs to resource class 1 due to its low energy content (6.4 kW/m). The k-means clustering method has been applied to extract 10 sea states that represent this wave climate as shown in Figure 2 and listed in

Table 1. Ten irregular sea states that represent the Marettimo deployment site.

Sea State	${\mathit{T}}_{\mathit{p}}$, s	${\mathit{H}}_{\mathit{s}}$, m	Probability O, %
1	3.82	0.24	8.06
2	5.13	0.44	14.62
3	6.20	0.61	17.80
4	7.18	0.90	18.01
5	8.30	0.73	12.10
6	8.43	1.92	9.58
7	9.68	1.08	8.68
8	10.24	2.76	5.78
9	11.56	1.46	3.30
10	12.99	3.69	2.07

. A weighted aggregation of these 10 irregular sea states are used to calculate the annual average power production of the WEC. It is assumed that all waves are unidirectional and propagate in the positive x-direction.

2.3. Equations of Motion

The following time-domain model describes the WEC response under the wave and PTO loads:

$$\mathbf{M}\ddot{\mathbf{x}}\left(t\right)={\mathbf{F}}_{exc}\left(t\right)+{\mathbf{F}}_{rad}\left(t\right)+{\mathbf{F}}_{visc}\left(t\right)+{\mathbf{F}}_{buoy}\left(t\right)+{\mathbf{F}}_{tens}\left(t\right),$$

(1)

where the $\mathbf{x}\in {\mathbb{R}}^{6\times 1}$ is the buoy position vector in $Oxyz$ coordinate system, $\mathbf{M}$ is a mass matrix, ${\mathbf{F}}_{exc}$ is the wave excitation force, ${\mathbf{F}}_{rad}$ is the wave radiation force, ${\mathbf{F}}_{visc}$ is the viscous drag force, ${\mathbf{F}}_{buoy}$ is the buoyancy force, ${\mathbf{F}}_{tens}$ is the tether tension force expressed in the Cartesian space that includes the pre-tension force and control (PTO) forces. The force acting along the k-th tether can be modelled as F_t,k$=F_{t0}+K_{pto}\Delta {\ell }_k+B_{pto}\Delta {\dot{\ell }}_k(k=1\dots 3)$ being proportional to the tether extension $\Delta \ell$, the rate of change of the tether length $\Delta \dot{\ell }$ and includes the initial tension $F_{t0}$. The PTO stiffness $K_{pto}$ and damping $B_{pto}$ coefficients take the same values for all three tethers. The transformation between the buoy velocity $\dot{\mathbf{x}}$ and the tether velocity vector $\dot{\mathbf{q}}={\left[\Delta {\dot{\ell }}_1\Delta {\dot{\ell }}_2\Delta {\dot{\ell }}_3\right]}^{\mathsf{T}}$ has a form of $\dot{\mathbf{q}}\left(t\right)={\mathbf{J}}^{-1}\left(\mathbf{x}\right)\dot{\mathbf{x}}\left(t\right)$, where ${\mathbf{J}}^{-1}\left(\mathbf{x}\right)\in {\mathbb{R}}^{3\times 6}$ is the inverse kinematic Jacobian that depends on the buoy position at each time instance [34]. So the tether force vector can be converted to the Cartesian space according to ${\mathbf{F}}_{tens}=-{\mathbf{J}}^{-\mathsf{T}}{\mathbf{F}}_t$.

The time-domain model in Equation (1) has a relatively high computation time and may not be suitable for optimisation purposes when a large number of evaluations are required. If to assume that all processes are Gaussian, it is possible to derive a spectral-domain model that can capture all required nonlinear forces using statistical linearisation technique [38,39]. The spectral-domain model approximates the system dynamics in the frequency domain by replacing all nonlinear terms with equivalent linear matrices [40]. The dynamic model in Equation (1) has two sources of nonlinearity: the viscous drag force F_visc and the generalised tether tension force F_tens. Due to the fact that geometric nonlinearity contained within Ftens is much weaker than the quadratic nonlinearity in F_visc, F_tens can be linearised around the zero position without loss of accuracy for the proposed configuration. If nonlinear effects from tethers become relevant, the equivalent terms can be derived as shown in [38,41,42]. Moreover, it should be noted that other nonlinear forces can be included in the model but omitted in this study, e.g., nonlinear Froude–Krylov force that becomes relevant when the buoy experiences large motion amplitudes [43]. As a result, a nonlinear dynamic Equation (1) is replaced by the equivalent frequency domain model:

$$\begin{array}{c}\hfill [-{\omega }^2(\mathbf{M}+\mathbf{A}(\omega ))+i\omega (\mathbf{B}(\omega )+{\mathbf{B}}_{pto}+{\mathbf{B}}_{eq})+{\mathbf{K}}_{pto}]\widehat{\mathbf{x}}(\omega )={\widehat{\mathbf{F}}}_{exc}(\omega ),\end{array}$$

(2)

where $\mathbf{x}\left(t\right)=\mathcal{R}\mathcal{e}\{\widehat{\mathbf{x}} e^{i\omega t}\}$, the radiation force is expressed using the frequency dependent added mass A(ω) and radiation damping matrix $\mathbf{B}\left(\omega \right)$, ${\widehat{\mathbf{F}}}_{rad}(\omega )=-(-{\omega }^2\mathbf{A}(\omega )+i\omega \mathbf{B}(\omega )) \widehat{\mathbf{x}}(\omega )$, the tether tension force is linearised as ${\widehat{\mathbf{F}}}_{tens}\left(\omega \right)=-(i\omega {\mathbf{B}}_{pto}+{\mathbf{K}}_{pto})\widehat{\mathbf{x}}\left(\omega \right)$ (see [44] for more details), and the viscous drag force is replaced by ${\widehat{\mathbf{F}}}_{visc}\left(\omega \right)=-i\omega {\mathbf{B}}_{eq}\widehat{\mathbf{x}}\left(\omega \right)$. The equivalent damping term B_eq is unknown and determined iteratively (for each wave condition separately) using the procedure explained in [38]:

$$\begin{array}{c}\hfill {\mathbf{B}}_{eq}=-{\displaystyle \langle }\frac{\partial {\mathbf{F}}_{visc}}{\partial \dot{\mathbf{x}}}{\displaystyle \rangle },\end{array}$$

(3)

where 〈·〉 indicates mathematical expectation, and the viscous force is interpreted as:

$$\begin{array}{c}\hfill {\mathbf{F}}_{visc}=-\frac{1}{2}{\rho }_w{\mathbf{C}}_d{\mathbf{A}}_d\left(\right||\dot{\mathbf{x}}\left|\right|\odot \dot{\mathbf{x}}),\end{array}$$

(4)

ρw is the density of water, C_d and A_d are the matrices of the drag coefficients and the cross-section areas of the buoy perpendicular to the direction of motion respectively, and ⊙ represents the Hadamard product (element-wise multiplication). Note that only the body velocity (not the relative fluid/body velocity) has been considered in the drag force formulation. A detailed methodology of how to incorporate the wave-particle velocity into the spectral-domain model is demonstrated in [45].

The following iterative procedure is used to estimate ${\mathbf{B}}_{eq}$ and approximate the response of the WEC in irregular waves:

Step 1.

Define the sea state and corresponding incident wave spectrum $S_{\eta }\left(\omega \right)$.

Step 2.

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

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

(5)

where ${\widehat{\mathbf{f}}}_{exc}$ is the vector of excitation force coefficients, and ${\left(\right)}^{*}$ denotes the conjugate transpose of a vector/matrix.

Step 3.

Calculate the WEC response matrix assuming ${\mathbf{B}}_{eq}={\mathbf{0}}_{6\times 6}$ in the first iteration:

$$\mathbf{H}\left(\omega \right)=[-{\omega }^2(\mathbf{M}+\mathbf{A}\left(\omega \right))+i\omega (\mathbf{B}\left(\omega \right)+{\mathbf{B}}_{pto}+{\mathbf{B}}_{eq})+{\mathbf{K}}_{pto}]{}^{-1}.$$

(6)

Step 4.

Establish the power spectral density matrix of the buoy motion:

$${\mathbf{S}}_{\mathbf{x}}\left(\omega \right)=\mathbf{H}\left(\omega \right){\mathbf{S}}_{\mathbf{F}}\left(\omega \right){\mathbf{H}}^{*}\left(\omega \right).$$

(7)

Step 5.

Calculate the covariance matrix of the WEC velocity:

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

(8)

Step 6.

Estimate the equivalent damping matrix ${\mathbf{B}}_{eq}$ using the analytical expression from [38]:

$$\begin{array}{c}\hfill {\mathbf{B}}_{eq}=-{\displaystyle \langle }\frac{\partial {\mathbf{F}}_{visc}}{\partial \dot{\mathbf{x}}}{\displaystyle \rangle }=\frac{1}{2}\sqrt{\frac{8}{\pi }}{\rho }_w{\mathbf{C}}_d{\mathbf{A}}_d{\sigma }_{\dot{\mathbf{x}}}^2.\end{array}$$

(9)

Step 7.

Check the convergence criteria:

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

(10)

where n corresponds to the iteration number, and the threshold is set to $\delta =0.01$. If this condition is not satisfied, go to Step 3.

It can take up to 10 iterations to estimate B_eq and the WEC response in irregular waves. Once calculated, the average power absorbed by each PTO unit $k=1\dots 3$ is calculated as [38]:

$${\overline{P}}_k=B_{pto}{\sigma }_{{\dot{q}}_k}^2,$$

(11)

where ${\sigma }_{{\dot{q}}_k}^2$ is the variance of the tether length rate change $\dot{\mathbf{q}}$:

$${\sigma }_{{\dot{q}}_k}^2={\int }_0^{\infty }{\omega }^2{S}_{q_k}\left(\omega \right)d\omega ,$$

(12)

and the transformation between the Cartesian coordinate system and the tether space is obtained using ${\mathbf{S}}_{\mathbf{q}}\left(\omega \right)={\mathbf{J}}_0^{-1}{\mathbf{S}}_{\mathbf{x}}\left(\omega \right){\mathbf{J}}_0^{-T}$, where ${\mathbf{J}}_0^{-1}={\mathbf{J}}^{-1}\left({\mathbf{x}}_0\right)$ is linearised about the nominal operating position ${\mathbf{x}}_0={\mathbf{0}}_{6\times 1}$.

The total power generated by three PTO units in an irregular wave with the significant wave height $H_s$ and peak wave period $T_p$ is:

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

(13)

The expected average annual power production from the WEC for a specific deployment site is estimated as:

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

(14)

where the matrix $\mathbf{O}(H_s,T_p)$ contains the occurrence probability of each sea state within the wave climate.

To demonstrate that the spectral-domain model is an effective tool that can fully capture the nonlinear dynamics of the considered WEC while significantly decreasing the computation time, a comparison of average power estimated using three different models is shown in Figure 3. The frequency-domain model is implemented based on Equation (2) assuming B_eq = 0, the spectral-domain model is specified in Equation (2) where B_eq is estimated iteratively for each sea state, and the time-domain model is represented by Equation (1). Good agreement is achieved between the spectral-domain and time-domain models, while the frequency domain model significantly overestimates power generation potential of the WEC.

2.4. Economic Model

Levelised cost of energy (LCoE) is used to measure the economic attractiveness of the proposed energy project. Due to the lack of publicly available information of the detailed cost estimations for wave energy technology, [46] proposed to approximate LCoE by the following equation:

$$\mathrm{LCOE}{\displaystyle (}\frac{€}{\mathrm{kWh}}{\displaystyle )}=\mathrm{RDC}\times {\displaystyle (}\frac{\mathrm{Energy}\left(\mathrm{MWh}\right)}{\mathrm{Mass}\left(\mathrm{kg}\right)}{\displaystyle )}{}^{-0.5},$$

(15)

where RDC is a site-specific coefficient that is set to 1 in this study, the characteristic mass of the system includes the mass of the buoy and the anchoring system.

The characteristic mass of the WEC is calculated using the following assumptions:

-

The mass of the buoy is calculated based on a given geometry as $m_b=0.5{\rho }_w\pi a^{2}H$;

-

The needed mass of the anchoring system (three piles) relays on the tether tension associated with buoyancy and the wave force, and can be approximated by $m_{as}\approx 0.116F_t^{peak}$ using case presented in [47] as a reference. The tether peak force ($99\%=2.57{\sigma }_{F_t}$) is estimated from the spectral-domain model.

As a consequence, the LCoE model applied in this research is:

$$\mathrm{LCOE}={\displaystyle (}\frac{8760P_{AAP}}{m_b+m_{as}}{\displaystyle )}{}^{-0.5}.$$

(16)

2.5. Implementation

To estimate the power output and LCoE for any WEC geometry, Equation (2) is solved inMATLAB. The mass matrix has a diagonal form M = diag(m_b, m_b, m_b, I_xx, I_yy, I_zz) with moments of inertia calculated for the cylindrical body. Hydrodynamic parameters of the WEC, including the added mass A(ω), hydrodynamic damping B(ω), and excitation force vector ${\widehat{\mathbf{F}}}_{exc}\left(\omega \right)$are estimated using a semi-analytical model [48,49]. B_eq is calculated based on the iterative procedure explained in Section 2.3.

Even though only one geometric shape (vertical cylinder) is used in the study, the magnitude of the viscous drag force, and the corresponding B_eq, are highly dependent of the ratio between the cylinder height to its diameter, especially for the heave mode. Therefore, it order to develop an optimisation procedure that can accommodate WEC geometries with various aspect ratios ($H/a$), the drag coefficient in heave is expressed as a function $C_{d_3}=-0.12(H/a)+1.2$ based on published data [50] shown in Figure 4. Drag coefficients in other directions are not sensitive to the cylinder aspect ratio and are kept fixed $C_{d_1}=C_{d_2}=1$ for surge and sway, and $C_{d_4}=C_{d_5}=0.2$ for roll and pitch. The irregular waves from Table 1 are modelled using the Bretschneider (modified Pierson–Moskowitz) spectrum according to [51].

3. Optimisation Configuration Models

In this research, The optimisation decision variables of the cylinder are including the radius of the buoy a, the aspect ratio that is considered as the proportion of the height over the radius of the buoy ($H/a$), two tether angles (attachment ${\alpha }_{ap}$ and inclination angle ${\alpha }_t$), two vectors of power take-off parameters, damping and stiffness coefficients represented ${\mathbf{b}}_{pto}={[B_{pto}^{\left(1\right)},B_{pto}^{\left(2\right)},\dots ,B_{pto}^{\left(N\right)}]}^T$ and ${\mathbf{k}}_{pto}={[K_{pto}^{\left(1\right)},K_{pto}^{\left(2\right)},\dots ,K_{pto}^{\left(N\right)}]}^T$, respectively. The length of each PTO vector is $N=10$. The whole number of decision designs are 24 which should be optimised in the following:

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

(17)

$${\mathbf{z}}_{\mathbf{2}}=[a,(H/a),{\alpha }_t,{\alpha }_{ap},{\mathbf{k}}_{pto}\in {\mathbb{R}}^{N\times 1},{\mathbf{b}}_{pto}\in {\mathbb{R}}^{N\times 1}].$$

(18)

We apply two fitness functions in order to maximise the power output and minimise the LCoE.

(i)

The average annual produce power output computed utilising Equation (14), that is maximised as

$$f_{O1}=\underset{\mathbf{z}}{argmax}{P}_{AAP}\left(\mathbf{z}\right),\mathrm{subject}\mathrm{to}:{\mathbf{z}}_{\mathbf{1}}\in [{\mathbf{z}}_{min},{\mathbf{z}}_{max}]$$

(19)

(ii)

The LCoE is minimised using the below equation that is specified in Equation (16):

$$f_{O2}=\underset{\mathbf{z}}{argmin}\mathrm{LCOE}\left(\mathbf{z}\right),\mathrm{subject}\mathrm{to}:{\mathbf{z}}_{\mathbf{2}}\in [{\mathbf{z}}_{min},{\mathbf{z}}_{max}]$$

(20)

Table 2. Boundary constraints of the cylinder parameters.

Parameter	Unit	Min	Max	Length
radius, a	m	1	20	1
height, H	m	1	30	1
aspect ratio, $(H/a)$		0.4	2	1
Tether inclination angle, ${\alpha }_t$	deg	10	80	1
Tether attachment angle, ${\alpha }_{ap}$	deg	10	80	1
PTO stiffness, $K_{pto}$	N/m	${10}^3$	${10}^8$	10
PTO damping, $B_{pto}$	N/(m/s)	${10}^3$	${10}^8$	10

shows the ranges of all design variables which are involved in the optimisation process.

4. Optimisation Algorithms

In this paper, we focus on two widespread optimisation strategies in order to maximise harnessed power and minimise the levelised cost of energy (LCoE) of a fully-submerged three-tether WEC. The first approach applies optimisation algorithms to all decision variables simultaneously. These design variables consist of the buoy geometry parameters (radius a, height H and aspect ratio ($H/a$)), the tether angles (inclination angle ${\alpha }_t$ and the tether attachment angle ${\alpha }_{ap}$), and the PTO parameters (spring stiffness ${\mathbf{k}}_{pto}$ and damping coefficients ${\mathbf{k}}_{pto}$). In total, there are 24 parameters that are optimised all-at-once.

The second strategy is to apply bi-level optimisation methods [52], which solve the problem using a two-level optimisation procedure, where one optimisation problem is nested within the other. The outer optimisation task is generally regarded as the upper-level optimisation problem, and the interior one is recognised as the lower-level optimisation problem. A significant characteristic of the bi-level optimisation problem is that the fitness functions of each level may be partly defined by variables advised by other levels. Following this strategy, we propose two bi-level optimisation methods and compare their performance with seven other well-known global search methods. The details of the optimisation algorithms performed for each strategy are outlined in

Table 3. The details of the optimisation methods settings. All approaches are restricted to the same evaluation number.

Methods	Settings
Nelder–Mead [53]	Nelder–Mead simplex direct search (NM)
1+1EA [54]	mutation step sizes are ${\sigma }_a={\xi }_1\times (U_a-L_a)$, ${\sigma }_H={\xi }_1\times (U_H-L_H)$,${\sigma }_{{\alpha }_t}={\sigma }_{{\alpha }_{ap}}={\xi }_1\times (U_{{\alpha }_t}-L_{{\alpha }_t})$, ${\sigma }_{K_{pto}}={\sigma }_{B_{pto}}={\xi }_2\times (U_{K_{pto}}-L_{K_{pto}})$, and Probability mutation rate=$\frac{1}{N}$, ${\xi }_1=0.3,{\xi }_2=0.01$
CMA-ES [55]	with the default settings and $\lambda =13$;
PSO [56]	with $\lambda =25$, $c_1=1.5$, $c_2=2$, $\omega =1$ ( decreased with a damping ratio $w_f=0.99$ exponentially);
GWO [35]	with $\lambda$= 25, $\alpha =2$ (linearly decreased to zero)
DE [57]	with $\lambda =25$, $F=0.5$, $P_{cr}=0.8$
SaDE [58]	with $\lambda =25$, $LP=50$, $NumSt=4$
LSHADE-EpSin [36]	$\lambda =25$, historical memory size $H=5$, $Nu{m}_{LS}=10$
Bi-level-1	SaDE +NM, WEC’s dimensions and tether angles are optimised in the lower-level, default settings of SaDE
Bi-level-2	LSHADE-EpSin + NM, WEC’s dimensions and tether angles are optimised in the lower-level, default settings of LSHADE-EpSin

.

4.1. All-at-Once Optimisation

Various factors associated with WEC design, tether angles and PTO parameters combined to form a non-convex, dynamic, constrained and large-scale optimisation problem. These challenges serve as our primary motivation for applying the meta-heuristics like evolutionary and swarm optimisation algorithms. We apply and compare the performance of seven well-known meta-heuristics that reliably optimise all decision variables of WECs all-at-once. This optimisation process leads to maximise the produced power and minimise the levelised cost of energy. The optimisation methods applied in this research include 1+1EA [59]; Differential Evolution (DE) [57], Covariance matrix adaptation evolution strategy (CMA-ES) [55], Particle Swarm Optimisation (PSO) [56], Grey Wolf Optimiser (GWO) [35] and two state-of-the-art self-adaptive optimisation methods including SaDE [58] and LSHADE-EpSin [36].

4.1.1. L-SHADE with an Ensemble Pool of Sinusoidal Parameter Adaptation (LSHADE-EpSin)

The Differential Evolution (DE) algorithm, and its adaptive and self-adaptive variants, are simple and robust evolutionary algorithms. Researchers from various fields of science and engineering have applied DE algorithms to various optimisation problems, notwithstanding problems with characteristic such being continuous, multi-modal, combinatorial or mixed variable. DE is able to obtain superior optimisation results across widely encountered real-world engineering problems [60,61]. Among a wide range of self-adaptive DE algorithms, LSHADE-EpSin performs outstandingly in solving different benchmarks and real-world problems [36]. LSHADE-EpSin is a modified version of the L-SHADE algorithm [62] with linear population size reduction and an ensemble pool of sinusoidal parameter adaptations. L-SHADE is a developed version of the SHADE algorithm [63] that practices a history-based parameter adaptation trajectory based on the JADE algorithm [64] which proposed the novel mutation strategy (current/to/pbest).

Mutation Strategy with External Archive

In LSHADE-EpSin, one of the best-performing mutation strategies for generating promising mutant vectors during the optimisation process is current-to-pbest/1 which is initially proposed by JADE. This mutation strategy can be seen in Equation (21).

$$v_{i,g}=x_{i,g}+F_{i,g}(x_{pbest,g}-x_{i,g})+F_{i,g}(x_{r_1,g}-x_{r_2,g})$$

(21)

where $x_{pbest,g}$ is chosen from the best solutions $N\times p(p\in [0,1\left]\right)$ of the current parent population (g). $x_{r_1,g}$ is randomly taken from the population and $x_{r_2,g}$ is randomly chosen from a combination of the current population and the external archive (A). The external archive keeps a record of the lower-ranking parents recently replaced by offspring.

Ensemble of Parameter Adaptation

An ensemble of parameter configurations is used in LSHADE-EpSin to control the adaptation of parameters. The adaptive parameters are associated with a combination of two sinusoidal formulas to adjust the scaling factor. Firstly, a non-adaptive sinusoidal adjustment technique is used to adjust the scale factor ($F_{i,g}$) which decreases during the optimisation process. Equation (22) shows this non-adaptive technique.

$$F_{i,g}=\frac{1}{2}\times (sin(2\pi \times freq\times g_{s_1}+\pi )\times \frac{ite{r}_{max}-g_{s_1}}{ite{r}_{max}}+1)$$

(22)

where $freq$ describes a pre-defined frequency for the sinusoidal function and $iter$ denotes the current generation number ($g_{s_1}<=\frac{ite{r}_{max}}{2}$). The second strategy for the adjustment of the scale factor is an adaptive sinusoidal adjustment method. This formulation can be seen in Equation (23).

$$F_{i,g}=\frac{1}{2}\times (sin(2\pi \times freq\times g_{s_1})\times \frac{g_{s_1}}{ite{r}_{max}}+1)$$

(23)

where $freq$ is an adaptive frequency based on a Cauchy distribution and a successful history-based of settings. $iter$ denotes the current generation number. One of the most effective DE parameter adaptation techniques is recording an archive of both mutation factors and probabilities of crossover based on their success during the optimisation process. The control parameters history-based was proposed by Zhang et al. [64] in JADE. In each generation of JADE, in order to generate an offspring, we have an array of the crossover probability rate that is produced based on a normal distribution of the mean μ_CR and variance at 0.1. The successful crossover probabilities S_CR are recorded and updated at each generation. The μ_CR is initialised by 0.5 and in the next generation it is updated by Equation (24).

$$\begin{array}{c}\hfill {\mu }_{CR}=(1-c)\times {\mu }_{CR}+c\times mea{n}_A\left(S_{CR}\right)\end{array}$$

(24)

where c is a constant generated between 0 and 1 randomly and mean_A is a simple arithmetic mean. Likewise, the mutation factor F_i of each x_i is separately generated at each generation, as stated in a Cauchy distribution with the mean μ_F and scale parameter 0.1. (Equation (25))

$$F_i={randc}_i({\mu }_F, 0.1)$$

(25)

where the ${randc}_i$ is the Cauchy distribution. All successful mutation factors are archived and point out as a set of S_F at the end of each generation. The value of μ_F is updated using Equation (26).

$$\begin{array}{c}\hfill {\mu }_F=(1-c)\times +c\times meanL\left(SF\right)\end{array}$$

(26)

where mean_L is the Lehmer mean [65] and computed as follows:

$$\begin{array}{c}\hfill meanL\left(SF\right)=\frac{\sum F\in SF{F}^2}{{\sum }_{F\in S_F}F}\end{array}$$

(27)

Linear Population Size Reduction

The LSHADE-EpSin algorithm benefits from a linear reduction in population size to fit the population size (N) iteratively at each generation as exposed in the following equation:

$$\begin{array}{c}\hfill N_{g+1}=Round\left[\right(\frac{N_{min}-N_{max}}{ite{r}_{max}})\times iter+N_{max}]\end{array}$$

(28)

where N_min is the minimum population size, and initialised at 4 that is required to make the current-to-pbest mutation strategy. The four required solutions are $x_i,x_{best}^p,x_{r_1}\mathrm{and}{x}_{r_2}.$ The mutant vector of this strategy is generated using Equation (29).

$$\begin{array}{c}\hfill V_{i,g}=x_{i,g}+F_i\times (x_{best,g}^p-x_{i,g})+F_i(x_{r_1,g}-x_{r_2,g})\end{array}$$

(29)

Local Search

In order to extend the exploitation capability of LSHADE-EpSin, a stochastic local search $is$ proposed that works based on Gaussian Walks. The local search is activated when the population size is less than 20 ($N_{ini}=25$), and 25 random samples are evaluated to exploit the neighbourhood of the best-found design among the current population. The Gaussian walks applied can be seen in Equation (30).

$$y_i=\mathcal{N}({\mu }_b,\sigma )+(r_1\times x_{best}-r_2\times x_i)$$

(30)

where $x_{best}$ is the best-found solution in the local search and ${\mu }_b$ is equal to $x_{best}$. $r_1$ and $r_2$ are two uniform random numbers from the range of $[0,1]$. Besides, the standard deviation ($\sigma$) of this Gaussian Walks is calculated using Equation (31).

$$\sigma =\frac{log\left(iter\right)}{iter}\times (x_i-x_{best})$$

(31)

4.2. Bi-Level Optimisation

In this paper, we propose two bi-level optimisation methods, including Bi-level-1 (SaDE+NM) and Bi-level-2 (LSHADE-EpSin+NM). We also provide a general formulation in order to maximise the harnessed power and minimise the LCoE of a cylindrical wave energy converter. These proposed approaches comprise two levels of optimisation tasks where one optimisation process is nested within the other. The exterior optimisation method (which is a global search method) is referred to as the leader’s (upper level) optimisation process. In the upper level, we apply two self-adaptive meta-heuristics, including Self-adaptive DE (SaDE) and LSHADE-EpSin. Both methods improve the ability of an adaptive learning strategy to fine-tune the control parameters and mutation strategy and demonstrate a considerable performance in optimising real engineering problems [66,67].

In the second level, the internal method is recognised as the follower’s (lower level) optimisation process. In the current study, the inner method is a Nelder–Mead (NM) simplex search method [68]. NM simplex is a downhill local search method, and it is straightforward to hybridise combine with other meta-heuristic methods. The primary reason for such hybridisation (or for using NM as the lower-level in a bi-level method) is to tune a more suitable trade-off between global optimality and computational budgets [69,70].

Figure 5 shows that the proposed bi-level optimisation framework consists of a global search method designed to optimise all decision variables in the upper-level, and both geometry parameters (radius and height) that given from upper-level decision vector are optimising in the lower-level. To adjust the geometry parameters of the cylinder, we use a local search method. The best-found geometry configuration in the lower-level will be replaced in the upper-level decision variables.

The pseudo-code of the proposed Bi-level-2 algorithm is shown in Algorithm 1. It can be seen that the algorithm is divided into two primary sections. At the top level, we have a self-adaptive DE (LSHADE-EpSin) employing two strategies to adjust the control parameters. These strategies are (1) adaptive sinusoidal increasing adjustment and (2) non-adaptive sinusoidal decreasing adjustment. The benefit of this ensemble approach is that it allows the algorithm to converge to a sufficient balance [36] between searching the neighbourhood of current bet-found solutions and the exploration of non-visited search space zones. In the lower-level, there are two nested inner local search methods. The initial local search is used to explore the search space of the cylinder dimensions (radius and height) where other decision variables are fixed. Next, both tether angles (inclination and attachment) are optimised using the second local search. In order to save computational budget, we define a performance criterion for both local search methods. This condition evaluates the local search performance; if the obtained power improvement cannot satisfy the criterion, Bi-level-2 will withdraw the local optimisation process and allocate the remaining budget to the global search method.

Tuning the Local Search

One of the significant parameters of the bi-level optimisation method is the maximum evaluation number ($Ma{x}_{eval}$) of the local search (NM). Tuning this variable plays an important role in obtaining a greater balance between saving on the computational budget and converging to the local optimum as much as possible. In order to tune the $Ma{x}_{eval}$, we perform the local search to optimise the WEC geometry parameters ($a,H$) and keep the other decision variables fixed. This experiment iterates ten times with different initial solutions. Meanwhile, the same tuning process runs to optimise both tether angles. Figure 6 shows the convergence curves of these experiments. We observe that the local search converges rapidly to a local optimum in the geometry and tether angles optimisation processes after 20 and 40 iterations, respectively on average. Therefore, we set the $Ma{x}_{eval}$ of the local search to 20 and 40 iterations.

5. Optimisation Results and Discussions

5.1. Multi-Modality of Search Space

In order to characterise the search space, we perform an experiment using a Nelder–Mead (NM) search method. Twenty random initial configurations are generated and NM is applied to optimise the absorbed power output. Figure 7 shows the trajectory of the NM performance during the optimisation process. It can be seen that the majority of the trajectories in the cylinder dimension (subplot (a)) converged to a specific area of the search space as expected. This is because large WECs can harness more power than small ones. The second observation is that the PTO search space is not uni-modal and each trajectory converged to different configurations (subfigure (c,d,e)).

5.2. Power Landscape Analysis

With regard to evaluating the impact of each buoy design variable on the level of produced power, we perform a sensitivity analysis experiment. Here, we assume both tether angles are kept fixed at 45°; note that this size is not optimal, because tether angles should be adjusted based on the buoy’s dimensions, as recommended by prior works [34]. Moreover, the search space of the $K_{pto}$ and $B_{pto}$ parameters are discretised, where each interval is ${10}^6$. In the next step, for each discrete configuration of PTO parameters, we evaluate the importance of the cylinder dimensions ($a,H$) using a grid search technique where the discretisation step size is 1 (m).

The results are shown in Figure 8, which includes 400 sub-figures. Each sub-figure represents the relationship of the cylinder radius and height sizes with the absorbed power, where the $K_{pto}$ and $B_{pto}$ are fixed. It is important to note that a variation in the size of the radius has a more substantial effect on the power output than a variation in the cylinder height. In this wide power landscape, we can see that the maximum produced powers are achieved when the PTO parameters are assigned around ${10}^7$, and the buoy radius and height sizes are large. However, it should be noted that the effect of PTO parameters on the absorbed power is more significant than the size of the cylinder dimensions.

Algorithm 1 Bi-level Optimisation method (LSHADE-EpSin+NM)
procedureBi-level Optimisation method
Initialization     $\mathit{P}=\{\langle a_1,H_1,{\alpha }_{t_1},{\alpha }_{{ap}_1},K_1^1,...,K_1^{10},B_1^1,...,B_1^{10}\rangle ,\dots$
$\dots ,\langle a_N,H_N,{\alpha }_{t_N},{\alpha }_{{ap}_N},K_N^1,...,K_N^{10},B_N^1,...,B_N^{10}\rangle \}$	▹ initial population
M:$\mu$F =$\mu$CR=0.5	▹ initialise memory of first control settings
M${}_{freq}$:$\mu$freq = 0.5,$Imp-rat{e}_d=Imp-rat{e}_{\alpha }=1$	▹ initialise memory of second control settings
Upper-Level (Global search method)
for $iter$ in $ite{r}_{max}$do	▹ termination criteria
if $iter>\frac{ite{r}_{max}}{2}$ then
Call second control parameter settings
$S_F=S_{CR}=\varnothing$	▹ Reset successful mean vectors
$r_i=rand(1,H)$	▹ Generate a random index, H is memory size
$F_i=randc(\mu F_{r_i},0.1)$,$C{R}_i=randn(\mu C{R}_{r_i},0.1)$
end if
if $iter\le \frac{ite{r}_{max}}{2}$ then
Call first control parameter settings
$c=rand(0,1)$
if $c<0.5$ then
$F_i=\frac{1}{2}\times (sin(2\pi \times freq\times iter+\pi )\times \frac{ite{r}_{max}-iter}{ite{r}_{max}}+1)$
else
$F_i=\frac{1}{2}\times (sin(2\pi \times freq\times iter)\times \frac{iter}{ite{r}_{max}}+1)$
end if
Generate $C{R}_i$ same as first control parameters (Equation 23)
end if
for $i=1$ to N do
Generate $p=rand(0,1)\times n$, $n=0.1\times N$
$v_i=x_i+F_i\times (x_{pbest}-x_i)+F_i\times (x_{r_1}-x_{r_2})$	▹ Mutation $current$-$to$-$pbest/1$
${\mathit{u}}_{\mathit{i},\mathit{iter}}^{\mathit{j}}=\{\begin{array}{cc}{\mathit{v}}_{\mathit{i},\mathit{iter}}^{\mathit{j}},\hfill & \mathrm{if}(\mathit{rand}<{\mathit{CR}}_{\mathit{i}})\mathrm{or}(j==j_{rand})\hfill \\ P_{i,iter}^j,\hfill & \mathrm{Otherwise}\hfill \end{array}$	▹ Binomial Crossover
${\mathit{P}}_{\mathit{i},\mathit{iter}+\mathit{1}}=\{\begin{array}{cc}{\mathit{u}}_{\mathit{i},\mathit{iter}},\hfill & \mathrm{if}(\mathit{f}\left({\mathit{u}}_{\mathit{i},\mathit{iter}}\right)>\mathit{f}\left({\mathit{P}}_{\mathit{i},\mathit{iter}}\right))Maximisation\hfill \\ P_{i,iter},\hfill & \mathrm{Otherwise}\hfill \end{array}$	▹ Selection
Store successful $F_i$ and $C{R}_i$
end for     Update the memory according to used settings     Update the population size by Equation (28)     $N_{diff}=N_g-N_{g+1}$     Sort $P_{iter}$ based on the fitness function     Remove worst solutions $N_{diff}$ from $P_{iter}$ AND Select the best solution $P_{best}$
Lower-Level (Local search method)
if $Imp-rat{e}_d$> 0.001% then	▹ Optimise Cylinder dimension
$P_{best}(a,H)=Nelder-Mead(P_{best}(a,H),Ma{x}_{eval})$      Compute improvement rate $Imp-rat{e}_d$     end if
if $Imp-rat{e}_{\alpha }$> 0.001% then	▹ Optimise tether angles
$P_{best}({\alpha }_t,{\alpha }_{ap})=Nelder-Mead(P_{best}({\alpha }_t,{\alpha }_{ap}),Ma{x}_{eval})$
Compute improvement rate $Imp-rat{e}_{\alpha }$
end if
Update $P_{iter}^{best}$ by the best-found NM configurations
end for
end procedure

5.3. The Annual Average Power Output Maximisation

In this section, we describe the optimisation results of our cylinder design experiments in order to maximise the annual average power output. Furthermore, we compare the performance of the optimisation algorithms outlined above in terms of best-found designs and speed of convergence.

Table 4. Best-found design parameters in order to maximise the average annual absorbed power.

Parameter	1+1EA	CMA-ES	PSO	GWO	DE	SaDE	LSHADE-EpSin	Bi-Level-1	Bi-Level-2
a [m]	16.62	16.10	19.99	16.68	15.46	15.50	15.49	15.61	14.51
H [m]	30	30	14.80	30	30	30	30	30	30
${\alpha }_t$ [deg]	70	26	60	14	48	26	39	50	10
${\alpha }_{ap}$ [deg]	10	13	63	28	10	11	29	40	67
${}_{i=1}^{N_k}{K}_{pto}(\times {10}^7)$	0.665	0.863	3.796	1.51	1.894	2.883	0.882	0.665	0.514
${}_{i=1}^{N_B}{B}_{pto}(\times {10}^7)$	2.765	9.928	4.676	1.51	3.775	4.036	2.479	2.095	1.129
$P_{AAP}\left[\mathrm{K}\mathrm{W}\right]$	259	248	239	261	259	261	262	265	279

reports the best-found cylinder designs using seven meta-heuristics and two new bi-level optimisation methods that produced the highest power output among all ten runs. Furthermore, it can be seen that Bi-level-2 performs better than other applied optimisation methods and that it can produce a considerable amount of power of 279 kW. The second observation is that almost all (8 out of 9) optimisation methods converged to the cylinder of 15 m radius with the largest possible height of 30 m. However, it should be noted that producing electricity using such large WECs can be expensive, due to the high manufacturing costs. In terms of the angles and PTO settings, a large range of values is proposed by all optimisation methods even though the maximised power output is not dramatically different. This fact proves that it is not straightforward to optimise a multi-mode WEC due to the strong dependencies between angles, PTO parameters, and the hydrodynamic model which dominates the power absorption (heave, surge or pitch).

Table 5. Performance comparison of various optimisation methods based on the maximum, minimum and average power output and LCoE of the best-found design per each experiment.

Power [MW]
	1+1EA	CMA-ES	PSO	GWO DE	SaDE	LSHADE-EpSin	Bi-Level-1	Bi-Level-2
Mean	0.2325	0.2329	0.2208	0.2537 0.2501	0.2537	0.2541	0.2551	0.2612
Min	0.1941	0.2121	0.1934	0.2467 0.2327	0.2498	0.2473	0.2526	0.2544
Max	0.2590	0.2476	0.2392	0.2615 0.2589	0.261	0.2621	0.261	0.2792
STD	0.0234	0.0117	0.0181	0.0049 0.0087	0.0036	0.0046	0.0032	0.0088
LCoE
	1+1EA	CMA-ES	PSO	GWO DE	SaDE	LSHADE-EpSin	Bi-Level-1	Bi-Level-2
Mean	0.0443	0.0303	0.0678	0.0315 0.0334	0.0309	0.028	0.0295	0.0268
Min	0.0316	0.0284	0.0556	0.0297 0.0282	0.0277	0.0248	0.0267	0.0243
Max	0.0599	0.0382	0.0794	0.0335 0.0514	0.0329	0.0361	0.0324	0.0285
STD	0.0109	0.0036	0.0071	0.0014 0.0079	0.0019	0.0041	0.0019	0.0012

presents the average best-found power output per each run for all optimisation methods. Bi-level-2 is not only capable of finding the best design configuration; it also performs the best average power output (Figure 9a) compared with other meta-heuristics. In terms of the convergence rate, Figure 10a depicts the applied optimisation method experiments during the 5000 evaluations. As we can see, GWO and LSHADE-SeSin rapidly converge to considerable settings; however, they could not sustain this upward trajectory and converge near locally optimal designs. Obviously, the fastest convergence rate is allocated to Bi-level-2.

5.4. LCoE Minimisation

In this section, we describe the second applied objective function related to LCoE and approximated as a ratio of the generated energy to the significant mass of the system. The best-found LCoE values and their relevant cylinder configurations which are obtained using nine meta-heuristic approaches are shown in

Table 6. Best-found design parameters in order to minimise the LCoE.

Parameter	1+1EA	CMA-ES	PSO	GWO	DE	SaDE	LSHADE-EpSin	Bi-Level-1	Bi-Level-2
a [m]	7.31	6.40	14.32	7.00	7.38	6.57	5.00	6.15	5.00
H/a	0.40	0.40	0.40	0.4	0.40	0.40	0.40	0.40	0.40
αt [deg]	28	29	10	10	31	25	35	31	34
αap [deg]	10	11	10	31	14	11	10	12	10
${}_{i=1}^{N_k}{K}_{pto}(\times {10}^7)$	0.647	0.919	3.90	0.651	3.50	0.383	2.094	0.77	2.071
${}_{i=1}^{N_B}{B}_{pto}(\times {10}^7)$	0.577	0.332	3.52	0.847	1.15	1.481	1.350	0.256	1.914
LCoE	0.0316	0.0284	0.0556	0.0297	0.0287	0.0277	0.0248	0.0267	0.0243
PAAP [kW]	53.1	43.6	131	51.4	64.8	50.6	27.1	43.5	28.3

. Interestingly, all optimisation methods (except PSO) converged to a narrow range of radii between 5 and 7.3 m, with the smallest possible aspect ratio of 0.4. This geometry leads to the fact that the power generation will be dominated by the heave mode rather than surge. Moreover, this is clearly seen from the optimised values of the tether angles as to absorb power from the vertical motion, the tether angles should be closer to vertical leading to α_t < 35 °. Another important finding is that the power production of WECs optimised for LCoE is relatively low leading to 28.3 kW.

Figure 9b shows the box-and-whiskers plot for the best configurations of the WEC which deliver the minimum LCoE for each run for nine search heuristics. It can be seen that the performance of Bi-level-2 is more reliable than that of the other meta-heuristic algorithms we applied. Both LSHADE-EpSin and Bi-level-1 show the next best average performances by 0.028 and 0.0295, respectively.

Investigating the convergence trajectories (Figure 10) from this experiment in the real wave model, it is clear that Bi-level-2 converges faster than other optimisation methods. It is noteworthy that among the seven optimisation methods in the all-at-once strategy, the LSHADE-EpSin convergence speed is substantially better than the others due to both adaptive and non-adaptive strategies in order to adjust the control parameters as well as to conduct an embedded local search in the initial iterations. However, it can be seen that the convergence rate of GWO is considerable in the initial 1000 evaluations.

In order to see the convergence performance of Bi-level optimisation algorithms, the search trajectory of the best agent in each generation for all decision variables is shown in Figure 11. Initially, we can see the high convergence ability of Bi-level-2 compared with DE in order to find and converge to the optimal range of both radius and height. Meanwhile, It can be observed that Bi-level-2 tends to explore promising areas of the tether angle search space broadly, and finally, to exploit the best values.

6. Conclusions

In this paper, two new bi-level optimisation methods are proposed with the aim of maximising the harnessed power output. These methods are also designed to minimise the levelised cost of energy of a fully-submerged, cylindrical WEC with three tethers for the wave climate of a Mediterranean sea site in the west of Sicily, Italy (featuring unidirectional irregular waves). The optimisation of a combination of WEC radius, height, tether inclination and attachment angles, and power take-off parameters is a relatively computationally expensive (5000 evaluations take around 15 h), multi-modal, large-scale and complex problem. These characteristics provided the principal motivation for investigating and proposing a faster and more reliable optimisation technique. With this in mind, we applied a bi-level strategy to optimise the design variables at various levels. A global search method was used at the upper level to optimise the parameters of the whole WEC’s. Furthermore, in the lower level, a Nelder–Mead (NM) simplex search method was applied to adjust the geometry settings and tether angles. To systematically compare the effectiveness of the proposed optimisation method, we considered seven state-of-the-art evolutionary and swarm algorithms. The experimental results showed that the bi-level method can outperform other meta-heuristics in terms of both convergence rate and the quality of WEC’s configuration. Moreover, according to the best-found configurations, if we focus on maximising the harnessed power output without considering the costs, a large cylindrical buoy is recommended. However, the cheapest energy can be delivered by a relatively small WEC with a radius of 5 m and a height of 2 m.