Symmetry-Breaking for Airflow Control Optimization of an Oscillating-Water-Column System

Global optimization problems are mostly solved using search methods. Therefore, decreasing the search space can increase the efficiency of their solving. A widely exploited technique to reduce the search space is symmetry-breaking, which helps impose constraints on breaking existing symmetries. The present article deals with the airflow control optimization problem in an oscillating-water-column using the Particle Swarm Optimization (PSO). In an effort to ameliorate the efficiency of the PSO search, a symmetry-breaking technique has been implemented. The results of optimization showed that shrinking the search space helped to reduce the search time and ameliorate the efficiency of the PSO algorithm.


Introduction
Natural resources are the best alternative to limited fossil fuel for the European energy market, especially that EU countries consume one-fifth of the global oil and gas supplies. This leads to the annual purchase of 350 billion Euros worth of petroleum and gas products from other countries [1]. Nowadays, an upward trend to develop renewable energy projects in European countries has been at the forefront of European policies in order to reduce fossil fuel consumption [2]. The conclusion of the Conference of Parties (COP21) has put in place further motivated goals for a 2 • C decrease [3], hence the need for Renewables to take a big part in the future of energy systems and markets. This calls for the importance of a stable and reliable energy mix development based on all renewable resources [4]. In fact, renewable energy sources are estimated to have provided 33% of total power generation in the UK and 40% in Germany and Spain in 2018 [5].
Even though Ocean energy is more reliable, predictable, and denser than other Renewables, Marine Renewable Energy (MRE) still remains underdeveloped and the least tapped energy industry [6]. The majority of MRE-based industries are up to the present time in the early-phase of development varying between theory and design, until the demonstration phase [1]. Although decades of research and development were invested, a mere 529 MW of MRE capacity has been recorded as operational at the end of 2017 and more than 93% of them are summarized in two main facilities [7]: Sihwa plant in South Korea with 254 MW [8] and La Rance tidal power station in France with 240 MW [9]. The main obstacles hindering the development of MRE are geographical installation that can affect nearby coastal and fishing areas and maritime routes, required and existing electricity infrastructure, potential environmental impacts, and legal and regulatory framework [7].

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

Wave Surface Dynamics
A monochromatic unidirectional wave has been considered to be input into the implemented numerical model of the OWC system. There exist 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 [46,47]. 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 neglecting turbulence, friction losses, and other energy losses [47].
The parameters of a wave are detailed in Figure 1, with SW L represents the "Still-Water-Level", 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 λ represents the interval between successive crests known as "wavelength" [47,48]. Therefore, the surface elevation for a sea wave is given as [49,50]: 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 [49]: k tanh(kh) = ω 2 /g (2) where g represents the acceleration gravity.

Capture Chamber Model
The volume of the air within the Oscillating Water Column's chamber is defined in [38,48] 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 [38,48]: with c = w/k. Once the chamber's geometry has been taken into account with Equation (4), the airflow velocity can be described as [38,48]: 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, is a self-rectifying axial-flow air turbine [50,51]. Self-rectifying air-turbines possess blades with special geometry allowing a unidirectional rotating motion regardless of the airflow direction [52][53][54]. The Wells turbine under study can be mathematically defined by Expressions (6)-(10) given in [18,55]: where dp represents the pressure drop, C a and C t represent the "power coefficients" and "torque coefficient", φ stands for the "flow coefficient", T t , K, and r represent the turbine's torque, constant, and mean radius, l, b, and n represent the blade's chord length, height, and number, ω 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 φ are 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 Expressions (11)- (16) given in [56,57]. Thus, the voltages of the stator and rotor in the dq frame can be defined as: with R s , R r represent the stator and rotor resistances, ω s , ω r represent the stator and rotor angular velocity, 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) with L ss , L rr , and L m representing the stator, rotor, and magnetizing inductances, respectively. The generated electromagnetic torque and its interaction with the turbine may be expressed as: with p representing the pair pole number and J representing the inertia of the system.

Stalling Behavior Problem
The stalling behavior in Wells turbines is a phenomenon which 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 that 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 flow coefficient φ defined in Section 2 (9) relies on the airflow velocity v x . In addition, as expressed in Section 2 (5), v x relies on both the wave amplitude and period. Indeed, the bigger the wave amplitude and the smaller the period are, the higher the pneumatic power and the airflow velocity are as shown in Figure 4.  The stall effect is demonstrated by investigating the uncontrolled OWC system with different sea states. The first sea condition is waves with a 10 s period and 0.9 m wave amplitude beginning at 0 s until 22.5 s. The second sea condition is waves with a 10 s period and a 1.2 m wave amplitude beginning at 22.5 s until 50 s. It may be observed in Figure 5 that during the first case the OWC offers a flow coefficient below the threshold value 0.3, but, during the second case, it surpasses 0.3. As a result, the attained turbine torque illustrated in Figure 6 shows no stalling during the first case and a significant decline at the crest of the second case because of the stall effect which reduces the obtained torque in terms of average value. The Wells turbine's stalling behavior can be evaded if the flow coefficient is constantly regulated [39,45]. From Expression (9), the flow coefficient relies on the airflow velocity in the turbine duct. Thus, adjusting the airflow speed v x will aid with 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 within the capture chamber, and 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 7. Tuning the PID controller in a complex system such as the OWC often is hard and tedious when using conventional methods and lacks 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 [35][36][37].

Control Problem Formulation
The PID tuning optimization problem for the airflow control scheme's objective is to compute the control design parameters T ∈ R 3 that represents the PID controller gains, i.e., This is achieved while minimizing the cost function. The Integral of Absolute Error (IAE) has been adapted as the cost function for this problem [38,39]: with e (.) being the error between the reference and the controlled variable, as detailed in Figure 8. The IAE cost function is minimized by considering some time-domain constraints, associated with 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 [38,39].
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 [38,39]: Here f : R 3 → R represents the cost function, S = x x x ∈ R 3 + , x x x low ≤ x x x ≤ x x x upper stands for the bounded search space for the control variables, and g i : R 3 → R, (i = 1, 2, 3) represent the constraints.
Many methods were proposed solve constraints such as Problem (18). The most often method relies on the application of penalties to the cost function. Thus, the external static penalty method has been adopted and may be written as [36,37]: where Λ j represents the "scaling penalty parameters" and n con represents the "number of constraints".
In literature, for optimization problems, the scaling parameters Λ j , utilized by (19), are gradually increased with each iteration in a linear way for the aim of enforcing the constraints. Generally, the quality of the computed solutions depends on the values of the scaling parameters. However, in this study, in order to make the proposed technique simpler, big, and constant scaling penalty parameters have been taken into consideration [38,39]. From the second part of (19), the objective function will increase with static penalties. These penalties are proportional to the degree of infeasibility of the constraints g j [38,39].

Symmetry-Breaking Constraints
Symmetries in combinatorial optimization enlarge the dimensions of the search space; hence, extra time is invested to visit new solutions that are symmetric to the already visited solutions. The search time of a combinatorial problem "P" can be reduced by adding new set of constraints "C", referred to as Symmetry-Breaking Constraints (SBC), such that some of the symmetric solutions are eliminated from the search space while preserving at least one solution [58][59][60].
The method of symmetry-breaking constraints can be used to take advantage of symmetries in many constrainted optimization problems, by adding constraints that eliminate symmetries and reduce the search space size. This method uses three types of symmetries [59,60]: • Variable symmetry where permuting variables are solution invariant defined as: • Value symmetry where permuting solution values are solution invariant defined as: • Variable/value symmetry where both variables and values permutation is solution invariant defined as: In the problem of PID tuning optimization, only value symmetries are considered which are broken by imposing further constraints on the original problem defined by Equation (18). As described by Equation (21), a value symmetry σ is a bijection on values which conserves solutions. That is, is a solution as well. The lex leader approach may be applied to break value symmetries [60]. Let ∑ be the set of all value symmetry permutations, these symmetries are broken by imposing: where X 1 to X n is any fixed ordering on the variables, n is the number of variables, and σ is a symmetry.

Particle Swarm Optimization
The conventional PSO algorithm uses a group of individuals formed of n p particles, marked as x 1 , x 2 , . . . , x n p , randomly scattered within an initially bounded search space, in order to find a global solution for a generic optimization problem like that of Equation (18) [38]. Each particle in the swarm has a position x x x i At the k th iteration, the position of the i th particle, with x i ∈ R d , changes depending on the future position and velocity equations: with w standing for the inertia factor, c 1 and c 2 representing the cognitive and the social scaling factors, r i 1,k and r i 2,k representing two random numbers uniformly scatted within the range [[0, 1]], and p p p i k and p p p g k are the best positions previously obtained for the ith particle and the entire population at the kth iteration [38].
The exploration and exploitation capabilities of the particle swarm optimization algorithm may be supplementary enhanced when using a linear evolution mechanism of the inertia factor [38,39], all while considering the iterations as follows: with k max representing the maximum number of iterations and w max and w min representing the maximum and minimum inertia factor typically set to 0.9 and 0.4 [61].
In conclusion, a pseudocode of the PSO algorithm for a minimization problem has been detailed in Algorithm 1.

Assess the associated fitness values
where pbest i k represents the best previous fitness of the i th particle and gbest k represents the best fitness in the whole swarm. 6. If the stopping criterion is met, the program finishes and stops searching with the obtained solution x x x * = arg min 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 Results
Due to the stochastic and irreproducible nature of PSO algorithms, validating its 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 iterations k max =100 and a population size n p =30 while running on a CPU Core i5, 3.30 GHz. Feasible solutions have been obtained in 75% of trials and in acceptable CPU calculation time. The attained optimization results have been introduced in Table 2, and the results show that, with the SBCs, the PSO algorithm performs faster and better.  Figure 9 illustrates the box-and-whisker plot of the results of the optimization with and without SBCs for Problem (18). The figure shows that the median line is lower for the PSO algorithm with SBCs than that for the PSO algorithm without SBCs, which indicates that the use of SBCs allows for obtaining better cost function in terms of mean value. Moreover, the box of the algorithm with SBCs is narrower than that the one of the box of the PSO algorithm without SBCs, which means that the algorithm is able to converge to the same region of interesting search space when using SBCs. Both indicators increase the reproductivity of the SBC-based PSO algorithm leading to the optimal solution. The behavior and performance of the PSO algorithm with and without the SBCs is further observed in the curves of the convergence history shown in Figure 10.
The convergence histories of the best case from optimization trials of Problem (18) using the PSO with and without SBCs show that the PSO algorithm in both cases converges to the same region which proves that it is able to find an interesting search space to explore. In addition, the curves show that with SBCs the algorithm converges to the interesting region in fewer iterations (64 iterations), whereas without SBCs it takes more iterations to reach the interesting region (91 iterations).
From the convergence curves of Figure 10, it is clear that the PSO algorithm manages to converge with fewer iterations to the region of interesting search space when using SBC, which explains the reduced simulation time of Table 2. This is due to the fact that the added symmetry-breaking constraints help narrow the area of search space and eliminates the need to visit symmetric solutions that have already been visited with other sets of variables.

OWC Performance
For the assessment of the performance of the suggested airflow-control on the OWC system, the optimized PID controller uses the parameters found of the mean case of the PSO optimization results. The evaluation will compare 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 using PSO (PSO-PID).
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 (PSO-PID) are illustrated in Figure 11. 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 PSO-PID controllers. However, when zooming to the curves, it is observed that the PSO-PID manages to provide a slightly closer flow coefficient to the threshold value than that of the ZN-PID. The obtained turbine torques of the PTO are shown in Figure 12. 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. On the other hand, the airflow control manages to avoid the stall effect and increases the torques in terms of average values.
The outputs of the OWC system, which are the generated power, are presented in Figure 13. As a result from the obtained torques of Figure 12, the generated power in the uncontrolled case is the lowest and the highest power is generated when using the PSO-PID.

Conclusions
The paper deals with the stall effect in the OWC system which is known to occur in Wells turbines when a certain airflow speed threshold is exceeded. The research work proposes an airflow control strategy to govern the valve in the air-chamber and avoid the stalling behavior. The PID controller is configured thanks to the Particle Swarm Optimization algorithm.
The PSO algorithm has been simulated to solve the airflow control problem and search for suitable control parameters, but, to further enhance its exploration capabilities, the use of symmetry-breaking constraints was introduced to shrink the search space by eliminating the symmetric solutions. The optimization and computational results indicate amelioration of the PSO performance.
The obtained control parameters from the PSO algorithm were used in the performance evaluation of the suggested airflow control. The operation of the uncontrolled OWC, the OWC using a traditionally tuned PID with the well known Ziegler-Nichols method (ZN-PID), and the OWC using the optimized PID with PSO algorithm (PSO-PID) were compared with the same sea conditions. The results show that the proposed airflow control scheme manages to successfully evade the stall effect in the Wells turbine. Moreover, the airflow control using the PSO-PID demonstrates a superior performance to that of the airflow control using the ZN-PID. by the MCIU/MINECO through RTI2018-094902-B-C21/RTI2018-094902-B-C22 (MCIU/AEI/FEDER, UE).

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

Nomenclature
The following 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 P Optimization problem x x x * Optimal solution for the problem f Cost function ϕ Penalty function Λ j Scaling penalty parameter for the jth constraint Best positions formerly found for the ith particle and the entire swarm at the kth iteration