# Spatial-Temporal ARX Modeling and Optimization for Polymer Flooding

Authors: Yulei Ge; Shurong Li; Songlin Lu; Peng Chang; Yang Lei
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A new polymer flooding model based on spatial-temporal decomposition and autoregressive model with external input (ARX) (STDARX model) is proposed. Karhunen-Loeve (K-L) decomposition is used to model the two-dimensional state parameters of reservoir (such as water saturation, pressure, and grid concentration). The polymer injection concentration and time coefficient got from the decomposition are taken as the input and output information. After being identified by least square method, the time iterative ARX models of all state variables are obtained, we build the ARX model among pressure, water saturation, grid concentration, and moisture content of production well, and identify it with recursive least-squares (RLS) method. After combining the above two models, we get the STDARX model of polymer flooding. The accuracy is proved by model with four injection wells and nine production wells through data which is obtained from mechanism model. In order to enhance the polymer flooding oil recovery when oil price is changing, iterative dynamic programming (IDP) is applied to optimize the STDARX model, to get the optimal injection of production scheme.


---

## Body

## 1. Introduction

It is of increasing necessity to produce oil fields more efficiently and economically because of the ever-increasing demand for petroleum worldwide. Chemical flooding is an effective measure to enhance oil recovery, and it has been widely used in practice. Polymer flooding is one of the most important parts of chemical flooding [1]. The basic principle is improving the fluid properties and reducing the moisture content by adding macromolecule polymer into water. Thus the oil recovery is enhanced. But the polymer flooding is difficult to exploit and needs much investment; moreover the polymer is often too expensive. It is very important to make a reasonable injection plan to obtain the biggest economic benefits.In recent years, some scholars have studied the optimal control approaches of polymer flooding. In 2004, Brouwer [2] gave a detailed description of dynamic water flooding, such as the mechanism modeling, smart wells, and the optimal control theory. In 2005, Zhang and Li [3] proposed the control vector parameterization with time parameters to solve the optimal control of polymer flooding, and the time nodes are improved. In 2008, Li et al. [4] adopted the SWIFT method to combine simplex method with penalty function method, used nonlinear programming method to solve the optimization problem, and optimized the injection concentration and slug size at the same time. In 2012, Lei et al. [5] presented a mix-integer iterative dynamic programming algorithm which incorporated a time normalization strategy and a special truncation procedure and solved the multiobjective problem of polymer flooding.But all of the above researches are studied on the basis of polymer flooding mechanism model. The model consists of a series of partial differential equations, which are the fluid seepage equation and the polymer absorption equation [6]. It has high dimension and heavy calculated amount and is very difficult to solve. Some scholars have attempted to model polymer flooding.Liu et al. [7] considered threshold pressure gradient, wellbore storage, and skin effect for polymer flooding and got the new model. Cheng et al. [8] presented a new model of polymer flooding by considering the viscoelasticity of polymer and the effect of shear rate on resistance coefficient. Rai et al. [9] built a model of surfactant and surfactant-polymer flooding for enhanced oil recovery using STARS (CMG) software, through a series of flooding experiments.Though many people have done researches on modeling of polymer flooding, nearly all the works are improving or enriching the primary model. What they have done is only to consider more influencing factors or inductive research; the model is still complex and difficult to apply. Since the polymer flooding is distributed parameter system, it is hard to model in common method. To find out a simple but accurate relatively model, we must introduce new method to build the distributed parameter model.Spatial-temporal decomposition technology is a very important method when distributing parameter modeling. It can reflect the information of time region and space region and gets the reasonable model. K-L (Karhunen-Loeve) decomposition [10] is one spatial-temporal decomposition method. In this method, the original data set is converted into the principal component space, and the cross correlation of single sample is reduced to the minimum. So under the condition of guaranteeing the accuracy, it can cut the system dimension as much as possible [11]. The initial system can be divided into spatial basis functions and temporal coefficients through K-L decomposition, which reflect the spatial distribution and the time variation, respectively. In addition, autoregressive model with external input (ARX) is a wide-used method to build the time model. It can model the time coefficients considering external input. Moreover, the identification accuracy can be determined artificially.Dynamic programming (DP) is an optimal method with reverse acquisition [12]. The basic idea is Bellman optimality principle. In this method, the decision-making process is divided into many stages. Computation is carried out from the last stage. Every optimal decision is got from the current stage. Then it is finished until the first stage. To overcome the shortcomings of low computation efficiency and curse of dimensionality for traditional DP, iterative dynamic programming (IDP) is proposed in 1990 by Luus [13, 14]. In terms of iteration, IDP avoids solving Hamilton-Jacobi-Bellman (HJB) equation and decreases the computing dimension. It does not need to solve complex nonlinear programming, has good robustness, and is easy to be carried out. So IDP is often applied in optimization of polymer flooding.In this paper, an ARX model of polymer flooding is obtained in terms of spatial-temporal decomposition. K-L decomposition is used to model the state parameters of reservoir (such as water saturation, pressure, and grid concentration) separately, and a series of time coefficients and spatial basis function can be obtained. The polymer injection concentration and spatial coefficient are taken as the input and output information. After being identified by RLS, the time iterative ARX models of all state variables are obtained. Then, we build the multivariable ARX model among pressure, water saturation, grid concentration, and moisture content of production well and identify it with RLS. After combing the above two models, we will get the whole polymer flooding ARX model with spatial-temporal decomposition (STDARX model). After being tested by mechanism model’s data, in order to enhance the polymer flooding oil recovery when oil price is changing, iterative dynamic programming is applied to get the optimal injection to production scheme for STDARX model.


## 2. Mechanism Model Description of Polymer Flooding

LetΩ
∈
R
2 denote the plain domain of reservoir with boundary ∂
Ω, n the unit outward normal on ∂
Ω, and x
,
y
∈
Ω the coordinate of a point in the reservoir. Given a terminal time t
f, there are N
w injection wells and N
o production wells. The injection and production wells are located at L
w
=
{
x
w
i
,
y
w
i
∣
i
=
1,2
,
…
,
N
w
} and L
o
=
{
x
o
j
,
y
o
j
∣
j
=
1,2
,
…
,
N
o
}, respectively. Then the mathematical description of polymer flooding [15] is as follows.The flow equation for oil phase(1)
∂
∂
x
r
o
h
∂
p
∂
x
+
∂
∂
y
r
o
h
∂
p
∂
y
-
1
-
f
w
q
out
=
h
∂
∂
t
ϕ
1
-
S
w
B
o
.The flow equation for water phase(2)
∂
∂
x
r
w
h
∂
p
∂
x
+
∂
∂
y
r
w
h
∂
p
∂
y
+
q
in
-
f
w
q
out
=
h
∂
∂
t
ϕ
S
w
B
w
.The flow equation for polymer component(3)
∂
∂
x
r
d
h
∂
c
p
∂
x
+
∂
∂
y
r
d
h
∂
c
p
∂
y
+
∂
∂
x
r
c
h
∂
p
∂
x
-
f
w
q
out
c
p
+
∂
∂
y
r
c
h
∂
p
∂
y
+
q
in
c
p
in
=
h
∂
∂
t
ϕ
p
S
w
c
p
B
w
+
ρ
r
1
-
ϕ
C
r
p
.The initial conditions and boundary conditions(4)
p
x
,
y
,
t
|
t
=
0
=
p
0
x
,
y
,
S
w
x
,
y
,
t
|
t
=
0
=
S
w
0
x
,
y
,
c
p
x
,
y
,
t
|
t
=
0
=
c
p
0
x
,
y
,
x
,
y
∈
Ω
,
∂
p
∂
n
∂
Ω
=
0
,
∂
S
w
∂
n
∂
Ω
=
0
,
∂
c
p
∂
n
∂
Ω
=
0
,
where r
o
=
K
k
r
o
/
B
o
μ
o, r
w
=
K
k
r
w
/
B
w
R
k
μ
w, r
d
=
D
ϕ
p
S
w
/
B
w, r
c
=
K
k
r
w
c
p
/
B
w
R
k
μ
p, K
x
,
y is the absolute permeability (μm2), D is the diffusion coefficient of polymer (m2/s), h is the thickness of reservoir bed (m), ρ
r is the rock density (kg/m3), and μ
o is the oil viscosity (mPa·s).The oil and water volume factorsB
o, B
w, the rock porosity ϕ, and the effective porosity to polymer ϕ
p are related to reservoir pressure p. Consider
(5)
B
o
=
B
o
r
1
+
C
o
p
-
p
r
,
B
w
=
B
w
r
1
+
C
w
p
-
p
r
,
ϕ
=
ϕ
r
1
+
C
r
p
-
p
r
,
ϕ
p
=
f
a
ϕ
,
where p
r is the reference pressure (MPa) and ϕ
r, B
o
r, and B
w
r denote the porosity, the oil, and water volume factors under the condition of p
r, respectively. f
a is the effective pore volume coefficient and C
o, C
w, and C
r denote the compressibility factors of oil, water, and rock, respectively.The oil and water relative permeabilityk
r
o and k
r
w are the functions of water saturation S
w and can be described as
(6)
k
r
w
=
k
r
w
r
o
S
w
-
S
w
c
1
-
S
w
c
-
S
o
r
n
w
,
k
r
o
=
k
r
o
c
w
1
-
S
w
-
S
o
r
1
-
S
w
c
-
S
o
r
n
o
,
where k
r
w
r
o denotes oil relative permeability at the irreducible water saturation S
w
c and k
r
o
c
w denotes the water relative permeability at the residual oil saturation S
o
r.The polymer solution viscosityμ
p (mPa·s), the permeability reduction factor R
k, and the amount absorbed per unit mass of the rock C
r
p (mg/g) are the functions of polymer concentration c
p. Consider
(7)
μ
p
=
μ
w
1
+
a
p
1
·
c
p
+
a
p
2
·
c
p
2
+
a
p
3
·
c
p
3
,
R
k
=
1
+
b
r
k
c
p
R
k
max
⁡
-
1
1
+
b
r
k
c
p
,
C
r
p
=
a
·
c
p
1
+
b
·
c
p
,
where μ
w denotes the viscosity of water with no polymer (mPa·s) and a
p
1
,
a
p
2
,
a
p
3
,
R
k
max
⁡
,
b
r
k
,
a
,
b are the constant parameters.The moisture contentf
w of production wells depends on p, S
w, and c
p, and it can be defined as
(8)
f
w
=
r
w
r
o
+
r
w
.When applied in computing the polymer flooding, the model can obtain thef
w of production wells and get the oil production.


## 3. Polymer Flooding Model with Spatial-Temporal ARX Method

### 3.1. The Two-Dimensional Spatial-Temporal Decomposition of State Parameters

In this paper, the K-L decomposition [16] method is adopted to build the two-dimensional distributed parameter spatial-temporal model of polymer flooding STDARX model [17]. As to this model, the input variables are pressure, water saturation, and grid concentration. The output variable is moisture content.Suppose the output of system isY
x
i
,
y
j
,
t
i
=
1
,
j
=
1
,
t
=
1
N
x
,
N
y
,
L, t
=
1,2
,
…
,
L denotes L individuals discrete points of output in time domain, i
=
1,2
,
…
,
N
x, j
=
1,2
,
…
,
N
y is N
x
×
N
y individuals state sample points in space.The output can be decomposed as follows:(9)
Y
x
,
y
,
t
=
∑
i
=
1
∞
φ
i
x
,
y
T
i
t
,
where φ
i
x
,
y
i
=
1
∞ denotes the spatial basis function, T
i
t
i
=
1
∞ denotes the time coefficient, and Y
x
,
y
,
t is the system output.The output is usually described as the finite form in specific engineering application:(10)
Y
n
x
,
y
,
t
=
∑
i
=
1
n
φ
i
x
,
y
T
i
t
,
where n mirrors the system energy; it depends the number of spatial basis function.The spatial basis function is unit orthogonal, which is defined as below:(11)
φ
i
x
,
y
,
φ
j
x
,
y
=
∬
Ω
φ
i
x
,
y
φ
j
x
,
y
d
x
d
y
=
0
,
i
≠
j
1
,
i
=
j
,
where f
x
,
g
x
=
∬
Ω
f
x
g
x
d
x
d
y is the inner product.If we know the spatial basis function, we can get the time coefficient:(12)
T
i
t
=
φ
i
x
,
y
,
Y
x
,
y
,
t
,
i
=
1,2
,
…
,
n
.Convert the problem of solving spatial basis function to the below optimization problem:(13)
min
⁡
φ
i
x
,
y
Y
x
,
y
,
t
-
Y
n
x
,
y
,
t
2
s
.
t
.
φ
i
,
φ
i
=
1
,
φ
i
∈
L
2
Ω
,
i
=
1,2
,
…
,
n
,
where f
x
,
y
=
f
x
,
y
,
f
x
,
y
1
/
2 is the norm and f
x
,
y
,
t
=
(
1
/
L
)
∑
t
=
1
L
f
x
,
y
,
t is the set average.The necessary condition of extremum existing to (13) is
(14)
∬
Ω
R
x
,
y
,
ξ
,
ζ
φ
i
ξ
,
ζ
d
ξ
d
ζ
=
λ
i
φ
i
x
,
y
φ
i
,
φ
i
=
1
,
i
=
1,2
,
…
,
n
,
where R
x
,
y
,
ξ
,
ζ
=
Y
x
,
y
,
t
Y
ξ
,
ζ
,
t is the correlation function between two points in the same plane. Thus the optimization problem is transformed to solve the problem described as (14).In normal conditions, (14) is very difficult to solve. But the snapshots [18] method can be applied when time nodes are less than space nodes; the calculated amount is cut down at the same time. The meaning of snapshots is as follows. When we choose a time node in time domain, all the output sampling values in two-dimensional plane at this time make up one snapshot. For example, Y
x
,
y
,
3, x
=
1,2
,
…
,
N
x, y
=
1,2
,
…
,
N
y, and N
=
N
x
·
N
y denote one snapshot of time node t
=
3.On the condition that time nodes are less than space nodes, let spatial basis function be expressed by a series of snapshots linear combination. Consider(15)
φ
i
x
,
y
=
∑
t
=
1
L
γ
i
t
Y
x
,
y
,
t
.In terms of (14), we can get
(16)
∬
Ω
R
x
,
y
,
ξ
,
ζ
φ
i
ξ
,
ζ
d
ξ
d
ζ
=
∬
Ω
1
L
Y
-
x
,
y
,
t
T
Y
-
ξ
,
ζ
,
t
×
Y
ξ
,
ζ
,
1
⋯
Y
ξ
,
ζ
,
L
T
γ
i
1
γ
i
2
⋯
γ
i
L
d
ξ
d
ζ
=
∬
Ω
1
L
Y
-
x
,
y
,
t
T
Y
-
ξ
,
ζ
,
t
Y
-
ξ
,
ζ
,
k
T
γ
i
d
ξ
d
ζ
=
λ
i
Y
-
x
,
y
,
t
T
γ
i
,
where Y
-
x
,
y
,
t
T
=
Y
x
,
y
,
1
,
Y
x
,
y
,
2
,
…
,
Y
x
,
y
,
L, Y
-
ξ
,
ζ
,
t
=
Y
ξ
,
ζ
,
1
,
Y
ξ
,
ζ
,
2
,
…
,
Y
ξ
,
ζ
,
L
T, and γ
i
=
γ
i
1
,
γ
i
2
,
…
,
γ
i
L
T.LetC
=
(
1
/
L
)
∬
Ω
Y
-
ξ
,
ζ
,
t
Y
-
ξ
,
ζ
,
k
T
d
ξ
d
ζ; (16) is transformed into
(17)
Y
-
x
,
y
,
t
T
C
γ
i
-
λ
i
γ
i
=
0
.We can get(18)
C
γ
i
=
λ
i
γ
i
.Thus, the problem of solving integral equation (14) is converted to solve the eigenvalue of (18). In this equation, γ
i is the number i eigenvector, C can be obtained from the experimental data, and λ
i is the eigenvalue of matrix C; then γ
i is determined in terms of the above parameters. Since matrix C is symmetric positive semidefinite and the spatial basis function which is obtained from (18) is orthogonal, the final spatial basis function can be obtained after normalizing.When determining the spatial basis function number ofn, we can sort the outputs as the dominant vector and select the appropriate dominant vector to reduce dimensions. In general, the bigger n is, the higher accuracy the system has. But the degree of state model will get higher at the same time; this will add difficulty to calculation. As long as we select the appropriate degree, the accuracy and calculated amount can be compromised, and the best model would be acquired.


### 3.2. STDARX Modeling of Polymer Flooding

The problem needs to be transformed to time iterative model when solved by DP. While ARX [19] can be used to build the time model; here the multivariable ARX is adopted to model the time coefficients. The specific description is as follows:
(19)
T
t
=
A
q
-
1
T
t
+
B
q
-
1
u
t
A
q
-
1
=
A
1
q
-
1
+
A
2
q
-
2
+
⋯
+
A
n
T
q
-
n
T
B
q
-
1
=
B
1
q
-
1
+
B
2
q
-
2
+
⋯
+
B
n
u
q
-
n
u
,
where T
t
=
T
1
t
,
T
2
t
,
…
,
T
n
T
t
T, q
-
1 is the backward shift operator which means backing up one sampling period, A
,
B are the matrix polynomial of q
-
1, u
t
=
u
1
t
,
u
2
t
,
…
,
u
n
u
t
T, and n
u and n
T are referred to as maximum lag with respect to the injection concentration and time coefficient, respectively.Convert the above problem to the linear regression model(20)
T
t
=
θ
·
Φ
t
,
(21)
Φ
t
=
T
t
-
1
T
,
…
,
T
t
-
n
T
T
T
u
t
-
1
T
,
…
,
u
t
-
n
u
T
T
,
(22)
θ
=
A
1
,
A
2
,
…
,
A
n
T
,
B
1
,
B
2
,
…
,
B
n
u
∈
R
n
×
n
n
T
+
m
n
u
.The recursive least-squares (RLS) shown in (23) are adopted to identify parameter θ in (22):
(23)
θ
^
t
=
θ
^
t
-
1
+
K
t
T
t
-
Φ
t
T
θ
^
t
-
1
,
K
t
=
P
t
-
1
Φ
t
Φ
t
T
P
t
-
1
Φ
t
+
μ
-
1
,
P
t
=
1
μ
I
-
K
t
Φ
t
T
P
t
-
1
,
where K
t denotes the weight matrix, 0
<
μ
<
1 is the forgetting factors, and P
t denotes the positive definite covariance matrix.After identifying the time coefficients, combining the spatial basis function solved in Section3.1, we can get the complete state parameter model:
(24)
T
^
t
=
∑
i
=
1
n
T
A
^
i
T
^
t
-
i
+
∑
j
=
1
n
u
B
^
j
u
t
-
j
,
Y
^
n
x
,
y
,
t
=
∑
i
=
1
n
φ
i
x
,
y
T
^
i
t
.The ARX models of state pressure, water saturation, and grid concentration can be obtained as (24). Then we take all states (pressure, water saturation, and grid concentration) as input variables and moisture content of production wells as the output variables. Do multivariable identification and get the relations between states and moisture content. Consider
(25)
f
w
t
=
U
q
-
1
c
p
t
+
V
q
-
1
S
w
t
+
W
q
-
1
p
t
,
U
q
-
1
=
∑
i
=
1
n
c
U
i
q
-
i
,
V
q
-
1
=
∑
i
=
1
n
s
V
i
q
-
i
,
W
q
-
1
=
∑
i
=
1
n
p
p
i
q
-
i
,
where f
w denotes the moisture content, c
p denotes the grid concentration, S
w denotes the water saturation, p denotes pressure, n
p
,
n
s
,
n
c and W
,
V
,
U denote the degrees and identification coefficients of pressure, water saturation, and grid concentration, respectively.Do parameter identification with RLS to (25), and get the ARX model of f
w. Then combine it with state models (pressure, water saturation, and grid concentration), the rounded STDARX model of polymer flooding is obtained. The model is shown in
(26)
f
w
l
t
=
U
q
-
1
c
p
t
+
V
q
-
1
S
w
t
+
W
q
-
1
p
t
,
l
=
1,2
,
…
,
9
(27)
T
^
m
t
=
∑
i
=
1
n
T
m
A
^
i
T
^
m
t
-
i
+
∑
j
=
1
n
u
m
B
^
j
u
t
-
j
Y
^
n
m
x
,
y
,
t
=
∑
i
=
1
n
m
φ
i
m
x
,
y
T
^
i
m
t
,
where n
T
m and n
u
m denote the degree of model, N
w and N
p denote the quantity of injection well and production well, respectively, Δ
T is the time step, m is the state identification parameter, m
=
1,2
,
3, m
=
1 denotes the model of state water saturation, m
=
2 denotes the model of state grid concentration, m
=
3 denotes the model of state pressure, and the rest parameters are the coefficients.


### 3.3. Model Verification for Polymer Flooding

#### 3.3.1. Reservoir Description

Suppose the oil field consists of four injections and nine production wells. All wells distribute uniformly; there is one injection well at the center of every four production wells. The specific distribution is shown in Figure1.Figure 1
The distribution diagram of well position.As to the oil field, the length is 630 m, the width is 630 m, and the thickness is 19.990 m. There are 7 layers in all, the thickness of each layer is 2.857 m, and the net thickness is 1.4286 m, the depth of upper surface is 2420 m, the porosity of every layer is 0.3, and the pore volume is 1.1097 × 106 m3. In order to simplify the modeling process, we only choose the first layer to research. The permeability distribution of the first layer in direction x
,
y that has been revised by RBF-HGRA [20] is shown in Figure 2. The initial grid concentration is 0 g/L. The grids of oil field are divided by three directions x
,
y
,
z, and the number is 21 × 21 × 7. The data of flow property is shown in Table 1.Table 1
Flow data.
—
Sign
Mean
Value
Sign
Mean
Value
Oil
ρ
o
Density(kg/
m
3
)
950
μ
o
Viscosity(mPa·s)
50
C
o
Compression coefficient(1/MPa)
5
×
1
0
-
6
B
o
Volume coefficient
1
Water
ρ
w
Density(kg/
m
3)
1000
μ
w
Viscosity(mPa·s)
0.458
C
w
Compression coefficient (1/MPa)
4
×
1
0
-
6
B
w
Volume coefficient
1
Polymer
a
Adsorption constant(c
m
3
/g)
0.03
b
Adsorption constant (L/g)
3.8
D
Diffusion coefficient (m
2
/s)
1
×
1
0
-
5
f
a
Attainable pore volume coefficient
1
R
k
max
⁡
Permeability reduction factor
1.15
β
b
r
k
p
Permeability reduction factor
1.2
a
p
1
Viscosity constant
15.426
a
p
2
Viscosity constant
0.4228
a
p
3
Viscosity constant
0.2749
Others
NaCl
Concentration(g/
m
3)
0.017
C
a
2+
M
g
2+
Concentration(g/
m
3)
0.00045Figure 2
Permeability of the first layer inx
,
y directions.The water injection rate of each injection well is 83 m3/d. The injection to production scheme is water flooding first and polymer flooding later. The whole process is from December 31, 2005, to December 31, 2013, sampled by month, and the total time nodes are 97 and the total space nodes are 441. The injection process of polymer is from October 1, 2006, to September 1, 2010. Three-slug injection scheme is adopted here. Four injection wells are I
4
-
251
,
I
4
-
256
,
I
5
-
234
,
I
5
-
263, and the unit is kg/m3. The injection concentration and the corresponding moisture content of production well are shown in Figures 3 and 4, respectively.Figure 3
Three-slug injection scheme.Figure 4
The moisture content of production wells.


#### 3.3.2. Modeling and Verification for Polymer Flooding

Given the polymer flooding above, we build the STDARX model with the method shown in Section3.2. In order to simplify the calculation, we only model for the first layer. The total time nodes are 97, and the total space nodes are 441.When determining the quantity of the spatial basis function, we define the mean square error (mse)
(28)
mse
=
∫
∑
e
x
,
y
,
t
2
d
x
∫
d
x
∑
Δ
t
.
Generally speaking, the more quantity the spatial basis function is, the more information of system the model reflects, and the higher accuracy the model has. But at the same time, the model will become too sensitive and the generalization capability will decrease. The dimension of model is high. On the other hand, the less quantity the spatial basis function is, the less the computation is; the model will reflect less information and the error can be very big. We do some tests to water saturation at different numbers of spatial basis function. The results are shown in Table 2; here the degrees of ARX model n
T
m, n
u
m are all equal to four.Table 2
mse for different number of spatial basis functions.
n
1
2
3
4
5
mse
0.669
0.478
0.225
0.098
0.087According to Table2, when n is less than 4, the error varies seriously with n. But when n is greater than 4, the error varies slightly with n. Four spatial basis functions can already reflect the information of system, so we let n
=
4 when modeling.To test the effect of STDARX model, we letu
I
4
-
251
*
=
2.3
1.9
0.8, u
I
4
-
256
*
=
2.1
2.4
1.2, u
I
5
-
234
*
=
3
2.1
1.6, and u
I
5
-
263
*
=
2.8
2.6
2, and the unit is kg/m3. Under the condition of fixed time node, we test the moisture content of production wells. The error between STDARX model output and mechanism model output for f
w is shown in Figure 5. The error between STDARX model and mechanism model for every state is shown in Figures 6, 7, and 8.Figure 5
The error of moisture content.Figure 6
The error diagram of water saturation.Figure 7
The error diagram of pressure.Figure 8
The error diagram of state polymer grid concentration.According to Figure5, the output of STDARX model is very close to the theoretical output. The error is very small. We define the below index st to evaluate the absolute value of mean value between STDARX model output and theoretical output:
(29)
st
=
1
N
∑
i
=
1
N
Y
i
-
Y
^
i
.After computing, the mean absolute errorst
=
0.0113
%, and the unit of f
w is %. The maximum st
=
2.161
%, mse
=
0.125. According to above data, the generalization of STDARX model is good; its effect is very close to the theoretical model.Since it is difficult to display the relation between two-dimensional plane error and time in three-dimensional space, here we only choose the middle line ofy direction to show for simplicity. According to Figures 6, 7, and 8, the errors are all very small for every state parameter. This illustrates the accuracy of STDARX model further. In addition, the injecting time is only from October 1, 2006, to September 1, 2010, 47 months in all. While the total time is 97 months, this leads to the error of f
w being smaller than real value. The other 50 months are water flooding, the model is not different from mechanism model, and the error is zero. There are two production wells in the middle line of y direction, so there could be big fluctuation at this place to state parameters.In conclusion, the spatial-temporal decomposition ARX (STDARX) model of polymer flooding is credible, and it has preferable accuracy and robustness.


## 3.1. The Two-Dimensional Spatial-Temporal Decomposition of State Parameters

In this paper, the K-L decomposition [16] method is adopted to build the two-dimensional distributed parameter spatial-temporal model of polymer flooding STDARX model [17]. As to this model, the input variables are pressure, water saturation, and grid concentration. The output variable is moisture content.Suppose the output of system isY
x
i
,
y
j
,
t
i
=
1
,
j
=
1
,
t
=
1
N
x
,
N
y
,
L, t
=
1,2
,
…
,
L denotes L individuals discrete points of output in time domain, i
=
1,2
,
…
,
N
x, j
=
1,2
,
…
,
N
y is N
x
×
N
y individuals state sample points in space.The output can be decomposed as follows:(9)
Y
x
,
y
,
t
=
∑
i
=
1
∞
φ
i
x
,
y
T
i
t
,
where φ
i
x
,
y
i
=
1
∞ denotes the spatial basis function, T
i
t
i
=
1
∞ denotes the time coefficient, and Y
x
,
y
,
t is the system output.The output is usually described as the finite form in specific engineering application:(10)
Y
n
x
,
y
,
t
=
∑
i
=
1
n
φ
i
x
,
y
T
i
t
,
where n mirrors the system energy; it depends the number of spatial basis function.The spatial basis function is unit orthogonal, which is defined as below:(11)
φ
i
x
,
y
,
φ
j
x
,
y
=
∬
Ω
φ
i
x
,
y
φ
j
x
,
y
d
x
d
y
=
0
,
i
≠
j
1
,
i
=
j
,
where f
x
,
g
x
=
∬
Ω
f
x
g
x
d
x
d
y is the inner product.If we know the spatial basis function, we can get the time coefficient:(12)
T
i
t
=
φ
i
x
,
y
,
Y
x
,
y
,
t
,
i
=
1,2
,
…
,
n
.Convert the problem of solving spatial basis function to the below optimization problem:(13)
min
⁡
φ
i
x
,
y
Y
x
,
y
,
t
-
Y
n
x
,
y
,
t
2
s
.
t
.
φ
i
,
φ
i
=
1
,
φ
i
∈
L
2
Ω
,
i
=
1,2
,
…
,
n
,
where f
x
,
y
=
f
x
,
y
,
f
x
,
y
1
/
2 is the norm and f
x
,
y
,
t
=
(
1
/
L
)
∑
t
=
1
L
f
x
,
y
,
t is the set average.The necessary condition of extremum existing to (13) is
(14)
∬
Ω
R
x
,
y
,
ξ
,
ζ
φ
i
ξ
,
ζ
d
ξ
d
ζ
=
λ
i
φ
i
x
,
y
φ
i
,
φ
i
=
1
,
i
=
1,2
,
…
,
n
,
where R
x
,
y
,
ξ
,
ζ
=
Y
x
,
y
,
t
Y
ξ
,
ζ
,
t is the correlation function between two points in the same plane. Thus the optimization problem is transformed to solve the problem described as (14).In normal conditions, (14) is very difficult to solve. But the snapshots [18] method can be applied when time nodes are less than space nodes; the calculated amount is cut down at the same time. The meaning of snapshots is as follows. When we choose a time node in time domain, all the output sampling values in two-dimensional plane at this time make up one snapshot. For example, Y
x
,
y
,
3, x
=
1,2
,
…
,
N
x, y
=
1,2
,
…
,
N
y, and N
=
N
x
·
N
y denote one snapshot of time node t
=
3.On the condition that time nodes are less than space nodes, let spatial basis function be expressed by a series of snapshots linear combination. Consider(15)
φ
i
x
,
y
=
∑
t
=
1
L
γ
i
t
Y
x
,
y
,
t
.In terms of (14), we can get
(16)
∬
Ω
R
x
,
y
,
ξ
,
ζ
φ
i
ξ
,
ζ
d
ξ
d
ζ
=
∬
Ω
1
L
Y
-
x
,
y
,
t
T
Y
-
ξ
,
ζ
,
t
×
Y
ξ
,
ζ
,
1
⋯
Y
ξ
,
ζ
,
L
T
γ
i
1
γ
i
2
⋯
γ
i
L
d
ξ
d
ζ
=
∬
Ω
1
L
Y
-
x
,
y
,
t
T
Y
-
ξ
,
ζ
,
t
Y
-
ξ
,
ζ
,
k
T
γ
i
d
ξ
d
ζ
=
λ
i
Y
-
x
,
y
,
t
T
γ
i
,
where Y
-
x
,
y
,
t
T
=
Y
x
,
y
,
1
,
Y
x
,
y
,
2
,
…
,
Y
x
,
y
,
L, Y
-
ξ
,
ζ
,
t
=
Y
ξ
,
ζ
,
1
,
Y
ξ
,
ζ
,
2
,
…
,
Y
ξ
,
ζ
,
L
T, and γ
i
=
γ
i
1
,
γ
i
2
,
…
,
γ
i
L
T.LetC
=
(
1
/
L
)
∬
Ω
Y
-
ξ
,
ζ
,
t
Y
-
ξ
,
ζ
,
k
T
d
ξ
d
ζ; (16) is transformed into
(17)
Y
-
x
,
y
,
t
T
C
γ
i
-
λ
i
γ
i
=
0
.We can get(18)
C
γ
i
=
λ
i
γ
i
.Thus, the problem of solving integral equation (14) is converted to solve the eigenvalue of (18). In this equation, γ
i is the number i eigenvector, C can be obtained from the experimental data, and λ
i is the eigenvalue of matrix C; then γ
i is determined in terms of the above parameters. Since matrix C is symmetric positive semidefinite and the spatial basis function which is obtained from (18) is orthogonal, the final spatial basis function can be obtained after normalizing.When determining the spatial basis function number ofn, we can sort the outputs as the dominant vector and select the appropriate dominant vector to reduce dimensions. In general, the bigger n is, the higher accuracy the system has. But the degree of state model will get higher at the same time; this will add difficulty to calculation. As long as we select the appropriate degree, the accuracy and calculated amount can be compromised, and the best model would be acquired.


## 3.2. STDARX Modeling of Polymer Flooding

The problem needs to be transformed to time iterative model when solved by DP. While ARX [19] can be used to build the time model; here the multivariable ARX is adopted to model the time coefficients. The specific description is as follows:
(19)
T
t
=
A
q
-
1
T
t
+
B
q
-
1
u
t
A
q
-
1
=
A
1
q
-
1
+
A
2
q
-
2
+
⋯
+
A
n
T
q
-
n
T
B
q
-
1
=
B
1
q
-
1
+
B
2
q
-
2
+
⋯
+
B
n
u
q
-
n
u
,
where T
t
=
T
1
t
,
T
2
t
,
…
,
T
n
T
t
T, q
-
1 is the backward shift operator which means backing up one sampling period, A
,
B are the matrix polynomial of q
-
1, u
t
=
u
1
t
,
u
2
t
,
…
,
u
n
u
t
T, and n
u and n
T are referred to as maximum lag with respect to the injection concentration and time coefficient, respectively.Convert the above problem to the linear regression model(20)
T
t
=
θ
·
Φ
t
,
(21)
Φ
t
=
T
t
-
1
T
,
…
,
T
t
-
n
T
T
T
u
t
-
1
T
,
…
,
u
t
-
n
u
T
T
,
(22)
θ
=
A
1
,
A
2
,
…
,
A
n
T
,
B
1
,
B
2
,
…
,
B
n
u
∈
R
n
×
n
n
T
+
m
n
u
.The recursive least-squares (RLS) shown in (23) are adopted to identify parameter θ in (22):
(23)
θ
^
t
=
θ
^
t
-
1
+
K
t
T
t
-
Φ
t
T
θ
^
t
-
1
,
K
t
=
P
t
-
1
Φ
t
Φ
t
T
P
t
-
1
Φ
t
+
μ
-
1
,
P
t
=
1
μ
I
-
K
t
Φ
t
T
P
t
-
1
,
where K
t denotes the weight matrix, 0
<
μ
<
1 is the forgetting factors, and P
t denotes the positive definite covariance matrix.After identifying the time coefficients, combining the spatial basis function solved in Section3.1, we can get the complete state parameter model:
(24)
T
^
t
=
∑
i
=
1
n
T
A
^
i
T
^
t
-
i
+
∑
j
=
1
n
u
B
^
j
u
t
-
j
,
Y
^
n
x
,
y
,
t
=
∑
i
=
1
n
φ
i
x
,
y
T
^
i
t
.The ARX models of state pressure, water saturation, and grid concentration can be obtained as (24). Then we take all states (pressure, water saturation, and grid concentration) as input variables and moisture content of production wells as the output variables. Do multivariable identification and get the relations between states and moisture content. Consider
(25)
f
w
t
=
U
q
-
1
c
p
t
+
V
q
-
1
S
w
t
+
W
q
-
1
p
t
,
U
q
-
1
=
∑
i
=
1
n
c
U
i
q
-
i
,
V
q
-
1
=
∑
i
=
1
n
s
V
i
q
-
i
,
W
q
-
1
=
∑
i
=
1
n
p
p
i
q
-
i
,
where f
w denotes the moisture content, c
p denotes the grid concentration, S
w denotes the water saturation, p denotes pressure, n
p
,
n
s
,
n
c and W
,
V
,
U denote the degrees and identification coefficients of pressure, water saturation, and grid concentration, respectively.Do parameter identification with RLS to (25), and get the ARX model of f
w. Then combine it with state models (pressure, water saturation, and grid concentration), the rounded STDARX model of polymer flooding is obtained. The model is shown in
(26)
f
w
l
t
=
U
q
-
1
c
p
t
+
V
q
-
1
S
w
t
+
W
q
-
1
p
t
,
l
=
1,2
,
…
,
9
(27)
T
^
m
t
=
∑
i
=
1
n
T
m
A
^
i
T
^
m
t
-
i
+
∑
j
=
1
n
u
m
B
^
j
u
t
-
j
Y
^
n
m
x
,
y
,
t
=
∑
i
=
1
n
m
φ
i
m
x
,
y
T
^
i
m
t
,
where n
T
m and n
u
m denote the degree of model, N
w and N
p denote the quantity of injection well and production well, respectively, Δ
T is the time step, m is the state identification parameter, m
=
1,2
,
3, m
=
1 denotes the model of state water saturation, m
=
2 denotes the model of state grid concentration, m
=
3 denotes the model of state pressure, and the rest parameters are the coefficients.


## 3.3. Model Verification for Polymer Flooding

### 3.3.1. Reservoir Description

Suppose the oil field consists of four injections and nine production wells. All wells distribute uniformly; there is one injection well at the center of every four production wells. The specific distribution is shown in Figure1.Figure 1
The distribution diagram of well position.As to the oil field, the length is 630 m, the width is 630 m, and the thickness is 19.990 m. There are 7 layers in all, the thickness of each layer is 2.857 m, and the net thickness is 1.4286 m, the depth of upper surface is 2420 m, the porosity of every layer is 0.3, and the pore volume is 1.1097 × 106 m3. In order to simplify the modeling process, we only choose the first layer to research. The permeability distribution of the first layer in direction x
,
y that has been revised by RBF-HGRA [20] is shown in Figure 2. The initial grid concentration is 0 g/L. The grids of oil field are divided by three directions x
,
y
,
z, and the number is 21 × 21 × 7. The data of flow property is shown in Table 1.Table 1
Flow data.
—
Sign
Mean
Value
Sign
Mean
Value
Oil
ρ
o
Density(kg/
m
3
)
950
μ
o
Viscosity(mPa·s)
50
C
o
Compression coefficient(1/MPa)
5
×
1
0
-
6
B
o
Volume coefficient
1
Water
ρ
w
Density(kg/
m
3)
1000
μ
w
Viscosity(mPa·s)
0.458
C
w
Compression coefficient (1/MPa)
4
×
1
0
-
6
B
w
Volume coefficient
1
Polymer
a
Adsorption constant(c
m
3
/g)
0.03
b
Adsorption constant (L/g)
3.8
D
Diffusion coefficient (m
2
/s)
1
×
1
0
-
5
f
a
Attainable pore volume coefficient
1
R
k
max
⁡
Permeability reduction factor
1.15
β
b
r
k
p
Permeability reduction factor
1.2
a
p
1
Viscosity constant
15.426
a
p
2
Viscosity constant
0.4228
a
p
3
Viscosity constant
0.2749
Others
NaCl
Concentration(g/
m
3)
0.017
C
a
2+
M
g
2+
Concentration(g/
m
3)
0.00045Figure 2
Permeability of the first layer inx
,
y directions.The water injection rate of each injection well is 83 m3/d. The injection to production scheme is water flooding first and polymer flooding later. The whole process is from December 31, 2005, to December 31, 2013, sampled by month, and the total time nodes are 97 and the total space nodes are 441. The injection process of polymer is from October 1, 2006, to September 1, 2010. Three-slug injection scheme is adopted here. Four injection wells are I
4
-
251
,
I
4
-
256
,
I
5
-
234
,
I
5
-
263, and the unit is kg/m3. The injection concentration and the corresponding moisture content of production well are shown in Figures 3 and 4, respectively.Figure 3
Three-slug injection scheme.Figure 4
The moisture content of production wells.


### 3.3.2. Modeling and Verification for Polymer Flooding

Given the polymer flooding above, we build the STDARX model with the method shown in Section3.2. In order to simplify the calculation, we only model for the first layer. The total time nodes are 97, and the total space nodes are 441.When determining the quantity of the spatial basis function, we define the mean square error (mse)
(28)
mse
=
∫
∑
e
x
,
y
,
t
2
d
x
∫
d
x
∑
Δ
t
.
Generally speaking, the more quantity the spatial basis function is, the more information of system the model reflects, and the higher accuracy the model has. But at the same time, the model will become too sensitive and the generalization capability will decrease. The dimension of model is high. On the other hand, the less quantity the spatial basis function is, the less the computation is; the model will reflect less information and the error can be very big. We do some tests to water saturation at different numbers of spatial basis function. The results are shown in Table 2; here the degrees of ARX model n
T
m, n
u
m are all equal to four.Table 2
mse for different number of spatial basis functions.
n
1
2
3
4
5
mse
0.669
0.478
0.225
0.098
0.087According to Table2, when n is less than 4, the error varies seriously with n. But when n is greater than 4, the error varies slightly with n. Four spatial basis functions can already reflect the information of system, so we let n
=
4 when modeling.To test the effect of STDARX model, we letu
I
4
-
251
*
=
2.3
1.9
0.8, u
I
4
-
256
*
=
2.1
2.4
1.2, u
I
5
-
234
*
=
3
2.1
1.6, and u
I
5
-
263
*
=
2.8
2.6
2, and the unit is kg/m3. Under the condition of fixed time node, we test the moisture content of production wells. The error between STDARX model output and mechanism model output for f
w is shown in Figure 5. The error between STDARX model and mechanism model for every state is shown in Figures 6, 7, and 8.Figure 5
The error of moisture content.Figure 6
The error diagram of water saturation.Figure 7
The error diagram of pressure.Figure 8
The error diagram of state polymer grid concentration.According to Figure5, the output of STDARX model is very close to the theoretical output. The error is very small. We define the below index st to evaluate the absolute value of mean value between STDARX model output and theoretical output:
(29)
st
=
1
N
∑
i
=
1
N
Y
i
-
Y
^
i
.After computing, the mean absolute errorst
=
0.0113
%, and the unit of f
w is %. The maximum st
=
2.161
%, mse
=
0.125. According to above data, the generalization of STDARX model is good; its effect is very close to the theoretical model.Since it is difficult to display the relation between two-dimensional plane error and time in three-dimensional space, here we only choose the middle line ofy direction to show for simplicity. According to Figures 6, 7, and 8, the errors are all very small for every state parameter. This illustrates the accuracy of STDARX model further. In addition, the injecting time is only from October 1, 2006, to September 1, 2010, 47 months in all. While the total time is 97 months, this leads to the error of f
w being smaller than real value. The other 50 months are water flooding, the model is not different from mechanism model, and the error is zero. There are two production wells in the middle line of y direction, so there could be big fluctuation at this place to state parameters.In conclusion, the spatial-temporal decomposition ARX (STDARX) model of polymer flooding is credible, and it has preferable accuracy and robustness.


## 3.3.1. Reservoir Description

Suppose the oil field consists of four injections and nine production wells. All wells distribute uniformly; there is one injection well at the center of every four production wells. The specific distribution is shown in Figure1.Figure 1
The distribution diagram of well position.As to the oil field, the length is 630 m, the width is 630 m, and the thickness is 19.990 m. There are 7 layers in all, the thickness of each layer is 2.857 m, and the net thickness is 1.4286 m, the depth of upper surface is 2420 m, the porosity of every layer is 0.3, and the pore volume is 1.1097 × 106 m3. In order to simplify the modeling process, we only choose the first layer to research. The permeability distribution of the first layer in direction x
,
y that has been revised by RBF-HGRA [20] is shown in Figure 2. The initial grid concentration is 0 g/L. The grids of oil field are divided by three directions x
,
y
,
z, and the number is 21 × 21 × 7. The data of flow property is shown in Table 1.Table 1
Flow data.
—
Sign
Mean
Value
Sign
Mean
Value
Oil
ρ
o
Density(kg/
m
3
)
950
μ
o
Viscosity(mPa·s)
50
C
o
Compression coefficient(1/MPa)
5
×
1
0
-
6
B
o
Volume coefficient
1
Water
ρ
w
Density(kg/
m
3)
1000
μ
w
Viscosity(mPa·s)
0.458
C
w
Compression coefficient (1/MPa)
4
×
1
0
-
6
B
w
Volume coefficient
1
Polymer
a
Adsorption constant(c
m
3
/g)
0.03
b
Adsorption constant (L/g)
3.8
D
Diffusion coefficient (m
2
/s)
1
×
1
0
-
5
f
a
Attainable pore volume coefficient
1
R
k
max
⁡
Permeability reduction factor
1.15
β
b
r
k
p
Permeability reduction factor
1.2
a
p
1
Viscosity constant
15.426
a
p
2
Viscosity constant
0.4228
a
p
3
Viscosity constant
0.2749
Others
NaCl
Concentration(g/
m
3)
0.017
C
a
2+
M
g
2+
Concentration(g/
m
3)
0.00045Figure 2
Permeability of the first layer inx
,
y directions.The water injection rate of each injection well is 83 m3/d. The injection to production scheme is water flooding first and polymer flooding later. The whole process is from December 31, 2005, to December 31, 2013, sampled by month, and the total time nodes are 97 and the total space nodes are 441. The injection process of polymer is from October 1, 2006, to September 1, 2010. Three-slug injection scheme is adopted here. Four injection wells are I
4
-
251
,
I
4
-
256
,
I
5
-
234
,
I
5
-
263, and the unit is kg/m3. The injection concentration and the corresponding moisture content of production well are shown in Figures 3 and 4, respectively.Figure 3
Three-slug injection scheme.Figure 4
The moisture content of production wells.


## 3.3.2. Modeling and Verification for Polymer Flooding

Given the polymer flooding above, we build the STDARX model with the method shown in Section3.2. In order to simplify the calculation, we only model for the first layer. The total time nodes are 97, and the total space nodes are 441.When determining the quantity of the spatial basis function, we define the mean square error (mse)
(28)
mse
=
∫
∑
e
x
,
y
,
t
2
d
x
∫
d
x
∑
Δ
t
.
Generally speaking, the more quantity the spatial basis function is, the more information of system the model reflects, and the higher accuracy the model has. But at the same time, the model will become too sensitive and the generalization capability will decrease. The dimension of model is high. On the other hand, the less quantity the spatial basis function is, the less the computation is; the model will reflect less information and the error can be very big. We do some tests to water saturation at different numbers of spatial basis function. The results are shown in Table 2; here the degrees of ARX model n
T
m, n
u
m are all equal to four.Table 2
mse for different number of spatial basis functions.
n
1
2
3
4
5
mse
0.669
0.478
0.225
0.098
0.087According to Table2, when n is less than 4, the error varies seriously with n. But when n is greater than 4, the error varies slightly with n. Four spatial basis functions can already reflect the information of system, so we let n
=
4 when modeling.To test the effect of STDARX model, we letu
I
4
-
251
*
=
2.3
1.9
0.8, u
I
4
-
256
*
=
2.1
2.4
1.2, u
I
5
-
234
*
=
3
2.1
1.6, and u
I
5
-
263
*
=
2.8
2.6
2, and the unit is kg/m3. Under the condition of fixed time node, we test the moisture content of production wells. The error between STDARX model output and mechanism model output for f
w is shown in Figure 5. The error between STDARX model and mechanism model for every state is shown in Figures 6, 7, and 8.Figure 5
The error of moisture content.Figure 6
The error diagram of water saturation.Figure 7
The error diagram of pressure.Figure 8
The error diagram of state polymer grid concentration.According to Figure5, the output of STDARX model is very close to the theoretical output. The error is very small. We define the below index st to evaluate the absolute value of mean value between STDARX model output and theoretical output:
(29)
st
=
1
N
∑
i
=
1
N
Y
i
-
Y
^
i
.After computing, the mean absolute errorst
=
0.0113
%, and the unit of f
w is %. The maximum st
=
2.161
%, mse
=
0.125. According to above data, the generalization of STDARX model is good; its effect is very close to the theoretical model.Since it is difficult to display the relation between two-dimensional plane error and time in three-dimensional space, here we only choose the middle line ofy direction to show for simplicity. According to Figures 6, 7, and 8, the errors are all very small for every state parameter. This illustrates the accuracy of STDARX model further. In addition, the injecting time is only from October 1, 2006, to September 1, 2010, 47 months in all. While the total time is 97 months, this leads to the error of f
w being smaller than real value. The other 50 months are water flooding, the model is not different from mechanism model, and the error is zero. There are two production wells in the middle line of y direction, so there could be big fluctuation at this place to state parameters.In conclusion, the spatial-temporal decomposition ARX (STDARX) model of polymer flooding is credible, and it has preferable accuracy and robustness.


## 4. Iterative Dynamic Programming Optimization for Polymer Flooding Model

### 4.1. Dynamic Programming Model Based on STDARX

Net present value (NPV) is adopted to evaluate the profit of polymer flooding [21]. The bigger NPV is, the better the profitability is. NPV is defined as
(30)
max
⁡
⁡
J
=
∫
0
t
f
C
in
t
-
C
out
t
1
+
r
-
t
d
t
,
where [
0
,
t
f
] denotes the running time period, r denotes the discount rate, 1
+
r
-
t denotes the discount coefficient, C
in
t denotes the income of crude oil, and C
out
t denotes the price of total polymer.The dynamic programming model based on STDARX is shown as follows.Performance index(31)
min
⁡
⁡
J
=
-
J
NPV
=
∑
t
=
0
T
∑
j
=
1
N
p
P
oil
t
q
out
1
-
f
w
t
1
+
r
-
t
Δ
T
/
t
0
Δ
T
·
∑
j
=
1
N
p
P
oil
t
q
out
1
-
f
w
t
∑
i
=
1
N
w
P
poly
q
in
c
p
in
t
-
∑
j
=
1
N
p
P
oil
t
q
out
1
-
f
w
t
.
State equation (STDARX model)(32)
f
w
l
t
=
U
q
-
1
c
p
t
+
V
q
-
1
S
w
t
+
W
q
-
1
p
t
,
l
=
1,2
,
…
,
9
,
T
^
m
t
=
∑
i
=
1
n
T
m
A
^
i
T
^
m
t
-
i
+
∑
j
=
1
n
u
m
B
^
j
u
t
-
j
,
Y
^
n
m
x
,
y
,
t
=
∑
i
=
1
n
m
φ
i
m
x
,
y
T
^
i
m
t
.
Control constraint(33)
0
≤
c
p
in
t
≤
c
max
⁡
,
t
=
1
,
…
,
T
0
≤
∑
t
=
1
T
c
i
t
q
i
≤
m
p
max
⁡
,where P
poly and P
oil
(
t
) are the oil price and polymer price, q
in and q
out are the water injection volume and oil production volume of one well, and c
p
in is the injection concentration of polymer. c
max
⁡ and m
p
max
⁡ are the upper bound of polymer concentration and polymer’s total amount, respectively. q
i denotes the migration velocity. The other parameters are in line with the before mentioned description.


### 4.2. Iterative Dynamic Programming Optimization with Varying Oil Price

LetN
w
=
4, N
p
=
9, P
poly
=
20 RMB/kg, c
max
⁡
=
3 kg/m3, and r
=
0.12. The crude oil price is predicted by the method of dynamic correcting ε
-SVR [22]. The unit is RMB
/
t, and the predicting result is shown in Figure 9.Figure 9
The price of crude oil.Utilize IDP to optimize the above model; the specific process is as follows.(1)
Divide the whole polymer injecting time (47 months) intop
=
3 period.
(2)
Parameter initialization. The initial constriction factorr
=
0.85, the number of discrete control variable R
=
5, choose the initial control strategy, and the reducing factor η
=
0.7. Iteration time is m and the circulation time is n.
(3)
Let the current circulation time beq, and do q
=
1. Let the current iteration time be j, and do j
=
1.
(4)
Do the current control region to ber
=
η
q
-
1
r
in.
(5)
Compute the state value of every period on the basis of current control strategy and initial time state. The state isx
k, k
=
1
,
…
,
p; do k
=
p
-
1.
(6)
Ifk
=
0, go to (6). If k
≠
0, get the current control strategy set as control strategy even generated method. Compute the state and its performance index J on the basis of the initial state x
(
k
) and the later period optimal control strategy which is obtained from the last calculation. Sort J and get its minimum value; its control strategy is the optimal control strategy. Do k
=
k
-
1, and go back to (5).
(7)
Zoom out the control region asr
=
γ
r, do j
=
j
+
1 until j
=
m, end this cycle, and do j
=
1.
(8)
Doq
=
q
+
1 until q
=
n. End the whole process.The injection is three-slug form, the time node is fixed, and the injection rate remains the same; we only optimize the injection concentration. Here we give an initial concentration to STDARX model,u
=
1.8
1.8
1.8
1.8 kg/m3 (all injection wells have the same concentration). After being solved by IDP, the results are shown in Figures 10, 11, 12, and 13.Figure 10
The comparison diagram ofc
p
in.Figure 11
The unoptimizedf
w and oil production.Figure 12
The optimizedf
w and oil production of production well.Figure 13
The averagef
w comparison of production well.The optimal injection concentration is shown in Figure10, and the optimized result of c
p
in is u
I
4
-
251
*
=
3
1.5
2.8 kg/m3,  u
I
4
-
256
*
=
2.4
2.1
1.2 kg/m3, u
I
5
-
234
*
=
1.5
2.2
2.9 kg/m3, and u
I
5
-
263
*
=
1.9
2.6
3 kg/m3. Figures 11 and 12 are the unoptimized and optimized for both moisture content and oil production, respectively.To strengthen the comparison effect, we take the mean value of moisture content at all time node and display them in Figure13. We can conclude that the f
w of optimized scheme decreases much less than the initial concentration obviously.On the condition of initial injection concentration, the whole consumed polymer is 842616 kg, oil production is 69018.91 m3, and net profit is 22.113 million RMB. While on the condition of optimal injection concentration, the whole consumed polymer is 1024500 kg, oil production is 76016.11 m3, and net profit is 24.176 million RMB. On the whole, increase 181884 kg polymer, but gain additional 6997.19 m3 crude oil and extra 2.063 million RMB profit. The oil recovery is enhanced by 9.2%.


## 4.1. Dynamic Programming Model Based on STDARX

Net present value (NPV) is adopted to evaluate the profit of polymer flooding [21]. The bigger NPV is, the better the profitability is. NPV is defined as
(30)
max
⁡
⁡
J
=
∫
0
t
f
C
in
t
-
C
out
t
1
+
r
-
t
d
t
,
where [
0
,
t
f
] denotes the running time period, r denotes the discount rate, 1
+
r
-
t denotes the discount coefficient, C
in
t denotes the income of crude oil, and C
out
t denotes the price of total polymer.The dynamic programming model based on STDARX is shown as follows.Performance index(31)
min
⁡
⁡
J
=
-
J
NPV
=
∑
t
=
0
T
∑
j
=
1
N
p
P
oil
t
q
out
1
-
f
w
t
1
+
r
-
t
Δ
T
/
t
0
Δ
T
·
∑
j
=
1
N
p
P
oil
t
q
out
1
-
f
w
t
∑
i
=
1
N
w
P
poly
q
in
c
p
in
t
-
∑
j
=
1
N
p
P
oil
t
q
out
1
-
f
w
t
.
State equation (STDARX model)(32)
f
w
l
t
=
U
q
-
1
c
p
t
+
V
q
-
1
S
w
t
+
W
q
-
1
p
t
,
l
=
1,2
,
…
,
9
,
T
^
m
t
=
∑
i
=
1
n
T
m
A
^
i
T
^
m
t
-
i
+
∑
j
=
1
n
u
m
B
^
j
u
t
-
j
,
Y
^
n
m
x
,
y
,
t
=
∑
i
=
1
n
m
φ
i
m
x
,
y
T
^
i
m
t
.
Control constraint(33)
0
≤
c
p
in
t
≤
c
max
⁡
,
t
=
1
,
…
,
T
0
≤
∑
t
=
1
T
c
i
t
q
i
≤
m
p
max
⁡
,where P
poly and P
oil
(
t
) are the oil price and polymer price, q
in and q
out are the water injection volume and oil production volume of one well, and c
p
in is the injection concentration of polymer. c
max
⁡ and m
p
max
⁡ are the upper bound of polymer concentration and polymer’s total amount, respectively. q
i denotes the migration velocity. The other parameters are in line with the before mentioned description.


## 4.2. Iterative Dynamic Programming Optimization with Varying Oil Price

LetN
w
=
4, N
p
=
9, P
poly
=
20 RMB/kg, c
max
⁡
=
3 kg/m3, and r
=
0.12. The crude oil price is predicted by the method of dynamic correcting ε
-SVR [22]. The unit is RMB
/
t, and the predicting result is shown in Figure 9.Figure 9
The price of crude oil.Utilize IDP to optimize the above model; the specific process is as follows.(1)
Divide the whole polymer injecting time (47 months) intop
=
3 period.
(2)
Parameter initialization. The initial constriction factorr
=
0.85, the number of discrete control variable R
=
5, choose the initial control strategy, and the reducing factor η
=
0.7. Iteration time is m and the circulation time is n.
(3)
Let the current circulation time beq, and do q
=
1. Let the current iteration time be j, and do j
=
1.
(4)
Do the current control region to ber
=
η
q
-
1
r
in.
(5)
Compute the state value of every period on the basis of current control strategy and initial time state. The state isx
k, k
=
1
,
…
,
p; do k
=
p
-
1.
(6)
Ifk
=
0, go to (6). If k
≠
0, get the current control strategy set as control strategy even generated method. Compute the state and its performance index J on the basis of the initial state x
(
k
) and the later period optimal control strategy which is obtained from the last calculation. Sort J and get its minimum value; its control strategy is the optimal control strategy. Do k
=
k
-
1, and go back to (5).
(7)
Zoom out the control region asr
=
γ
r, do j
=
j
+
1 until j
=
m, end this cycle, and do j
=
1.
(8)
Doq
=
q
+
1 until q
=
n. End the whole process.The injection is three-slug form, the time node is fixed, and the injection rate remains the same; we only optimize the injection concentration. Here we give an initial concentration to STDARX model,u
=
1.8
1.8
1.8
1.8 kg/m3 (all injection wells have the same concentration). After being solved by IDP, the results are shown in Figures 10, 11, 12, and 13.Figure 10
The comparison diagram ofc
p
in.Figure 11
The unoptimizedf
w and oil production.Figure 12
The optimizedf
w and oil production of production well.Figure 13
The averagef
w comparison of production well.The optimal injection concentration is shown in Figure10, and the optimized result of c
p
in is u
I
4
-
251
*
=
3
1.5
2.8 kg/m3,  u
I
4
-
256
*
=
2.4
2.1
1.2 kg/m3, u
I
5
-
234
*
=
1.5
2.2
2.9 kg/m3, and u
I
5
-
263
*
=
1.9
2.6
3 kg/m3. Figures 11 and 12 are the unoptimized and optimized for both moisture content and oil production, respectively.To strengthen the comparison effect, we take the mean value of moisture content at all time node and display them in Figure13. We can conclude that the f
w of optimized scheme decreases much less than the initial concentration obviously.On the condition of initial injection concentration, the whole consumed polymer is 842616 kg, oil production is 69018.91 m3, and net profit is 22.113 million RMB. While on the condition of optimal injection concentration, the whole consumed polymer is 1024500 kg, oil production is 76016.11 m3, and net profit is 24.176 million RMB. On the whole, increase 181884 kg polymer, but gain additional 6997.19 m3 crude oil and extra 2.063 million RMB profit. The oil recovery is enhanced by 9.2%.


## 5. Conclusion

In this paper, a new modeling method for polymer flooding based on spatial-temporal decomposition and ARX is proposed. The main conclusions obtained from this study are as follows.(1)
The STDARX model of polymer flooding is established. This is the first time to build the distributed polymer flooding model without considering mechanism equations in this area. It avoids the solution for complex mechanism model and simplifies the calculation. Through verifying, the new model has good accuracy and generalization.
(2)
A new two-dimensional spatial-temporal decomposition idea is put forward. By use of ARX and RLS, this idea has good effect in practice. The performance can be adjusted artificially by changing the number of spatial basis functions and the degree of ARX.
(3)
IDP is introduced to optimize STDARX model of polymer flooding. NPV is treated as the evaluation index, and the variation of oil price with time is considered at the same time. Through optimization, we get the optimal injection concentration and improve the moisture content and oil production obviously.
(4)
In this paper, we only build two-dimensional model for simplifying. But it is not so reasonable, because reservoir modeling is a distribute-state-spatial model, the three-dimensional model can also be established in future research.



---

*Source: 713091-2014-10-30.xml*