FIG. 5 is a comparison diagram of the wellbore temperature after simulated annealing optimization and the actually-measured temperature in the present application; and

FIG. 6 is a comparison diagram of the wellbore temperature after simulated annealing optimization and the wellbore temperature calculated by the numerical model using the optimized parameters in the present application.

DETAILED DESCRIPTION OF EMBODIMENTS

The present application will be further described in detail below in conjunction with the embodiments and drawings.

In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are some, but not all of the embodiments of the present application. Based on the embodiments of the present application, all other embodiments obtained by a person ordinarily skilled in the art without creative work fall within the scope of protection of the present application. Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope claimed in the present application, but merely represents selected embodiments of the present application.

Example

A wellbore temperature optimization and prediction method integrating machine learning and numerical models includes the following steps:

Step S1: establishing a wellbore-formation transient heat transfer model, based on the principle of energy conservation combined with the heat transfer mechanism of each control area of a wellbore-formation, in consideration of the influence of heat generated by fluid circulation frictional resistance, and a complex heat source term on the wellbore temperature,

- where the established wellbore-formation transient heat transfer model is solved discretely by using the fully implicit finite difference method, which mainly includes the following contents:
- the difference between the incoming heat and outgoing heat of the heat transfer microelement plus the heat released by the internal heat source should be equal to the heat change of the microelement. The heat transfer control equation of each area is as shown in formula (1).

ρ
⁢

C
⁡
(

∂
T

∂
t

+

v
r

⁢

∂
T

∂
r

+

v
z

⁢

∂
T

∂
z

)

=

k
r

⁢

∂
T

∂
r

+

k
⁢

∂
2

T

∂

r
2

+

k
⁢

∂
2

T

∂

z
2

+
S

(
1
)

In the formula (1), T represents a temperature, ° C.; p represents a fluid density, kg/m³; C represents fluid specific heat capacity, J/(kg·° C.); k represents fluid thermal conductivity, W/(m·° C.); t represents time, s; S represents the complex heat source term; r and z represent a radial direction and an axial direction respectively; and v represents a flow velocity, m/s.

In the process of the fluid in the drill string flowing downward, it is assumed that there is an axial velocity. Ignoring the velocity gradient of the fluid in the radial direction and the heat conduction of the drill string in the axial direction, the in-drill-string heat transfer model may be expressed as formula (2):

ρ
m

⁢

C
m

(

∂

T
pi

∂
t

+

V
zp

⁢

∂

T
pi

∂
z

)

=

k
m

r

⁢

∂

T
pi

∂
r

+

S
pf

(
2
)

Convective heat transfer occurs between the fluid in the drill string and the wall of the drill string, and the boundary condition may be expressed as formula (3):

-

k
m

(

∂

T
pi

∂
r

)

r
=

r
pi

=

h
pi

(

T
p

-

T
pi

)

(
3
)

in the formulas (2) and (3), T_piand T_prepresent the temperatures of the fluid in the drill string and the wall of the drill string respectively, ° C.; μ_mrepresents the density of the drilling fluid, kg/m³; C_mrepresents the specific heat capacity of the drilling fluid, J/(kg·° C.); k_mrepresents the thermal conductivity of the drilling fluid, W/(m·° C.); h_pirepresents the convection heat transfer coefficient of the inner wall of the drill string, W/(m²·° C.); r_pirepresents the inner radius of the drill string, m; v_zprepresents the flow velocity of the fluid in the drill string, m/s; and S_pfrepresents the heat generated by the frictional resistance of the fluid circulation in the drill string, J.

The influence factor affecting the temperature of the wall of the drill string is the heat exchanged between the inner and outer wall surfaces thereof and the fluid through convective heat transfer. The heat transfer model of the wall of the drill string may be expressed as formula (4):

ρ
p

⁢

C
p

⁢

∂

T
p

∂
t

=

k
p

r

⁢

∂

T
p

∂
r

(
4
)

For the convective heat transfer generated between the fluid and the inner and outer walls of the drill string, the boundary conditions may be expressed by formulas (5) and (6):

-

k
p

(

∂

T
p

∂
r

)

r
=

r
pi

=

h
pi

(

T
pi

-

T
p

)

(
5
)

-

k
ρ

(

∂

T
p

∂
r

)

r
=

r
po

=

h
po

(

T
a

-

T
p

)

(
6
)

in the formulas (4) to (6), pp represents the density of drill string, kg/m³; C_prepresents the specific heat capacity of the drill string, J/(kg·° C.); k_prepresents the thermal conductivity of the drill string, W/(m·° C.); r_porepresents the inner radius of the drill string, m; and h_porepresents the convective heat transfer coefficient of the outer wall of the drill string, W/(m²·° C.).

During the upward flowing of the annular fluid, convective heat transfer occurs between the annular fluid and the outer wall surface of the drill string and between the annular fluid and the wall surface of the well wall. Ignoring the velocity gradient of the fluid in the radial direction and the heat conduction of the drill string in the axial direction, the annular heat transfer model may be expressed as formula (7):

ρ
m

⁢

C
m

(

∂

T
a

∂
t

+

V
za

⁢

∂

T
a

∂
z

)

=

k
m

r

⁢

∂

T
a

∂
r

+

S
a

(
7
)

When heat exchange occurs in the form of heat convection for the outer wall surface of the casing and the wall surface of the well wall, the boundary conditions of the two wall surfaces may be expressed by formulas (8) and (9):

-

k
m

(

∂

T
a

∂
r

)

r
=

r
po

=

h
po

(

T
W

-

T
a

)

(
8
)

-

k
f

(

∂

T
a

∂
r

)

r
=

r
w

=

h
w

(

T
f

-

T
a

)

(
9
)

In the formulas (7) to (9), T_a, T_wand T_frepresent the temperatures of the annular fluid, the wall surface of the well wall and the formation, respectively, ° C.; v_zarepresents the flow velocity of the annular fluid, m/s; h_wrepresents the convective heat transfer coefficient of the well wall, W/(m²·° C.); k_frepresents the thermal conductivity of the formation, W/(m·° C.); r_wrepresents the radius of the well wall, m; and S_arepresents the heat generated by the complex heat source term, J.

The complex heat source term S_amainly includes heat generated by the fluid circulation frictional resistance, and heat generated by drill string rotation, rock breaking by the drill bit, and nozzle pressure drop.

{circle around (1)} To calculate the heat generated by the annular fluid circulation frictional resistance, the τ_wequation for the flow rate Q of the given drilling fluid if firstly solved as shown in the formula (10):

Q
=

wh
2

2
⁢

K

1
2

⁢

τ
W
2

⁢

(

m

1
+

2
⁢
m

)

⁢

(

τ
w

-

τ
y

)

1
+
m

m

⁢

(

τ
w

+

m

1
+
m

⁢

τ
y

)

(
10
)

In the above:

h
=

D
2

-

D
1

2

;

w
=

π
2

⁢

(

D
2

+

D
1

)

.

Then, the generalized flow index Nis calculated according to the above formula. The calculation formula is as shown in formula (11):

3
⁢
N

1
+

2
⁢
N

=

3
⁢
m

1
+

2
⁢
m

[

1
-

(

1

1
+
m

)

⁢

x

-

(

m

1
+
m

)

⁢

x
2

]

(
11
)

In the above:

x
=

τ
y

τ
w

.

In the formulas (10)-(11), D₂represents the annulus outer diameter, and D₁represents the annulus inner diameter, mm; K represents the consistency coefficient; m represents the flow index; τ_wrepresents the shear stress of the wall surface, Pa; τ_yrepresents the yield stress, Pa; and w represents the flowing width of the fluid.

Then the annular fluid Reynolds number is calculated, and the fluid state (laminar flow, turbulent flow, transitional flow) is judged through the relationship between it with the critical Reynolds number. The calculation formula of the Reynolds number is as shown in formula (12):

R
⁢
e

=

ρ
⁢

V
a

(

D
2

-

D
1

)

u

(
12
)

In the formula (12), V_arepresents the flow velocity of annular fluid, m/s; and u represents the fluid viscosity, Pa·s.

The calculation formula of the critical Reynolds number is as shown in the formula (13) and the formula (14):

Re

1

=

2

,
TagBox[",", "NumberComma", Rule[SyntaxForm, "0"]]

100
·

[

N
0.331

(

1
+

1
.
4

⁢
02
⁢
κ

-

0
.
9

⁢
77
⁢

κ
2

)

-

0.019

ε
⁢

N
κ

-
0.868

]

(
13
)

Re
2

=

2

,
TagBox[",", "NumberComma", Rule[SyntaxForm, "0"]]

900
·

[

N

-
0.039

⁢

Re
1
0.307

]

(
14
)

In the formulas (13)-(14): N represents the generalized flow index, K represents the ratio of the annulus outer diameter to the annulus inner diameter, ε represents the eccentricity.

If Re<Re₁, the fluid state is the laminar flow; if Re<Re₂, the fluid state is the turbulent flow; and if Re₁<Re<Re₂, the fluid state is the transitional flow.

When the annular fluid is the laminar flow, the calculation formulas of frictional resistance coefficient when the annulus is non-eccentric and when it is eccentric are as shown in the formula (15) and the formula (16) respectively:

f
L

=

2
⁢
4

Re

(
15
)

f

L
2

=

f
L

·

[

1
-

0.072

ε
N

⁢

κ
0.8454

-

1.5

ε
2

⁢

N

⁢

κ
0.1852

+

0.96

ε
3

⁢

N

⁢

κ
0.2527

]

(
16
)

When the annular fluid is the turbulent flow, the calculation formulas of frictional resistance coefficient when the annulus is non-eccentric and when it is eccentric are as shown in formula (17) and formula (18) respectively:

f
T
0.5

=

4

N
0.75

⁢

log
⁡
(

Re
⁢

f
T

1
-

N
/
2

)

-

0
.
3

⁢
9

N
1.2

(
17
)

f

T
2

=

f
T

·

[

1
-

0.048

ε
N

⁢

κ
0.8454

-

2
3

⁢

ε
2

⁢

N

⁢

κ
0.1852

+

0
.
2

⁢
58
⁢

ε
3

⁢

N

⁢

κ
0.2527

]

(
18
)

When the annular fluid is the transitional flow, the calculation formulas of frictional resistance coefficient is as shown in formula (19):

f
tr

=

f
L

+

(

Re
-

Re
1

)

⁢

(

f
T

-

f
L

)

Re
2

-

Re
1

(
19
)

The pressure drop caused by circulation frictional resistance may be expressed as shown in formula (20):

dp
dz

=

2
⁢

ρ
m

⁢

fv
2

D
2

-

D
1

(
20
)

The heat generated by the circulation frictional resistance may be expressed as shown in formula (21):

S
f

=

Q
·

(

dp
dz

)

(
21
)

{circle around (2)} The heat generated by the friction of the drill bit breaking the rock at the bottom of the well may be expressed by equation (22):

S
d

=

ξ
⁢

fWR
⁡
(

D
2

+
Dd
+

d
2

)

7
4.3

(

D
+
d

)

(
22
)

In the formula (22): S_drepresents the heat generated by the drill bit breaking the rock, J; ξ=2.5939, represents the correction factor; f represents the frictional coefficient between the drill bit and the formation; w represents the drilling pressure, N; R represents the rotation speed, rpm; D represents the outer diameter of the drill bit, m; and d represents the inner diameter of the drill bit, m.

{circle around (3)} The nozzle pressure drop of the drill bit refers to the pressure drop value of the drilling fluid after the fluid passes through the nozzle of the drill bit. When the fluid flow rate and the size of the nozzle are constant, the expression of the drill bit pressure loss may be deduced based on the energy balance equation in hydrodynamics as shown in formula (23):

Δ
⁢

p
b

=

5
·

10
3

⁢

ρ
m

⁢

Q
2

C
2

⁢

A
0
2

(
23
)

In the formula (23): A₀²represents the equivalent diameter of the nozzle, m²; Q represents the flow rate of the drilling fluid, m³/s; and C represents the nozzle flow coefficient, which is usually taken as 0.95.

The heat generated by the nozzle pressure drop may be expressed through equation (24):

S
p

=

Δ
⁢

p
b

·
Q

(
24
)

{circle around (4)} The rotation of the drill string causes the fluid to generate spiral flow and work, generating heat. The Reynolds number Re and Taylor number Ta are often used to describe the influence of the rotation on the pressure drop caused by frictional resistance. The pressure drop caused by rotation frictional resistance may be expressed by the following formula (25), and the heat generated by the rotation of the drill string is expressed by formula (26):

(

dp
a

dz

)

ω

=

0
.
3

⁢
6
⁢

(

1
⁢

3
.
5

+

τ
y

ρ
⁢

v
2

)

0.428

⁢

ε
0.158

⁢

n
0.054

⁢

Ta
0.0319

⁢

Re
0.042

⁢

κ
⁡
(

1
κ

-
1

)

-
0.0152

(
25
)

S
rl

=

Q
·

(

dp
a

dz

)

ω

(
26
)

in the above:

Ta
=

r
1

(

r
2

-

r
1

)

3

⁢

(

ρ
m

⁢
w

μ

)

2

;

κ
=

r
1

r
2

;

In the formulas (25)-(26), τ_yrepresents the yield value of the Hershel-Bulkley fluid, Pa; ε represents the eccentricity; n represents the flow index; Re represents the Reynolds number; Ta represents the Taylor number; k represents the ratio of the wellbore radii; and u represents the fluid viscosity, Pa·s.

In summary, the complex heat source term S_aof the annular fluid may be expressed as formula (27):

S
a

=

S
p

+

S
d

+

S
f

+

S
rl

(
27
)

Due to the area around the well wall is composed of rocks, cement ring, casing and static drilling fluid, the heat flowing thereof is mainly in the form of heat conduction. Therefore, the temperature model of the formation near the well wall may be expressed as formula (28):

ρ
f

⁢

C
f

⁢

∂

T
i

∂
t

=

k
f

r
i

⁢

∂

T
i

∂
r

+

k
f

⁢

∂
2

T
i

∂

r
2

(
28
)

In the formula (28), Ti represents the temperature of each area near the well wall, ° C.; r_irepresents the radial distance of each area, m; i represents the number of radial units, usually greater than 6; ρ_frepresents the density of the rock near the well wall, kg/m³; and C_frepresents the specific heat capacity of the rock near the well wall, J/(kg·° C.).

For the established mathematical model, the initial conditions for solving the model are given as follows:

{circle around (1)} Before drilling, the wellbore-formation has reached a thermodynamic equilibrium state. Therefore, the temperature of each control area (the fluid in the drill string, drill string wall, annular fluid and formation near the well wall, etc.) of the wellbore-formation is the temperature of the original formation. The mathematical model thereof may be expressed as:

T

k
,
j

=

T
s

+

g
f

⁢
z
⁢
cos
⁢
θ

(
29
)

In the formula (29): T_kjrepresents the temperature of the fluid in the drill string, drill string wall, annular fluid and formation near the well wall, ° C.; θ represents the well inclination angle, °; k represents the number of grid cells in the radial direction of the wellbore-formation, 1≤k≤11; j represents the number of grid cells in the axial direction; T_srepresents the land surface temperature, ° C.; gf represents the geothermal gradient, ° C./100 m; and Z represents the well depth, m.

{circle around (2)} On the ground, the fluid temperature in the drill string is the inlet temperature of the fluid:

T
p

(

z
=
0

,
 
t

)

=

T
in

(
30
)

In the formula (30): T_inrepresents the inlet temperature of the drilling fluid, ° C.

The boundary conditions for solving the model include:

{circle around (1)} The temperatures of the fluid in the drill string, the drill string wall and the annular fluid are all equal at the bottom depth of the well, that is:

T
ρ

(

z
=
H

,
 
t

)

=

T
pi

(

z
=
H

,
 
t

)

=

T
a

(

z
=
H

,

t

)

(
31
)

{circle around (2)} Far away from the well wall, the formation temperature is not disturbed, that is, the temperature far away from the well wall is the original temperature:

∂

T
fk

(

r
,
z
,
t

)

∂
r

❘
"[LeftBracketingBar]"

r
→
∞

=

0
⁢

or
⁢

T
fk

(

r
→
∞

,
z
,
t

)

=

T
s

+

g
f

⁢
z
⁢
cos
⁢
θ

(
32
)

{circle around (3)} On the land surface, there is no heat exchange between each control area of the wellbore-formation and the atmosphere:

∂

T
fk

(

r
,
z
,
t

)

∂
z

❘
"[LeftBracketingBar]"

z
→
0

=

0
⁢

or
⁢

T
fk

(

r
,

z
=
0

)

=

T
s

(
33
)

In order to accurately solve the transient temperature change of the wellbore, the finite difference method is used to discretize each heat transfer model of the wellbore-formation in time and space; then, the temperature data of the same time and space nodes at these discrete points are integrated to construct a linear equation system. According to equations (7)-(9), the annular heat transfer model may be discretely expressed as:

ρ
m

⁢

qC
m

⁢

∂

T
a

∂
z

+

2
⁢
π
⁢

r
po

⁢

h
po

(

T
p

-

T
a

)

-

2
⁢
π
⁢

r
w

⁢

h
w

(

T
a

-

T
w

)

+

Q
Sa

=

πρ
m

⁢

C
m

(

r
w
2

-

r
po
2

)

⁢

∂

T
a

∂
t

(
34
)

Then, the formula (34) is discretized in space and time, and may be expressed as:

ρ
m

⁢
q
⁢

C
m

⁢

T

a
,

j
+
1

n
+
1

-

T

a
,
j

n
+
1

Δ
⁢

z
j

+

2
⁢
π
⁢

r

p
⁢
o

⁢

h

p
⁢
o

(

T

p
,
j

n
+
1

-

T

a
,
j

n
+
1

)

-

2
⁢
π
⁢

r
w

⁢

h
w

(

T

a
,
j

n
+
1

-

T

w
,
j

n
+
1

)

+

Q

S
⁢
a

=

π
⁢

ρ
m

⁢

C
m

(

r
w
2

-

r
po
2

)

⁢

T

a
,
j

n
+
1

-

T

a
,
j

n

Δ
⁢
t

(
35
)

By combining the same terms in the formula (35), the model may be expressed as the following linear equation:

(
36
)

2
⁢
π
⁢

r

p
⁢
o

⁢

h

p
⁢
o

·

T

ρ
,
j

n
+
1

-

[

ρ
m

⁢
q
⁢

C
m

Δ
⁢

z
j

+

2
⁢
π
⁢

r

p
⁢
o

⁢

h

p
⁢
o

+

2
⁢
π
⁢

r
w

⁢

h
w

+

ρ
m

⁢

C
m

⁢

π
⁡
(

r
w
2

-

r
po
2

)

Δ
⁢
t

]

·

T

a
,
j

n
+
1

+

ρ
m

⁢
q
⁢

C
m

Δ
⁢

z
j

·

T

a
,

j
+
1

n
+
1

+

2
⁢
π
⁢

r
w

⁢

h
w

·

T

w
,
j

n
+
1

=

Q

S
⁢
a

-

ρ
m

⁢

C
m

⁢

π
⁡
(

r
w
2

-

r
po
2

)

Δ
⁢
t

·

T

a
,
j

n

In the formulas (34) to (36), q represents the flow rate, m³/s; j represents the number of well depth nodes; n represents the number of time nodes; Tⁿ⁺¹_pjrepresents the temperature of the drill string wall at well depth j and time n+1, ° C.; Tⁿ⁺¹_ajrepresents the annular temperature at the well depth j and time n+1, ° C.; Tⁿ⁺¹_wjrepresents the temperature of the well wall at the well depth j and time n+1, ° C.; Tⁿ_ajrepresents the annular temperature at the well depth j and time n, ° C.; Δz_jrepresents the depth step; Δt represents the time step; and Q_sarepresents the heat generated by the complex heat source term, J.

The above is the wellbore-formation transient heat transfer numerical module established based on the principle of energy conservation. By comprehensively considering the heat transfer mechanisms of different regions and the influence of the complex heat source term (including heat generated by fluid circulation frictional resistance, and heat generated by drill string rotation, rock breaking by the drill bit, and nozzle pressure drop) on the temperature, the corresponding transient heat transfer control equations are established for the inside of the drill string, the wall of the drill string, annulus and area near the well wall from inside to outside, the corresponding initial conditions and boundary conditions are given, and the finite difference method is used to discretely solve each equation. The established numerical models have high accuracy in predicting and calculating the wellbore temperature.

Step S2: obtaining the initial data set composed of relevant parameters by combining the actually-measured data of the on-site well group;

By sorting out the on-site actually-measured data of multiple wells, the well depth, drilling time, circulation time, inlet temperature, mechanical drilling speed, drilling pressure, rotation speed, flow rate and wellbore temperature are selected to form a data set for machine learning, where the bottom hole temperature is the prediction label and the others are feature labels. Some data sets are as shown in Table 1:

TABLE 1

Some datasets used for machine learning

Well
Drilling

Inlet
Mechanical
Drilling
Rotation
Flow
Bottom hole

depth
time
Circulation
temperature
drilling speed
pressure
speed
rate
temperature

(m)
(h)
time (h)
(° C.)
(m/h)
(KN)
(r/min)
(L/s)
(° C.)

5782.00
1.82
2.90
22
7.69
13
80
30.5
133.2

5806.00
2.28
3.00
22
10.53
13
80
30.5
133.2

5820.00
1.29
2.19
22
10.85
13
80
30.5
131.7

5836.00
2.08
3.40
22
7.69
13.5
80
30.5
131.7

5855.00
2.02
2.86
22
9.41
13.5
80
30.5
127.3

5873.00
1.98
2.80
20
9.09
13.5
80
30.5
133.2

5892.00
1.96
2.83
20
9.69
13.5
80
30.5
134.0

5921.00
2.39
3.00
21
12.13
13.5
80
30.5
134.0

5938.00
1.60
2.72
22
10.63
13.5
80
30.5
134.0

5958.00
1.72
2.76
21
11.63
13.5
80
30.5
134.0

5980.00
1.70
2.91
22
12.94
13.5
80
30.5
136.2

6008.00
1.97
2.84
22
14.21
13.5
60
30.5
134.4

6020.00
2.39
2.83
19
5.02
13.5
60
30.5
134.7

6039.00
2.04
2.97
21
9.31
13.5
60
30.5
134.7

6060.00
2.32
2.86
22
9.05
13.5
60
29.5
136.2

6081.00
2.32
3.01
22
9.05
14.5
70
29.5
133.9

6098.00
1.69
2.79
21
10.06
14.5
65
29.5
131.7

6123.00
2.41
3.17
22
10.37
14.5
65
29.5
131.7

6143.00
1.92
2.61
22
10.42
14.5
65
29.5
133.9

6164.00
2.30
2.82
24
9.13
14.5
65
29.5
131.7

6185.00
2.24
2.85
25
9.38
15
65
29.5
131.0

6209.00
2.50
2.97
25
9.60
15
65
29.5
131.0

6236.00
2.66
2.86
25
10.15
15
65
29.5
131.0

6257.00
2.22
2.89
25
9.46
15
65
29.5
131.7

6276.00
1.96
2.75
25
9.69
15
65
29.5
131.7

6291.00
2.38
2.98
25
6.30
15
65
29.5
131.7

6305.00
2.87
2.90
24
4.88
15
65
29.5
133.2

6307.00
1.98
3.08
24
1.01
15
65
29.5
133.2

Step S3: normalizing the initial data set;

Since there are dimensional inconsistency among the parameters in the data set, data normalization is applied to proportionally transform the original data into the [0,1] interval to ensure that the data have uniform magnitude.

Some normalized data is as shown in Table 2:

TABLE 2

Part of the data set after normalization

Well
Drilling
Circulation
Inlet
Mechanical
Drilling
Rotation

depth
time
time
temperature
drilling
pressure
speed
Flow rate

(m)
(h)
(h)
(° C.)
speed (m/h)
(KN)
(r/min)
(L/s)

0.716
0.547
0.099
0.125
0.229
0.500
0.800
0.444

0.726
0.746
0.113
0.125
0.326
0.500
0.800
0.444

0.732
0.319
0.000
0.125
0.337
0.500
0.800
0.444

0.739
0.659
0.168
0.125
0.229
0.625
0.800
0.444

0.747
0.634
0.093
0.125
0.287
0.625
0.800
0.444

0.754
0.616
0.085
0.063
0.277
0.625
0.800
0.444

0.762
0.608
0.089
0.063
0.297
0.625
0.800
0.444

0.775
0.793
0.113
0.094
0.381
0.625
0.800
0.444

0.782
0.453
0.074
0.125
0.329
0.625
0.800
0.444

0.790
0.504
0.079
0.094
0.364
0.625
0.800
0.444

0.800
0.496
0.100
0.125
0.408
0.625
0.800
0.444

0.812
0.612
0.090
0.125
0.452
0.625
0.000
0.444

0.817
0.793
0.089
0.031
0.137
0.625
0.000
0.444

0.825
0.642
0.108
0.094
0.284
0.625
0.000
0.444

0.834
0.763
0.093
0.125
0.275
0.625
0.000
0.222

0.842
0.763
0.114
0.125
0.275
0.875
0.400
0.222

0.850
0.491
0.083
0.094
0.310
0.875
0.200
0.222

0.860
0.802
0.136
0.125
0.321
0.875
0.200
0.222

0.869
0.591
0.058
0.125
0.322
0.875
0.200
0.222

0.878
0.754
0.088
0.188
0.278
0.875
0.200
0.222

0.886
0.728
0.092
0.219
0.286
1.000
0.200
0.222

0.897
0.841
0.108
0.219
0.294
1.000
0.200
0.222

0.908
0.909
0.093
0.219
0.313
1.000
0.200
0.222

0.917
0.720
0.097
0.219
0.289
1.000
0.200
0.222

0.925
0.608
0.078
0.219
0.297
1.000
0.200
0.222

0.931
0.789
0.110
0.219
0.181
1.000
0.200
0.222

0.937
1.000
0.099
0.188
0.132
1.000
0.200
0.222

0.938
0.616
0.124
0.188
0.000
1.000
0.200
0.222

Step S4: training the wellbore temperature prediction model by using the random forest algorithm, and then optimizing hyperparameters of the random forest algorithm by using the genetic algorithm.

The random forest algorithm is an integrated learning method that improves the performance and stability of the model mainly by constructing multiple decision trees. The principle of random forest algorithm is as follows.

First, the Bagging method is used to perform random sampling, that is, the training data set is randomly sampled with replacement:

Bagging: a data set containing N samples is random sampled with replacement to generate K subsample sets, each of which also contains N samples. For the k-th subsample set, the definition is as shown in formula (37):

(
37
)

D
k

=

{

(

x
i

,

y
i

)

|

i

∈

random
⁢

sampling
⁢

with
⁢

replacement
⁢

from
⁢

{

1
,
2
,
…

,
N

}

}

Secondly, a decision tree is constructed for each subsample, and at the node of each decision tree, a part of features are randomly selected to determine the best splitting, thereby increasing the diversity of the model and reducing the risk of overfitting:

- node splitting: at the node of each decision tree, m features (m<M, where M is the total number of features) are randomly selected, the splitting point of each feature is calculated, and the feature and its splitting point that can maximally reduce the Gini index for splitting are selected. The Gini index is defined as shown in formula (38):

Gini
⁡
(
D
)

=

1
-

∑

k
=
1

K

p
k
2

(
38
)

In the formula (38): p_krepresents the proportion of class-k samples in the data set, and K represents the total number of classes.

Finally, the prediction results of all decision trees are integrated and the average value is taken as the prediction result of the random forest algorithm, as shown in formula (39):

y
^

=

1
B

⁢

∑

b
=
1

B

y
^

b

(
39
)

In the formula (39): ŷ represents the final prediction result, ŷ_brepresents the prediction result of the b-th decision tree, and B represents the total number of decision trees.

For the preprocessed data set, the training samples and test samples are divided in a ratio of 8:2.

Meanwhile, in order to further improve the accuracy of the random forest prediction model, the random forest hyperparameters n_estimators (number of decision trees), max_depth (maximum depth of the tree), min_samples_leaf (minimum number of samples) and min_samples_split (minimum number of split samples) are optimized by using the genetic algorithm (GA) combined with the five-fold cross-validation method. The genetic algorithm population size is set to 40, the mutation probability is 0.05, and the maximum number of iterations is 50. Different hyperparameter combinations are iteratively generated and evaluated to find the optimal parameter combination. After optimization, the optimal hyperparameter combination is as shown in Table 3 below:

TABLE 3

GA optimization result

Hyperparameter
Set Value

n_estimators
150

max_depth
5

min_samples_leaf
2

min_samples_split
4

After optimization by the genetic algorithm, the random forest prediction model is evaluated using the determination coefficient (R²). The calculation formula of the determination coefficient is as shown in formula (40):

R
2

=

1
-

∑

i
=
1

n

(

y
i

-

x
i

)

2

∑

i
=
1

n

(

y
i

-
y

)

2

(
40
)

In the formula (40): n represents the number of test sets, x_irepresents the prediction value, y_irepresents the true value, and y represents the average value of the true values.

The determination coefficient R²is used to measure the correlation between the true value and the prediction value, and the value thereof is between 0 and 1. The closer the value is to 1, the better the prediction performance of the model is. The determination coefficient of the optimized wellbore temperature prediction model is about 0.978, which proves that this model has a very high prediction accuracy. The comparison between the model prediction value and the true value of the wellbore temperature is as shown in FIG. 2 and FIG. 3. Based on this prediction model, the feature importance analysis of each feature parameter on the bottom hole temperature is performed, which can reflect the contribution of each feature to the prediction performance of the entire model. The result of the feature importance analysis is as shown in FIG. 4. The feature importance of the feature X_jto the model is calculated as shown in formula (41):

FI
⁡
(

X
j

)

=

1
B

⁢

∑

t
=
1

B

Δ
⁢

l
⁡
(

t
,

X
j

)

(
41
)

In the formula (41): B represents the total number of decision trees, and Δl(t, X_j) represents the impurity reduction caused by the feature X_jat the t-th tree.

Step S5: globally optimizing the algorithm by using an annealing algorithm to obtain an optimized wellbore temperature.

After determining the feature parameters to be optimized, the simulated annealing algorithm is used for global optimization, to find the optimal parameters to minimize the wellbore temperature. The simulated annealing algorithm, which is a general probabilistic algorithm based on the principle of solid annealing, is mainly used to find the optimal solution of a proposition in a large search space.

Its core formula is mainly based on the Metropolis criterion: that is, certain inferior solutions are allowed to be accepted, instead of being blindly rejected. That is, even if a local optimal solution is found, the iteration continues until the global optimal solution is found. Assuming that the energy state of the solid here is E_old, by slightly disturbing the internal structure of the solid to obtain a new state, whose energy is E_new, the probability of accepting the state after the disturbance is as shown in formula (42):

P
=

{

1

E
⁢

(

x

n
⁢
e
⁢
w

)

<

E
⁢

(

x

o
⁢
l
⁢
d

)

exp
⁢

(

-

E
⁡
(

x

n
⁢
e
⁢
w

)

-

E
⁡
(

x

o
⁢
l
⁢
d

)

T

)

E
⁢

(

x

n
⁢
e
⁢
w

)

≥

E
⁢

(

x

o
⁢
l
⁢
d

)

(
42
)

In the formula (42): T represents the current temperature.

The attenuation function of the temperature is controlled to be as shown in formula (43):

T

k
+
1

=

α
⁢

T
k

,

k
=
0

,
1
,

2
⁢

…

(
43
)

In the formula (43), a represents a constant close to 1, usually taken as 0.5˜0.99.

The basic process of simulated annealing is divided into the following four steps: {circle around (1)} giving an initial temperature T₀(i.e., the current optimal solution), and randomly generating an initial solution x₀(i.e., x_old) at the same time, and calculating E(x₀) of the current objective function; {circle around (2)} executing the third step if the current temperature reaches that meeting the internal circulation stop condition, otherwise, giving a disturbance to the current solution to generate a new solution x₁(x_new), and calculating ΔE=E(X₁)−E(x₀), where if ΔE<0, the new solution is accepted, and if ΔE≥0, the new solution is accepted with the probability

exp
⁡
(

-

E
⁡
(

x

n
⁢
e
⁢
w

)

-

E
⁡
(

x
old

)

T

)

by using the Metropolis criterion, and the new solution is used as the current solution when the new solution is accepted, otherwise, the second step is repeated; {circle around (3)} k=k+1, executing the temperature attenuation function, and executing the fourth step if the termination condition is met, otherwise, executing the second step; and {circle around (4)} ending the operation, where T reaches the termination temperature, i.e. the global optimal solution is found.

First, the original parameter values of the data set are used as the initial values of the optimization parameters, and the values are stored in a NumPy array initial_parameters. Next, the search space (i.e., parameter range) of five optimization parameters is given. Here, the range of each parameter in the data set is used for assignment and stored in the parameter_bounds list. Then the dual_annealing method in the simulated annealing algorithm is called, and the above initial values and search space are assigned to this function to minimize the objective function (wellbore temperature), and finally the optimal parameter value is returned. Finally, the optimal parameters are applied to a copy of the original data row (the original parameters are replaced with the optimal parameters), the new data row is converted to DataFrame format, and the loaded normalization class is used for normalization, and finally, the loaded prediction model is used to predict the optimized new data row, to finally obtain the optimized wellbore temperature. The optimized parameter combination and the optimized temperatures are saved as Excel files for further data analysis. The comparison between the actually-measured temperatures and the optimized wellbore temperatures are as shown in FIG. 5. It can be seen from the figure that in the well depth range of 4,000 m to 5,800 m, the difference between the optimized temperature and the initial temperature is small, and the optimization range is small. After the well depth reaches 5,800 m, the optimized temperature is 3-4° C. lower than the initial temperature, indicating that the optimization algorithm has played a certain role in temperature control. The recommended parameters of some well sections after optimization are as shown in Table 4:

TABLE 4

Optimized parameters and results of some well sections