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

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.


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

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 SW L represents the "Still-Water-Level" and h, called "sea depth", represents the interval from the sea floor to SW L. H marks the interval from wave trough to wave crest called "wave height". A measures the distance between SW L and the wave crest known as "wave amplitude" and λ 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) = A sin (ωt − kxθ) = H/2 sin (ωt − kxθ) (1) where ω is the wave frequency, x is the wave horizontal coordinate, θ marks the the angular opening from the x-axis to waves' direction and k represents the wave number linked to ω with relation (2) as described in [29]: k tanh(kh) = ω 2 /g, where g represents the acceleration gravity.

Capture Chamber Model
The volume of the air within the Oscillating Water Column's chamber is defined in [23,28] as: 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]: 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]: where T w is the wave period and D is the duct diameter.

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]: where dp represents the pressure drop; C a and C t represent the "power coefficients" and "torque coefficient", respectively; φ 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; ω r represents the "angular velocity"; a stands for the "cross-sectional area"; and ρ 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 φ as shown in Figure 3.

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: where R s and R r represent the stator and rotor resistances, ω s and ω 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: ψ ds = L ss i ds + L m i dr ψ qs = L ss i qs + L m i qr (13) ψ dr = L rr i dr + L m i ds ψ qr = L rr i qr + L m i qs (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: where p represents the pair pole number and J represents the inertia of the system.

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 ω 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 φ surpasses a critical value 0.3, the torque coefficient C t declines considerably because the rotational speed ω 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-s period and 0.8-m wave amplitude (Figure 4a,b). The second sea condition examines waves with a 10-s period and a 1.3-m wave amplitude (Figure 4c,d).
As shown in Figure 4, when the waves are low (i.e., A = 0.8 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 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].

Airflow Control Problem Formulation
The PID tuning optimization problem for the airflow control scheme's objective is to compute the best control design parameters X X X + . 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]: where e (.) is the error between the reference and the controlled variable, as detailed in Figure 6.

PID controller
Capture chamber

Kp Ki Kd
ThroƩle valve 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 δ (%) 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]: where f : R 3 → R represents the cost function, S = X X X ∈ R 3 + , LB ≤ X X X ≤ UB stands for the bounded search space for the control variables, and g i : R 3 → R, (i = 1, 2, 3) represents the constraints.

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.

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(j) are represented by musical instruments and the fantastic harmony is the desired optimal solution, where x(j) ∈ [LB(j), UB(j)], j = 1, 2, ..., n and n is the number of variables. Every musical practice represents another iteration limited to a maximum Number of Improvisations (N I) and the quality of the results are evaluated based on the pitches of the instruments.
First an initial group of harmony vectors X X X i = {x i (1), x i (2), . . . , x i (n)} are randomly generated as: where UB(j) and LB(j) are the upper and lower bonds for the j 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: where HMS is the Harmony Memory Size.
To improvise new harmony vectors X X 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 X X X new is generated, the harmony memory HM will be updated based on the fitness of the new vector and the worst existing vector X X X w . The new vector X X X new will then replace the worst vector X X 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: If the maximum number of improvisations N I is reached then stop the program and return the best harmony vector X X 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.

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: where k is the iteration or the current improvisation and PAR max and PAR min are the maximum and minimum adjustment rate, respectively. At the same time, the degree of adjustment is decreased by decreasing BW exponentially as: where BW max and BW 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(k) and BW(k) as detailed in pseudocode of Algorithm 2: Algorithm 2: Improved Harmony Search Algorithm 1. Define parameters of IHS algorithm: n, HMS, HMCR, PAR max , PAR min , BW max , BW min , N I. 2. Initialize HM with random harmony vectors using (19) and calculate the cost function of each vector. 3. Calculate the k th pitch adjustment rate PAR(k) and bandwidth distance BW(k). 4. Improvise a new Harmony vector X X X new based on the three rules of Figure 8. 5. Update the HM with X X X new if ( f (X X X new ) < f (X X X w )) as X X X w = X X X new . 6. If the maximum number of improvisations N I is reached then stop the program and return the best harmony vector X X 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., PAR max , PAR min , BW max and BW min ). Moreover, the fact that the parameter PAR continues to increase without settling even when an interesting region is reached made this method questionable.

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 th element x B (l) from the best harmony vector X X X B to the j th decision variable x new (j) in the new harmony vector X X 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, N I. 2. Initialize HM with random harmony vectors using (19) and calculate the cost function of each vector. 3. Improvise a new Harmony vector X X X new based on the three rules of Figure 9. 4. Update the HM with X X X new if ( f (X X X new ) < f (X X X w )) as X X X w = X X X new . 5. If the maximum number of improvisations N I is reached then stop the program and return the best harmony vector X X X B . Otherwise, go back to step 3.
This method eliminates the problem of BW selection unlike in previous variants.

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.

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 X X 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. with mean values HMCRm(PARm) and a Standard Deviation of (SD) 0.01(0.05) [42]. HMCR and PAR are then recalculated for every certain Learning Period of (LP) to adapt them at every phase of the search process.

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: where BW min and BW max are the minimum and maximum bandwidth distances.

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 (j) of the best harmony vector X X X B in the HM to the corresponding decision variable x new (j) of the new harmony vector X X 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.

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, HMCR max , HMCR min , PAR max , PAR min , SD, LP, N I. 2. Initialize BW max , BW min , HMCR m and PAR m . Set learning counter l p 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 HMCR m and PAR m . Calculate BW(k) from BW max , BW min . 5. Improvise a new Harmony vector X X X new based on the three rules of Figure 10. 6. Update the HM with X X X new if ( f (X X X new ) < f (X X X w )) as X X X w = X X X new and record HMCR and PAR. 7. If (l p = LP) then recalculate HMCR m and PAR m and reset l p to 1. Otherwise, increment l p by 1. 8. If the maximum number of improvisations N I is reached, then stop the program and return the best harmony vector X X X B . Otherwise, go back to step 4.

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.

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) N I = 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. 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. 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.

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.

Performance With Regular Waves
During the simulations, two regular wave conditions were considered the first wave is weak with a wave period T = 10 s and a wave amplitude A = 0.8 m from 0 s to 22.5 s. From 22.5 s to 50 s, we considered a stronger wave with a wave period T = 10 s and a wave amplitude A = 1.3 m as shown in Figure 13. 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.

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.  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.

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.

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

Nomenclature
The folowing symbols are used in this manuscript: Error between the reference variable and measured variable u Control signal obtained from the PID controller K p , K i , K d Proportional, integral and derivative gains of the PID controller x x x * Optimal solution for the problem f Cost function