Application of a particle swarm optimization for shape optimization in hydraulic machinery

A study of shape optimization has become increasingly popular in academia and industry. A typical problem is to find an optimal shape, which minimizes (or maximizes) a certain cost function and satisfies given constraints. Particle Swarm Optimization (PSO) has received a lot of attention in past years and is inspired by social behaviour of some animals such as flocking behaviour of birds. This paper focuses on a possibility of a diffuser shape optimization using particle swarm optimization (PSO), which is coupled with CFD simulation. Influence of main parameters of PSO-algorithm and later diffuser shapes obtained with this method are discussed and advantages/disadvantages summarized.


Introduction
A study of size and shape optimization is very popular in industry and academia.There are a variety of methods to get the desired size/shape.These methods include nowadays well-liked particle swarm optimization (PSO) method.
This paper builds upon knowledge of an article [1] and extends the possibilities of optimizing a hydraulic diffuser, which is located behind a runner of a swirl turbine and where a recovery of kinetic energy into pressure energy takes place (Figure 1).-metres) [1,13] Main goal in [1] and also in this paper is to modify the shape of the diffuser for the purpose of maximizing a coefficient of pressure recovery c p (1): where p 1 (Pa) is static pressure at inlet of the computational domain, p 2 (Pa) is static pressure at outlet of the domain; ȡ (kg/m 3 ) is density of water, v 1 (m/s) is bulk velocity at inlet of the domain.Note: The value of velocity at inlet is constant and does not match with the real velocity profile behind the runner of the swirl turbine [1].

Particle swarm optimization
Fig. 2. Particle swarm optimization -concept [14] The PSO algorithm (PSOA) [2] was introduced in 1995 by J. Kennedy and R. Eberhart and was discovered through a simulation of a simplified social model.This algorithm is based on a bird flocking, fish schooling and swarming theory in particular.
The PSO algorithm uses a swarm of particles (subjects/individuals) randomly distributed in computational area.Individual particles exchange information about their velocity (= step size), position and fitness (fitness = solution value of the examined function).Fig. 3. GBEST model (left) and LBEST model [7] Each particle in PSO algorithm is treated as a point in n-dimensional space.The ith particle is represented as ܺ ூ ൌ ሺ‫ݔ‬ ଵ ǡ ‫ݔ‬ ଶ ǡ ǥ ǡ ‫ݔ‬ ሻ.The best previous position of any particle is recorded and represented as ܲ ூ ൌ ሺ‫‬ ଵ ǡ ‫‬ ଶ ǡ ǥ ǡ ‫‬ ሻ.The index of the best particle in the swarm is represented by index g, and velocity for ith particle is ܸ ூ ൌ ሺ‫ݒ‬ ଵ ǡ ‫ݒ‬ ଶ ǡ ǥ ǡ ‫ݒ‬ ሻ [3,4].
In GBEST model (Figure 3) all particles are directly connected to the best solution (of the examined function) in the swarm, but on the other hand in LBEST model (Figure 3) each particle is influenced by a smaller number of neighbouring particles of the swarm array.Typically LBEST neighbourhoods consist of two neighbours, one on each side -ring lattice [7].

Global Best PSO
The particles of the swarm are moving according to the following equation [3,4]: where c 1 and c 2 are positive constants; rand() and Rand() are two random vectors from range (0, 1); p gn is the best value of the examined function obtained by the whole swarm (GBEST).
The positions of particles are afterwards computed as follows [3,4]:

Local Best PSO
In this concept, velocity is computed similar as in Global Best PSO [3,4]: only difference is p ln , which is the best value of the examined function in the particle neighbourhood (LBEST).
Note: Random vectors rand() and Rand() were later in this study (in MATLAB code) created by function rand.

Improvements and modifications of Global Best PSO algorithm (PSOA)
Some references prefer Global Best PSO (pros: fast convergence, for simple problems; cons: liable to being trapped in local minimum than Local Best PSO) [3,9,12].On the other hand, some references favour Local Best PSO (pros: for complex problems; cons: slower algorithm) [3,9,12].In the middle stands [9] -both of these concepts show similar performance therefore GBEST or LBEST cannot be preferred.This all leads to conclusion that the advantage of the chosen PSO mainly depends on experience of user and on the type of examined problem.Global Best PSO algorithm is used in this paper.It can be modified in several ways for example [4,5,6,7,8].In this study were considered only four alternations, which are mentioned in the following text.

Maximal velocity v max
In just few iterations of PSO algorithm particle velocity can increase its value to the large magnitudes.This problem can be solved with velocity constraint -the maximal velocity v max [8]:

Inertia weight w
In [4] was introduced "A modified Particle Swarm Optimizer" and the inertia weight w was implemented into the equations (2,3).Equation ( 2) changed to: The inertia weight w is used to balance the global and local search abilities in the swarm.Parameter w could be treated in two possible ways: constant or adaptive.

Constant inertia weight w
In [4] authors suggested to choose inertia weight from interval 0.9<w<1.2 to ensure, that PSOA will have the best chance to find the global optimum (maximum of coefficient c p ) in moderate number of iterations.
Authors in [8] claimed that inertia weight w=0.8 is good "starting" choice for the algorithm, but it depends mainly on the magnitude of maximal velocity v max .

Adaptive inertia weight w
The inertia weight can be also computed adaptively.Many researchers have recommended that w should be large in the exploration state (the beginning of the algorithm) and small in the exploitation state (the end of the algorithm).In [5] authors proposed that w decreases linearly: where iteration represents current iteration of PSOA; max_iteration represents maximal number of iterations of PSOA; w max and w min are maximal/minimal inertia weights and they are usually set to 0.9 for w max and 0.4 for w min [4,5].

Parameters c 1 , c 2
Parameter c 1 represents "self-cognition".It can maintain diversity of the swarm (helps particle to explore computational area) [2].Parameter c 2 represents "social influence".This parameter forces the swarm to converge to the current global best solution of the examined function.It can also assist to the faster convergence [2].

Constant c 1 and c 2
In [2] authors recommended constant values of the parameters c 1 and c 2 equal to integer 2. This setting ensures that PSOA will have the best chance to find global optimum (maximum of the coefficient c p in our case).

Adaptive c 1 and c 2
Parameters c 1 and c 2 (c 1,2 shortly) can be treated in the same ways as inertia weight in -linearly [10]: and where iteration represents current iteration of PSOA; max_iteration represents maximal number of iterations of PSOA; c 1i and c 1f are maximal/minimal values of c 1 parameter and are usually set to 2.5 and 0.5 respectively; c 2i and c 2f are minimal/maximal values of c 2 parameter and are usually set to 0.5 and 2.5 respectively [10].

Specific computational area
In this study it was necessary to define computational area (the swarm territory) for the swarm due to possibility of creation of degenerated and inadmissible designs (geometries) -red rhomboid "FGHI" in Figure 4.If any particle emerged outside this area during the algorithm, it was immediately assigned to the relevant area border.

PSOA -Overview
In MATLAB software a master code was created, which calls and handles particle swarm optimization algorithm (PSOA -also programmed in MATLAB), mesh creation in Ansys ICEM CFD and CFD simulation in Ansys FLUENT.

Particle swarm optimization algorithm (PSOA)
Shape optimization using PSOA consists of several main steps (Figure 5).For this study, the ending criterion of the algorithm was set to the maximal number of iterations (max_iteration).Initial velocities (velocity initialization -Figure 5) were set as zeros [12].Every particle of the swarm contains information about locations of control points and value of c p and looks like this: (x 1 ; y 1 ; x 2 ; y 2 ; c p ).

Computational mesh for CFD
Fig. 6.Initial 2D domain mesh [1] During the shape optimization the computational 2D mesh was created in Ansys ICEM CFD.The mesh creation was handled with the text script, which ICEM can easily read and be controlled with.Mesh always met necessary quality requirements -quad type of mesh, skewness, and higher mesh resolution in near wall region (Figure 6).

CFD
Ansys FLUENT was also controlled via text scriptjournal.Total number of iterations was set to 1500 to ensure the adequate computational convergence.After each simulation, FLUENT saves essential variables (p 1 , p 2 -to the transcript file) for coefficient of pressure recovery c p -equation (1).

Evaluation procedure
Shape change was confined to section in between points E and D (Figure 1).This line could be parameterized in many different ways.This study focuses only on parameterization with Bézier curves, which have two free control points (i.e.four optimization parameters).Several test cases of chosen parameters were made (Figure 7).Every test case was repeated 5 times to capture the basic behaviour of the swarm.The ending criterion was set to the maximal number of iterations max_iteration=50.Two main attributes of the swarm were examined: the magnitudes of c p of the GBEST particle (Figure 8

Initial shape
The value of coefficient of pressure recovery c p for initial shape (Figure 1) was 0.736824 [13].

Swarm of 5 particles
Five particles in the swarm were chosen as the smallest reasonable number for the first six test cases.Average computational time for this choice was 28107 (s).
Free control points of Bézier curve were randomly generated in computational area at every beginning of PSOA (Figure 5).Every particle consists of two control points -one control point is represented as "x" and the other one as "+" (Figure 8).This selection of parameters shows the widest spread of final values of the coefficient of pressure recovery (Figure 11).Despite this disadvantage examined function has improved compared to the initial shape.C p value of the GBEST particle in the swarm limited by v max =0.05 reached up to 0.810670.

Adaptive w, constant c 1,2
In these test cases only inertia weight w is adaptivelinearly decreases its value through the optimization process -equation ( 7).This fact should in theory help the swarm more precisely locate global maximum of c p at the end of the PSO algorithm than the swarm with constant w.
Figures 13 -16 capture the examined attributes of the swarm of five particles with adaptive parameter w and constant parameters c 1,2 .The only difference between these two test cases was again in maximal velocities v max =0.025 and v max =0.05.
Note: Only one combination of w max and w min was tested and they were set to 0.9 and 0.4 respectively.This combination was based on [4].

Adaptive parameters w, c 1,2
The last two test cases contain swarms with adaptive parameters w and c 1,2 .The modification of these parameters was realized via linear decrease -equations (7,8,9).This combination of parameters forces the swarm to quick convergence to the global maximum of c p , but one danger could emerge -convergence to some local minimum without proper "search" of the computational area.
Figures 17 -20 capture the examined attributes of the swarm of five particles with adaptive parameter w and adaptive parameters c 1,2 .The only difference between these two test cases was once again in maximal velocities v max =0.025 and v max =0.05.
Note: Only one combination of w max , w min [4] and one combination of maximal/minimal values of c 1 /c 2 [10] was tested.

Partial results
Every test case improved the value of the coefficient c p compared to the initial shape -Table 2. The best diffuser design is shown in Figure 21.
For constant parameters w and c 1,2 is better choice smaller maximal velocity v max =0.025.This fact strongly depends on the size of computational area for the swarm.On the other hand -when at least one of the parameter is adaptive, the bigger maximal velocity (v max =0.05) seems to be a good choice to achieve better results.
Test case with adaptive w and c 1,2 shows the most stable values of the coefficient c p (Figure 19), on the other hand the most unstable test case was with parameters w=0.8; c 1,2 =2; v max =0.05 (Figure 10).
The disadvantage of each test case is in final positions of the control points of Bézier curves -in every single test they were spread out across the whole computational area and they never converged to a specific location.Note: The coefficient of pressure recovery c p increases as a result of weaker velocity gradients and shear layers (it was confirmed by theoretical analysis in [11]).

Swarm of 10 particles
The last six test cases contain swarms of ten particles.Average computational time of this choice was 56925 (s).
Once again: Free control points of Bézier curve were randomly generated in computational area at every beginning of PSOA (Figure 5).Every particle consists of two control points -one control point is represented as "x" and the other one as "+" (Figure 22).Note: Swarms with larger number of particles were not tested due to enormous size of computational time.
Test procedure was the same as in case of swarm of five -evaluation technique was executed according to the diagram in Figure 7.

Partial results
Every test case improved the value of c p coefficient compared to the initial shape -Table 3. The design of the best attempt is shown in Figure 35.Trend of smaller maximal velocity for constant parameters continues also for the swarms of ten particles -Figure 23.Ten particles should generally search the computational area thoroughly (compared to five particles), but the highest value of c p falls behind the swarms of five (for comparison Tables 2 and 3).
The disadvantage of almost every test case is in final positions of the control points of Bézier curves, but compared to the swarms of 5 particles, there is a promising combination of parameters w=0.9/0.4;c 1,2 =2.5/0.5;v max =0.025 -see Figures 31 and 32.This combination shows stable c p values and stable positions of control points in the end of the test; on the other hand this test case converged to some local maximum and did not bring the highest value of the examined function (coefficient of c p ). Note: The coefficient of pressure recovery c p increases as a result of weaker velocity gradients and shear layers (it was confirmed by theoretical analysis in [11]).Application of this optimization method strongly depends on examined problem.Present paper builds upon knowledge of [1] and focuses its energy on optimization of the hydraulic turbine diffuser.The shape optimization of the turbine diffuser became a basic model prototype of various methods [1,13].
Two types of swarms were constructed -with five and ten particles and with different combinations of parameters (Figure 7).Every time PSO managed to improve value of c p compared to the initial shape (Figure 1) -c p increases as a result of weaker velocity gradients and shear layers.Final diffuser designs (Figure 21 and 35) had similar shape and achieved similar c p -0.812915 for swarm of five and 0.812857 for the swarm of ten particles.
The main advantage of PSO: high value of coefficient of pressure recovery, because particles can "search" the computational area quite effectively.But on the other hand there are also some disadvantages: a need of creation of the computational area (Figure 4) to avert of production of degenerated designs and also a computational time was an issue -every iteration a great number of CFD simulation had to be performed.
Overall the PSO shows promising results -Table 4.  The comparison of computational time for all mentioned methods (Table 4) has not been done due to use of different computers (clusters).
There are some ways for improvement -future work will focus on some possible modifications of PSO e.g.: 1. Particle reduction during the algorithm -exploit a benefit of a large number of particles at the beginning of the algorithm and then if a certain criterion is met, algorithm will start to reduce the size of the swarm.This alternation brings potential savings in the computational time.2. Combination of two optimization methods -at the beginning of the optimization process usage of PSO for its robustness and later in the end of the process an application of a straight forward method as Nelder-Mead to ensure good convergence.Also this modification brings potential savings in computational time.

Fig. 5 .
Fig. 5. PSOA value of velocity at inlet is constant and does not match with the real velocity profile behind the runner of the swirl turbine[1].

Figure 7 .
Figure 7. Selections of the main parameters

Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .
Fig. 8. Random sets of control points -5 particles Note: Randomness of point generation was ensured in MATLAB by function rand.Rand returns random number from interval (0, 1).Constant parameters w, c 1,2 Figures 9 -12 capture the examined attributes of the swarm of five particles with constant parameters w and c 1,2 .The only difference between these two test cases was in maximal velocity v max =0.025 (v max =0.05).

Fig. 21 .
Fig. 21.The best diffuser design achieved by swarm of five

Fig. 22 .
Fig. 22. Random sets of control points in computational area -10 particles

Fig. 35 .
Fig. 35.Best diffuser design achieved by swarm of ten

Table 4 .
Comparison shows that the PSO algorithm produced slightly better values of c p than both Nelder-Mead algorithm and Adjoint solver.Control points positions do not match, when PSO and Nelder-Mead are compared -Figure36.But in general -all three methods show similar foundation of the final design -the tapered shape at the inlet of the diffuser and thus the crucial difference lies in the rest of the parametrized curve (Figure36).