Enhanced Oil Recovery for ASP Flooding Based on Biorthogonal Spatial-Temporal Wiener Modeling and Iterative Dynamic Programming

Because of the mechanism complexity, coupling, and time-space characteristic of alkali-surfactant-polymer (ASP) flooding, common methods are very hard to be implemented directly. In this paper, an iterative dynamic programming (IDP) based on a biorthogonal spatial-temporal Wiener modeling method is developed to solve the enhanced oil recovery for ASP flooding. At first, a comprehensive mechanism model for the enhanced oil recovery of ASP flooding is introduced. Then the biorthogonal spatial-temporal Wiener model is presented to build the relation between inputs and states, in which the Wiener model is expanded on a set of spatial basis functions and temporal basis functions. After inferring the necessary condition of solutions, these basis functions are determined by the snapshot method. Furthermore, a theorem is proved to identify parameters in the Wiener model. Combined with the least square estimation (LSE), all unknown parameters are determined. In addition, the ARMA model is applied to build the model between states and outputs, whose parameters are identified by recursive least squares (RLS). Thus, the whole modeling process for ASP flooding is finished. At last, the IDP algorithm is applied to solve the enhanced oil recovery problem for ASP flooding based on the identification model to obtain the optimal injection strategy. Simulations on the ASP flooding with four injection wells and nine production wells show the accuracy and effectiveness of the proposed method.


Introduction
With the old oil fields entering the later period of development, the moisture content of a reservoir is increasing, and the oil production is reducing [1].How to update the technique to ensure oil recovery is an important measure to stabilize oil production.ASP flooding is an important tertiary oil recovery technique which is widely studied [2].It can enhance the oil production evidently by use of the interaction among alkali, surfactant, and polymer which can improve the physicochemical property of a reservoir.The mechanism model of ASP flooding is a complex distributed parameter system.It is hard to get the optimal injection strategy by a common method, because of the features of infinite dimensions, spatial-temporal coupling, and complex nonlinear behavior.In addition, there is a lack of a uniform mathematical model description because of the uncertainty of alkali reaction.It is important to build an accurate and easily to be applied model in the research of ASP flooding.
In an enhanced oil recovery problem of ASP flooding, the inputs are the injection concentrations of alkali, surfactant, and polymer; the output is the moisture content of production wells; the state variables contain water saturation, pressure, and grid concentrations; and the performance index is usually the net present value (NPV).In applications, many different kinds of methods are used to search the optimal injection concentration.One kind is to solve the mechanism equations [3]: fluid equations, seepage equations, and so on.However, this method usually involves a lot of math operations.The process of mathematic treatment always needs a huge amount of calculation.Though many scholars have done researches on modeling of ASP flooding, nearly all the works are about improving or enriching the primary model.The main work they have done is to consider more influencing factors or inductive research; the model becomes more complex and difficult.So the identification method is considered to approximate the ASP flooding system so that the mathematical operation process can be simplified.It will be helpful to improve the computational efficiency and reduce the computational complexity of the algorithm [4].
In the current research, the block-oriented nonlinear models have been widely used because of their simple structures and abilities to approximate nonlinear systems and considering the dynamic characteristics.This kind of method consists of two parts: the static nonlinear block and the dynamical linear block.The Hammerstein model and Wiener model are two common structures.The only difference between them is the sort order of two blocks.For the Wiener model, the dynamical linear block is in the front, while for the Hammerstein model, it is in the back [5].In [6], an optimal two-stage identification algorithm was presented for Hammerstein-Wiener systems where two static nonlinear elements surrounded a linear block.Reference [7] introduced one Wiener modeling method, in which the output was expanded on a set of spatial basis functions and temporal coefficients.After identifying the temporal coefficients, the final model is built.At present, most researches about the block-oriented models focus on a lumped parameter system mainly.Very few of them involve the DPS.Though the above methods can get a relatively good model, their essence is a lumped parameter method.In the case of model matching and modeling accuracy and the spatialtemporal characteristics, unmodeled dynamics need be further studied.
Spatial-temporal decomposition technology [8,9], which comes from Fourier series expansion, is a very useful modeling method for DPS.A spatial-temporal variable can be expanded into an infinite number of spatial basis functions and corresponding temporal coefficients.Generally speaking, the first few primary spatial basis functions can reflect the maximal inner information of the system, which provides a good approximation because of their separation properties [10].In this way, the spatial-temporal method can obtain a finite-dimensional model.In spite of these advances, a conventional spatial-temporal modeling method only considers the spatial information, while the temporal information which is important for analyzing system performance is ignored.The biorthogonal decomposition [11], which has wide applications in image processing field, can expand the distributed parameter on the spatial basis function and temporal basis function.It can reflect both space information and time information.
An important condition of modeling is to guarantee the accuracy which is highly dependent on the choice of spatial basis functions.In particular, Karhunen-Loeve (KL) decomposition, which is called proper orthogonal decomposition and principal component analysis [12,13], is a popular spatial-temporal decomposition approach to find principal spatial structures and reduce the dimension.Among all linear expansions, KL expansion is the most efficient and, for a given approximation error, the number of KL bases required is minimal.As a result, KL decomposition can help to reduce the model dimension and the number of estimated parameters.
Dynamic programming (DP) is a very good optimization procedure [14].The basic idea is Bellman's principle of optimality.In DP, the decision-making process is divided into many stages, and optimization computations are done sequentially from the last stage to the first stage.To overcome the shortcomings of requiring a very large number of grid points (the menace of expanding grid) to get reasonable accuracy, Luus [15,16] suggested the use of a small number of grid points, and to get sufficient accuracy to use DP in an iterative way, where after every iteration, the grid points are brought closer together.The details of this method called iterative dynamic programming (IDP) and a computer program written in FORTRAN for a typical optimal control problem are given by Luus [17].Here, IDP is adopted to solve the enhanced oil recovery of ASP flooding.
From the above, a new IDP algorithm based on biorthogonal spatial-temporal modeling is developed to solve the enhanced oil recovery of ASP flooding in this paper.At first, the mechanism model for the enhanced oil recovery of ASP flooding is given.Then the biorthogonal spatial-temporal modeling is presented to build the relation between the input (injection concentration of displacing agents) and states (water saturation).After giving the necessary condition for the existence of basis function, the spatial basis function and temporal basis function are obtained by the snapshot method based on the state samples.To identify the parameter in the Wiener model, the LSE is adopted.During this process, a theorem is proved to help the identification.In addition, the relation between states and output (moisture content) is established by ARMA which is identified by RLS.Furthermore, combining the established model with the index NPV, the approximation optimization problem for ASP flooding is solved by IDP.At last, the proposed method is applied to solve the enhanced oil recovery of ASP flooding with four injection wells and nine production wells.The simulation result shows the accuracy and generalization ability of the proposed method.
This article is organized as follows: The mechanism model description of ASP flooding is described in Section 2. Section 3 presents the biorthogonal spatial-temporal Wiener modeling method based on biorthogonal decomposition.In Section 4, the optimization procedure with IDP is given.In Section 5, the simulation process and result with our spatial-temporal modeling method and IDP are shown.Finally, a few conclusions are summarized in Section 6.

Mechanism Model Description for the Enhanced Oil Recovery of ASP Flooding
The following are the basic assumptions for ASP flooding: (1) The stratum satisfies heterogeneity, the whole reservoir is isothermal, and the adsorption 2 Complexity process complies with the Langmuir isothermal adsorption equation.
(2) The oil displacement system is composed of alkali, surfactant, and polymer.All displacing agents exist in the water phase, while only the surfactant exists in the oil phase.
(3) The Darcy law is fit for the flow of the oil phase and water phase.
(4) The phase equilibrium is set up instantly for all kinds of adsorptions which satisfy the generalized Fick's law.
For a reservoir with the region x, y, z ∈ Ω, Ω ∈ R 3 , ASP flooding consists of five main components: oil, water, polymer, surfactant, and alkali.Based on [18,19], the below mathematical model of ASP flooding can be drawn: The seepage continuity equation for the oil phase is The seepage continuity equation for the water phase is The adsorption diffusion equation for the polymer is The adsorption diffusion equation for the surfactant is The adsorption diffusion equation for the alkali is Furthermore, the dosage limit of displacing agents is defined as The slug injection strategy is usually used in ASP flooding, in which the slug only refers to the stage with displacing agents injected.Suppose that there are P slugs; then P = 3 denotes the three-slug injection which is usually adopted in industrial production.Figure 1 gives the details.
The injection concentration limit and slug size limit satisfy below conditions: where c Θ max denotes the maximum injection concentration of displacing agents and T i denotes the time length of each slug which has the integer restriction, T i = t i − t i−1 .Note that P + 1 th stage denotes the water flooding, in which c Θin ≡ 0.
The initial conditions and the terminal constraint of water content are

Complexity
The boundary conditions are selected as Among all the above equations, K is the absolute permeability of rock; k ro and k rw are the relative permeability of oil and water, respectively; p o and p w are the pressure of oil and water, respectively; S o and S w are the oil saturation and water saturation, respectively; c Θ , Θ = p, s, OH are the concentrations of polymer, surfactant, and alkali, respectively; c ij is the mass concentration of component i in solution j; μ o and μ w are the viscosity of oil and water, respectively; ϕ, ϕ p , and ϕ s are the rock porosity, reachable porosity of polymer, and reachable porosity of surfactant, respectively; ϕ p,s = ϕf a,s , f a , f s are the reachable porosity factors; R k is descending coefficient of relative permeability; K a is the speed factor of the ion exchange and adsorption capacity; K b is the adsorption constant of surfactant; v w is the seepage velocity; R OH is the alkali consumption; C rp and C rs are the adsorption quality of unit mass of rock of polymer and surfactant, respectively; q o and q w are the flow rate of oil and water in the standard state; q c , q d , and q e are the transport velocity of shaft flooding agents; D i , i ∈ w, o, OH is the diffusion coefficient; D ij , i ∈ w, o , j ∈ s, p is the diffusion coefficient of component j in solution i; M Θp denotes the maximum usage of displacing agents Θ; and ∇ is the Hamilton operator, which can be written as ∇ = ∂/ ∂x i + ∂/∂y j + ∂/∂z k in a rectangular coordinate system.
The flow coefficients of the oil phase and water phase are r o = Kk ro /B o μ o and r w = Kk rw /B w R k μ w , respectively.The moisture content f w x, y, z, t is the function of all spatial states p, S w , c p , c s , and c OH and can be defined as The concentration of surfactant is defined as c s = q w c ws + q o c os q w + q o 12 Define q in x, y, z, t ≥ 0 and q out x, y, z, t ≥ 0 as the injection term and production term, respectively, at the well location x, y, z .In view of all the above descriptions, q o , q w , q c , q in , q out , and the control variables c Θin only work at the location of injection wells and production wells.
Suppose that the location sets of injection wells and production wells are ψ w = x, y, z j,k , j ∈ κ w , k ∈ ϑ wj and ψ p = x, y, z j,k , j ∈ κ p , k ∈ ϑ p .The definitions of flow terms (q o , q w , q c , q in , and q out ) in (1), ( 2), (3), (4), and (5) can be described as x, y, z ∉ ψ p ⋃ ψ w , x, y, z ∉ ψ p ⋃ ψ w , In the optimization model of enhanced oil recovery for ASP flooding, the performance index is functional.In this paper, the maximum of NPV is adopted as the performance index, which means the profit maximization [20].

14
where t f denotes the terminal flooding time, χ denotes the discount rate, t a denotes the total discounting period, P oil denotes the price of oil, P Θ denotes the price of three displacing agents, and P cost is the production cost per day.
The other physicochemical algebraic equations are presented in Appendix.

Biorthogonal Spatial-Temporal Wiener Modeling Method
The ASP flooding system is a complex distributed parameter system which is very difficult to be modeled with a general method.In this paper, a biorthogonal spatial-temporal Wiener modeling method based on Karhunen-Loeve decomposition is presented, which is shown in Figure 2. The spatial information and temporal information of samples are reflected on the spatial basis function and temporal basis function.During the modeling process for ASP flooding, the main problem is to identify an appropriate spatialtemporal model according to the input (injection concentration of displacing agents c Θin , Θ = OH, s, p ), output (moisture content f w ), and system state (pressure P, grid concentration c Θ , and water saturation S w ).For easy description, a set υ is introduced to represent all state variables.The whole model for ASP flooding can be divided into three stages: (1) Determine the spatial basis function and temporal basis function according to system state υ x, y, z, t .
(3) Model for the moisture content f w t of production wells and corresponding system states υ x, y, z, t .
3.1.Spatial-Temporal Wiener System.The structure of traditional Wiener model is presented in [21].According to the specialty of ASP flooding, the entire system is transformed into the DPS and the output identification is added.The specific structure is shown in Figure 3.It mainly consists of two parts: the distributed parameter identification for states and the lumped parameter identification for moisture content.For simplicity, only a one-dimensional method is shown in this paper.The expression of the three-dimensional method is similarly to the one-dimensional situation.
The distributed parameter identification contains a linear dynamical block followed by a static nonlinear block N • : ℝ⟶ ℝ, the description is as follows.
where h x, z, q 1 × 1 is the transfer function of dynamical block, τ and t are the time variables, x and z denote spatial variables defined on the domain Ω, while q is the time forward shift operator, v x, t denotes the system intermediate variable, c Θin z, t ∈ ℝ and υ x, t ∈ ℝ stand for the input and state output of the system at time t separately, and d x, t ∈ ℝ includes the unmodeled dynamics and the stochastic disturbance.In addition, the integral operator is used for spatial operation and sum operator for temporal operation [22,23].Here only the single-input single-output system is considered for simplicity.It is obvious that this method is also applicable to the multi-input multi-output system [24].This part can realize the modeling between inputs and system states.
The output identification, which is also the lumped parameter identification, contains a static block F • : ℝ⟶ ℝ.Its aim is to build the relation between system states and target outputs.The formula description is shown in (16).
where f w t ∈ ℝ denotes the target output of the system and d ′ x, t ∈ ℝ stands for the unmodeled dynamics and the stochastic disturbance.
From the above, the whole model for ASP flooding is given.Then the problem is to identify N and h according to the input-state data and F on the basis of state-output data.

Basis Function Determination Method
3.2.1.Biorthogonal Spatial-Temporal Decomposition.KL decomposition can acquire a lower dimension model than other spatial-temporal methods which are expanded by linear basis functions [4].That is to say, these basis functions obtained by KL decomposition can reflect the system information at the most extent when the number of basis function is constant.Here the KL decomposition is adopted to find the system dominant characteristics from sampled data.For a general system, suppose that the system output y x i , t N,L i=1,t=1 is uniformly distributed on the time and space domain, the orthogonal spatial basis function is φ i x ∞ i=1 , the temporal basis function is ψ i t ∞ i=1 , and the coefficient is α i ∞ i=1 .Then the output can be expressed as follows [25][26][27]: The basis function is orthogonal; that is, where •, • denotes the inner product.
In practical application, only a few dominant bases which can reflect most system information are chosen.Then

Complexity
To obtain the dominant bases, we can get the below minimization problem: Bring in a Lagrangian multiplier and convert (20) into an unconstrained optimization problem: Take the variation of (21-a) with respect to φ i x and Equation ( 22) can be rewritten as The necessary condition of the extreme value for this functional problem is δJ = 0. Since φ i x and ψ i t can be an arbitrary function, the below equation can be developed: For all sample points in the time domain and space domain, T f x, t dt = ∑ L t=1 f x, t and Ω f x, t dx = ∑ N z=1 f x, t can be drawn approximately.Here, z = 1, 2, … , N correspond to x 1 , x 2 , … , x N .Then (24) can be written as In (25-a), α i ψ i t can be regarded as the projection of y x, t on φ i x ; that is, α i ψ i t = y x, t , φ i x .Thus, Considering the mutual independence for i = 1, … , n, (27) is true if and only if the below condition is satisfied: f w (t) F(g) c in (z,t) Figure 3: Spatial-temporal Wiener system.
From the above, the necessary condition for the optimization problem in (20) Therefore, the problem of solving ( 20) is transformed into the problem to find the solution of necessary condition (30).The snapshot method is introduced to obtain the basis function.Since the sample data always distributes discretely, the numerical solution method is usually adopted to solve (30).An eigenvalue problem of the matrix can be obtained by discretizing (30).

Snapshot Method.
The snapshot method [28] can get the solution of (30) effectively.Assume that the spatial basis function φ i x and temporal basis function ψ i t are the linear combination of a series of snapshots: Define the correlation function of two points as follows: Then the optimization problem in ( 30) is converted into the below matrix eigenvalue problem: where According to (34), a series of γ i and ϖ i can be obtained.Then substitute them into (31); the spatial basis function and temporal basis function can be determined.Since matrix C is a symmetric positive semidefinite, the basis functions corresponding to γ i and ϖ i are orthogonal.After standardization treatment, we can get the final spatial basis function φ i x and temporal basis function ψ i t .

Dimension Selection.
In the above problem, suppose that the maximal number of the nonzero eigenvalue is K, which satisfies K ≤ min N, L .Sort these eigenvalues in a descending order, then For the determination of basis function, the more number the dominant basis function is, the better accuracy the system has.But this will lead to the complexity of modeling.So only a few domain bases which can describe most energy are selected here.Define that the total system energy is equal to the sum of all eigenvalues.Then the bigger the value is, the more energy the basis function reflects.For every basis function ϕ i • with eigenvalue μ i , which is obtained by biorthogonal spatial-temporal decomposition, the proportion of energy is Generally speaking, when the proportion is more than 99%, it can be concluded that these bases can reflect the most energy; the corresponding n is the dimension finally.
Experience shows that for most spatial-temporal systems, the front several basis functions can represent nearly all the energy.For an arbitrary basis function θ i • n i=1 , the below equation can be obtained [29]: 2  36 This shows that the KL method can get a lower order than other methods in modeling.

Biorthogonal Spatial-Temporal Model.
For ASP flooding, the spatial basis function φ i x n i=1 and temporal basis function ψ i t l i=1 can be obtained by decomposing the state sample information with the method in Section 3.2.And the input basis function ϕ i x m i=1 can be obtained by decomposing the injection concentration of displacing agents c Θin 7 Complexity combined with ψ i t l i=1 .Then the expanded biorthogonal spatial-temporal model is presented as follows: The input c Θin can be expressed as where c i = c Θin z, t , ϕ i z ψ i t denotes the coefficient.
The unmodeled dynamics and error d x, t can be expressed as where d i t = d x, t , φ i x ψ i t denotes the coefficient.Suppose that the block h x, z, τ in ( 15) is absolutely integrable on time domain 0, ∞ at any spatial point x and z.That is to say, the dynamical block is stable.Then the mathematical description of h x, z, τ is where α i,j,k ∈ ℝ i = 1, … , n, j = 1, … , m, k = 1, … , l denotes the constant coefficients corresponding to the basis functions φ i x , ϕ j z , and ψ k τ , respectively.Since the unmodeled dynamics and error are added between the two blocks of the Wiener model, we cannot get the relation between the input and the state directly.The method of solving the inverse function is usually used.To illustrate the mathematical description more precisely, a theorem is presented in this section.
Lemma 1 (see [30]).Assume that y = f u , u ∈ D , and If both y = f u and u = φ x exist inverse functions and they are u = f −1 y and x = φ −1 u , respectively, the composite function y = f φ x must exist an inverse function; the function is x = φ −1 f −1 y .Theorem 1.For a biorthogonal spatial-temporal Wiener system in Figure 3, suppose that ξ i x is the spatial basis function, ς i t is the temporal basis function, ρ i is the coefficients, and P • : ℝ⟶ ℝ and P −1 • : ℝ⟶ ℝ are the inverse functions of each other.If υ x, t = N v x, t = ∑ n i=1 ρ i ξ i x ς i t P v x, t , for an arbitrary given set of spatial basis function ϖ i x and temporal basis function w i t , there must be a set of coefficients μ i which can make v x, t = N −1 υ x, t = ∑ n i=1 μ i ϖ i x w i t P −1 υ x, t true.
Proof.According to the property of the inverse function that every implicit function has an inverse function, there must be an inverse function for N • .
As υ x, t = ∑ n i=1 ξ i x ρ i t P v x, t , P • and P −1 • are the inverse functions of each other, according to Lemma 1: Take the derivative of (40) with respect to y x, t ; then where u = 1/ ∑ n i=1 ρ i ξ i x ς i t υ x, t .For v x, t = N −1 υ x, t , take the derivative with respect to υ x, t ; then where u = υ x, t .
To make sure that v x, t = N −1 υ x, t and υ x, t = N v x, t are the inverse function of each other, according to the property of derivatives of inverse functions, (41) has to be equal to (42).Then μ i ϖ i x w i t satisfies the below relation: From our experience, there are infinite pairs of w i t , ϖ i x and μ i .Let χ = 1/n∑ n i=1 ρ i ξ i x ς i t ; then If μ i ϖ i x w i t ≡ 1/n 2 χ, then for an arbitrary given set of w i t and ϖ i x , there is the corresponding set of μ i = 1/ n 2 χw i t ϖ i x , which can make (41) true.Therefore, Theorem 1 is proved.
In practical application, it is no necessary that (43) is true precisely.The appropriate values of w i t , ϖ i x , and μ i can be selected approximately under a given error.Among all the combinations, one group of w i t , ϖ i x , and μ i by the KL theory can be found.
Assume that N • is invertible; then v x, t can be expressed as (45) according to Theorem 1.
where B i is the identification parameter.

Complexity
In addition, v x, t = h x, z, q c Θin z, t + d x, t .That is, Uniting ( 38), (45), and (46), then Project (48) on the temporal basis function ψ s t and spatial basis function φ h x with the Galerkin method.Considering the orthogonal property of temporal basis function, (48) can be inferred.
Since the spatial basis function is orthogonal, n equations can be drawn from (49).
Without the loss of generality, take α 1,1 = I.Equation (51) can be rewritten as where The only unknown parameter in (52) is A which can be obtained by a quadratic criterion on the prediction errors [31] as follows: Parameter A can be identified by using the LSE method [32]: So far the value of A has been determined, that is, α and B. Parameter c i is decomposed from the input variable; it is the optimization variable actually.Thus, all the unknown parameters have been determined.According to this method, the relation between the input (c Θin ) and the state (υ) can be built.

Output Identification.
To acquire the mathematical description between system state υ and output f w t , the multivariable ARMA identification [31] is introduced which is identified by the recursive least squares (RLS) [33].Suppose that there are n s states in ASP flooding.The ARMA model is , and U n s q −1 = ∑ N n s i=1 U n s ,i q −i ; q −1 is the backward shift operator which means backing up one sampling period; U 0 , U 1 , … , U n s are the matrix polynomial of q −1 ; N 0 denote the time lag for output; and N 1 , … , N n s denote the time lag for all states.

Complexity
Convert the above problem into the linear regression model, where Identify parameter θ with RLS in (58).

58
where K t denotes the weight matrix, 0 < μ < 1 is the forgetting factor, and P t denotes the positive definite covariance matrix.
Then the modeling for output identification is finished.In conclusion, the whole model for ASP flooding is determined, which is composed of ( 15) and (55).Adding the net present value as the performance index, an optimization problem for ASP flooding can be constructed.Then the iterative dynamic programming (IDP) algorithm is applied to solve the optimal injection strategy.

Optimization Procedures
The biorthogonal spatial-temporal Wiener modeling method in Section 3 is also applicable to a three-dimensional space.For a practical reservoir given in Section 2, which is a three-dimensional system, since the moisture content f w of production wells is related to the water saturation S w of corresponding wells directly, only the relation between S w and f w for the same production well is considered here.Then build the model for ASP flooding with proposed method.Equations ( 15) and ( 55) can be transformed into (59-b) and (59-c).According to the discussions above, the below optimization problem can be drawn:

c
where x, y, z is the location of the production well, N 0 , N 1 denote the time lag of f w , Ŝw , and ξ, ζ, ς denotes the coordinate of an arbitrary point in the three-dimensional space.c Θin , Θ = OH, s, p is the injection concentration for three displacing agents, which is also the optimization variable in (59-a)-( 59-c).
The next work is to get the optimal injection strategy for the enhanced oil recovery of ASP flooding based on the approximation optimization problem.
In this section, the IDP [34] is adopted to solve this problem based on the biorthogonal spatial-temporal Wiener model.For simplicity, only the injection concentration is selected as the control variable.The specific solution procedure for ASP flooding is presented as follows: (1) Divide the whole process of ASP flooding into P + 1 stages.The optimization variable is c Θin k , k = 1, … , P. Note that c Θin P + 1 ≡ 0.
(2) In initialization, give the slug length and select the initial control feasible region r Θin and the initial control c 0 Θin .The number of the initial discrete state grid is M. Select R values for every grid.Set contraction factor γ, maximal iteration l max , and convergence accuracy ε.Start from generation l = 1.
(3) Set the current control domain r l Θ = r Θin .(4) In the current feasible region, generate M − 1 injection concentrations with a uniform strategy [35] in the region is the optimal control strategy obtained from the previous generation.Specially, P denotes the optimal injection concentration of the previous iteration.When generating injection concentrations which not conform to (8), we deal with the control variables as follows: Calculate the system state from t P−1 to t P in (8).Then select the optimal injection concentration c * l Θin P + 1 which is closest to the value of current t P and calculate the system state from t P to t P+1 .Compute the total performance index from t P−1 to t P+1 .Compare and record the injection concentration c * l Θin P with the maximal NPV for every grid at stage P.
Take the optimal injection concentration obtained from step (7) as the center of the feasible domain in the next iteration.
(9) Add the iteration time l = l + 1.If J new − J old ≥ ε, go to step (4), and continue computing and iterating.If J new − J old < ε, then save the values of optimal injection concentration, and finish the optimization process.
As for the physicochemical equation constraints, the penalty function method [36] is used to deal with these constraints.The other process keeps the same.

Solutions of Enhanced Oil Recovery for ASP
Flooding Based on Spatial-Temporal Modeling and IDP 5.1.Reservoir Description.Suppose that the reservoir of ASP flooding consists of four injection wells and nine production wells.All wells distribute uniformly; there is one injection well at the center of every four production wells.The distribution of wells is shown in Figure 4.

Reservoir Parameters.
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 the upper surface is 2420 m, the porosity of every layer is 0.3, and the pore volume is 1.1097 × 10 6 m 3 .The initial grid concentration of ASP is 0 g/L.The initial permeability, initial pressure, and initial water saturation are shown in Figures 5-7.The grids of the reservoir are divided into three directions x, y, and z.The grids in x and y are divided into 21 units, and the grids in z are divided into 7 units.The total number is 21 × 21 × 7.
The water injection rate of each injection well is 83 m 3 /d.The recovery rate of production wells is defined in Table 1.Using reservoir simulation software CMG 2010, the whole oil production process can be obtained by simulation based on the mechanism model of ASP flooding in Section 2. The whole production time lasts for 96 months, and the sampled output can be obtained by the snapshot method.The injecting time of ASP flooding is served as the initial time, so the total number of time nodes is L = 97 and the total number of space nodes is N = 21 × 21 × 7. The other reservoir parameters for ASP flooding keep the same with [37].

Modeling and
Verification for ASP Flooding.Given the ASP flooding above, we build the approximate model with the method proposed in this paper.In order to motivate the system sufficiently, simulations are executed 50 times on    12 Complexity CMG 2010.The software randomly generates different injection strategies for every month, and the process of injection totally lasts for 48 months.After obtaining the samples, the modeling process with the proposed method is carried out on MATLAB R2016b.In order to improve the generalization of the identification model, the below basis function selection method is used here.(a) Apply the biorthogonal spatial-temporal decomposition method to get the spatial basis function, the temporal basis function, and coefficients of the 50 sets of grid water saturation samples.(b) Use every group of basis function to reconstitute the system states with other 49 groups of coefficients, and calculate the average mean square errors of all sampling points for these 49 times.(c) After executing this process for 50 times, select the group of basis function with the smallest mean square error to serve as the final basis functions for modeling.Then build the spatial-temporal model for ASP flooding with the proposed method.The specific description is given in (59-a)-( 59-c).
Regarding the determination of the number of the spatial basis functions, the error index is defined as follows: where e x, y, z, t = Ŝw x, y, z, t − S w x, y, z, t .Generally speaking, the more number the basis function is, the more information of the system model reflects, and the higher accuracy the model has.But at the same time, the model will be sensitive to the external disturbance and the generalization ability will decrease.Besides, the dimension of the model increases.On the other hand, the less quantity the basis function has, the model will reflect less information and the error may be very big.It is important to choose the proper amounts of the spatial basis function and temporal basis function.In order to test the influence on the accuracy of modeling, we choose the different numbers of basis function to reconstruct the system with 13 Complexity the coefficients and calculate the average RMSE which is shown in Table 2.
Note that the number of spatial basis function is equal to that of temporal basis function.This is decided by the matrix operations.According to Table 2, the accuracy of the model increases with the number of basis function.However, when the number is over 4, the RMSE of the model increases slightly.So the number of both spatial basis function and temporal basis function adopted in this paper is 4.Besides that, the time lag of moisture content is N 0 = 5; the time lag of water saturation is N 1 = 6.
With regard to the whole modeling process for ASP flooding, the errors are as follows: The mean RMSE for grid water saturation is 0.0285, and the mean RMSE for the moisture content of production wells is 1.1378%.This demonstrates the good modeling ability.
In order to verify the generalization ability, the injection concentrations of ASP flooding are given randomly: The displacing agent injection totally lasts for 48 months which are divided into 3 slugs uniformly, and the rest time is water flooding.The specific injection can be found in Figure 8(a).Because the whole system is five-dimensional considering the time and value of water saturation, it cannot be plotted in the picture directly.We choose the 21 grids ranging from 1, 11, 3 to 21, 11, 3 at the third layer; then the result and error of modeling are shown Figures 8 and 9.
To analyze the generalization ability further, we compute the error with (63).The mean absolute errors for the moisture content of production wells and for the grid water saturation are 1.2970% and 0.0318%, respectively.Figure 10 compares the moisture content results.The dotted line denotes the moisture content of numerical reservoir simulation software CMG 2010.It represents the real result, which is called the original curve in Figure 10.And the continuous line is the moisture content of the proposed modeling method.It can be known that the errors of nine production wells are very small, which verifies that the whole model has better generalization ability.
From the above, the biorthogonal spatial-temporal Wiener modeling method can model for ASP flooding well.It has good modeling accuracy and generalization ability.

Solutions of Enhanced
Oil Recovery for ASP Flooding with IDP.For now, the specific spatial-temporal identification model for ASP flooding has been acquired.Then apply the IDP algorithm according to Section 4 to optimize the optimal injection strategy for alkali, surfactant, and polymer.5.3.1.Parameter Setting.P = 3, P p = 6 85$/kg, P s = 8 31$/kg, P a = 2 86$/kg, the production cost per day is P cost = 500 $/day, the discount rate χ = 1 5 × 10 −3 , and the oil price is referenced to the WTI price [18], that is, 55$/barrel, which   In order to illustrate the accuracy of the proposed method, the mechanism model is also solved with the IDP algorithm in Section 4 as a comparison.Before computation, the mechanism model is processed by the finite difference method.The results are presented in Figures 11-15.For simplicity, we use "PM" to denote the proposed method in this paper and "MM" to denote the optimization based on the mechanism model.After finishing the optimization, the optimal injection strategy of the proposed method is c OHin = 3 5, 2 59, 1 8 kg/m 3 , c pin = 2 26, 1 17, 0 69 kg/m 3 , c sin = 2 42, 1 5, 0 49 kg/m 3

65
In Figure 11, the green line and black line denote the control strategies of MM and PM corresponding to the mechanism model and identification model, respectively.It is clear that the result of PM is very close to that of MM, which can explain the effectiveness of PM partly.
The moisture content and oil production are shown in Figures 12 and 13.The full line denotes the moisture of nine production wells, and the dotted line denotes the oil production of nine production wells.In order to show the oil flooding effect clearly, the average moisture content and average oil production for nine production wells are computed and showed in Figures 14 and 15.From the comparative analysis, the moisture content of the proposed method is very close to that of the mechanism model for both moisture content and oil production.For the result of the proposed method, the moisture content descends obviously at the early injection stage.This is because the more production is gained at the early stage, the more NPV is discounted to the present time.While the injection of displacing agents can improve the oil production obviously, the injection 15 Complexity concentration at the early stage has a larger value.But after a period of development time, since the oil reserves are decreasing and oil flooding efficiency is attenuating, the moisture content is rising gradually.
To illustrate the optimization effect, some essential indexes are calculated by the MM and proposed method.The specific result is shown in Table 3.
By comparison, the results of the PM and MM are very close for nearly all the indexes.Although the production and NPV have little differences, these can be ignored compared with their big cardinal number, since only a tiny difference on injection strategy can accumulate a certain influence on NPV.After being substituted by the identification model in this paper, the oil production is decreased by 311t, and the NPV is decreased by 0.08 million $ than the optimization result of the mechanism model.The relative error of oil production is only 0.22%, and that of NPV is only 0.11%.
In a word, the proposed method is effective for the modeling and optimization.It has good accuracy and generalization and can be used to solve the enhanced oil recovery problem of ASP flooding in oil exploitation.

Conclusion
In this paper, an iterative dynamic programming (IDP) method based on biorthogonal spatial-temporal Wiener modeling is developed to solve the enhanced oil recovery for alkali-surfactant-polymer (ASP) flooding.The main contributions are the following: (1) a biorthogonal spatialtemporal Wiener modeling method is developed; (2) in order to identify the static nonlinear part, a theorem is given and proved; (3) combined with IDP, the mechanism model is replaced by identification which is built by the proposed method in solving the optimization problem of ASP flooding.The simulation result shows the accuracy and generalization ability of the proposed method.This paper presents an effective method to solve the enhanced oil recovery of ASP   16 Complexity flooding which can avoid solving the complex mechanism equations.It is also suitable for solving other complex optimization problems in engineering.

Appendix
The relevant physicochemical algebraic equations are shown below.

A. Alkali Consumption
Many kinds of conditions can influence the alkali consumption, while only the adsorption consumption, acidoid consumption, and heavy metal ion consumption are considered in this paper to describe the mechanism well and simply.The alkali consumption is given through the chemical reaction term r i , which is defined as follows: where r i is the alkali consumption per unit volume of different effects; the specific forms are as follows [20]: (a) The rapid alkali consumption r 1 is caused by the ion exchange between Na + in solution and H + on the rock surface.This effect can be expressed as the Langmuir isothermal adsorption equation: where r 1 is the maximal alkali consumption of this effect and a 1 is the coefficient; they are set by experimental data.
(b) The alkali consumption r 2 caused by acidoid is The curve of r 2 under different acidity and alkalinity can be obtained through an experiment.
(c) The alkali consumption r 3 caused by Ca 2+ , Mg 2+ , and acidoid in oil is where K Ca 2+ and K Mg 2+ denote the solubility product.

B. Surface Tension
Surface tension varies with the concentration of alkali, surfactant, and polymer.The form is [20] σ = σ c OH , c p , c s B 1

C. Capillary Pressure
The capillary pressure of the compound system can be described as the function of oil/water capillary pressure and surface tension; that is [20], where C PC and N PC denote the constant and S n is the wet phase saturation.where C 0 rs and a s are related to the ion concentration, pH max is the pH value of the injection fluid, and b s is the coefficient [37].

E. Polymer Adsorption
Considering that the relation of polymer adsorption is reversible with salinity and is irreversible with concentration, the adsorption can be described as follows according to the Langmuir isothermal adsorption law [36]: where C max rp is the maximal adsorption capacity to the rock surface of the polymer under different salinity and a 1 and b 1 are the adsorption equilibrium constants.

F. Relative Permeability
The relative permeability of oil phase k ro and water phase k rw is usually given by the interpolation method in the process of model solution.Here, the formula proposed by Ge et al. in [38] is adopted.
where A, B, C, and D denote the coefficients which are identified by the data of k ro and k rw .Since the addition of the displacing agent, the permeability reduction factor has to be considered.
where q θ and R k denote the adsorption capacity of displacing agent and permeability reduction factor of the water phase under different salinity, respectively.17 Complexity

G. Liquid Viscosity
For the displacing liquid of ASP flooding, the relation between viscosity μ and shear rate γ can be expressed by the Meter equation as [38] μ ASP = μ w + μ 0 p − μ w 1 + γ/γ 1/2 P a −1 , G 1 where γ 1/2 is the shear rate when viscosity μ is equal to the mean value of μ 0 p and μ w and P a is the economic coefficient.

H. Residual Saturation
The residual saturation is relevant to the capillary number N c ; it can be obtained through the experiments [37].
where S rj denotes the residual saturation of phase j.

Figure 1 :
Figure 1: Three-slug injection scheme.The slug size is the injection concentration, and the slug length is the injection time of displacing agents.

10 Complexity l = 1 ( 5 )( 6 )
corresponds to the initial condition.Calculate the system as (59-b), and get M state trajectories.M state values S w k − 1 , k = 1, … , P + 1 at the beginning of each stage constitute M state grids of the corresponding stage.Start from P + 1 th stage, and calculate the state and performance index as (59-a)-(59-c).Record the result.Start from the Pth stage, and generate R injection concentration for every state grid S w k − 2 as (60).a diagonal matrix in which all elements on diagonal are random numbers belonging to −1, 1 and c * l−1 θin

( 7 )( 8 )
Shift the time step forward, and repeat step (6) until the first slug.Compute the optimal index of the whole time domain as (59-a), and save the corresponding injection concentration c * l Θin k with the maximal NPV.Shrink the feasible domain of decision variables:

Figure 4 :Figure 7 :
Figure 4: The distribution diagram of well position.

Table 1 :Figure 6 :
Figure 6: The distribution diagram of initial pressure.The unit is MPa.

Figure 5 :
Figure 5: The distribution diagram of initial permeability.

Figure 9 :
Figure 9: Modeling error of water saturation.

Figure 10 :Figure 11 :
Figure 10: Comparison for the moisture content of nine production wells.

Figure 12 :Figure 13 :
Figure 12: The moisture content and oil production of production wells for PM.

Figure 15 :
Figure 15: The average oil production comparison of 9 production wells.

Figure 14 :
Figure 14: The average moisture content comparison of production wells.

Table 2 :
RMSE for different numbers of basis function.