Self-Adaptive Global-Best Harmony Search Algorithm-Based Airflow Control of a Wells-Turbine-Based Oscillating-Water Column

Abstract

The Harmony Search algorithm has attracted a lot of interest in the past years because of its simplicity and efficiency. This led many scientists to develop various variants for many applications. In this paper, four variants of the Harmony search algorithm were implemented and tested to optimize the control design of the Proportional-Integral-derivative (PID) controller in a proposed airflow control scheme. The airflow control strategy has been proposed to deal with the undesired stalling phenomenon of the Wells turbine in an Oscillating Water Column (OWC). To showcase the effectiveness of the Self-Adaptive Global Harmony Search (SGHS) algorithm over traditional tuning methods, a comparative study has been carried out between the optimized PID, the traditionally tuned PID and the uncontrolled OWC system. The results of optimization showed that the Self-Adaptive Global Harmony Search (SGHS) algorithm adapted the best to the problem of the airflow control within the wave energy converter. Moreover, the OWC performance is superior when using the SGHS-tuned PID.

1. Introduction

Marine Renewable Energy (MRE) is trending after solar and wind energy in the R&D sector and energy markets these years. This is due to the fact that solar and wind energy industries have reached the point of maturity and reliability. Moreover, MRE is an abundant source of untapped energy in many forms, in fact, it is estimated that harnessing merely 0.2% of the unused global ocean energy may provide sufficient power to meet power demands [1]. Additionally, a 529 MW of MRE installed capacity has been recorded as operational by the end of 2017 [2]. So, in the efforts of reducing dependency on depleting fossil fuel resources and utilizing an environmentally friendly resource of energy, it was inevitable to turn to ocean energy for countries with low solar and wind energy. This led authorities and policymakers to assess and invest more in MRE in places like Hawaii, India, Thailand, Brazil, and many others [3,4,5,6].

Wave energy is considered the most exploited resource of MRE thanks to its availability and predictability. Many Wave Energy Converters (WECs) were developed, yet no particular concept has reached the commercial maturity [7]. The three major obstacles of the industrialization of WECs are associated with reliability, efficiency, and economic viability, but extensive efforts are being spent on the improvement of these aspects [8]. Reliability is achieved by refining the WEC designs and optimizing control strategies [9,10], while economic viability and efficiency require the optimal selection of the site, the system, and materials in order to reduce costs and environmental impact while guaranteeing maximum energy absorption [11,12].

Many ideas were investigated in order to conquer some of the obstacles, such as WEC integration in ongoing projects of maritime structure and breakwaters [13,14,15]. Other ideas proposed the multiple-use of offshore platforms by combining different energies such as the wind and wave energy converters in [16] or by combining different converters such as the Overtopping device and the Oscillating Water Column (OWC) in [17].

This paper discusses the airflow control of an example of WEC integrated into a breakwater which is the NEREIDA Multiple Oscillating Water Column (MOWC) [18]. This facility is installed in the breakwater of Mutriku in the north of Spain and it is an OWC system based on Wells turbine for wave energy conversion. The main drawback of Power-Take Off (PTO) systems equipped with Wells turbine is their power-limitation due to the occurrence of the stall effect [19,20]. In this context, an airflow control based on a Proportional-Integral-Derivative (PID) controller has been proposed and the recourse to the Harmony Search algorithms and its variants has been proposed in order to facilitate the control design of the parameters.

PID controller tuning is a delicate and complex task when lacking a systematic approach. To solve the problem design of PID controllers, the optimization theory has been proven to be an effective method to tune and optimize the controller parameters [21,22]. Numerous optimization algorithms were investigated and tested with the PID controller in different applications, for instance, the Particle Swarm Optimization (PSO) [23,24], the Water Cycle Algorithm (WCA) [21,24], and the Harmony Search Algorithm (HSA) [22,24].

The remainder of the article has been arranged as the following: Section 2 introduces the modelling section, describing all parts of the OWC system. Section 3 introduces the stalling behavior and the problem formulation of the proposed airflow control. Section 4 presents variants of the Harmony Search algorithm. Section 5 presents the Self-Adaptive Global-Best Harmony Search (SGHS) algorithm adopted to optimally tune the parameters of the Proportional-Integral-Derivative (PID) controller. Section 6 details the tests and simulations carried out to demonstrate the efficiency of the SGHS algorithm, then the performance of the proposed SGHS-PID airflow control in two different wave conditions versus the uncontrolled case. Finally, we finish the article with some concluding remarks in Section 7.

2. Model Statement

This section describes the modelling of the different subsystems of the Oscillating Water Column (OWC) shown in Figure 1, which includes would be including the mathematical models of the wave input, capture chamber, Wells turbine and the Doubly-Fed Induction Generator (DFIG).

2.1. Wave Surface Dynamics

A monochromatic unidirectional wave has been considered as to be the input to the implemented numerical model of the OWC system. There exists numerous wave theories in the literature to express the surface dynamics of ocean wave like the Cnoidal wave theory, second- and higher -order Stokes theory and Airy linear theory [26,27]. In this paper, the Airy wave theory has been adopted because it presents the simplest description and it is the most widely used thanks to it neglecting turbulence, friction losses and other energy losses [27].

The parameters of a wave are detailed in Figure 1, where $SWL$ represents the “Still-Water-Level” and h, called “sea depth”, represents the interval from the sea floor to $SWL$. H marks the interval from wave trough to wave crest called “wave height”. A measures the distance between $SWL$ and the wave crest known as “wave amplitude” and $\lambda$ representing the interval between successive crests known as “wavelength” [27,28]. Therefore, the surface elevation for a sea wave is given as [29,30]:

$$z(x,t)=Asin\left(\omega t-kx\theta \right)=H/2sin\left(\omega t-kx\theta \right)$$

(1)

where $\omega$ is the wave frequency, x is the wave horizontal coordinate, $\theta$ marks the the angular opening from the x-axis to waves’ direction and k represents the wave number linked to $\omega$ with relation (2) as described in [29]:

$$ktanh\left(kh\right)={\omega }^2/g,$$

(2)

where g represents the acceleration gravity.

2.2. Capture Chamber Model

The volume of the air within the Oscillating Water Column’s chamber is defined in [23,28] as:

$$V\left(t\right)=V_c+\frac{w_{c}H}{k}sin\left(k{l}_c/2\right)sin\left(\omega t\right),$$

(3)

where $V_c$, $w_c$ and $l_c$ represent the chamber’s volume, inner width and length, respectively.

The volume flow rate in the chamber can be obtained from Equation (3) and defined as [23,28]:

$$Q\left(t\right)=w_{c}cHsin\left(\frac{k{l}_c}{2}\right)cos\left(\omega t\right),$$

(4)

where $c=w/k$.

Once the chamber’s geometry has been taken into account with Equation (4), the airflow velocity can be described as [23,28]:

$$v_x\left(t\right)=\frac{Q\left(t\right)}{S}=\frac{8Ac{w}_c}{\pi D^2}sin\left(\frac{\pi l_c}{c{T}_w}\right)cos\left(\frac{2\pi }{T_w}t\right)$$

(5)

where $T_w$ is the wave period and D is the duct diameter.

2.3. Wells Turbine Model

The OWC is fitted with a Wells turbine, shown in Figure 2, which is a self-rectifying axial-flow air turbine [30,31]. Self-rectifying air-turbines possess blades with special geometry allowing a unidirectional rotating motion regardless of the airflow direction [32,33,34].

The Wells turbine under study can be mathematically defined by the expressions (6)–(10) given in [20,35]:

$$\begin{array}{ccc}\hfill dp& =& C_{a}K\left(1/\phantom{1a}a\right)\left[v_x^2+{\left(r{\omega }_r\right)}^2\right]\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill K& =& \rho lbn/2\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill T_t& =& r{C}_{t}K\left[v_x^2+{\left(r{\omega }_r\right)}^2\right]\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill \varphi & =& v_x{\left(r{\omega }_r\right)}^{-1}\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill Q& =& a{v}_x\hfill \end{array}$$

(10)

where $dp$ represents the pressure drop; $C_a$ and $C_t$ represent the “power coefficients” and “torque coefficient”, respectively; $\varphi$ stands for the “flow coefficient”; $T_t$, K and r represent the turbine’s torque, constant and mean radius, respectively; l, b, and n represent the blade’s chord length, height and number, respectively; ${\omega }_r$ represents the “angular velocity”; a stands for the “cross-sectional area”; and $\rho$ represents the air density.

The characteristic curves of the Wells turbine under study are formed by the power coefficient $C_a$ and the torque coefficient $C_t$ versus the flow coefficient $\varphi$ as shown in Figure 3.

2.4. Doubly-Fed Induction Generator Model

In the OWC system under study, the Wells turbine drives a Doubly Fed Induction Generator to deliver electrical power to the grid. In a $dq$ diphase frame, the DFIG generator can be defined with the expression (11)–(16) given in [36,37]. Thus, the voltages of the stator and rotor in the $dq$ frame can be defined as:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle v_{ds}=R_s{i}_{ds}+\frac{d{\psi }_{ds}}{dt}-{\omega }_s{\psi }_{qs}}\\ {\displaystyle v_{qs}=R_s{i}_{qs}+\frac{d{\psi }_{qs}}{dt}+{\omega }_s{\psi }_{ds}}\end{array}\right.\end{array}$$

(11)

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle v_{dr}=R_r{i}_{dr}+\frac{d{\psi }_{dr}}{dt}-{\omega }_r{\psi }_{qr}}\\ {\displaystyle v_{qr}=R_r{i}_{qr}+\frac{d{\psi }_{qr}}{dt}+{\omega }_r{\psi }_{dr}}\end{array}\right.\end{array}$$

(12)

where $R_s$ and $R_r$ represent the stator and rotor resistances, ${\omega }_s$ and ${\omega }_r$ represent the stator and rotor angular velocity, and $i_{ds}$, $i_{qs}$ $i_{dr}$ and $i_{qr}$ represent the d-q stator and rotor currents.

The flux linkage at the stator and the rotor can be described by:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle {\psi }_{ds}=L_{ss}{i}_{ds}+L_m{i}_{dr}}\\ {\displaystyle {\psi }_{qs}=L_{ss}{i}_{qs}+L_m{i}_{qr}}\end{array}\right.\end{array}$$

(13)

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle {\psi }_{dr}=L_{rr}{i}_{dr}+L_m{i}_{ds}}\\ {\displaystyle {\psi }_{qr}=L_{rr}{i}_{qr}+L_m{i}_{qs}}\end{array}\right.\end{array}$$

(14)

where $L_{ss}$, $L_{rr}$ and $L_m$ represent the stator, rotor and magnetizing inductances, respectively.

The generated electromagnetic torque and its interaction with the turbine may be expressed as:

$${\displaystyle T_e=\frac{3}{2}p\left({\psi }_{ds}{i}_{qs}-{\psi }_{qs}{i}_{ds}\right)}$$

(15)

$$\frac{J}{p}\frac{d{\omega }_r}{dt}=T_e-T_t$$

(16)

where p represents the pair pole number and J represents the inertia of the system.

3. Control Statement

3.1. Stalling Behavior of the Wells Turbine

The stalling behavior in Wells turbines is a phenomenon that restricts the produced power. It happens in the event that the airflow speed $v_x$ rises; however, the rotational velocity ${\omega }_r$ is slow because the generator is unable to spin quick enough to match the incoming airflow of strong waves. This behavior is visible in Figure 3b which demonstrates when the flow coefficient $\varphi$ surpasses a critical value $0.3$, the torque coefficient $C_t$ declines considerably because the rotational speed ${\omega }_r$ is unable to match the airflow velocity $v_x$.

The stall effect is explained by operating the uncontrolled OWC plant with different sea states. The first sea condition examines waves with a $10\text{-}\mathrm{s}$ period and $0.8\text{-}\mathrm{m}$ wave amplitude (Figure 4a,b). The second sea condition examines waves with a $10\text{-}\mathrm{s}$ period and a $1.3\text{-}\mathrm{m}$ wave amplitude (Figure 4c,d).

As shown in Figure 4, when the waves are low (i.e., $A=0.8 \mathrm{m}$) the Wells turbine will have a low flow coefficient, which in this case does not exceed the threshold value $0.3$ (see Figure 4a). Hence, the resulting turbine torque is not affected by the stalling behavior (see Figure 4b). However, when the waves are high (i.e., $A=1.3 \mathrm{m}$) the Wells turbine will have a higher flow coefficient that exceeds the threshold value $0.3$ (see Figure 4c). Hence, the resulting turbine torque is affected by the stalling behavior (see Figure 4d).

The Wells turbine’s stalling behavior can be evaded if the flow coefficient is constantly regulated [24,25]. From the expression (9), the flow coefficient relies on the airflow velocity in the turbine duct. Thus, adjusting the airflow speed $v_x$ will aid in evading the stall effect; therefore, an airflow control strategy has been suggested.

The implementation of the airflow control puts to use the air valve set to use within the capture chamber, this device can be used to vary the pressure and airflow in the OWC system. The actuator of the air valve is controlled using a PID controller, as explained by the scheme of Figure 5.

Tuning the PID controller in a complex system such as the OWC often is hard and tedious when using conventional methods and lacking an appropriate systematic design approach. In order to tune the PID controller, the use of optimization theory has been suggested as a promising recourse to easily calculate and optimize all PID gains [21,22].

3.2. Airflow Control Problem Formulation

The PID tuning optimization problem for the airflow control scheme’s objective is to compute the best control design parameters ${\mathit{X}}_B=\left(x_B\left(1\right),x_B\left(2\right),x_B\left(3\right)\right)\in {\mathbb{R}}^3$ that represents the PID controller gains, i.e., $\mathit{X}=\left(K_p,K_i,K_d\right)\in {\mathbb{R}}_{+}^3$. This is achieved while minimizing the cost function. The Integral of Absolute Error (IAE) has been adapted as the cost function for this problem [23,24]:

$$f_{IAE}\left(\mathit{X}\right)={\int }_0^{+\infty }\left|e\left(\mathit{X},t\right)\right|dt,$$

(17)

where $e\left(.\right)$ is the error between the reference and the controlled variable, as detailed in Figure 6.

The IAE cost function is minimized by considering some time-domain constraints, associated to the rise and settling times ($t_r$ and $t_s$), the steady-state error $E_{ss}$, and the overshoot $\delta \left(\%\right)$ criteria of the closed-loop step response [23,24].

The PID tuning problem formulation for the airflow control has been formulated as a constrained and nonlinear optimization problem with expression (18). The tuning problem can be solved by the PSO algorithm [23,24]:

$$\left\{\begin{array}{c}minimizef\left(\mathit{X}\right)\hfill \\ \mathit{X}=\left(K_p,K_i,K_d\right)\in \mathcal{S}\subseteq {\mathbb{R}}_{+}^3\hfill \\ \mathrm{subject}\mathrm{to}:\hfill \\ g_1\left(\mathit{X}\right)=E_{ss}-E_{ss}^{\max }\le 0\hfill \\ g_2\left(\mathit{X}\right)=t_r-t_r^{\max }\le 0\hfill \\ g_3\left(\mathit{X}\right)=t_s-t_s^{\max }\le 0\hfill \end{array}\right.$$

(18)

where $f:{\mathbb{R}}^3\to \mathbb{R}$ represents the cost function, $\mathcal{S}=\left\{\mathit{X}\in {\mathbb{R}}_{+}^3,LB\le \mathit{X}\le UB\right\}$ stands for the bounded search space for the control variables, and $g_i:{\mathbb{R}}^3\to \mathbb{R},\left(i=1,2,3\right)$ represents the constraints.

4. Harmony Search Algorithm and Its Variants

Before introducing the Self-Adaptive Global-Best Harmony Search (SGHS) algorithm, we will briefly introduce the developed variant algorithms leading up to it.

4.1. Harmony Search Algorithm

The Harmony Search (HS) algorithm was first introduced in early 2000 by Z.W. Geem et al. in [38]. The HS algorithm is inspired by the musical process of musicians in the search for a fantastic harmony by aesthetic estimation. The variables $x\left(j\right)$ are represented by musical instruments and the fantastic harmony is the desired optimal solution, where $x\left(j\right)\in \left[LB\left(j\right),UB\left(j\right)\right]$, j = $1,2,\dots ,n$ and n is the number of variables. Every musical practice represents another iteration limited to a maximum Number of Improvisations ($NI$) and the quality of the results are evaluated based on the pitches of the instruments.

First an initial group of harmony vectors ${\mathit{X}}_i$ = $\left\{x_i\left(1\right),x_i\left(2\right),\dots ,x_i\left(n\right)\right\}$ are randomly generated as:

$$x_i\left(j\right)=\left(UB\left(j\right)-LB\left(j\right)\right)\times r+LB\left(j\right)$$

(19)

where $UB\left(j\right)$ and $LB\left(j\right)$ are the upper and lower bonds for the $j^{\mathrm{th}}$ variable and r is a random uniform number between 0 and 1.

The harmony vectors are then combined to form the Harmony Memory (HM) with a total of $HMS$ vectors to store the best harmony improvisation vectors based on their cost function which is stored in HM as well as:

$$HM=\left[\begin{array}{c}{\mathit{X}}_1\\ {\mathit{X}}_2\\ ⋮\\ {\mathit{X}}_i\\ ⋮\\ {\mathit{X}}_{HMS}\end{array}\right]=\left[\begin{array}{cc}\left.\begin{array}{ccccc}{x}_1\left(1\right)& \dots & x_1\left(j\right)& \dots & x_1\left(n\right)\\ x_2\left(1\right)& \dots & x_2\left(j\right)& \dots & x_2\left(n\right)\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ x_i\left(1\right)& \dots & x_i\left(j\right)& \dots & x_i\left(n\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ x_{HMS}\left(1\right)& \dots & x_{HMS}\left(j\right)& \dots & x_{HMS}\left(n\right)\end{array}\right|& \begin{array}{c}f\left({\mathit{X}}_1\right)\\ f\left({\mathit{X}}_2\right)\\ ⋮\\ f\left({\mathit{X}}_i\right)\\ ⋮\\ f\left({\mathit{X}}_{HMS}\right)\end{array}\end{array}\right]$$

(20)

where $HMS$ is the Harmony Memory Size.

To improvise new harmony vectors ${\mathit{X}}_{new}$, three rules are considered; the memory consideration, the pitch adjustment, and the random selection. The new vector could be selected from the vectors of the HM by testing a random number $r_1$ for memory consideration based on the Harmony Memory Consideration Rate ($HMCR$), and further will be pitch adjusted using a predefined Bandwidth ($BW$) based on the Pitch Adjustment Rate ($PAR$), otherwise it is randomly selected as explained in Figure 7.

Once the new harmony vector ${\mathit{X}}_{new}$ is generated, the harmony memory HM will be updated based on the fitness of the new vector and the worst existing vector ${\mathit{X}}_w$. The new vector ${\mathit{X}}_{new}$ will then replace the worst vector ${\mathit{X}}_w$ in HM if its fitness value is better.

Finally, the HS algorithm can be summarized by the steps detailed in the pseudocode of Algorithm 1:

	Algorithm 1: Harmony Search Algorithm
1.	Define parameters of HS algorithm: n, $HMS$, $HMCR$, $PAR$, $BW$, $NI$.
2.	Initialize HM with random harmony vectors using (19) and calculate the cost function of each vector.
3.	Improvise a new Harmony vector ${\mathit{X}}_{new}$ based on the three rules of Figure 7.
4.	Update the HM with ${\mathit{X}}_{new}$ if $\left(f\left({\mathit{X}}_{new}\right)<f\left({\mathit{X}}_w\right)\right)$ as ${\mathit{X}}_w$ = ${\mathit{X}}_{new}$.
5.	If the maximum number of improvisations $NI$ is reached then stop the program and return the best
	harmony vector ${\mathit{X}}_B$. Otherwise, go back to step 3.

$HMCR$ maintains the balance between the exploration and exploitation; on the other hand, $PAR$ is responsible for the refinement of the solutions by a distance BW. Therefore, the setting of these three parameters greatly influences the efficiency of the algorithm.

4.2. Improved Harmony Search Algorithm

The Improved Harmony Search (IHS) algorithm was introduced in 2007 by M. Mahdavi et al. in [39]. IHS was developed in an effort to improve the convergence of the solutions by dynamically varying $PAR$ and $BW$ at every new improvisation k as explained by Figure 8.

The probability of pith adjustment is increased by increasing $PAR$ linearly as:

$$PAR\left(k\right)=PA{R}_{min}+\left(\frac{PA{R}_{max}-PA{R}_{min}}{NI}\times k\right)$$

(21)

where k is the iteration or the current improvisation and $PA{R}_{max}$ and $PA{R}_{min}$ are the maximum and minimum adjustment rate, respectively.

At the same time, the degree of adjustment is decreased by decreasing $BW$ exponentially as:

$$BW\left(k\right)=B{W}_{max}exp\left(\frac{1}{NI}\times ln\left(\frac{B{W}_{min}}{B{W}_{max}}\right)\times k\right)$$

(22)

where $B{W}_{max}$ and $B{W}_{min}$ are the maximum and minimum bandwidth, respectively.

Similarly to the HS algorithm, the IHS algorithm follows the same steps with the inclusion of the calculation of $PAR\left(k\right)$ and $BW\left(k\right)$ as detailed in pseudocode of Algorithm 2:

	Algorithm 2: Improved Harmony Search Algorithm
1.	Define parameters of IHS algorithm: n, $HMS$, $HMCR$, $PA{R}_{max}$, $PA{R}_{min}$, $B{W}_{max}$, $B{W}_{min}$, $NI$.
2.	Initialize HM with random harmony vectors using (19) and calculate the cost function of each vector.
3.	Calculate the $k^{\mathrm{th}}$ pitch adjustment rate $PAR\left(k\right)$ and bandwidth distance $BW\left(k\right)$.
4.	Improvise a new Harmony vector ${\mathit{X}}_{new}$ based on the three rules of Figure 8.
5.	Update the HM with ${\mathit{X}}_{new}$ if $\left(f\left({\mathit{X}}_{new}\right)<f\left({\mathit{X}}_w\right)\right)$ as ${\mathit{X}}_w$ = ${\mathit{X}}_{new}$.
6.	If the maximum number of improvisations $NI$ is reached then stop the program and return the best
	harmony vector ${\mathit{X}}_B$. Otherwise, go back to step 3.

This method strongly enhances the capabilities of the IHS algorithm in terms of precision once the algorithm has converged to an interesting region of the search space. However, this method introduces the problem of bounds selection (i.e., $PA{R}_{max}$, $PA{R}_{min}$, $B{W}_{max}$ and $B{W}_{min}$). Moreover, the fact that the parameter $PAR$ continues to increase without settling even when an interesting region is reached made this method questionable.

4.3. Global-Best Harmony Search Algorithm

The Global-Best Harmony Search algorithm (GHS) has been developed and introduced in 2008 by Omran and Mahdavi in [40] based on the concept of the Particle Swarm Optimization algorithm [41]. The GHS has been proposed to deal with the limitations of the HS algorithm as a neighborhood metaheuristic which does not make use of its own past experience. The idea is to directly consider the best harmony vector in HM and simplify the pitch adjustment process as explained in Figure 9.

The new pitch adjustment rule randomly selects the $l^{\mathrm{th}}$ element $x_B\left(l\right)$ from the best harmony vector ${\mathit{X}}_B$ to the $j^{\mathrm{th}}$ decision variable $x_{new}\left(j\right)$ in the new harmony vector ${\mathit{X}}_{new}$.

The GHS algorithm can be summarized by the steps of Algorithm 3:

	Algorithm 3: Global-Best Harmony Search Algorithm
1.	Define parameters of GHS algorithm: n, $HMS$, $HMCR$, $PAR$, $NI$.
2.	Initialize HM with random harmony vectors using (19) and calculate the cost function of each vector.
3.	Improvise a new Harmony vector ${\mathit{X}}_{new}$ based on the three rules of Figure 9.
4.	Update the HM with ${\mathit{X}}_{new}$ if $\left(f\left({\mathit{X}}_{new}\right)<f\left({\mathit{X}}_w\right)\right)$ as ${\mathit{X}}_w$ = ${\mathit{X}}_{new}$.
5.	If the maximum number of improvisations $NI$ is reached then stop the program and return the best
	harmony vector ${\mathit{X}}_B$. Otherwise, go back to step 3.

This method eliminates the problem of $BW$ selection unlike in previous variants.

5. Self-Adaptive Global-Best Harmony Search Algorithm

The Self-Adaptive Global-Best Harmony Search (SGHS) algorithm has been developed and introduced in 2010 by Q.K. Pan et al. in [42]. The SGHS algorithm is based on the GHS algorithm with a few modifications in an effort to enhance its capabilities. The SGHS presents three major changes to the algorithm; (1) self-adaptation of $HMCR$ and $PAR$, (2) dynamic evolution of $BW$ and (3) new improvisation scheme.

5.1. Self-Adaptation of $HMCR$ and $PAR$

First of all, the parameters $HMCR$ and $PAR$ are no longer fixed values and are instead learned to adapt to the problem and the evolution of the search process. In fact, a very large $HMCR$ favors local search which increases the convergence rate, whereas a small value will favor the exploration which will diversify the HM. On the other hand, if $PAR$ is big, it favors the exploitation of the best harmony vector ${\mathit{X}}_B$ and passing it’s information to the next improvisation; but a small value will favor the perturbation of the values in HM to diversify, hence increasing the exploration. Therefore, Q.K. Pan et al. considered $HMCR$ and $PAR$ as normally distributed values in the interval $[0.9,1.0]\left(\right[0.0,1.0\left]\right)$ with mean values $HMCRm\left(PARm\right)$ and a Standard Deviation of ($SD$) $0.01\left(0.05\right)$ [42]. $HMCR$ and $PAR$ are then recalculated for every certain Learning Period of ($LP$) to adapt them at every phase of the search process.

5.2. Dynamic Evolution of $BW$

Secondly, the bandwidth distance $BW$ parameter was initially kept fixed in HS algorithm which resulted in neighboring solutions to the actual optimal solution. This was than dynamically increased right up to the end of the search in the IHS algorithm without settling down once the algorithm reaches an interesting region that could diverge the search.

For the SGHS algorithm $BW$ should be dynamically varying at the beginning but settles down in mid search to favor local search once an interesting region is reached. Therefore $BW$ is defined as:

$$BW\left(k\right)=\left\{\begin{array}{cc}B{W}_{max}-\frac{B{W}_{max}-B{W}_{min}}{NI}\times 2k\hfill & ift<NI/\phantom{NI2}2\\ B{W}_{min}\hfill & ift\ge NI/\phantom{NI2}2\end{array}\right.$$

(23)

where $B{W}_{\min }$ and $B{W}_{\max }$ are the minimum and maximum bandwidth distances.

5.3. New Improvisation Process

The improvisation process of the SGHS algorithm reinstates the use of bandwidth adjustment $BW$ but in the memory consideration rule. Moreover, the pitch adjustment rule is modified to assign the elements $x_B\left(j\right)$ of the best harmony vector ${\mathit{X}}_B$ in the HM to the corresponding decision variable $x_{new}\left(j\right)$ of the new harmony vector ${\mathit{X}}_{new}$, unlike the GHS algorithm, which assigns them randomly. This scheme ensures the use of the features of the best vector and offers the possibility to refine the solutions through the bandwidth parameter, as explained by the scheme of Figure 10.

5.4. SGHS Computational Procedure

The computational process SGHS algorithm combines all three new changes in the algorithm which is summarized by the pseudocode of Algorithm 4.

	Algorithm 4: Self-Adaptive Global-Best Harmony Search Algorithm
1.	Define parameters of SGHS algorithm: n, $HMS$, $HMC{R}_{max}$, $HMC{R}_{min}$, $PA{R}_{max}$, $PA{R}_{min}$, $SD$, $LP$, $NI$.
2.	Initialize $B{W}_{max}$, $B{W}_{min}$, $HMC{R}_m$ and $PA{R}_m$. Set learning counter $lp$ to 1.
3.	Initialize HM with random harmony vectors using (19) and calculate the cost function of each vector.
4.	Self-adapt $HMCR$ and $PAR$ according to $HMC{R}_m$ and $PA{R}_m$. Calculate $BW\left(k\right)$ from $B{W}_{max}$, $B{W}_{min}$.
5.	Improvise a new Harmony vector ${\mathit{X}}_{new}$ based on the three rules of Figure 10.
6.	Update the HM with ${\mathit{X}}_{new}$ if $\left(f\left({\mathit{X}}_{new}\right)<f\left({\mathit{X}}_w\right)\right)$ as ${\mathit{X}}_w$ = ${\mathit{X}}_{new}$ and record $HMCR$ and $PAR$.
7.	If $\left(lp=LP\right)$ then recalculate $HMC{R}_m$ and $PA{R}_m$ and reset $lp$ to 1. Otherwise, increment $lp$ by 1.
8.	If the maximum number of improvisations $NI$ is reached, then stop the program and return the best
	harmony vector ${\mathit{X}}_B$. Otherwise, go back to step 4.

6. Results and Discussion

The performance evaluation of the suggested optimization for the airflow control in the OWC has been carried out by numerical simulations using a numerical wave-to-wire model on Matlab/Simulink. The OWC wave-to-wire model is configured using the parameters of NEREIDA detailed in

Table 1. OWC parameters from the NEREIDA wave power plant.

Capture Chamber	Wells Turbine	DFIG Generator
$w_c$ = 4.5 m	n = 5	$R_s$ = 0.5968 $\mathsf{\Omega }$	$P_{rated}$ = 18.45 kW
$l_c$ = 4.3 m	b = 0.21 m	$R_r$ = 0.6258 $\mathsf{\Omega }$	$V{s}_{rated}$ = 400 V
${\rho }_a$ = 1.19 kg/m${}^3$	l = 0.165 m	$L_{ss}$ = 0.0003495 H	$f_{rated}$ = 50 Hz
${\rho }_w$ = 1029 kg/m${}^3$	r = 0.375 m	$L_{rr}$ = 0.324 H	$p=2$
	a = 0.4417 m${}^2$	$L_m$ = 0.324 H

.

6.1. Optimization and Computational Results

Due to the stochastic and irreproducible nature of optimization algorithms, validating their performance is supposed to be via statistical analysis on the goodness of the found solutions of several trials. Thus, the suggested particle swarm optimization algorithm has been simulated 20 times with a maximum Number of Improvisations (iterations) $NI$ = 100 and a Harmony Memory Size $HMS$ = 20 while running on an Intel Core i5, 3.30 GHz CPU. Feasible solutions have been obtained in 80% of trials and in acceptable CPU calculation time.

Figure 11 illustrates the box-and-whisker plot of the results of the optimization of all four variants of the Harmony Search algorithm (i.e., HS, IHS, GHS, SGHS) for Problem (18).

The figure shows that the obtained solutions from all four algorithms are in the same region of the search space. From a statistical point of view, we focus on the average values which are 4.015, 3.575, 3.875, and 3.450 for HS, HIS, GHS and SGHS, respectively. However, in general, it is obvious that, in terms of average value (red line) and of minimum value (bottom whisker), among 20 trials of every algorithm the SGHS algorithm presents superiority over the previous variants. However, it is to be noted that the box of the IHS algorithm is the narrowest which is thanks to its pitch adjustment rules that favor exploitation and precision. On the other hand, the GHS show less favorable results with wider box, meaning dispersed solutions, this may be due to the lack of solution refinement via the $BW$ parameter.

The significant outcomes of the 20 trials using all four algorithms are detailed in

Table 2. Details of optimization results from 20 trials of Problem (18).

Algorithm	Min	Max	Median	Mean	Standard Deviation
HS	3.900	4.900	4.015	4.400	0.530
IHS	3.100	4.500	3.575	3.800	0.540
GHS	3.402	4.950	3.875	4.176	0.825
SGHS	2.874	4.300	3.450	3.587	0.697

. The SGHS algorithm is the best in terms of average value but not significantly better. SGHS has the lowest minimum, maximum, median and mean values and the HS algorithm has the lowest standard deviation.

To further understand the behavior of the algorithms Figure 12 illustrates the most typical convergence curves of HS, IHS, GHS and SGHS along with the curves of previously tested algorithms—the Particle Swarm Optimization with decreasing inertia (PSO-In), the Fractional-Order Particle Swarm Optimization Memetic Algorithm (FPSOMA) and the Water Cycle Algorithm (WCA).

From the convergence histories, it can be noticed that all algorithms successfully converge to the same region of the search space (between 3 and 4) but the optimal is from the SGHS algorithm. The other variants of the HS algorithm managed to enhance the exploration and exploitation capabilities but with different outcomes. In the case of the IHS algorithm, the precision has been improved thanks to the increase of the exploitation with the new pitch adjustment rules (21) and (22) but because $PAR$ continues to increase until the end, somehow IHS diverged from the optimal region. In the case of the GHS algorithm, it converged quicker to the best solution thanks to the use of the best vector in the improvisation process, however, the GHS failed to further converge because of the lack of refinement rules using $BW$ this behavior is known as the premature convergence. Finally, in the case of the SGHS algorithm, it can be noticed that a gradual convergence took place thanks to the self-adaptation of $HMCR$ and $PAR$ every $LP=10$ improvisations. Moreover, the reinstatement of the refinement rules using the bandwidth parameter $BW$ allowed SGHS to obtain more precise solutions than those of the IHS algorithm. The Water Cycle Algorithm presents close results to the HS algorithm, but the new variants of the HS algorithm are better in comparison to other previous algorithms tested on the OWC system.

The parameters obtained from the mean case of feasible results from 20 trials of optimization using the HS, IHS, GHS and SGHS algorithms are given in

Table 3. PID parameters obtained by optimization algorithms.

Algorithm	Proportional Gain (${\mathit{K}}_{\mathit{p}}$)	Integral Gain (${\mathit{K}}_{\mathit{i}}$)	Derivative Gain (${\mathit{K}}_{\mathit{d}}$)
HS	2088	103.50	77.33
IHS	1806	112.87	66.88
GHS	2115	132.18	78.33
SGHS	2364	147.75	87.55

. The $K_p$, $K_i$ and $K_d$ values are within a close range, proving that the algorithms converge to the same interesting search region.

6.2. OWC Performance Evaluation

For the assessment of the performance of the suggested airflow-control on the OWC system, the optimized PID controller using the parameters found of the mean case of the SGHS optimization results. The evaluation compares the performance of the uncontrolled OWC, the OWC using a traditionally tuned PID using the well known Ziegler–Nichols method (ZN-PID) and the OWC using the optimized PID with SGHS algorithm (SGHS-PID) with the parameters given in Table 3.

6.2.1. Performance With Regular Waves

During the simulations, two regular wave conditions were considered the first wave is weak with a wave period T = $10 \mathrm{s}$ and a wave amplitude A = $0.8 \mathrm{m}$ from $0 \mathrm{s}$ to $22.5 \mathrm{s}$. From $22.5 \mathrm{s}$ to $50 \mathrm{s}$, we considered a stronger wave with a wave period T = $10 \mathrm{s}$ and a wave amplitude A = $1.3 \mathrm{m}$ as shown in Figure 13.

The flow coefficients of the OWC in the uncontrolled case and with the airflow control using the traditionally tuned PID (ZN-PID) and the optimized PID (SGHS-PID) are illustrated in Figure 14.

The figure shows that in the uncontrolled case the flow coefficient surpasses the threshold value 0.3 which will provoke the stall effect, whereas in both controlled cases the flow coefficients were regulated below 0.3 thanks to both ZN-PID and SGHS-PID controllers. However, when zooming to the curves it is observed that the SGHS-PID manages to provide a closer flow coefficient to the threshold value than that of the ZN-PID.

Figure 15 shows that in the uncontrolled case, the airflow speed continues to increase; however, in both controlled cases, the airflow speed decreased thanks to both ZN-PID and SGHS-PID controllers, but with a slight superiority for the SGHS-PID.

The obtained turbine torques of the PTO in the uncontrolled case, ZN-PID controller and SGHS-PID controller are shown in Figure 16a, Figure 16c,e, respectively. Furthermore, the generated powers for the uncontrolled case, ZN-PID controller and SGHS-PID controller are shown in Figure 16b, Figure 16d,f, respectively.

It may be observed that during the uncontrolled case the torque has been affected by the stall effect which reduced it in terms of average value to $61.16$$N.m$. On the other hand, the airflow control manages to avoid the stall effect and increase the torques in terms of average values to $64.22$$N.m$ and $66.71$$N.m$ for ZN-PID and SGHS-PID, respectively. As a result from the obtained torques, the generated powers in the uncontrolled case is the lowest with $-15.52$ $\mathrm{kW}$ and the highest power is generated when using the SGHS-PID with $-17.38$ $\mathrm{kW}$ followed by $18.55$ $\mathrm{kW}$ when using a ZN-PID in the proposed airflow control scheme.

6.2.2. Performance with Real Wave Data

For this study case, real surface elevation measurements of waves in Mutriku obtained by the Acoustic Doppler Current Profiler on 12 May 2014, from 00:00:00 a.m. to 00:00:50 a.m., shown in Figure 17, were considered.

The flow coefficients of the OWC with real measured wave data input in the uncontrolled case, and with airflow control using the traditionally tuned PID (ZN-PID) and the optimized PID (SGHS-PID), are illustrated in Figure 18.

Figure 18 shows that in the uncontrolled case, the flow coefficient exceeds the threshold value 0.3, which will provoke the stall effect in two time periods with strong waves. The first wave train is from 2 s to 11 s and the second strong wave train is from 37 s to 50 s. In both controlled cases, the flow coefficients were regulated below 0.3 thanks to both ZN-PID and SGHS-PID controllers, with a slight superiority for SGHS-PID.

Figure 19 shows that in the uncontrolled case, the airflow speed continues to increase uncontrollably in both periods of the strong waves; however, in both controlled cases, the airflow speed was decreased thanks to both ZN-PID and SGHS-PID controllers, but with a slight superiority for SGHS-PID.

The obtained turbine torques of the PTO in the uncontrolled case, ZN-PID controller and SGHS-PID controller are shown in Figure 20a,c,e, respectively. And the generated powers for the uncontrolled case, ZN-PID controller and SGHS-PID controller are shown in Figure 20b, Figure 20d,f, respectively.

As with regular waves, it may be observed that during the uncontrolled case, the torque has been affected by the stall effect, which reduces it in terms of average value. On the other hand, the airflow control manages to avoid the stall effect and increase the torques in terms of average values in both high wave regions. The resulting generated powers in the uncontrolled case is the lower than the powers generated when using the SGHS-PID and the ZN-PID in the proposed airflow control scheme.

7. Conclusions

The paper proposes an airflow control strategy to deal with the undesired stalling phenomenon of the Wells-turbine-based oscillating water columns. The proposed strategy takes into consideration the aerodynamic characteristics of the studied Wells turbine and implements a control using a PID controller to govern an air valve situated within the capture chamber of the OWC system.

To efficiently tune the parameters of the PID with such a complex system, a recourse to optimization theory has been proposed, specifically by using developed variants of the Harmony Search algorithm. The algorithms were implemented and tested to get to the final variant (i.e., SGHS). The parameters obtained by the SGHS algorithm were then used in the PID controller to simulate the OWC system with the proposed airflow control scheme.

Finally, in order to show the effectiveness of the proposed strategy and tuning technique, two comparative study cases have been carried out between the uncontrolled OWC and the controlled OWC using the airflow control with a Zeigler–Nichols-tuned PID (ZN-PID) and a SGHS algorithm-tuned PID (SGHS-PID). The first case study considers two different sea conditions of regular waves, whereas the second case study considers real measured wave data. The results of the studies demonstrate successful avoidance of the stall effect with both airflow control cases. However, the control using the SGHS algorithm shows superior performance to the traditional tuning method.