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

Authors: Jiacheng Tan; Liqun Xu; Kailai Zhang; Chao Yang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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 modified 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 effectiveness and practicability of IQPSO were evaluated by numerical experiments. In this paper, taking one concrete face rock-fill dam (CFRD) as a case, a back analysis for seepage parameters was accomplished by utilizing the proposed optimization algorithm and the steady seepage field of the dam was analysed by the finite 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 coefficient of CFRD with high accuracy.


---

## Body

## 1. 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. The 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. The back analysis method, based on the measured data and numerical simulation results, can cost-effectively obtain rational seepage parameters. Therefore, it is necessary to develop a back analysis method to supplement and improve the work of determining seepage parameters.The 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–9], have been widely used in back analysis for mechanical parameters. The 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.The 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–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–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. The 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. Thus, 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. The 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.


## 2. 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. The 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. This 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. The position of the ith particle in the tth iteration is the vector Xit=xi1t,xi2t,…,xiDt,i=1,2,…,N without a velocity vector. The ith particle position corresponding to the historical optimal fitness in the tth iteration is stored as Pbestt=Pbesti1t,Pbesti2t,…,PbestiDt, and the particle position corresponding to the global optimal fitness in the previous tth iteration is also stored and updated as Gbestt=Gbesti1t,Gbesti2t,…,GbestiDt. Each particle is updated according to the following equations [33]:(1)Mbestt=1N∑i=1NPbestit,(2)Pit=r1⋅Pbestit+1−r1⋅Gbestit,(3)Xit+1=Pit−α⋅Mbestt−Xit⋅ln1/r2,r2≥0.5,(4)Xit+1=Pit+α⋅Mbestt−Xit⋅ln1/r2,r2<0.5,where Mbestt is the mean best position which is defined as the mean of all the best positions of the particles in the tth iteration. Pit is the random position between Pbestit and Gbestit in the tth iteration, and r1 and r2 are random numbers distributed uniformly on [0, 1]. The 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. The value is adaptively allocated per equation:(5)α=αmin+αmax−αmin⋅T−tT,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:(6)xidt=xdmax,xidt≥xdmax,xidt,xdmin≤xidt≤xdmax,xdmin,xidt≤xdmin,where xidt represents the dth parameter value of the ith particle in the tth iteration. xdmax and xdmin are the maximum and the minimum of the dth parameter, respectively.


### 2.2. Three 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. Therefore, three improvement strategies are introduced to QPSO to improve the algorithm performance.


#### 2.2.1. 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. The 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. Then, 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. This regulation strategy effectively improves the population diversity and can effectively avoid local convergence during evolution.(7)ρXi=∑j=1N∑d=1D1,0.95≤xidxjd≤1.05,0,others,i=1,2,⋯,N,(8)dXi=ρXiN⋅D,i=1,2,⋯,N,where ρXi 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. dXi represents the concentration of the ith particle. The 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.


#### 2.2.2. 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. The global optimal solutionGbest 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 Pm to suffer genetic mutation operations. The location of mutation particles is shown in equation (9). This 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:(9)xidt=xidt+xdmax−xidt×r×1−tT,r>0.5,xidt+xidt−xdmin×r×1−tT,r≤0.5,where r is a random number distributed uniformly on [0, 1].A schematic diagram of the immune mechanism is shown in Figure1. The 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. Then, 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.Figure 1
Schematic diagram of the immune mechanism.


#### 2.2.3. 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. The 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. The 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. The MATLAB worker is a slave process, which is mainly responsible for the calculation of fitness values. The communications between the master and slave processes are shown in Figure2. 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.Figure 2
Communications between the master process and the slave processes.


### 2.3. IQPSO Algorithm Procedures

The steps of the IQPSO algorithm are as follows and the flow chart of the IQPSO algorithm is shown in Figure3:Step 1: set initial parameters. Population sizeN, dimension D, maximum iteration number T, population diversity threshold ξ, and mutation probability Pm 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. The 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. EliminateM high-concentration particles and keep N-M low-concentration particles.Step 4: implement vaccination. Half ofM 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. The particles in the second pattern may suffer genetic mutations based on the mutation probability Pm. The 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. The 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.Figure 3
Flow chart of the IQPSO 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. The 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. This 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. The position of the ith particle in the tth iteration is the vector Xit=xi1t,xi2t,…,xiDt,i=1,2,…,N without a velocity vector. The ith particle position corresponding to the historical optimal fitness in the tth iteration is stored as Pbestt=Pbesti1t,Pbesti2t,…,PbestiDt, and the particle position corresponding to the global optimal fitness in the previous tth iteration is also stored and updated as Gbestt=Gbesti1t,Gbesti2t,…,GbestiDt. Each particle is updated according to the following equations [33]:(1)Mbestt=1N∑i=1NPbestit,(2)Pit=r1⋅Pbestit+1−r1⋅Gbestit,(3)Xit+1=Pit−α⋅Mbestt−Xit⋅ln1/r2,r2≥0.5,(4)Xit+1=Pit+α⋅Mbestt−Xit⋅ln1/r2,r2<0.5,where Mbestt is the mean best position which is defined as the mean of all the best positions of the particles in the tth iteration. Pit is the random position between Pbestit and Gbestit in the tth iteration, and r1 and r2 are random numbers distributed uniformly on [0, 1]. The 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. The value is adaptively allocated per equation:(5)α=αmin+αmax−αmin⋅T−tT,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:(6)xidt=xdmax,xidt≥xdmax,xidt,xdmin≤xidt≤xdmax,xdmin,xidt≤xdmin,where xidt represents the dth parameter value of the ith particle in the tth iteration. xdmax and xdmin are the maximum and the minimum of the dth parameter, respectively.


## 2.2. Three 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. Therefore, three improvement strategies are introduced to QPSO to improve the algorithm performance.


### 2.2.1. 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. The 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. Then, 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. This regulation strategy effectively improves the population diversity and can effectively avoid local convergence during evolution.(7)ρXi=∑j=1N∑d=1D1,0.95≤xidxjd≤1.05,0,others,i=1,2,⋯,N,(8)dXi=ρXiN⋅D,i=1,2,⋯,N,where ρXi 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. dXi represents the concentration of the ith particle. The 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.


### 2.2.2. 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. The global optimal solutionGbest 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 Pm to suffer genetic mutation operations. The location of mutation particles is shown in equation (9). This 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:(9)xidt=xidt+xdmax−xidt×r×1−tT,r>0.5,xidt+xidt−xdmin×r×1−tT,r≤0.5,where r is a random number distributed uniformly on [0, 1].A schematic diagram of the immune mechanism is shown in Figure1. The 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. Then, 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.Figure 1
Schematic diagram of the immune mechanism.


### 2.2.3. 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. The 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. The 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. The MATLAB worker is a slave process, which is mainly responsible for the calculation of fitness values. The communications between the master and slave processes are shown in Figure2. 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.Figure 2
Communications between the master process and the slave processes.


## 2.2.1. 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. The 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. Then, 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. This regulation strategy effectively improves the population diversity and can effectively avoid local convergence during evolution.(7)ρXi=∑j=1N∑d=1D1,0.95≤xidxjd≤1.05,0,others,i=1,2,⋯,N,(8)dXi=ρXiN⋅D,i=1,2,⋯,N,where ρXi 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. dXi represents the concentration of the ith particle. The 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.


## 2.2.2. 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. The global optimal solutionGbest 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 Pm to suffer genetic mutation operations. The location of mutation particles is shown in equation (9). This 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:(9)xidt=xidt+xdmax−xidt×r×1−tT,r>0.5,xidt+xidt−xdmin×r×1−tT,r≤0.5,where r is a random number distributed uniformly on [0, 1].A schematic diagram of the immune mechanism is shown in Figure1. The 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. Then, 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.Figure 1
Schematic diagram of the immune mechanism.


## 2.2.3. 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. The 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. The 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. The MATLAB worker is a slave process, which is mainly responsible for the calculation of fitness values. The communications between the master and slave processes are shown in Figure2. 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.Figure 2
Communications between the master process and the slave processes.


## 2.3. IQPSO Algorithm Procedures

The steps of the IQPSO algorithm are as follows and the flow chart of the IQPSO algorithm is shown in Figure3:Step 1: set initial parameters. Population sizeN, dimension D, maximum iteration number T, population diversity threshold ξ, and mutation probability Pm 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. The 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. EliminateM high-concentration particles and keep N-M low-concentration particles.Step 4: implement vaccination. Half ofM 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. The particles in the second pattern may suffer genetic mutations based on the mutation probability Pm. The 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. The 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.Figure 3
Flow chart of the IQPSO algorithm.


## 3. Numerical Experiments

In order to verify the feasibility and validity of the proposed IQPSO algorithm, three typical single-peak functions, Sphere, Step, and Rosenbrock, and three complex multipeak functions, Ackley, Rastrigin, and Griewank, are selected for numerical experiments. The mathematical formula, the searching space, and the optimum solution of six test functions are shown specifically in Table1. The experimental results of the IQPSO algorithm are compared with PSO, QPSO, and IPSO.Table 1
Six typical test functions.
Function nameMathematical formulaSearching spaceOptimal solutionSpherefx=∑i=1Dxi2[−5.12, 5.12]0Stepfx=∑i=1Dxi+0.52[−30, 30]0Rosenbrockfx=∑i=1D−1100xi+1−xi22+xi−12[−30, 30]0Ackleyfx=−20exp−0.21/n∑i=1nxi2−exp1/n∑i=1ncos2πxi+20+e[−30, 30]0Rastriginfx=10n+∑i=1nxi2−10cos2πxi[−5.12, 5.12]0Griewankfx=1/4000∑i=1nxi2−∏i=1ncosxi/i+1[−600, 600]0In the numerical experiment, the particle swarmN, the maximum iteration number T, the population diversity threshold ξ, the mutation probability Pm, the acceleration constants c1 and c2, 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. The 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. The smaller the value is, the higher the solution accuracy is.Table 2
Comparison of average optimum values on test functions.
Test functionDimensionPSOQPSOIPSOIQPSOSphere305.50E − 055.56E − 112.64E − 068.84E − 16607.44E − 019.88E − 049.00E − 036.08E − 07903.66E + 001.50E − 011.16E − 015.18E − 04Step303.34E − 021.75E − 093.60E − 036.50E − 04601.16E + 004.70E − 025.10E − 023.37E − 02906.35E + 005.33E + 004.39E − 012.46E − 01Rosenbrock309.39E + 016.99E + 015.54E + 012.58E + 01606.15E + 025.02E + 023.38E + 021.67E + 02902.37E + 031.22E + 038.08E + 024.92E + 02Ackley304.62E + 001.99E − 058.38E − 037.69E − 08606.40E + 003.06E − 012.78E + 002.40E − 03907.61E + 002.26E + 004.64E + 005.42E − 02Rastrigin304.61E + 011.92E + 014.43E + 011.42E + 01601.61E + 028.19E + 011.01E + 024.99E + 01903.25E + 021.97E + 021.56E + 021.20E + 02Griewank304.21E + 005.92E − 059.91E − 011.18E − 09601.45E + 012.50E − 012.46E + 005.08E − 04902.48E + 011.49E + 004.10E + 005.11E − 02Figure 4
Convergence process curves of average fitness for four algorithms (60 dimensions): (a) Sphere function; (b) Step function; (c) Rosenbrock function; (d) Ackley function; (e) Rastrigin function; (f) Griewank function.
(a)(b)(c)(d)(e)(f)As shown in Table2, 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. The 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.


## 4. Case Study

### 4.1. 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. The normal storage level and the dead water level are 900.00 m and 860.00 m, respectively. The dam has the standard CFRD structure, including an upstream concrete face slab, bedding material, transition material, primary rock-fill zone, secondary rock-fill zone, and impervious curtain. The bottom boundary of the impervious curtain enters the water-proof layer (q < 3 Lu) and its deepest depth is 45 m. The maximum cross section of the dam body is shown in Figure 5.Figure 5
The maximum cross section of the dam body.The permeability coefficients of five materials need to be optimized. These 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 Table3.Table 3
Permeability coefficient of the dam body and dam foundation (cm/s).
Inversion materialRange of permeability coefficientDetermined materialValue of permeability coefficientRelative impermeable layer3 × 10−6∼3 × 10−4Upstream face slab1 × 10−7Aquitard layer6 × 10−6∼6 × 10−4Bedding material5 × 10−4Weak weathered layer8 × 10−5∼8 × 10−3Transition material3 × 10−1Primary rock-fill1 × 10−1∼1 × 101Secondary rock-fill1 × 10−1Impervious curtain3 × 10−6∼3 × 10−4Upstream blanket1 × 10−2Upstream arbitrary fill material1 × 100Downstream slope5 × 10−2Dam crest material3.5 × 10−7


### 4.2. Analysis of Measured Osmotic Pressure Data

The osmometers in the project are marked in red in Figure5 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.Figure 6
Time curve of upstream water level and measured hydraulic head of typical osmometer (P11).Comparing the time curve of the upstream water level and the measured hydraulic head of the typical osmometer in Figure6, we find out that the hydraulic head of the measured point (P11) changes periodically during the reservoir stable storage period. That 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). The upstream water level is 887.36 m, and the downstream water level is 800.59 m. The measured hydraulic head values (listed in Table 4) during this calculation period are taken as the basic data in back analysis.Table 4
Measured hydraulic head of osmometers in back analysis (m).
Osmometer numberP11P12P13P15P16P17P18P19Measured hydraulic head31.0621.2612.829.437.565.245.765.96


### 4.3. Finite Element Analysis

The finite element method [38, 39] is a common numerical method to solve the seepage field of earth-rock dams, and the numerical 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. The three-dimensional finite element model is established (shown in Figure 7), and 7582 nodes and 3686 elements are built by using the superelement division method. The 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. The 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.Figure 7
Finite element mesh of unit width dam section.


### 4.4. Optimized Back Analysis for Permeability Coefficient

#### 4.4.1. Mathematical Model

The 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:(10)minFk=∑i=18wihic−him2,(11)s.t.aj≤kj≤bj,j=1,2,…,5.In the above model, equation (10) is the objective function based on the least-squares criterion. k=k1,k2,⋯,k5T is a group of permeability coefficients to be inverted. hic and him are the ith computed and measured hydraulic head, respectively, and wi is the ith weight which adds up to 1 (wi (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. aj and bj are the jth minimum and maximum of kj, respectively.Next, based on the IQPSO algorithm mentioned in Section2, 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. The particle swarm N, the maximum iteration number T, the population diversity threshold ξ, the mutation probability Pm, the acceleration constants c1 and c2, and the initial and final inertial weights are set as 30, 200, 0.35, 0.5, 2.0, 2.0, 0.9, and 0.4, respectively. The 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. Then, 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.


#### 4.4.2. Back Analysis Results

The optimal fitness curves for the three algorithms are shown in Figure8, 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 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.Figure 8
The optimal fitness curves for three algorithms.Table 5
Results of the simulated and measured hydraulic head for three algorithms.
AlgorithmOsmometer numberP11P12P13P15P16P17P18P19IPSOMeasured value (m)808.06806.26805.82805.43805.56804.24804.76802.96Simulated value (m)808.65807.18805.62804.30805.19806.02804.95803.56Absolute error (m)0.590.92−0.20−1.13−0.361.780.190.60Relative error (%)3.275.63−1.29−7.30−2.3412.511.264.66QPSOSimulated value (m)808.54807.06805.50804.20805.09805.92804.86803.56Absolute error (m)0.480.80−0.32−1.23−0.471.680.100.60Relative error (%)2.664.94−2.02−7.98−2.9911.810.674.66IQPSOSimulated value (m)808.41806.67805.43804.17805.06805.66804.84803.55Absolute error (m)0.350.41−0.39−1.26−0.501.420.080.59Relative error (%)1.932.52−2.49−8.16−3.199.970.544.56Compared 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. This 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.The permeability coefficients acquired from IQPSO are listed in Table6. The 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. The hydraulic head is reduced by 78.10 m, accounting for 90.01% of the total water head. The 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.Table 6
Inversion results of permeability coefficients by IQPSO (cm/s).
Inversion materialRange of permeability coefficientValue of permeability coefficientRelative impermeable layer3 × 10−6∼3 × 10−42.58 × 10−4Aquitard layer6 × 10−6∼6 × 10−44.05 × 10−4Weak weathered layer8 × 10−5∼8 × 10−37.91 × 10−3Primary rock-fill1 × 10−1∼1 × 1012.03 × 10−1Impervious curtain3 × 10−6∼3 × 10−43.39 × 10−5Figure 9
Contour map of the maximum section of the dam body and dam foundation.


## 4.1. 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. The normal storage level and the dead water level are 900.00 m and 860.00 m, respectively. The dam has the standard CFRD structure, including an upstream concrete face slab, bedding material, transition material, primary rock-fill zone, secondary rock-fill zone, and impervious curtain. The bottom boundary of the impervious curtain enters the water-proof layer (q < 3 Lu) and its deepest depth is 45 m. The maximum cross section of the dam body is shown in Figure 5.Figure 5
The maximum cross section of the dam body.The permeability coefficients of five materials need to be optimized. These 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 Table3.Table 3
Permeability coefficient of the dam body and dam foundation (cm/s).
Inversion materialRange of permeability coefficientDetermined materialValue of permeability coefficientRelative impermeable layer3 × 10−6∼3 × 10−4Upstream face slab1 × 10−7Aquitard layer6 × 10−6∼6 × 10−4Bedding material5 × 10−4Weak weathered layer8 × 10−5∼8 × 10−3Transition material3 × 10−1Primary rock-fill1 × 10−1∼1 × 101Secondary rock-fill1 × 10−1Impervious curtain3 × 10−6∼3 × 10−4Upstream blanket1 × 10−2Upstream arbitrary fill material1 × 100Downstream slope5 × 10−2Dam crest material3.5 × 10−7


## 4.2. Analysis of Measured Osmotic Pressure Data

The osmometers in the project are marked in red in Figure5 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.Figure 6
Time curve of upstream water level and measured hydraulic head of typical osmometer (P11).Comparing the time curve of the upstream water level and the measured hydraulic head of the typical osmometer in Figure6, we find out that the hydraulic head of the measured point (P11) changes periodically during the reservoir stable storage period. That 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). The upstream water level is 887.36 m, and the downstream water level is 800.59 m. The measured hydraulic head values (listed in Table 4) during this calculation period are taken as the basic data in back analysis.Table 4
Measured hydraulic head of osmometers in back analysis (m).
Osmometer numberP11P12P13P15P16P17P18P19Measured hydraulic head31.0621.2612.829.437.565.245.765.96


## 4.3. Finite Element Analysis

The finite element method [38, 39] is a common numerical method to solve the seepage field of earth-rock dams, and the numerical 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. The three-dimensional finite element model is established (shown in Figure 7), and 7582 nodes and 3686 elements are built by using the superelement division method. The 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. The 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.Figure 7
Finite element mesh of unit width dam section.


## 4.4. Optimized Back Analysis for Permeability Coefficient

### 4.4.1. Mathematical Model

The 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:(10)minFk=∑i=18wihic−him2,(11)s.t.aj≤kj≤bj,j=1,2,…,5.In the above model, equation (10) is the objective function based on the least-squares criterion. k=k1,k2,⋯,k5T is a group of permeability coefficients to be inverted. hic and him are the ith computed and measured hydraulic head, respectively, and wi is the ith weight which adds up to 1 (wi (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. aj and bj are the jth minimum and maximum of kj, respectively.Next, based on the IQPSO algorithm mentioned in Section2, 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. The particle swarm N, the maximum iteration number T, the population diversity threshold ξ, the mutation probability Pm, the acceleration constants c1 and c2, and the initial and final inertial weights are set as 30, 200, 0.35, 0.5, 2.0, 2.0, 0.9, and 0.4, respectively. The 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. Then, 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.


### 4.4.2. Back Analysis Results

The optimal fitness curves for the three algorithms are shown in Figure8, 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 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.Figure 8
The optimal fitness curves for three algorithms.Table 5
Results of the simulated and measured hydraulic head for three algorithms.
AlgorithmOsmometer numberP11P12P13P15P16P17P18P19IPSOMeasured value (m)808.06806.26805.82805.43805.56804.24804.76802.96Simulated value (m)808.65807.18805.62804.30805.19806.02804.95803.56Absolute error (m)0.590.92−0.20−1.13−0.361.780.190.60Relative error (%)3.275.63−1.29−7.30−2.3412.511.264.66QPSOSimulated value (m)808.54807.06805.50804.20805.09805.92804.86803.56Absolute error (m)0.480.80−0.32−1.23−0.471.680.100.60Relative error (%)2.664.94−2.02−7.98−2.9911.810.674.66IQPSOSimulated value (m)808.41806.67805.43804.17805.06805.66804.84803.55Absolute error (m)0.350.41−0.39−1.26−0.501.420.080.59Relative error (%)1.932.52−2.49−8.16−3.199.970.544.56Compared 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. This 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.The permeability coefficients acquired from IQPSO are listed in Table6. The 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. The hydraulic head is reduced by 78.10 m, accounting for 90.01% of the total water head. The 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.Table 6
Inversion results of permeability coefficients by IQPSO (cm/s).
Inversion materialRange of permeability coefficientValue of permeability coefficientRelative impermeable layer3 × 10−6∼3 × 10−42.58 × 10−4Aquitard layer6 × 10−6∼6 × 10−44.05 × 10−4Weak weathered layer8 × 10−5∼8 × 10−37.91 × 10−3Primary rock-fill1 × 10−1∼1 × 1012.03 × 10−1Impervious curtain3 × 10−6∼3 × 10−43.39 × 10−5Figure 9
Contour map of the maximum section of the dam body and dam foundation.


## 4.4.1. Mathematical Model

The 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:(10)minFk=∑i=18wihic−him2,(11)s.t.aj≤kj≤bj,j=1,2,…,5.In the above model, equation (10) is the objective function based on the least-squares criterion. k=k1,k2,⋯,k5T is a group of permeability coefficients to be inverted. hic and him are the ith computed and measured hydraulic head, respectively, and wi is the ith weight which adds up to 1 (wi (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. aj and bj are the jth minimum and maximum of kj, respectively.Next, based on the IQPSO algorithm mentioned in Section2, 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. The particle swarm N, the maximum iteration number T, the population diversity threshold ξ, the mutation probability Pm, the acceleration constants c1 and c2, and the initial and final inertial weights are set as 30, 200, 0.35, 0.5, 2.0, 2.0, 0.9, and 0.4, respectively. The 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. Then, 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.


## 4.4.2. Back Analysis Results

The optimal fitness curves for the three algorithms are shown in Figure8, 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 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.Figure 8
The optimal fitness curves for three algorithms.Table 5
Results of the simulated and measured hydraulic head for three algorithms.
AlgorithmOsmometer numberP11P12P13P15P16P17P18P19IPSOMeasured value (m)808.06806.26805.82805.43805.56804.24804.76802.96Simulated value (m)808.65807.18805.62804.30805.19806.02804.95803.56Absolute error (m)0.590.92−0.20−1.13−0.361.780.190.60Relative error (%)3.275.63−1.29−7.30−2.3412.511.264.66QPSOSimulated value (m)808.54807.06805.50804.20805.09805.92804.86803.56Absolute error (m)0.480.80−0.32−1.23−0.471.680.100.60Relative error (%)2.664.94−2.02−7.98−2.9911.810.674.66IQPSOSimulated value (m)808.41806.67805.43804.17805.06805.66804.84803.55Absolute error (m)0.350.41−0.39−1.26−0.501.420.080.59Relative error (%)1.932.52−2.49−8.16−3.199.970.544.56Compared 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. This 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.The permeability coefficients acquired from IQPSO are listed in Table6. The 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. The hydraulic head is reduced by 78.10 m, accounting for 90.01% of the total water head. The 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.Table 6
Inversion results of permeability coefficients by IQPSO (cm/s).
Inversion materialRange of permeability coefficientValue of permeability coefficientRelative impermeable layer3 × 10−6∼3 × 10−42.58 × 10−4Aquitard layer6 × 10−6∼6 × 10−44.05 × 10−4Weak weathered layer8 × 10−5∼8 × 10−37.91 × 10−3Primary rock-fill1 × 10−1∼1 × 1012.03 × 10−1Impervious curtain3 × 10−6∼3 × 10−43.39 × 10−5Figure 9
Contour map of the maximum section of the dam body and dam foundation.


## 5. 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 mechanism-based quantum particle swarm optimization algorithm was proposed in this paper to solve the inversion problem. IQPSO was tested by numerical experiments. Then, 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. The 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. The numerical experimental results of IQPSO were obviously better than the other three algorithms, which proved the effectiveness and practicability of the algorithm.(2)
The 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. The maximum absolute error by IQPSO was 1.42 m, and the maximum relative error was 9.97%. The analysis results of the steady seepage field indicated that during the normal operation of the reservoir, the antiseepage system composed of “Face slab-toe plate-impervious 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.



---

*Source: 2191079-2020-06-22.xml*