Research on Horizontal Directional Drilling (HDD) Trajectory Design and Optimization Using Improved Radial Movement Optimization

Abstract

In practice, the drilling path of horizontal directional drilling (HDD) projects is usually constructed by trial and error based on a preliminary designing trajectory. This study aimed to propose and test a method to predesign and optimize the drilling paths automatically, with the view of improving the efficiency of HDD design preparations. Alternating straight and curvilinear segments is a commonly used method for designing drilling paths, especially the “straight—curvilinear—horizontal straight—curvilinear—straight” five-segment arrangement. The catenary method was proposed to design the drilling path with the advantage of lower friction for the mechanical constraints. However, it is difficult to be implemented with technology limitations due to its continuously changing curvature. In this study, five-segment trajectories were combined with the catenary trajectory to utilize their advantages using the improved radial movement optimization (IRMO) algorithm. Drilling mud pressure was considered in the processes of the mechanical design to avoid collapse or possible instability. Two different examples were tested in different scenarios, theoretical and practical. The results show that the IRMO algorithm has a great potential for automatically designing and optimizing preliminary drilling paths with low time-consumption and high feasibility.

1. Introduction

Horizontal directional drilling (HDD) is a trenchless technology applied to install underground pipelines with minimal impacts on the environment or damage to existing infrastructure such as roadways and other surface structures [1]. According to the predesigned drilling trajectory, a small pilot hole is drilled first, and then the pilot hole is enlarged by a reamer that replaces the drilling head. The pipes or installations of the well are pulled back into the reamed hole from the exit toward the entry. Due to its lower cost, higher flexibility, and weaker surrounding influence, HDD has gradually become more popular than the conventional excavation methods adopted to install and replace or repair pipelines in a city [2]. The development of HDD in China has continued in the past two decades due to the broad market and stable economic growth of China. Up till November 2020, the longest known HDD in China had reached 5.2 km [3,4]. The geological conditions of the crossing area are much more complex in longer-distance crossing HDD practices, which results in more engineering problems and costs. In practice, most of the HDD drilling paths are determined by seeking or tracking the moving or static drill head to adjust it from the entry point to the exit simultaneously and iteratively [5]. Scholars have carried out various studies on improving the accuracy or efficiency of the seeking or tracking systems [6,7,8,9]. However, this method of iterative trial and correction is difficult to apply in complex practices, such as crossing rivers or buildings which may hinder the sensors. On the other hand, many scholars have focused on the predesign and optimization of HDD trajectories to ensure suitability and efficiency. For conventional trajectory design, the theoretical geometric model of trajectory was mainly developed to ascertain the accurate position of the wells, such as the “Improved Tangential Method” [10]. Moreover, some scholars added more factors in the model to realize more complex designs closed to actual projects [11,12,13]. In China, Zhou [14] took the drilling rate in meters per hour as the optimization object to establish a dynamic trajectory design optimization model, which could find the optimal drilling path with the minimum drilling time. Lu [15] focused on complex underground pipelines and obstacles and presented an optimization method for avoiding pipelines and obstacles in trajectory design. Niu [16] conducted research on the theory of the guide strength of the drilling path, which analyzed the contribution of the curved segment to the drilling path design. When it comes to artificial intelligence applications, an optimization algorithm is more intuitive and reliable than repetition tests and traditional human-computer interaction techniques [17]. The drilling path design and optimization were also developed further by several optimization algorithms. The catenary trajectory design method was introduced by Wiśniowski [18,19,20], who tried to fit the five-segment trajectory with the catenary trajectory to design the drilling path geometrically. The genetic algorithm (GA) was used to realize the fitting of the two trajectories efficiently. In another study, the ant colony optimization (ACO) algorithm was used to optimize HDD alignment by Patino-Ramirez, et al. [21], which ensured the minimal drill path length and the reduced cost associated with it. Notably, the ACO optimization method added consideration of mechanical constraints (drilling mud pressure and pipe integrity) and geometric constraints (remaining in the construction domain) simultaneously in the design.

The artificial intelligence algorithm is an excellent way to figure out a single or multidimensional nonlinear objective optimization problem which the design or optimization of trajectory can be transferred to. On the basis of previous studies, this paper tested a new numerical method, using improved radial movement optimization (IRMO) to preliminarily design and optimize the HDD trajectories and its parameters. According to the theory of catenary trajectory in previous literature, its shape is closed to the natural stress distribution of pipelines along the length so as to enable a lower drilling or pullback force during the process [18,19,20]. So, a five-segment trajectory was attempted to be realized by fitting with the catenary trajectory to get the largest similarity in shape. What is more, the drilling mud pressure is calculated in each iteration step simultaneously for consideration of the well-bore stability evaluation. The use of the IRMO algorithm can make the whole construction processing automatic with the intention to reduce the time-consuming manual adjustment. Section 2 introduces the concepts of the catenary trajectory design and five-segment trajectory design for HDD projects. Then we present the design method based on the geometrical constraints for the whole construction and each parameter. The objective functions for this method used in the IRMO algorithm are also summarized. Besides, the mechanical design is also introduced for the borehole stability consideration based on drilling mud pressure. Section 3 summarizes the framework and implementation of the IRMO algorithm applied to the trajectory design. Section 4 presents two design examples to compare and analyze the results determined by IRMO. The conclusion of this study is presented in Section 5.

2. HDD Trajectory Designs

2.1. Geometric Designs

The catenary trajectory is a chain curve with natural deflection, such as seen with ropes, cables or chains with uniform weight suspended between two points (Figure 1a). In petroleum engineering, the well path in the catenary profile has been proved to reduce the wellbore friction and torque compared to conventional trajectories [22,23]. The method for planning the catenary well path was proposed in later studies, which avoided a trial-and-error procedure and provided excellent maneuverability of planning requirements [24]. In 1985, McClendon [25] proposed an HDD trajectory design method based on a catenary trajectory. However, its gradient curvature made it difficult to be implemented in practice. Therefore, Wiśniowski et al. [20] proposed a new method that combined the catenary trajectory with a five-segment trajectory (Figure 1b) to design a drill path that was easier to be implemented using the available technology.

The catenary trajectory is actually a hyperbolic cosine, whose standard mathematical formula is written as the following:

z = a × cosh(x/a)

(1)

a = Np/q

(2)

where, a is a parameter determined by the unit weight of pipeline q (N/m) and pullback force Np (kgf) [20]. So, the shape of the catenary trajectory is mainly associated with the crossing pipelines and the capacity of the pulling back equipment. By Equation (3), the catenary trajectory can be described in a specific coordinate, where the left point of the catenary trajectory is set as the origin point. Therefore, the depth of each point on the catenary trajectory (zc_i) can be fully calculated with increments. As shown in Figure 2, A and H are the horizontal and vertical distances of the exit and the entry points, respectively. (x₀, z₀) is the catenary vertex coordinate.

$$\{\begin{array}{c}{x}_i=i \times \mathit{dx}\\ {\mathit{zc}}_i=(\mathit{Np}/q) \times \mathit{cos}h\left(\right(x_i - x_0)/(\mathit{Np}/q) - 1) - z_0\end{array}$$

(3)

Another concept for HDD trajectory design is the segment trajectory method, combining the most conventional methods of straight and curvilinear segments. Mostly, the geometric profile of segment trajectory (shown in Figure 1b) commonly consists of an entry straight segment (L₁), an entry curved segment (L₂) with bending radius (R₂), a central horizontal straight segment (L₃), an exit curved segment (L₄) with bending radius (R₄), and an exit straight segment (L₅) respectively. Based on this combination form, the point coordinates of each segment can be calculated clearly by geometric constraints. As shown in Figure 3, once the angle γ and radius R of the curved segments are set, the increments of the trajectory depth $\mathit{dz}$ and horizontal distance $\mathit{dx}$ can be determined according to Equation (4). Furthermore, combining the inclination angle θ and length L of the adjacent straight segments, the depth $z_i$ of another adjacent segment can be fully determined with increasing horizontal distance $\mathit{dx}$ (Equations (5) and (6)).

$$\{\begin{array}{c}\mathit{dx}=R \times (\mathit{sin}\gamma - \mathit{sin}(\gamma - d\gamma \left)\right)\\ \mathit{dz}=R \times (\mathit{cos}\gamma - \mathit{cos}(\gamma - d\gamma \left)\right)\end{array}$$

(4)

$$\alpha =\theta +\gamma$$

(5)

$$\{\begin{array}{c}{x}_i=i \times \mathit{dx}\\ z_i=L_1 \times \mathit{sin}\theta +R \times \left(\mathit{cos}\theta - 1\right)+\frac{x_{i }{- L_1}_{ }\times \mathit{sin}\theta - R \times (\mathit{cos}\theta +\mathit{cos}\alpha ) }{\mathit{tan} \alpha }\end{array}$$

(6)

According to Equations (4)–(6) above, the designed alternating straight and curvilinear segments can be detailed with tiny segments with coordinates. The whole geometry can be constrained only by four kinds of geometric parameters within the assigned domain, including the inclined angle of the straight segment θ (clockwise is positive and counterclockwise is negative), the length of segment L_i, the angle of curved segment γ, and the curved segment radius R. If the central segment of the trajectory is considered to be horizontal and straight in design, the geometric parameters curve angle γ and deviation angle θ will be united (θ₁ = γ₁, θ₂ = γ₂), so that the geometrical design parameters of the five-segment trajectory can be simplified and containing the following:

• Length of entry straight segment, L₁
• Inclined angle of entry inclined straight segment, θ₁
• Radius of entry curved segment, R₂
• Central straight segment length, L₃
• Inclined angle of exit inclined straight segment, R₄
• Radius of exit curved segment, θ₂

In practice, the straight segments of drill path, L₁, L₃, and L₅, are generally determined by the gyratory drilling processes of the drill head, while the curved segments are determined by the direct jacking [16]. The entry angle θ₁ (for the entry straight segment L₁) and exit angle θ₂ (and the exit straight segment L₅) of the trajectory are commonly restricted in valid ranges. Patino-Ramirez et al. [21] suggested the ranges could be defined between 5 and 18 degrees for trajectory design, while Wiśniowski et al. [19] recommended it could be set between 6 and 16 degrees. Another study gave higher acceptable values of deviation angles between 8 and 30 degrees [26]. Its determination is usually constrained by constructability limitations which impart conditions of geometric design. The central straight segment L₃ is commonly the core part of the trajectory, which must satisfy the requirements of specific burial depths and crossing distance for the engineering. To avoid surface disturbance and damage to installed products or pipelines, the minimum depth of cover along the trajectory should not be less than the minimum allowable H in the standards.

The bending radius (R₂ and R₄) of curved segments mostly depends on the pipeline strength which is limited by the material properties or diameter. To prevent excessive bending stress in the pipelines, the radius of the curved segments should be limited to the minimum bending radius. More details of the minimum bending radius can be found in the published literatures [16,19,21]. In this study, the bending radius adopted was between 2000 and 3000 m for pipelines with diameter over 1000 mm.

After preliminary investigation, the assigned domains of the locations of entry and exit points are decided upon for the predesigns. Project risks are evaluated as well, which might result in more stringent limitations in the predesigns. Based on these limitations, possible trajectories are designed (fitted) to ensure the maximum similarity with the catenary trajectory in unifying the plane coordinate system by the least squares method. The sum of squares SOS is determined as the similarity measurement function to find the best fit (Equation (7)). z_i and zc_i are the depths of the five-segment trajectory and catenary trajectory determined by Equations (3)–(6) above. The lower the difference between the depths, the higher is the similarity of the two trajectories.

Fitness_S = SOS = Σ(Δz_i)² = Σ(z_i − zc_i)²

(7)

According to the six parameters obtained above, the length of the preliminary design trajectory can be given in Equation (8) with specific geometric constraints.

Fitness_L = ΣL_i = L₁ + R₂ × θ₁ + L₃ + R₄ × θ₂ + (A −L₁ × sin θ₂ −L₃ −R₄ × cos θ₂)/sin θ₂

(8)

2.2. Mechanical Design Considering Drilling Mud Pressure

Taking into account the collapse and instability of the drill paths, we tested, step by step, the feasibility and rationality of each geometrically possible trajectory based on drilling mud pressure theory. The engineering geological conditions of the crossing area for calculation are simplified in this study. The drilling fluid is a non-Newtonian fluid mixed HDD process. It is necessary to exceed the minimum required mud pressure (MRP) to maintain the recirculation of the drilling mud. The minimum mud pressure consists of the static drilling mud pressure Ps and the drilling mud pressure loss Pl [27,28], as shown in Equation (9) [27,28].

MRP = Ps + Pl

(9)

Simultaneously, the injected pressure of the drilling fluid must be controlled below the maximum allowable drilling mud pressure (MAP) to avoid cavity shear failure (blowout) of the drill path [21]. The MAP is calculated by the Delft equation [27,29], which assumes that the soil surrounding the borehole experiences perfectly elastic deformation and that the far-field stress around the borehole is isotropic. Once the deformation reaches the plastic threshold defined by the Mohr–Coulomb criterion (the internal pressure reaches the limited pressure), shear failure of the borehole will occur as a blowout. Although the Delft equation has been continuously revised to obtain more accurate analytical solutions [30], Delft equation-based formulations are still the most commonly accepted formulations and are therefore adopted in this study, Equation (10) [29].

MAP = Plim

(10)

Once the geometrically possible trajectory is determined, the MRP and MAP of the whole construction can be calculated predictably for the mechanical considerations. For all of the steps, the maximum allowable pressure (MAP) of the fluid is predicted to avoid the triggered blowout.

3. Improved Radial Movement Optimization

3.1. Concepts of IRMO

The IRMO algorithm [31,32,33,34], which can quickly solve multidimensional objective functions of nonlinear constraint problems, is a global optimization algorithm that further improves the data structure of the traditional radial movement optimization (RMO) algorithm [34]. The IRMO algorithm simulates a group of particles [X_i,j] moving in a gradually shrinking solution space. In IRMO, one particle position, which contains the information of variables, represents a solution vector referring to a complete HDD trajectory in this study. All of the solution vectors form a solution matrix [X_i,j]. After the expression of the objective function and the range of the variables are determined, a group of initial particles are generated randomly in the IRMO algorithm. The position of each particle is evaluated by the value of the fitness function (the value of the sum of squares SOS or total length L); therefore, the local optimal solution Rbest (radial best) and global optimal solution Gbest (global best) can be obtained by one-by-one comparisons. The particle corresponding to the best value of the function is chosen as the initial center. To avoid the algorithm overreliance on the central particle and the loss of potential good particles of the previous generation, which contribute to deviation from the optimization direction of the global optimal solution, IRMO proposed generating prepositioned particles [Y_i,j]. The fitness values of [Y_i,j] are compared with those of [X_i,j], pre-Rbest and pre-Gbest to evaluate the optimal information that will be updated or persist. As the iteration progresses, the solution space (ranges of variables) will be narrowed so that the central particle, determined by Rbest and Gbest, will move toward the optimum. When the positions of the particles are narrowed to the smallest solution space, the optimum solution is output. The principle of updating the central particle is shown in Figure 4.

3.2. Implementation of IRMO

According to the studies above, one solution of the five-segment trajectory can be determined by an M dimensional vector containing six variables: L₁, θ₁, R₂, L₃, R₄, and θ₂. The N trajectories form a standard matrix [X_N,M] with N rows and M columns (M = 6). The standard matrix [X_N,M] realizes the connection between the IRMO algorithm and the HDD trajectory optimization (shown in Equation (11)).

$$\left[{\mathrm{X}}_{\mathrm{N},\mathrm{M}} \right]=\left[\begin{array}{ccc}\begin{array}{cc}{L}_{1,1}^1& {\theta }_{1,2}^1\end{array}& \begin{array}{cc}{R}_{1,3}^2& L_{1,4}^3\end{array}& \begin{array}{cc}{R}_{1,5}^4& {\theta }_{1,6}^2\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{L}_{N,1}^1& {\theta }_{N,2}^1\end{array}& \begin{array}{cc}{R}_{N,3}^2& L_{N,4}^3\end{array}& \begin{array}{cc}{R}_{N,5}^4& {\theta }_{N,6}^2\end{array}\end{array}\right]$$

(11)

After setting the ranges of each variable and the constraint conditions, optimization begins with the initial matrix [X_N,M] initiated by Equation (12).

X_i,j = X_{min j} + (X_{max j} − X_{min j}) × rand(0, 1)

(12)

where X_{max j} is the upper limit of the j-th variable in [X_i,j], and X_{min j} is the lower limit.

The initial central particle cp¹ is chosen from the initial [X_N,M] which has the best fitness value. The central particles cp^k will be updated by Rbest and Gbest, as shown in Equation (13). c₁ and c₂ are proportional coefficients that affect the convergence speed and accuracy of the calculation. The values of c₁ and c₂ are set as 0.5 and 0.4, respectively, in this study.

cp^k = cp^k−1 + c₁ × (Gbeat^k−1 − cp^k−1) + c₂ × (Rbeat^k−1 − cp^k−1)

(13)

Following Equation (14), prepositioned particles [Y_i,j] are generated by the previous central particles.

Y_i,j^k = cp^k + (X_{max j} − X_{min j}) × rand(−0.5,0.5) × w^k

(14)

w^k = 1.009 × exp(−(1/generation + 0.03307)²/0.2454)

(15)

To calculate the fitness values of [Y_i,j]: If the fitness values of [Y_i,j^k] are better than those of [X_i,j^k], the [X_i,j^k+1] will be updated to [Y_i,j^k], otherwise [X_i,j^k] will persist to the next generation. The determination of inertia weight w^k uses a Gaussian function (Equation (15)).

Based on the objective functions given above (Equations (7) and (8)), the design and optimization of the HDD trajectories in this study were implemented by IRMO on MATLAB 2021b. The implementation of the IRMO algorithm for HDD trajectory optimization is shown in detail in Figure 5.

4. Comparison and Analysis

This study aimed to assess the performance of IRMO in trajectory design by two examples originating from other literature. The first one gives the details of geometric design processes which are compared with the genetic algorithm (GA) to highlight the advantages of IRMO. To ensure the stability of this IRMO method for practicable consideration, the parameters are discussed based on twenty times calculation. Therefore, a further optimization, shortening the total length of trajectories, is proposed and tested to overcome the limitation of unstable parameter design. Another example tests the IRMO designs with the actual design using published data from constructed projects. To test the mechanical design part, predictions of the maximum allowable drilling mud pressure of the geometric design trajectories are compared with other studies.

4.1. Comparison and Analysis of Numerical

The first example [19] is used to test the fitting process of this method using the IRMO algorithm. The catenary trajectory parameters and five-segment trajectory parameters are shown in

Table 1. Parameters and range setting.

Parameters of Catenary Trajectory [19]	Parameters of Five-Segment Trajectory
Horizontal distance of exit and entry point: A = 1000 m Vertical distance of exit and entry point: H = −15 m Casing unit weight: q = 80 kg/m Pullback force: Np = 25,000 kgf	Length of the entry straight segment, L1∈[1, 300] m; Length of the central segment, L3∈[1, 300] m; Radius of curved segments, R2, R4∈[1, 300] m; Angle of the entry straight segment, θ1∈[1, 300] °; Angle of the entry straight segment, θ2∈[1, 300] °;

.

To determine the influence of different algorithm parameters on the fitting results, the value of R² is used to manifest the similarity of fit. As shown in

Table 2. Comparison of goodness of fit for different algorithm parameters.

Series	dx [m]	N	Generation	SOS	R2	Time [ms]
1	50	100	100	10.25	0.9991	765
2	50	50	100	21.01	0.9981	487
3	5	50	100	70.68	0.9991	532
4	5	50	250	57.27	0.9993	1275
5	1	50	250	287.22	0.9995	1680

Note: R² = 1 − SSE/SSR, R²∈(0,1); the closer the value of R² is to 1, the better the goodness of fit.

, the difference in the value of R² is quite minor and is greater than 0.998 in each case. The overlap of the fitted trajectory images for series 1, 3, and 5 is relatively high (Figure 6). Moreover, the time resource of each case increases as the calculation becomes more complex.

The optimal fitting result shows a lower SOS value searched by the IRMO algorithm (in

Table 3. Comparison of optimal fitting results.

Method	SOS	Time [ms]
I: Searching the entire solution space [19]	370.19	4183
II: GA [19]	366.53	212.55
III: IRMO (this study)	287.22	1680

), but more time costing. Compared with the results of the genetic algorithm (GA), the IRMO algorithm can give more accurate calculations with an acceptable extra time costing.

During the 100 iterations, the global optimal SOS tends to a stable convergence when approaching the approximate 50th generation (shown in Figure 7). The optimal result can be improved when the total number of iterations is increased.

As shown in

Table 4. Comparison of parameter optimization results.

Method	L1 [m]	θ1 [°]	R2 [m]	θ2 [°]	R4 [m]	L3 [m]
GA [12]	—	−9.3	2473.96	7.44	2378.51	100.78
IRMO (this study)	—	−8.83	2678.16	6.49	2133.25	100.00

, the optimization result of L₃ in this study is largely similar to that of the GA [19]. The other parameter optimization results have θ₁, θ_2, and R₂ that are all smaller than those of the GA, while R₂ is larger. The difference between the optimization results is attributed to the mutual influence and restriction of each parameter.

Taking Series 4 in Table 3. as an example, the SOS is calculated 20 times consecutively (Figure 8), and the 20 trajectories are shown in Figure 9.

The calculations completed by the IRMO showed great stability and accuracy, with the standard deviation of the SOS being 3.61244, consistent with the trajectory images. In Figure 9, different colors are used to distinguish each segment of the five-segment trajectory. The fitting results of each segment are highly similar as well.

The parameter collocation of the optimal fitting trajectory cannot be identical here. As the fitting similarity is higher, the mutual influences and restrictions of each parameter on the others will be enlarged causing different results. To achieve a high similarity of fitting, the results of the parameters, especially L₁, R₂, and R₄, fluctuate in a small range after 20 calculations (Figure 10, Figure 11 and Figure 12), while the fitting results are not much different (Figure 9). The results of L₃, θ₁, and θ₂ fluctuate relatively consistently, especially the length of the central segment L₃. The optimization of L₃ tends to the lower limit of its initial range (as shown in Table 1) with no difference in the 20 searches, which indicates that a shorter L₃ is beneficial to the trajectory fitting. The results of θ₁ and θ₂ fluctuate only in the range of 0.5° (Figure 10), with standard deviations of 0.29313 and 0.212049, respectively. However, the optimization results of L₁, R₂, and R₄ contrast obviously with the former findings. This shows that the fluctuations in the fitting are less correlated with L₁, R₂, and R₄.

Wiśniowski et al. [19] suggested fixing the length of the central segment L₃ and narrowing the ranges of the other parameters to further optimize the five-segment trajectory in a smaller solution space. The final optimization is clearly to be applied in practice. However, the optimization result that he obtained is one among numerous possible results. In this study, to achieve more precise and stable results of parameters, we narrowed the parameter limitations according to the volatility ranges and set the Fitness_L(Equation (8)) as the purpose to search the shortest trajectory to reduce the costs associated with the length of trajectory. The narrowed ranges are concluded in

Table 5. The narrowed optimization ranges of trajectory parameters.

The Narrowed Optimization Ranges of Trajectory Parameters
Length of the entry straight segment, L1∈[60, 110] m;
Length of the central straight segment, L3∈[100, 120] m;
Radius of curved segments, R2∈[2500, 3000] m, R4∈[1800, 2400] m;
Inclination angle of the entry straight segment, θ1∈[−10, −8] °;
Inclination angle of the exit straight segment, θ2∈[6, 8] °;

, and the parameters of the IRMO algorithm are: N = 50, M = 6, and generation = 250.

After searching 20 times by IRMO in the narrowed optimization ranges above, the total length of the trajectory is obviously shortened (Figure 13) compared with Figure 6 and Figure 9, although the fitting degree is lower.

The values of the shortest length searched 20 times are shown in Figure 14. The standard deviation is only 0.05813, in which the maximum decrease is 16.92 m, the minimum decrease is 10.62 m, and the average decrease is 14.53 m. The Figure 15 shows the convergence of the total length when the iteration of the IRMO is set as 250. The global optimum tends to be stable as well.

A set of trajectory parameters before and after optimization is compared in

Table 6. The trajectory design parameters before and after optimization.

	L1 [m]	θ1 [°]	R2 [m]	θ2 [°]	R4 [m]	L3 [m]	L5 [m]	Total Length [m]
Before	94.76	−8.87	2699.64	6.44	2063.20	100.00	289.06	1013.45
After	95.64	−8.00	3000	6.16	2330.27	104.96	131.18	1001.04

. The results after optimization of θ₁ and θ₂ are smaller than before, while the results after optimization of R₂, R₄, L₁, and L₃ are larger. Moreover, the length of the exit straight segment L₅ is exhibited, which demonstrates that the optimization effect of the total length is mainly contributed by the exit straight segment L₅.

Figure 16, Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 show the convergence of six parameters over 250 iterations in searching. The iterative convergence of the parameters carried by the central particle and the average of each generation of particles are given in the figures simultaneously. We note that the convergence of the central particle and average is relatively more stable (the fluctuation is smaller). This means that for the IRMO algorithm, if the searches for the global optimal value have unstable fluctuations (possibly falling into the local optimum), the central particle could maintain the optimal position to avoid the result becoming stuck in a local optimum (the optimum eventually coincides with the convergence of the central particle).

To highlight the optimization trends of parameters, we tagged the upper limit and lower limit (Table 5) in the following figures for each parameter. In Figure 16 and Figure 17, both θ₁ and θ₂ are trending to a lower value of their ranges, which means that the lower inclination of entry or exit straight segments is more conducive to optimizing the total length of the five-segment trajectory.

For the curved segments, R₂ and R₄, perform the different convergence trends in iteration. As shown in Figure 18, the final R₂ is neither close to the upper limit of its searching range nor the lower limit. Meanwhile, the values of R₄ converge to the upper limit to obtain the optimum total length (Figure 19). It could be considered that R₂ plays a more influential role in the length optimization.

For the length of entry straight segment L₁, an obvious convergence can be observed to its upper limit, especially after 100 generations. However, the length of the central straight segment L₃ finally trends to nearly 110 m.

It is hard to say what are the relations between the lower limit or upper limit with the optimization of R₂ and L₃ (Figure 18 and Figure 21), but it could be concluded that the determinations of R₂ and L₃ are not of as much importance as the others. To some extent, their influence on the optimal results is weaker. Therefore, it is advisable to pay more attention to optimizing other parameters, including L₁ R₄, θ₁, and θ₂, in actual projects. On the other hand, the optimization result of L₃ (104.96 m) is larger than the fitting result (100.00 m) (Table 6), which is conducive to the trajectory in practice if there needs to be a longer horizontal crossing.

This process of length optimization could certainly reduce the total length of the whole construction, especially when considering the cost associated with the length. As another possibility, it also proves that the trajectory can be designed and optimized with the purpose of minimal total length as in the ACO method [21]. However, for the IRMO method, the outputs of the design parameters are more specific with six parameters available which may help the engineers’ judgement. However, in this example, the maximum optimized length is only a fraction of the total length. Its practical value may need to have more practices to be proved in further study. However, it may provide valuable advice for determining certain parameters when those parameters are volatile during fitting.

4.2. Comparison of an Actual Drill Path Design

Another example containing mechanical design, the drill path of the Maxi-HDD Qin River crossing project in China [4], was tested to support this method of predesigning trajectory. The crossing project, including two parallel HDD crossing with 1016 mm diameter and 1750 m driven length, was located at an alluvial plain where the Yellow River and the Qin River intersect. Consisting of a 7° entry angle, 8° exit angle, 2290 m entry and exit radius of curved segments, 1015 m horizontal segment, 33.50 m maximum burial depth and 29.78 m minimum burial depth below the riverbed, the drill path was carefully designed manually for preventing blowout, as shown in Figure 22.

According to the engineering data provided in the literature [4], the geometrical design parameters are limited in the ranges given in

Table 7. Parameters and range setting.

Parameters of Catenary Trajectory	Parameters of Five-Segment Trajectory
Horizontal distance of exit and entry point: A = 1750 m Vertical distance of exit and entry point: H = −1 m Casing unit weight: q = 100 kg/m Pullback force: Np = 1,000,000 kgf	Length of the entry straight segment, L1∈[1, 100] m; Length of the central segment, L3∈[1000, 1500] m;
	Radius of curved segments,	case 1: R2&R4∈[2000, 3000] m; case 1: R2&R4∈[2300, 3000] m;
	Inclination angle of the entry straight segment, θ2∈[−18, 5] °; Inclination angle of the exist straight segment, θ4∈[5, 18] °;

. The minimum bending radius limitation was given differently for two cases, 2000 m and 2300 m, to illustrate its effects on trajectory design.

In terms of mechanical design, heavy loam is used to evaluate the maximum allowable mud pressure, including 1.57 g/cm³ dry density, 20.0 kPa cohesion, 15.0° friction angle, and 0.74 coefficient of lateral earth pressure. The MAP predicted at the location of minimum burial depth is 1.63 MPa, which is very close to the real pressure of about 1.7 MPa [4]. All segments of the trajectory were tested for the pressure predictions through the whole process of geometric design.

As shown in Figure 23 and

Table 8. Comparison of parameter design results.

	L1 [m]	θ1 [°]	R2 [m]	θ2 [°]	R4 [m]	L3 [m]	Minimum Burial Depth [m]	Total Length [m]
Actual design	—	−7	2290	8	2290	1015	29.78	—
Case 1 (IRMO)	100	−7.49	2000	6.04	2000	1000	31.90	1752.99
Case 2 (IRMO)	50.65	−8.00	2300	7.07	2300	1000	31.24	1752.26

, although the fitting similarity is largely reduced from the last example (mostly due to the difference of parameters ranges), both cases of trajectory design by IRMO show a great similarity with the actual design. The parameter design results arrive at similar trends as the previous conclusions. With the different minimum bending radius, the entry straight segment (L₁) is designed with nearly 50 m difference while the other parameters are not much different in design. It proves that the fitting similarity has less influence on the design results than we thought. The design parameters could be designed suitably by fitting with the catenary trajectory. What is more, there is no need for extra length optimization due to the tiny differences in the total length of the designs caused by strict limitations, especially for the central segment length L₃.

5. Conclusions

• With the application of intelligence algorithms, specifically the IRMO algorithm, the design and optimization of HDD trajectories could obviously improve the efficiency and accuracy in this study. All the calculations and optimizations were tested on a 2.30 GHz computer (CPU: Intel i5-6300 HQ), and the total time of fitting and length optimization was stabilized within 60 s, which shows great potential for shortening the manual adjustment time required in actual projects.
• This study proposed and tested a new method to design the HDD trajectory combining five-segment trajectory with catenary trajectory. With two processes, the five-segment trajectory could be designed with a great similarity to the catenary trajectory, and simultaneously with the shortest length, thereby reducing the associated costs. Both processes were verified through data from the literature or from actual projects. It was proved that the catenary trajectory has great potential to design a traditional drilling path with precise and achievable parameters.
• Six concise parameters (L₁, L₃, R₂, R₄, θ₁, and θ₂) were concluded herein to design trajectories and be optimized by the IRMO algorithm. According to the analysis of the parameters, it was found that the limits of L₃ and R₂ have fewer impacts on the optimization results. Therefore, the optimization ranges of the other four parameters should be set carefully during the trajectory design.
• By improving the structure data of the radial movement optimization, the obtained IRMO algorithm had a great ability to solve the extremum value of the multidimensional nonlinear objective function. Benefitting from the refined parameter matrix and efficient data structure, the authors believe IRMO has further potential in more complex trajectory design such as 3D construction or in avoiding obstacles.
• It should be emphasized that the importance of the accuracy and viability of the results by algorithm are mostly based on the set of the constraints and validated ranges of input parameters which largely rely on the experience of trained engineers. This method hopes to provide more advice in terms of geometrical parameters and drilling mud pressure for current practice preparations with a large consideration of site and soil property factors. There are more factors, such as geology, rock properties, and drilling rigs that should be discussed in further study to ensure the mechanical design. What is more, the mechanical advantages of the five-segments trajectory similar to the catenary trajectory also need more experimental studies to be verified.