Spatial-Temporal ARX Modeling and Optimization for Polymer Flooding

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, andmoisture 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.


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 2 Mathematical Problems in Engineering 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.

Mechanism Model Description of Polymer Flooding
Let Ω ∈  ] . ( The flow equation for water phase The flow equation for polymer component where   is the reference pressure (MPa) and   ,   , and   denote the porosity, the oil, and water volume factors under the condition of   , respectively.  is the effective pore volume coefficient and   ,   , and   denote the compressibility factors of oil, water, and rock, respectively.
The oil and water relative permeability   and   are the functions of water saturation   and can be described as where   denotes oil relative permeability at the irreducible water saturation   and   denotes the water relative permeability at the residual oil saturation   .The polymer solution viscosity   (mPa⋅s), the permeability reduction factor   , and the amount absorbed per unit mass of the rock   (mg/g) are the functions of polymer concentration   .Consider where   denotes the viscosity of water with no polymer (mPa⋅s) and  1 ,  2 ,  3 ,   max ,   , ,  are the constant parameters.
The moisture content   of production wells depends on ,   , and   , and it can be defined as When applied in computing the polymer flooding, the model can obtain the   of production wells and get the oil production.

Polymer Flooding Model with
Spatial-Temporal ARX Method

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 is {(  ,   , )} The output can be decomposed as follows: where {  (, )} ∞ =1 denotes the spatial basis function, {  ()} ∞ =1 denotes the time coefficient, and (, , ) is the system output.
The output is usually described as the finite form in specific engineering application: where  mirrors the system energy; it depends the number of spatial basis function.
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, (, , 3),  = 1, 2, . . .,   ,  = 1, 2, . . .,   , and  =   ⋅   denote one snapshot of time node  = 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 In terms of ( 14), we can get where (, , )  = ((, , 1), (, , 2), . . ., (, , )), ((, , )) = ((, , 1), (, , 2), . . ., (, , ))  , and We can get Thus, the problem of solving integral equation ( 14) is converted to solve the eigenvalue of (18).In this equation,   is the number  eigenvector,  can be obtained from the experimental data, and   is the eigenvalue of matrix ; then   is determined in terms of the above parameters.Since matrix  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 of , we can sort the outputs as the dominant vector and select the appropriate dominant vector to reduce dimensions.In general, the bigger  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.

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: where () = ( 1 (),  2 (), . . .,    ())  ,  −1 is the backward shift operator which means backing up one sampling period, ,  are the matrix polynomial of  −1 , () = ( 1 (),  2 (), . . .,    ())  , and   and   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 The recursive least-squares (RLS) shown in (23) are adopted to identify parameter  in (22): where () denotes the weight matrix, 0 <  < 1 is the forgetting factors, and () denotes the positive definite covariance matrix.
After identifying the time coefficients, combining the spatial basis function solved in Section 3.1, we can get the complete state parameter model: 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   () =  ( −1 )   () +  ( −1 )   () +  ( −1 )  () , where   denotes the moisture content,   denotes the grid concentration,   denotes the water saturation,  denotes pressure,   ,   ,   and , ,  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   .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 +  ( −1 )  () ,  = 1, 2, . . ., 9 where    and    denote the degree of model,   and   denote the quantity of injection well and production well, respectively, Δ is the time step,  is the state identification parameter,  = 1, 2, 3,  = 1 denotes the model of state water saturation,  = 2 denotes the model of state grid concentration,  = 3 denotes the model of state pressure, and the rest parameters are the coefficients.

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 Figure 1.
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 × 10 6 m 3 .In order to simplify the modeling process, we only choose the first layer to research.The permeability distribution of the first layer in direction ,  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 , , , and the number is 21 × 21 × 7. The data of flow property is shown in Table 1.
The water injection rate of each injection well is 83 m 3 /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 4-251, 4-256, 5-234, 5-263, and the unit is kg/m 3 .The injection concentration and  the corresponding moisture content of production well are shown in Figures 3 and 4, respectively.

Modeling and Verification for Polymer Flooding.
Given the polymer flooding above, we build the STDARX model with the method shown in Section 3.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) 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    ,    are all equal to four.According to Table 2, when  is less than 4, the error varies seriously with .But when  is greater than 4, the error varies slightly with .Four spatial basis functions can already reflect the information of system, so we let  = 4 when modeling.
To test the effect of STDARX model, we let  * 5-234 = (3 2.1 1.6), and  * 5-263 = (2.8 2.6 2), and the unit is kg/m 3 .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   is shown in Figure 5.The error between STDARX model and mechanism model for every state is shown in Figures 6, 7, and 8.
According to Figure 5, 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: After computing, the mean absolute error st = 0.0113%, and the unit of   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 twodimensional plane error and time in three-dimensional space, here we only choose the middle line of  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   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  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.

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 where [0,   ] denotes the running time period,  denotes the discount rate, (1 + ) − denotes the discount coefficient,  in () denotes the income of crude oil, and  out () denotes the price of total polymer.The dynamic programming model based on STDARX is shown as follows.
Performance index where  poly and  oil () are the oil price and polymer price,  in and  out are the water injection volume and oil production volume of one well, and  in is the injection concentration of polymer. max and   max are the upper bound of polymer concentration and polymer's total amount, respectively.  denotes the migration velocity.The other parameters are in line with the before mentioned description.

Iterative Dynamic Programming Optimization with Varying Oil Price.
Let   = 4,   = 9,  poly = 20 RMB/kg,  max = 3 kg/m 3 , and  = 0.12.The crude oil price is predicted by the method of dynamic correcting -SVR [22].The unit is RMB/, and the predicting result is shown in Figure 9.
Utilize IDP to optimize the above model; the specific process is as follows.
(2) Parameter initialization.The initial constriction factor  = 0.85, the number of discrete control variable  = 5, choose the initial control strategy, and the reducing factor  = 0.7.Iteration time is  and the circulation time is .
(3) Let the current circulation time be , and do  = 1.
Let the current iteration time be , and do  = 1.
(4) Do the current control region to be  =  −1  in .
(5) Compute the state value of every period on the basis of current control strategy and initial time state.The state is (),  = 1, . . ., ; do  =  − 1.
(6) If  = 0, go to (6).If  ̸ = 0, get the current control strategy set as control strategy even generated method.Compute the state and its performance index  on the basis of the initial state () and the later period optimal control strategy which is obtained from the last calculation.Sort  and get its minimum value; its control strategy is the optimal control strategy.Do  =  − 1, and go back to (5).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 strengthen the comparison effect, we take the mean value of moisture content at all time node and display them in  Figure 13.We can conclude that the   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.91m 3 , 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.11m 3 , and net profit is 24.176 million RMB.On the whole, increase 181884 kg polymer, but gain additional 6997.19m 3 crude oil and extra 2.063 million RMB profit.The oil recovery is enhanced by 9.2%.

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.

Figure 4 :
Figure 4: The moisture content of production wells.

Figure 5 :Figure 6 :
Figure 5: The error of moisture content.

Figure 9 :
Figure 9: The price of crude oil.

Figure 10 :Figure 11 :
Figure 10: The comparison diagram of  in .

Figure 12 :
Figure 12: The optimized   and oil production of production well.

Figure 13 :
Figure 13: The average   comparison of production well.
2denote the plain domain of reservoir with boundary Ω,  the unit outward normal on Ω, and (, ) ∈ Ω the coordinate of a point in the reservoir.Given a terminal time   , there are   injection wells and   production wells.The injection and production wells are located at   = {(  ,   ) |  = 1, 2, . . .,   } and   = {(  ,   ) =   /    ,   =   /      ,   =     /  ,   =     /      , (, ) is the absolute permeability (m 2 ),  is the diffusion coefficient of polymer (m 2 /s), ℎ is the thickness of reservoir bed (m),   is the rock density (kg/m 3 ), and   is the oil viscosity (mPa⋅s).The oil and water volume factors   ,   , the rock porosity , and the effective porosity to polymer   are related to reservoir pressure .Consider

Table 2 :
mse for different number of spatial basis functions.