A Biological Immune Mechanism-Based Quantum PSO Algorithm and Its Application in Back Analysis for Seepage Parameters

Back analysis for seepage parameters is a classic issue in hydraulic engineering seepage calculations. Considering the characteristics of inversion problems, including high dimensionality, numerous local optimal values, poor convergence performance, and excessive calculation time, a biological immune mechanism-based quantum particle swarm optimization (IQPSO) algorithm was proposed to solve the inversion problem. By introducing a concentration regulation strategy to improve the population diversity and a vaccination strategy to accelerate the convergence rate, the modiﬁed algorithm overcame the shortcomings of traditional PSO which can easily fall into a local optimum. Furthermore, a simple multicore parallel computation strategy was applied to reduce computation time. The eﬀectiveness and practicability of IQPSO were evaluated by numerical experiments. In this paper, taking one concrete face rock-ﬁll dam (CFRD) as a case, a back analysis for seepage parameters was accomplished by utilizing the proposed optimization algorithm and the steady seepage ﬁeld of the dam was analysed by the ﬁnite element method (FEM). Compared with immune PSO and quantum PSO, the proposed algorithm had better global search ability, convergence performance, and calculation rate. The optimized back analysis could obtain the permeability coeﬃcient of CFRD with high accuracy.


Introduction
In hydraulic engineering, a seepage calculation is used to obtain hydraulic factors such as head, discharge, and gradient by utilizing the basic seepage parameters for seepage stability analysis and operation management [1]. Determining correct seepage parameters is an important issue in seepage calculation, which directly influences the accuracy and rationality of seepage field analysis. Generally, there are three methods to determine seepage parameters including the empirical formula method, testing method, and back analysis method. Based on mathematical assumptions and empirical estimations, the results of the empirical formula method are inaccurate and not suitable in complex structures. e testing method, including indoor and in situ tests, can acquire accurate permeability coefficients. However, when the number of samples is small, the results are inaccurate for whole structures, and then, the number of samples is large and the cost is excessive. e back analysis method, based on the measured data and numerical simulation results, can cost-effectively obtain rational seepage parameters. erefore, it is necessary to develop a back analysis method to supplement and improve the work of determining seepage parameters. e essence of the back analysis method is to acquire material parameters while minimizing the error between the computed and measured values. Considering that the back analysis method requires repeated iterations and numerous calculations, many scholars have recently introduced various optimization algorithms into back analysis to solve the global optimal solution. Some intelligent algorithms, such as genetic algorithm (GA) [2], particle swarm optimization (PSO) [3], support vector machine (SVM) [4,5], artificial bee colony (ABC) [6], and artificial neural network (ANN) [7][8][9], have been widely used in back analysis for mechanical parameters. e work in the seepage field has also achieved some progress. For example, based on the concept of "nonlinear mapping," Chi et al. [10], Ni and Chi [11], and Liu et al. [12] built different neural network optimization models for the forward finite element analysis process and established a corresponding relationship between the measured values and simulated values. An improved GA was used by Wang and Liu [13,14] to solve the seepage parameter inversion problem in a fractured rock mass. Prasad and Rastogi [15] developed a groundwater-system inversion method based on the GA and built an inversion model of regional aquifers. By combining the FEM and adaptive GA, a back analysis method based on a stress-seepage coupling model was established by Deng et al. [16,17]. Although much work has been performed on back analysis algorithms, the optimization of the algorithm is not perfect. Due to the poor global search ability, limited convergence performance, and long FEM calculation time, these traditional algorithms cannot be applied practically in hydraulic engineering inversion problems. e PSO algorithm proposed by Kennedy and Eberhart [18] has attracted great attention in recent years, because of its simple structure and fast convergence speed. Many scholars [19][20][21] have realized the application of PSO in back analysis for seepage parameters. Some mechanisms have also been developed to improve the PSO algorithm, such as immune particle swarm optimization (IPSO) [22][23][24], simulated annealing particle swarm optimization (SAPSO) [25], and chaotic particle swarm optimization (CPSO) [26]. However, considering the PSO is not a global convergence algorithm, conventional improvement is limited. Inspired by Clerc and Kennedy's [27] analysis of the convergence of traditional PSO, Sun et al. [28,29] proposed a quantum particle swarm optimization algorithm (QPSO), which is a new algorithm model developed based on quantum mechanics theory. e particles with quantum behaviour have the capacity to search in the whole feasible solution place during iteration, so QPSO has a better global optimization ability. us, QPSO is utilized as an optimization algorithm in back analysis in this paper. Furthermore, to improve the population diversity and convergence performance, the concentration regulation strategy and the vaccination strategy based on a biological immune mechanism [30,31] are introduced to QPSO. In addition, a simple multicore parallel computing strategy is applied to reduce the FEM calculation time.
Based on the above three improvement strategies, we proposed an immune quantum particle swarm optimization algorithm in this paper and applied it in a back analysis for seepage parameter optimization. e effectiveness and practicability of IQPSO were assessed by numerical experiments. A case study was adopted to validate the effectiveness of the proposed algorithm for solving complex inversion problems in engineering seepage.

Biological Immune Mechanism-Based
Quantum PSO Algorithm 2.1. Basic QPSO Algorithm. QPSO is an improved PSO algorithm based on quantum mechanics theory, which has been proven to be a global convergence algorithm [32]. Assuming that the PSO system is similar to a quantum particle system, particles have quantum behaviour. e velocity and position of particles in quantum space cannot be determined at the same time, so the state of a particle is depicted by a wave function, instead of velocity. is method can make the particle appear at any point in the space with a certain probability, so that the particle can search in the whole feasible solution space and has a good global optimization ability.
According to the principle of QPSO, N is the size of the particle population, which acts as a search for better permeability coefficients. In addition, there are T iterations for D parameters. e position of the ith particle in the tth iteration is the vector , i � 1, 2, . . . , N without a velocity vector. e ith particle position corresponding to the historical optimal fitness in the tth iteration is stored as Pbest(t) � [Pbest i1 (t), Pbest i2 (t), . . . , Pbest i D (t)], and the particle position corresponding to the global optimal fitness in the previous tth iteration is also stored and updated as Each particle is updated according to the following equations [33]: where Mbest(t) is the mean best position which is defined as the mean of all the best positions of the particles in the tth iteration. P i (t) is the random position between Pbest i (t) and Gbest i (t) in the tth iteration, and r 1 and r 2 are random numbers distributed uniformly on [0, 1]. e parameter α in equations (3) and (4) is called the contraction expansion coefficient, which can be tuned to control the convergence rate of QPSO. In the early stage of evolution, a larger α value is beneficial to global search and can avoid falling into local convergence; with the increase of evolutionary algebra, a smaller α value is conducive to the convergence of the algorithm. e value is adaptively allocated per equation: where α max is the initial contraction expansion value, α min is the final contraction expansion value, and the values are 1.0 and 0.5, respectively. In addition, to prevent particles from exceeding the range of values, the following constraints are set: where x i d (t) represents the dth parameter value of the ith particle in the tth iteration. x dmax and x dmin are the maximum and the minimum of the dth parameter, respectively.

ree Improvement Strategies on QPSO.
Although QPSO can improve the global convergence of PSO effectively, it is found that the particles will still converge to the global optimal position during evolution, resulting in a loss of population diversity and a decrease of particle searching ability. erefore, three improvement strategies are introduced to QPSO to improve the algorithm performance.

Concentration Regulation Strategy.
In order to improve the particle searching ability and the population diversity, a concentration regulation strategy based on the particle similarity values is adopted. Based on the idea of the "Hamming distance," the particle similarity value is calculated by using equation (7) and the particle concentration is calculated by using equation (8). By inhibiting (or promoting) the particles with high (or low) concentration, the strategy can self-regulate to generate an appropriate number of new particles and maintain the particle concentration balance. e maximum concentration of particles is used to control the quality of the population. When the maximum concentration of particles reaches the threshold ξ, it is considered that the population diversity is poor and the evolutionary ability of the particle is not good. en, the particle concentration values will be sorted in decreasing turn and the subpopulation number M is determined by multiplying the maximum concentration dmax by the population size N. Finally, the first M high-concentration particles are eliminated, and the left N-M low-concentration particles are kept. is regulation strategy effectively improves the population diversity and can effectively avoid local convergence during evolution.
where ρ(X i ) represents the proximity degree of similarity between all particles and the ith particle. For example, if a particle in the population is close to the ith particle in the dth dimension, the value will be 1, otherwise, 0. d(X i ) represents the concentration of the ith particle. e larger the value is, the more similar the particles are and the worse the population diversity is. N and D represent the population size and the particle dimension, respectively.

Vaccination Strategy.
To accelerate the convergence rate and improve the algorithm performance, the immune memory operator is introduced to preserve the excellent particles in the evolution, which can prevent population degradation and guide the evolution process of the algorithm directionally. e global optimal solution Gbest in each iteration will be preserved and updated as a "vaccine." Furthermore, two different vaccination patterns are taken into account; one is that half of M high-concentration particles (mentioned in Section 2.2.1) will be replaced by new random particles. In another pattern, the other half of particles are selected to carry out vaccination by the "vaccine" particle. Particularly in the second pattern, these particles will be replaced by the memory particle first, and then, they will be chosen based on the mutation probability P m to suffer genetic mutation operations. e location of mutation particles is shown in equation (9). is vaccination strategy can not only ensure the existence of excellent particles (high fitness values) but also maintain the randomness of particles, which strongly improves the population quality and convergence performance: where r is a random number distributed uniformly on [0, 1]. A schematic diagram of the immune mechanism is shown in Figure 1. e particles in the initial population have low fitness values but a good evolutionary ability. As the iteration continues, the particles begin to aggregate, and the fitness values of some particles become higher but the evolutionary ability becomes worse. In this paper, the biological immune mechanism is used to improve particle performance. According to the two strategies mentioned above, M high-concentration particles are replaced and adjusted.
en, merged with the previously kept N-M particles, a new evolutionary population is produced. To be clear, in Figure 1, the red circle represents the global optimal particle Gbest in each generation, the green circle is the random particle, the yellow circle is the mutation particle based on Gbest, and the other grey circles are the initial population particles.

Multicore Parallel Computing Strategy.
With the popularity of multicore processors and the development of parallel computing research, parallel computing has gradually become an important way to improve computing efficiency. Computing tasks are assigned to multiple CPUs and performed simultaneously, which offers a good way to reduce calculation time.
e existing parallel programming models can support many programming languages and have powerful functions such as MPI [34,35] and OpenMP [36]. However, these models are too complex, abstract, and difficult to program. In this paper, the MATLAB platform is applied in parallel computing, which is a simple and easy way to develop parallel optimization programs to realize multicore parallel computation of the IQPSO algorithm.
During the optimization process, the calculation of the fitness value is the most time-consuming part. It is a reasonable choice to adopt the Master-Slave parallel structure. In the MATLAB platform, the "parpool" function is used to configure and open a parallel pool and the "client and worker" pattern is used to execute the "parfor" loop in parallel computing. e MATLAB client is the main process responsible for controlling the entire calculation process, including particle initialization, particle location updating, fitness value comparison, concentration regulation strategy, and vaccination strategy. e MATLAB worker is a slave process, which is mainly responsible for the calculation of fitness values. e communications between the master and slave processes are shown in Figure 2. Parameter n in Figure 2 represents the number of slave processes. In particular, to avoid slave process idling, the particle number of the particle swarm is an integer multiple of n.

IQPSO Algorithm Procedures.
e steps of the IQPSO algorithm are as follows and the flow chart of the IQPSO algorithm is shown in Figure 3: Step 1: set initial parameters. Population size N, dimension D, maximum iteration number T, population diversity threshold ξ, and mutation probability P m are determined. Initialize the position of each particle and calculate the fitness value. In addition, update the historical optimal position of the particle Pbest, the mean best position of the population Mbest and the global optimal position of the population Gbest which, is stored in the memory bank as a "vaccine." Step 2: initiate population diversity judgement. e concentration of each particle is calculated using equations (7) and (8), and the population diversity is judged. If the maximum concentration of particles reaches the population diversity threshold ξ, then go to Step 3; otherwise, jump to Step 5.
Step 3: realize concentration regulation. Sort by the particle concentration values and control subpopulation size with maximum concentration. Eliminate M high-concentration particles and keep N-M low-concentration particles.
Step 4: implement vaccination. Half of M high-concentration particles are replaced by new random particles that meet the constraint, and the other half of particles are selected to carry out vaccination with the memory particle. e particles in the second pattern may suffer genetic mutations based on the mutation probability P m . e N-M particles are kept in Step 3, and the M particles generated in Step 4 are merged to make up a new population.
Step 5: Update the particles. e position of particles is updated using equations (1)∼(6), and a new generation of N particles is produced. According to the objective function and a simple multicore parallel computing strategy, the fitness value of each particle is calculated and the positions of Pbest, Mbest, and Gbest are updated. Finally, the memory bank will be reset.
Step 6: Implement algorithm iterations. Judge whether the number of iterations reaches the predetermined value; if not, return to Step 2; otherwise, end the iteration and output the global optimal solution and fitness value.

Numerical Experiments
In order to verify the feasibility and validity of the proposed IQPSO algorithm, three typical single-peak functions, Sphere, Step  Figure 1: Schematic diagram of the immune mechanism.
In the numerical experiment, the particle swarm N, the maximum iteration number T, the population diversity threshold ξ, the mutation probability P m , the acceleration constants c 1 and c 2 , and the initial and final inertial weights are set as 40, 1000, 0.35, 0.5, 2.0, 2.0, 0.9, and 0.4, respectively. In addition, the searching hyperspace dimension D is considered to be 30, 60, and 90. To ensure the stability of the calculation results, each test function is performed 50 times independently using four different algorithms. e average optimum value of the results in different dimensions is also listed in Table 2, and the best calculation results between four algorithms are roughly marked. Figure 4 shows the convergence process curves of the logarithmic average fitness values for four algorithms in 60 dimensions. e smaller the value is, the higher the solution accuracy is.
As shown in Table 2, the results of the IQPSO algorithm are extremely close to the optimal value 0 and obviously better than those obtained from other three algorithms in almost all cases.
e proposed algorithm has better performance whether the test function is single-peak or multipeak. Figure 4 indicates that, in the early stage of evolution, the four algorithms can search for the optimal solution quickly, and the fitness value curve decreases continuously. However, as the evolution continues, the other three algorithms begin to fall into a local optimal solution and appear to converge prematurely, while the fitness value curve of the IQPSO algorithm (marked in red) continues to decline and is able to search for the better solution. Actually, if given enough iteration, the optimal solution can ultimately be found. In conclusion, the experimental results show that the proposed algorithm in this paper has higher solution accuracy, stronger global search ability, and better robustness.

Engineering Overview.
In order to verify the effectiveness of the proposed algorithm in hydraulic engineering optimized back analysis, a three-dimensional finite element model of a Chinese concrete-face rock-fill dam (CFRD) was established. Based on a geological survey and monitor data, a back analysis for seepage parameters was accomplished with the proposed IQPSO algorithm. More specifically, the dam is located in the Sinan River, Yunnan Province, China, and the maximum height of the dam is 115 m. e normal storage level and the dead water level are 900.00 m and 860.00 m, respectively. e dam has the standard CFRD structure, including an upstream concrete face slab, bedding material, transition material, primary rock-fill zone, secondary rockfill zone, and impervious curtain. e bottom boundary of the impervious curtain enters the water-proof layer (q < 3 Lu) and its deepest depth is 45 m. e maximum cross section of the dam body is shown in Figure 5. e permeability coefficients of five materials need to be optimized. ese materials are impervious curtain, primary rock-fill zones, and three strata of the dam foundation. Based on geological prospecting data and engineering experience, permeability coefficients of these five materials are limited to a reasonable range. Assuming that all the materials are isotropic, permeability coefficients are shown in Table 3.

Analysis of Measured Osmotic Pressure Data.
e osmometers in the project are marked in red in Figure 5 and are distributed at the downstream side of the impervious curtain and the dam foundation. In view of the measurement bias and the sensitivity of back analysis, the osmotic pressure data of points P11-P13 and P15-P19, a total of 8 observation points, are selected as the basic information for back analysis. Figure 6 shows the upstream water level of the dam and the measured hydraulic head of a typical osmometer (P11) from 19 January 2009 to 12 February 2014.
Comparing the time curve of the upstream water level and the measured hydraulic head of the typical osmometer in Figure 6, we find out that the hydraulic head of the measured point (P11) changes periodically during the reservoir stable storage period. at is, with the upstream water level rising and declining, the hydraulic head increases and decreases, which is consistent with the conclusion drawn by Wu [37] in his book. Considering that there is a seepage process during the upstream water level fluctuation, the measured hydraulic head lags behind the reservoir water level, which is called the "lag effect." Hence, one complete and stable storage period is selected as the calculation period (2009-8-20-2010-8-20). e upstream water level is 887.36 m, and the downstream water level is 800.59 m. e measured hydraulic head values (listed in Table 4) during this calculation period are taken as the basic data in back analysis.

Finite Element Analysis.
e finite element method [38,39] is a common numerical method to solve the seepage field of earth-rock dams, and the numerical Mathematical Problems in Engineering simulation technology is relatively mature. In this paper, a three-dimensional seepage "CNPM" finite element program compiled by the research group is used to calculate the steady seepage field. According to the similar engineering experience, the numerical model domain is set as follows. In the x-direction, the length of the dam foundation beyond the dam body is the height of the dam along the upstream and downstream. In the z-direction, the depth of the foundation is twice the height of the dam. e three-dimensional finite element model is established Start Multicore parallel computing t = 1 Step 1 Step 2 Step 3 Step 4 Step 5 Step 6 Multicore parallel computing (shown in Figure 7), and 7582 nodes and 3686 elements are built by using the superelement division method. e reservoir water level data are given as the boundary condition of the upstream and downstream water levels, and other parts of the model are considered as impermeable boundary in FEM calculation.
e measured values in Table 4 are selected as the water-head constraint condition. With the unit width dam section, the hydraulic head values are calculated under different combinations of permeability coefficients by FEM. Additionally, 13 colours represent 13 material regions, as shown in Figure 7.

Mathematical Model.
e mathematical model for permeability coefficients of the dam is based on the calculation results of the FEM and the measured hydraulic head values. Generally, according to the hydraulic head data, only the ratio of permeability coefficients in each material area could be inverted. However, if given some certain permeability coefficients of certain materials, permeability coefficients of other materials can be basically determined. In this case, the optimization mathematical model is established as follows: s.t. a j ≤ k j ≤ b j , (j � 1, 2, . . . , 5).
In the above model, equation (10) is the objective function based on the least-squares criterion. k � [k 1 , k 2 , · · · , k 5 ] T is a group of permeability coefficients to be inverted. h c i and h m i are the ith computed and measured hydraulic head, respectively, and w i is the ith weight which adds up to 1 (w i (i � 1, 2, · · ·, 8) � [0.2, 0.2, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1] in this case). Equation (11) is the constraint condition that gives the range of permeability coefficients. a j and b j are the jth minimum and maximum of k j , respectively.
Next, based on the IQPSO algorithm mentioned in Section 2, the optimization back analysis for permeability coefficients is carried out and the inversion results are compared with QPSO and IPSO. In addition, in order to allow for the time superiority of a multicore parallel computing strategy, a back analysis by serial IQPSO is also performed. e particle swarm N, the maximum iteration number T, the population diversity threshold ξ, the mutation probability P m , the acceleration constants c 1 and c 2 , and the initial and final inertial weights are set as 30, 200, 0.35,  Step , 600] 0 0.5, 2.0, 2.0, 0.9, and 0.4, respectively. e running environment is a 64-bit operating system, Intel (R) Core (TM) i5-8400 CPU, 6 cores, 2.80 GHz. Hence, 6 cores are applied in parallel computing. Based on the three-dimensional seepage "CNPM" finite element program, MATLAB is used to adjust the permeability coefficients by modifying the input file of the program and to read the hydraulic head of the nodes. en, based on equation (10), fitness function is established to optimize the permeability coefficients in a specific range. In particular, this case sets 200 iterations because with approximately 125 iterations, the algorithm can be considered to converge in a tentative calculation.
However, for different cases, iterations of the algorithm should be adjusted according to the actual results of the tentative calculation.

Back Analysis Results.
e optimal fitness curves for the three algorithms are shown in Figure 8, and the simulated and measured hydraulic head of the steady seepage field at the measuring points are listed in Table 5, which also contains the absolute error and relative error of the simulated and measured hydraulic head. From the above analysis, all three algorithms can achieve convergence under the 200- iteration condition. Compared with QPSO and IPSO, IQPSO has a smaller optimal objective function value and a smaller error between the measured values and the simulated values, which indicates that the proposed algorithm in this paper has higher accuracy, stronger global search ability, and better robustness. Based on the error    Mathematical Problems in Engineering analysis method [10], the maximum absolute error by IQPSO is 1.42 m and the maximum relative error is 9.97%, within a reasonable range, which strongly indicates that the permeability coefficients are reasonable.
Compared with serial IQPSO, the computation time of parallel IQPSO is decreased dramatically from 22.28 h to 4.52 h. According to the parallel performance analysis [40], the computing acceleration ratio and parallel   efficiency are 4.93 and 0.82, respectively, which indicates that the time superiority of the multicore parallel strategy is remarkable. is parallel strategy not only provides convenience for the debugging and operating of the model but also offers a good approach for applying various optimization algorithms to complex hydraulic engineering inversion problems. e permeability coefficients acquired from IQPSO are listed in Table 6. e contour map of the maximum section of the dam body and dam foundation under steady seepage field is shown in Figure 9. Given the permeability coefficients of upstream face slab and other materials, the ratio of permeability coefficients predicated from seepage hydraulic water data is the absolute value of materials. It can be seen from the contour map that the antiseepage system composed of "face slab-toe plate-impervious curtain" has a good antiseepage effect. e hydraulic head is reduced by 78.10 m, accounting for 90.01% of the total water head. e phreatic line of the dam (marked red in Figure 9) is at a low position and decreases continuously, which indicates that the dam body is waterless and has good seepage stability.

Conclusion
Aiming at the characteristics of seepage parameter inversion problems, including high dimensionality, numerous local optimal values, poor convergence performance, and excessive calculation time, a biological immune mechanismbased quantum particle swarm optimization algorithm was proposed in this paper to solve the inversion problem. IQPSO was tested by numerical experiments. en, an optimized back analysis for permeability coefficients of one concrete face rock-fill dam was realized and the steady seepage field of the dam was analysed. e main conclusions are as follows: (1) Based on the concentration regulation strategy and the vaccination strategy, the global search ability and the convergence performance of the novel algorithm had been improved significantly. e numerical    experimental results of IQPSO were obviously better than the other three algorithms, which proved the effectiveness and practicability of the algorithm. (2) e back analysis results showed that compared with QPSO and IPSO, IQPSO had a smaller objective function value and fewer errors between the calculated hydraulic head values and the measured values. e maximum absolute error by IQPSO was 1.42 m, and the maximum relative error was 9.97%. e analysis results of the steady seepage field indicated that during the normal operation of the reservoir, the antiseepage system composed of "Face slab-toe plateimpervious curtain" had a good antiseepage effect.
(3) Compared with serial IQPSO, the computing time was obviously shortened by parallel IQPSO and the acceleration ratio and parallel efficiency were 4.93 and 0.82, respectively, which indicated that the time superiority of the multicore parallel strategy was remarkable.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest.