Analysis of a Mathematical Model for Drilling System With Reverse Air Circulation by Using the ANN-BHCS Technique

In the present article, mathematical analysis of drilling system with reverse air circulation is presented by a novel hybrid technique of feedforward artificial neural network (ANN) and biogeography based cuckoo search (BHCS) algorithm. A series solution is constructed with unknown weights for the differential equations representing the drilling problem. Five numerical cases are analysed to show the effectiveness of our method for the solution of differential equations. From the experimental outcomes, it is investigated that our soft computing procedure has a better rate of convergence to the best solution as compared to state-of-the-art techniques. From solution graphs, it is established that our results are in agreement with the reference solutions. It is noted that our technique is easy to implement and can be used for any mathematical model containing nonlinear differential equations. The graphical abstract of this article is given in <xref ref-type="fig" rid="fig1">Figure (1)</xref>. <fig id="fig1" position="float" orientation="portrait"><label>FIGURE 1.</label><caption><p>Graphical illustration of the soft computing procedure followed in this paper.</p></caption> <graphic position="float" orientation="portrait" xlink:href="sulai1-3107405.eps"/> </fig>


I. INTRODUCTION
In the technology of drilling air or different gases are used evacuation fluids to bring the cutting particles formed as a result of drilling at the lowest end of the borehole out to the surface. The technology is in use since the 1950s [1]. It improves the penetration rate in hard formation [2] and points out the problems of loss during circulation in naturally fractured and depleted reservoirs [3]- [6]. There are two kinds of techniques that are used in air and gas drilling. In the field of petroleum engineering; direct and reverse circulation is used for drilling. In straight boreholes, the reverse circulation technique is used for drilling. Compressed air is utilized as The associate editor coordinating the review of this manuscript and approving it for publication was Aniruddha Datta. a mean of drilling fluid in this method, see Figure. (3). The compressed air is injected into the top of the annulus between the two walls of the drilling pipe, the air flows in annulus towards the lowest end of the borehole, the cutting particles at the lowest end of the borehole enter into compressed air and the compressed air brings the cutting particles up to the surface [7]. One of the advantages of this method is that the rock cuttings are continuously returned to the surface. It means that the system of reverse air circulation drilling is a closed system. Therefore, the volume of injected compressed air is proportional to the velocity of the air [8]. The efficiency of reverse circulation of air in drilling depends on whether the cutting particles can be shifted or not from the bottom to the central channel then up to the surface. The air's velocity must satisfy the demands of continuous returning of cutting particles to clean the lowest end of the borehole. Therefore, it is very useful to discuss the velocity of the cutting particles during the drilling process.
The motion of the cutting particles in direct air circulation and gas drilling has already been studied by researchers [9]- [13]. In this drilling technique, the cutting particles are brought up to the surface from the bottom through the annular region between the walls of the well bore and the drill string. The cutting particles weaken the wall of the well bore that causes the instability of well bore [14]- [19]. The demand of airflow rate is large to keep clean the lowest end of the borehole as the cross-sectional area of the annular region is large that increases the drilling cost. The method of reverse air circulation is free from these problems as in this method the cutting particles are carried out to the surface through the drill pipe of dual walls. In this method, the velocity of air is larger than the velocity in other methods with the same volume of air injection because the dual wall drill pipe has a smooth surface and small crosssectional area.
However, there still exist some problems. In the method of reverse air circulation, the drilling bit has some particular structure. Some researchers, like Bulroc and Numa, have primarily focussed on external jetting nozzles [20], [21]. Meanwhile, some studies have concentrated on bottom jetting nozzle bits [22]. The velocity of air is high because it flows through the nozzles of the bit to the lowest end of the borehole. Initially, the cutting particles sit still, and then with the help of aerodynamic drag force it starts to move. When the sliding frictional force controls the aerodynamic force, the particles remain motionless. In a situation, when the particles of compressed air are too large, its velocity is very low. It is impossible to change the bit's structure in order to obtain the smaller sizes of the particles and it is difficult to increase the volume of injected air to increase the velocity of the air, because the air compressor has limited workability. Thus, the breaking of the particle is repeated until its size is small enough to carry it out to the surface. Repeating the breaking of the particle wastes a lot of energy.
If the bits are of different structures then the size of the particles will be different [23]. If the particles have different sizes then the velocity of air, to keep clean the lowest part of the borehole, will also be different. Thus, the volume of air injection will be different. The compressibility of the air because of an increase in the depth of the borehole gradually decreases the volume of air injection [24]. In addition, the volume of the air also changes gradually as it returns. As the air velocity is changed, it affects the capacity of compressed air used to carry the cutting particles out to the surface. Hence, it is beneficial to analyze the motion of cutting particles to design the bit's structure and volume of air injection. These factors can help keep the lower end of the borehole clean and avoid the repetition in the breaking of particles.
The dilute phase pneumatic conveying process is involved in the motion of the cutting particles [8], [25], [26]. In the operations of reverse air circulation, single-phase flow involves the fluid flow from the top of the annular region to the lowest end of the borehole. Two-phase flow occurs when the fluid flow entrains the cuttings formed at the end of the borehole.
Hence, various soft computing approaches have been implemented for the solution of different real-life problems [27]- [38], [38]- [50]. Furthermore, analog active filter design, mathematical models of orthopedic implants, CMOS ring oscillation model, analysis of drilling parameters based on Pareto optimality are discussed and investigated in [51]- [55]. Mathematical models of drilling processes are analysed by many researchers by applying soft computing techniques such as Genetic Algorithm (GA), Ant Colony Optimization (ACO), Particle Swarm Optimization (PSO), Cuckoo Search (CS), Tabu Search (TS), Intelligent Water Drop (IWD), Swarm Intelligent (SI), Artificial intelligence (AI), Firefly Algorithm (FA) and some modified and hybrid approaches. A review of the soft computing approaches used for calculating solutions to drilling problems is given in Table (1). Exploration and exploitation are two essential characteristics of any metaheuristic for a balanced search in the search domain. Many of the researchers have proposed hybrids of these heuristics for a balanced search. Most of these hybrids fail to tackle problems involving complex objective functions with no prior information about their landscapes. Exploration and exploitation are two essential characteristics of any metaheuristic for a balanced search in the search domain. Many of the researchers have proposed hybrids of these heuristics for a balanced search. Most of these hybrids fail to tackle problems involving complex objective functions with no prior information about their landscapes [56].
Keeping in view the work done on the analysis of the drilling problem, the authors of this manuscript were interested in designing a better soft computing approach that can handle mathematical models representing real-life problems with ease of implementation and minimum effort with less arbitrary parameter settings. This paper has studied the mathematical modeling of the velocity of cutting particles in a radial direction. The effects of air velocity and the size of the particle are considered in this model. The model was solved analytically by Zhu et al. [24]. We have designed a hybrid of Artificial Neural Networks (ANNs) and BHCS algorithm, named the ANN-BHCS algorithm. The ANN-BHCS algorithm solves the model, and the results are compared with other state-ofthe-art approaches.

II. BASIC ASSUMPTIONS
The lowest end of the borehole is a confined circular space because the drill bit has a particular structure. Along the circumferential direction of the bottom of the borehole, various breaking points are found whose number depends on the number of the cutting teeth of the bit. The height from the bit's cutting face to the end of the borehole is not longer than the cutting tooth, and at the breaking point, a single particle of cutting rocks is accommodated. The Particle that is driven by the flow of air and slowed down due to sliding frictional force moves from the point where it breaks to the inlet of the central channel. As the bit operates, the particles of the cutting rocks are continuously stored in the inlet and move up from the central channel to the surface. That motion of the particles from a point where they break, to the inlet is called motion in the radial direction.
On the basis of technical features of the drilling with the reverse circulation of air, a single particle of the cutting rocks in the radial direction is considered as a research object. Perfect gas law can be used to approximate compressible air. Furthermore, in the drilling with straight boreholes, the cutting particles are uniformly distributed in the compressed air when they get into compressible air [26]. The rotation and the temperature of the drill pipe do not affect the flow field of the borehole. Moreover, the terminal effects of the particles and the drilling structures are not considered [57], [58].

A. EQUIVALENT DIAMETER
The particles of cutting rocks have different shapes on the basis of some factors such as the bit's structure, condition of formation and rotary drilling circulation medium, etc. It is difficult to derive the mathematical form of the particle's motion, so we assume all the particles to be spherical. The particle's diameter based on the equivalent volume is calculated as: where d sv and V s represent the diameter and volume of the particle, respectively.

B. MODELLING THE PARTICLE'S MOTION IN RADIAL DIRECTION
Considering single particle in the radial direction, we concentrate on the case where drag force controls the sliding frictional force and the particle moves in the radial direction from the breaking point into the central channel (Fig.2). In the radial direction, the aerodynamic drag force F d is where C D , A r , ρ g , v g , and v s are drag coefficients, the projected area of the particle, the density of air, velocities of air and particle, respectively. The particle weight G s is given as: where m s , ρ s , and g denote the mass of the particle, density of the particle, and gravitational acceleration respectively. The frictional coefficient and weight of the particle determine the sliding frictional force. Here, the mathematical form of frictional force f r is given as: where µ r is the coefficient of friction of cuttings at the lowest end of the borehole. Using Eqs. (2 and 4), we obtain the relationship of the critical diameter d cd of the particle and  the critical velocity v cv of compressed air at the end of the borehole as Rearranging Eq. (5), v cv can be obtained as Eq. (6) describes the critical condition for the motion of the cutting particles. This condition means that the particles with smaller diameters than critical diameters are brought by air and the particles with bigger diameters than the critical diameters are broken again and again until their diameters become smaller than the particle's critical diameter. We discuss the motion of the particles having diameters not bigger than the critical diameters of the particles. According to Newton's second law, the motion for the flow of entraining VOLUME 9, 2021 cutting is given by Substituting Eqs. (2)-(4) in Eq. (7).
We assume that the end of the borehole is big enough. The particle is accelerated by giving the aerodynamical drag force. At last, the particle shows suspension movement with a high velocity. Let the aerodynamical drag coefficient is C D and the particle's suspension velocity is v sv , then the result is The drag coefficient depends on Reynolds number and follows Newton's law of resistance. The Reynolds numbers Re r and Re sv correspond to C D and C D respectively, and the difference between their values is very small. Hence, according to Newton's law of resistance, C D and C D is be written as where η g d sv and η g is the air's coefficient of viscosity; k represents the exponent in Newton's law of resistance, whose value is in the range of 0 to 1.
From Eqs. (9) and (11), C D is substituting Eq. (12) into Eq. (8), the particle's motion in radial direction is given as The dimensionless variables and parameters that are used in Eq. (13) are given by Therefore, Eq.(13) becomes In drilling with reverse air circulation, the unit of diameter is in millimeter, so the value of k in Eq. (15) is 0, and Eq. (15) becomes Eq. (16) represents the mathematical model of the particle's motion in the radial direction.

IV. REFERENCE SOLUTIONS
The particles of the cutting rocks are static at the end of the borehole, hence the values of t and v s are 0. Using initial condition, Eq. (16) is solved as given below [24]:

V. THE ANN-BHCS APPROACH
In this paper, we have solved a problem that arises in the field of petroleum engineering. The mathematical model of the velocity of cutting particles in the radial direction during the process of reverse circulation of air in drilling is considered. An ANN based approach is developed and log-sigmoid is taken as activation function. The training of the weights is performed by a hybrid of biogeography based optimization and cuckoo search algorithm (BHCS). We have named our approach as ANN-BHCS algorithm.

A. CONSTRUCTION OF ANN MODEL
The general form of the approximate solution of ODEs of reverse circulation air drilling is considered in our research and the nth derivative of the solution is given as Eq.(18) [79]: here α i , ω i and β i are unknown weights whose values are real and are bounded, f is used for activation function and j is used for the number of neurons in our ANN model. Log-sigmoid is selected as an activation function for ANN which is given as: Approximate solution and its first derivative for our proposed ANN model is given as:

B. OBJECTIVE FUNCTION
The fitness function for the problem is a sum of the mean squared errors E1 and E2 which is given as: where E 1 is related to ODE and E 2 is related to the initial or boundary conditions. For Eq.(16), E 1 and E 2 are given as: In Eqs. (21,22), α, ω and β are the weights that can be adjusted in order to minimize E 1 and E 2 such that E approaches to 0. Hence, the approximate solutionφ of the problem will be near to the exact solution ϕ.

C. CUCKOO SEARCH
Inspired by the cuckoo bird's breeding behavior, a metaheuristic algorithm was developed which is called the cuckoo search algorithm [80]. The female bird lays eggs in the nests of other birds and they unintentionally raise her brood. When the host bird finds the cuckoo's egg in the nest, it either starts making her brood elsewhere or throw the egg out of her nest [81].
In the cuckoo search algorithm, the egg of the host bird shows a solution, and the egg of the cuckoo bird shows a new candidate solution. The CS algorithm works according to the three rules that are as follows [82]: (1) at a time, every cuckoo bird lays only one egg and put it in the host's nest; (2) those nests which have eggs of high quality which are better solutions will go to the next generation, and (3) the number of hosts' nests is fixed, and the host bird can find an alien egg with certain probability. Assuming x iD ) as the position for the ith egg (solution) then updated solution x i new is generated by Levy flights as given below: where the product is entry-wise multiplication; the Levy flight exponent is denoted by β; the step size for a cuckoo is determined by a positive parameter α; the best solutions in the current population are denoted by x g ; u and v are used for random numbers: where is used for Gamma function, and β controls the value of σ u . There is a discovery operator in CS which is used to replace the discovered nests according to the probability pa.
To update the solution, the following equation is used:  where x ij new is the jth element of the ith solution x i new ;x r1,j and x r2,j are the jth elements of the two solutions x r1 and x r2 , where r1 and r2 are two different integers in the interval [1,NP], where NP represents the population size, discovery probability is denoted by pa, P and rand are some random numbers belong to the interval [0, 1].

D. BIOGEOGRAPHY-BASED OPTIMIZATION
Biogeography based optimization (BBO) is an evolutionary algorithm that is inspired by different characteristics of species living in the islands [83]. In BBO, each habitat is considered as a candidate solution having some habitat's suitability index (HSI), which is employed for the measurement of the quality of a habitat. A habitat (solution) is represented by some suitability index variables (SIV). Two types of operators are used in BBO, migration, and mutation that are employed for the evolution of the population. In the migration process, the solutions with high HSI share their features with the solutions having low HSI and the solutions with low HSI accept new features from the solutions with high HSI.
In BBO, the population is randomly initialized with NP habitats (solutions). Each generation sorts the population from the best to the worst and each habitat is assigned with immigration and emigration rates λ and µ respectively: here I and E are used for immigration and emigration rates respectively where I = E = 1; S i represents the number of species in population and S i = NP − i. Accordingly, for the best solution, the S i value is NP − 1 and for the second best solution the S i value is NP − 2 and for the worst solution, the S i value is 0. The migration mixes the features within the population that modifies the solutions. After migration, to modify the solutions, BBO also uses the mutation operator.

E. BBO BASED HETEROGENEOUS CUCKOO SEARCH ALGORITHM
CS and BBO are hybridized because CS uses the Levy flights to modify the solution which is good at exploration and BBO employs the migration operator to modify the solution which is good at exploitation. Combining the exploration and exploitation, a hybrid metaheuristic is developed which is known as BBO based heterogeneous cuckoo search (BHCS) algorithm. The proposed BHCS algorithm has two main steps that are heterogeneous cuckoo search and biogeography based discovery. The two steps are explained in the next sections.

1) STRATEGY OF HETEROGENEOUS CUCKOO SEARCH
In first step, BHCS uses the Levy flights and quantum mechanism based heterogeneous cuckoo search. This strategy is inspired by quantum mechanism and was first presented in [81], [84]. The rules to update the solutions by the heterogeneous cuckoo search are given as follows [81], [84]: where L = δ ln(1/η), ε = δ exp(η), x g is used for solution that is best at current iteration;x = 1 NP NP i=1 x i represents the mean of the solutions; sr and η are randomly selected numbers belong to the interval [0, 1]. Eq. (31) shows that heterogeneous cuckoo search employs three equations to update the solutions with the same probabilities. The first equation is based on Levy flights in original CS and the second and third equations to update the solutions are based on quantum mechanism. Updating the solutions using heterogeneous rules diversify the search space and follows the direction towards the real global region.

2) BIOGEOGRAPHY-BASED DISCOVERY OPERATOR
In the second step, New solutions are generated using a discovery operator. When the host bird finds an alien egg with probability pa, it abandons the old nest and makes a new nest on the basis of the migration operator. Initially, solutions are sorted from best to worst, then an immigration rate µ is assigned to each solution: where E represents the emigration rate having maximum value of 1; S i = NP − i represents the number of species in solutions. In biogeography based discovery operator, the solutions having best fitness share their characteristics with other solutions which help to enhance the exploitation.

F. OVERALL BHCS ALGORITHM
The BHCS algorithm uses a cascading structure for the implementation of its two search steps. The coordination between the search strategies of the heterogeneous cuckoo search and biogeography based discovery operator can efficiently balance the exploitation and exploration.

VI. RESULTS AND DISCUSSION
In this paper, we have implemented the ANN-BHCS algorithm for the solution of the ODEs that model the motion of particles in the radial direction during the drilling process with the reverse circulation of air. The motion of the cutting particles in the radial direction is discussed. The solutions obtained by Zhu et al. [24] are considered as reference solutions. To compare the results of the ANN-BHCS algorithm, we have solved the problem by genetic algorithm (GA), particle swarm optimization (PSO), bat algorithm and multiverse optimizer (MVO). Four statistical operators are employed to test the efficiency of the designed algorithm. These operators are MAD, TIC and ENSE. Mathematical forms of these operators are given as follows: here n denotes the number of grid points in approximate solution.

A. PARTICLE'S VELOCITY IN RADIAL DIRECTION
In the model of the radial velocity of the particles, we have considered five cases on the basis of the values of µ r and β. The initial condition imposed on these cases is ϕ(0) = 0 because at the bottom of the borehole, the particles are at rest and their velocity is 0. Each hidden layer in the ANN architecture has 10 neurons with 30 unknown weights. The step size is taken as h = 0.2, and there are 31 grid points in the entire domain that is [0,6]. We have simulated the ANN-BHCS algorithm with 100 runs to assess the reliability and convergence of our approach.

1) CASE 1
In this case, the values of µ r = 0.2 and β = 0.2. Using the values, Eq.(16) becomes: using Eq.(37), the minimization objective function for this case is given by: We have taken ten neurons in our ANN model. The approximate solutions obtained using the ANN-BHCS algorithm contains 10 terms, one neuron corresponds to one term of solution. Convergence of the fitness values for 100 runs is given in Figure(6a).
The minimum fitness value obtained for this case is 1.1977E − 08. Weights obtained by the ANN-BHCS algorithm to minimize the fitness function are plotted in Figure(5a). Using these weights, the series solution for case 1 is given as follows:    The minimum values of the fitness function, MAD, TIC and ENSE are 1.2974E − 08, 7.8026E − 05, 4.4484E − 05 and 5.5261E − 08 respectively, which shows that the solutions obtained by the ANN-BHCS algorithm are very close to the analytical solution.
Convergence analysis of performance metrics for case 1 is given in Table (5). Histograms and box-plots for performance metrics and fitness function are given in Figure ( Histogram and box plots of MAD values show that more than 80% of the values are less than or equal to E − 03. The box plot of MAD values shows that more than 75% of the values are in the range E − 02 to E − 06. Similarly, the histogram of TIC values shows that more than 90% of the values are less than or equal to E − 03 and the box plot of the TIC values shows that more than 75% of the values are in the range E − 03 to E − 06.
Histogram of the ENSE values shows that more than 90% of the values are less than or equal to E − 03 and box plot shows that 75% of the values are in the range E − 05 to E − 08. The above discussion justifies the efficiency of the ANN-BHCS algorithm for the solution of ODEs.

2) CASE 2
In this case, the values of µ r = 0.2 and β = 0.6. Using these values, Eq.(16) becomes using Eq.(40), the minimization objective function for this case can be written as: The fitness value obtained for this case is 4.2369E − 10. Convergence of the fitness values for 100 runs is given in Figure(9a). Weights obtained to minimize the fitness function are given in Figure( Numerical solutions obtained for this case are given in Table (7). The solutions are also plotted in Figure (8b), which shows that the solution obtained by the ANN-BHCS algorithm is very close to the solution obtained by [24].
The accuracy of the solution obtained by the ANN-BHCS algorithm is also clear from the comparison of absolute errors in Table ( Figure (9c) shows the values of performance metrics and fitness function. Table (9) and Figure (9d) illustrate the values of performance metrics in terms of best, mean, and worst. Convergence analysis of performance metrics for 100 runs is given in Table (10).
Histograms and box plots of performance metrics and fitness function are given in Figure (10). Histogram of the fitness values shows that 100% of the values are less than or equal to E − 03 and box plot of the fitness values shows that more than 75% of the values are between E − 06 and E − 10.
Histogram of the MAD values shows that more than 95% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E −03 to E −06. Histogram of the TIC values shows that 100% of the values are less than or equal to E −03 and box plot shows that more than 75% of the values are in the range E −04 to E −06.   Similarly, a histogram of the ENSE values shows that more than 95% of the values are less than or equal to E − 03 and box plot shows that more than 75% of the values are between E − 06 to E − 10. The above results show that ANN-BHCS algorithm can efficiently solve the ODE of case 2.

3) CASE 3
In this case, the values of µ r = 0.2 and β = 1. Using these values in Eq.(16) we obtain: using Eq.(43), the minimization fitness function for this case is given by: There are 10 neurons in the proposed ANN model. The solution obtained using the ANN-BHCS algorithm contains 10 terms corresponding to the number of neurons. The fitness value obtained for this case is 2.9009E − 11.
The convergence of the fitness values for 100 simulations is plotted in Figure( Numerical solutions obtained by the ANN-BHCS algorithm and other algorithms for case 3 are given in Table (12). The solutions are also plotted in Figure (11b).
Absolute errors of the ANN-BHCS algorithm and other algorithms are compared in Table ( 13) and Figure (11c). From Table (13) and Figure (11c), we can see that ANN-BHCS algorithm gives a solution with minimum absolute errors.
Statistical analysis of the absolute errors is given in Table ( The values of performance metrics and fitness function in terms of best, mean and worst are presented in Table ( 14) and are also shown in Figure (12d). Convergence analysis of fitness values and performance metrics are given in Table (15).
Histograms and box plots of fitness values and performance metrics are presented in Figure (13). Histogram of the fitness values shows that 100% of the values are less than or equal to E − 03 and the box plot of the fitness values shows that more than 75% of the values are in the range E − 07 to E − 11. Histogram of the MAD values shows that more than 95% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 04 to E − 06. Histogram of TIC values shows that 95% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are between E − 05 and E − 07. Histogram of the ENSE values shows that more than 95% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 07 to E − 11. The results of case 3 show the accuracy and efficiency of the ANN-BHCS algorithm.

4) CASE 4
In this case, the values of µ r = 0.6 and β = 0.2. Using these values, Eq.(16) becomes using Eq.(46), the minimization fitness function for this case is written as: The fitness value obtained for this case is 5.2611E − 08. Convergence of the fitness values for 100 runs is given in Figure(15a). Weights obtained to minimize the fitness function Numerical solutions obtained by ANN-BHCS algorithm and other algorithms for case 4 are presented in Table( 17). The solutions are also plotted in Figure(14b). Absolute errors of the ANN-BHCS algorithm and other algorithms are compared in Table ( 18) and Figure (14c) which shows that the ANN-BHCS algorithm can obtain a solution with less absolute errors than other algorithms.
Statistical analysis of absolute errors of the ANN-BHCS algorithm is given in Table(  to E − 03 which shows the accuracy in the solutions obtained by the ANN-BHCS algorithm. The values of performance metrics for 100 runs are plotted in Figure(15c). The minimum, mean and maximum values of the performance metrics and fitness function are given in Table( 19) and are also plotted in Figure(15d). Convergence analysis of the fitness values and performance metrics is presented in Table ( 20).
Histograms and box plots of the values of the fitness function and performance metrics are given in Figure(16).
Histogram of the fitness values shows that more than 80% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 03 to E − 08. Histogram of the MAD values shows that over 90% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the E − 03 to E − 05. Histogram of TIC values shows that more than 90% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 04 to E − 05. Histogram for ENSE values shows that more than 95% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the VOLUME 9, 2021 values are between E − 05 and E − 08. These results for case 4 show the efficiency of the ANN-BHCS algorithm.

5) CASE 5
In this case, the values of µ r = 1 and β = 0.2. Using these values, Eq.(16) becomes using Eq.(49), the minimization objective function for this case can be written as: Numerical solutions obtained by the ANN-BHCS algorithm and other algorithms for case 5 are presented in Table( 22) and solutions are also plotted in Figure(17b).
The absolute errors of the ANN-BHCS algorithm are compared with other algorithms in Table ( 23) and Figure (17c). The absolute errors of the ANN-BHCS algorithm are less than that of other algorithms which shows that solution obtained by the ANN-BHCS algorithm is better than other algorithms. Statistical analysis of absolute errors of the ANN-BHCS algorithm is given in Table(     The values of performance metrics for 100 runs are plotted in Figure(18c). The minimum, mean and maximum values of the performance metrics and fitness function are presented  in Table( 24) and plotted in Figure(18d). The convergence of fitness values and performance metrics is given in Table ( 25).
Histograms and box plots of fitness values and the performance metrics are given in Figure(19). Histogram of the fitness values shows that more than 80% of the values are less than equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 04 to E − 08.
Histogram of the MAD values shows that 90% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E −04 to E −05. Histogram of the TIC values shows that more than 90% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 04 to E − 05. Histogram of the ENSE values shows that more than 90% of the values are less than or equal to E − 03 and the box plot shows that more than 75% of the values are in the range E − 06 to E − 08. The above results show the efficiency of the ANN-BHCS algorithm.
The fitness value obtained for this case is 5.2611E − 08. The convergence of the fitness values for 100 runs is given in Figure(18a). Weights obtained to minimize the fitness function are given in Figure(17a). Using these weights, the series solution for this case is given as:  Table ( 30) and Figure (21b). The results show that the solution is getting better as the population size increases.

VIII. CONCLUSION
We have used a hybrid of biogeography based heterogeneous cuckoo search algorithm and artificial neural networks to solve ODEs that model the particles' velocity in a radial direction in reverse circulation air drilling. The results obtained by the ANN-BHCS algorithm are compared with the analytical results given in [24] and other algorithms. We have considered five cases of the problem based on the different values of β and µ r . The ANN-BHCS algorithm's efficiency is checked by calculating the absolute errors, MAD, TIC, and ENSE, for five cases. The analytical solutions, which are in terms of log-sigmoid function for all the cases, are given in Eqs. (39,42,45,48) and (51). Comparison of numerical solutions for the five cases is given in Tables. (2,7,12,17) and (22). Absolute errors for all the cases are given in Tables (3,8,13,18) and (23). The absolute errors of ANN-BHCS algorithm are from E − 03 to E − 07 for case 1, from E − 05 to E − 07 for case 2, from E − 06 to E − 09 for case 3, from E − 04 to E − 06 for case 4 and from E − 03 to E − 07 which shows the accuracy of the solutions obtained by the ANN-BHCS algorithm. The accuracy of the results can also be seen from the good agreement of the ANN-BHCS algorithm solutions with analytical solutions which are given in Figs. (5b,8b,11b,14b) and (17b). The efficiency of the ANN-BHCS algorithm is obvious from convergence analysis of values of performance metrics which are given in Tables (5,10,15,20,25) and Figs. (7,10,13,16) and (19). The results show that the algorithm can efficiently  minimize the fitness function and obtains the best solution as more than 90% values of the fitness function and performance metrics are less than E − 03.  Velocity of air. v s Velocity of particle. G s Weight of the particle. m s Mass of the particle.      He is also the Director of the Computational and Applied Science for Smart Innovation Cluster (CLASSIC Research Cluster), KMUTT. He has authored or coauthored more than 800 international peerreviewed journals. His main research interests include fixed point theory and applications, computational fixed point algorithms, nonlinear optimization, control theory, and optimization algorithms.
MAHARANI ABU BAKAR received the Ph.D. degree in mathematical sciences from the University of Essex, in 2015. She is currently working as a Senior Lecturer with the Department of Applied Mathematics, Universiti Malaysia Terengganu (UMT). She is also a member of the Applied Mathematics and Scientific Computing Research Interest Group, UMT. She has authored and coauthored more than 20 research articles and conference articles. Her research interests include numerical analysis and its relations with the solutions of problems, which involve large scale of variables, parallel computing, and cloud computing. Recently, she has been working on solving the high-dimensions of partial differential equations solutions by using machine learning, particularly using support vector regression (SVR) and artificial neural networks (ANNs) models.
MIFTAHUDDIN received the bachelor's degree from the University of Hasanuddin, South Sulawesi, the master's degree from Bandung Institute of Technology, West Java, and the Ph.D. degree from the University of Essex, Colchester, U.K. He is currently working as a Lecturer, a Researcher, and an Entrepreneur with the Statistics Studies Program, Department of Statistics, Faculty of Mathematics and Sciences, Syiah Kuala University, Banda Aceh. His areas of expertise are statistical modeling, machine learning, data analysis, big data, data sciences, biometrics, time series and regression analysis, semiparametric models, stochastic processes, actuarial sciences, financial analysis, and geostatistics (time-spatial). His supervisor was Prof. Berthold Lausen, the Head of Mathematical and Statistical Sciences, University of Essex.