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

Authors: Shurong Li; Yulei Ge; Yuhuan Shi
Journal: Complexity
 (2018)

---

## Abstract

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.


---

## Body

## 1. 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 spatial-temporal 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 Section2. 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.


## 2. 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 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 regionx,y,z∈Ω,Ω∈R3, 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(1)∇⋅KkroBoμo∇po−ρogh+qo=∂∂tϕ1−SwBo.The seepage continuity equation for the water phase is(2)∇⋅KkrwBwRkμw∇pw−ρwgh+qw=∂∂tϕSwBw.The adsorption diffusion equation for the polymer is(3)∇⋅KkrwcpBwRkμw∇pw−ρwgh+∇⋅Dw+DwpϕpSwBw∇cp+qc=∂∂tρr1−ϕCrp+ϕpSwcpBw.The adsorption diffusion equation for the surfactant is(4)∇⋅KkrwcwsBwRkμw∇pw−ρwgh+∇⋅KkrocosBoRkμo∇po−ρogh+∇⋅Do+DosϕsSoBo∇cos+∇⋅Dw+DwsϕsSwBw∇cws+qd=∂∂tϕsSocosBo+ϕsSwcwsBw+ρr1−ϕCrs.The adsorption diffusion equation for the alkali is(5)∇⋅DOHϕSwBw∇cOH+∇⋅KkrwcOHBwRkμw∇pw−ρwgh+ROH+qe=∂∂tKa1+KbcOH21−ϕSwcOH+ϕSwcOHBw.Furthermore, the dosage limit of displacing agents is defined as(6)∑i=1P∭ΩqincΘindσ⋅Ti≤MΘmax.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 areP slugs; then
(7)cΘint=cΘini,ti−1≤t≤ti,i=1,2,…,P,0,ti−1≤t≤tf,i=P+1.P
=
3 denotes the three-slug injection which is usually adopted in industrial production. Figure 1 gives the details.Figure 1
Three-slug injection scheme. The slug size is the injection concentration, and the slug length is the injection time of displacing agents.The injection concentration limit and slug size limit satisfy below conditions:(8)0≤cΘin≤cΘmax,∑i=1PTi=tP,where cΘmax denotes the maximum injection concentration of displacing agents and Ti denotes the time length of each slug which has the integer restriction, Ti=ti−ti−1. Note that P+1th stage denotes the water flooding, in which cΘin≡0.The initial conditions and the terminal constraint of water content are(9)px,y,z,t∣t=0=p0,Swx,y,z,t∣t=0=Sw0,cΘx,y,z,t∣t=0=cΘ0,x,y,z∈Ω.The boundary conditions are selected as(10)∂p∂n∣∂Ω=0,∂Sw∂n∣∂Ω=0,∂cΘ∂n∣∂Ω=0.Among all the above equations,K is the absolute permeability of rock; kro and krw are the relative permeability of oil and water, respectively; po and pw are the pressure of oil and water, respectively; So and Sw are the oil saturation and water saturation, respectively; cΘ,Θ=p,s,OH are the concentrations of polymer, surfactant, and alkali, respectively; cij 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=ϕfa,s,fa,fs are the reachable porosity factors; Rk is descending coefficient of relative permeability; Ka is the speed factor of the ion exchange and adsorption capacity; Kb is the adsorption constant of surfactant; vw is the seepage velocity; ROH is the alkali consumption; Crp and Crs are the adsorption quality of unit mass of rock of polymer and surfactant, respectively; qo and qw are the flow rate of oil and water in the standard state; qc, qd, and qe are the transport velocity of shaft flooding agents; Di,i∈w,o,OH is the diffusion coefficient; Dij,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 ∇=∂/∂xi+∂/∂yj+∂/∂zk in a rectangular coordinate system.The flow coefficients of the oil phase and water phase arero=Kkro/Boμo and rw=Kkrw/BwRkμw, respectively. The moisture content fwx,y,z,t is the function of all spatial states p, Sw, cp, cs, and cOH and can be defined as
(11)fw=rwrw+ro.The concentration of surfactant is defined as(12)cs=qwcws+qocosqw+qo.Defineqinx,y,z,t≥0 and qoutx,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, qo, qw, qc, qin, qout, 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,zj,k,j∈κw,k∈ϑwj and ψp=x,y,zj,k,j∈κp,k∈ϑp. The definitions of flow terms (qo,qw,qc,qin,andqout) in (1), (2), (3), (4), and (5) can be described as
(13)qo=−1−fwqout,x,y,z∈ψp,0,x,y,z∉ψp,qw=−fwqout,x,y,z∈ψp,qin,x,y,z∈ψw,0,x,y,z∉ψp⋃ψw,qc=qwcp,x,y,z∈ψp,qwcpin,x,y,z∈ψw,0,x,y,z∉ψp⋃ψw,qd=qwcws+qocos,x,y,z∈ψp,qwcsin,x,y,z∈ψw,0,x,y,z∉ψp⋃ψw,qe=qwcOH,x,y,z∈ψp,qwcOHin,x,y,z∈ψw,0,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)maxJNPV=∫0tf1+χ−t/ta∭ΩPoil1−fwqoutt−∑ΘPΘqincΘintdσ−Pcostdt,where tf denotes the terminal flooding time, χ denotes the discount rate, ta denotes the total discounting period, Poil denotes the price of oil, PΘ denotes the price of three displacing agents, and Pcost is the production cost per day.The other physicochemical algebraic equations are presented in Appendix.


## 3. 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 Figure2. 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 spatial-temporal model according to the input (injection concentration of displacing agents cΘin,Θ=OH,s,p), output (moisture content fw), and system state (pressure P, grid concentration cΘ, and water saturation Sw). 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.(2)
Identify the biorthogonal spatial-temporal Wiener system.(3)
Model for the moisture contentfwt of production wells and corresponding system states υx,y,z,t.Figure 2
The schematic diagram for biorthogonal spatial-temporal Wiener modeling.


### 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.Figure 3
Spatial-temporal Wiener system.The distributed parameter identification contains a linear dynamical block followed by a static nonlinear blockN·:ℝ⟶ℝ, the description is as follows.
(15)υx,t=Nvx,t=N∑τ=0t∫Ωhx,z,τcΘinz,t−τdz+dx,t,where hx,z,q1×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, vx,t denotes the system intermediate variable, cΘinz,t∈ℝ and υx,t∈ℝ stand for the input and state output of the system at time t separately, and dx,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 blockF·:ℝ⟶ℝ. Its aim is to build the relation between system states and target outputs. The formula description is shown in (16).
(16)f⌢wt=Fυx,t+d′x,t,where f⌢wt∈ℝ 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 identifyN and h according to the input-state data and F on the basis of state-output data.


### 3.2. 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 outputyxi,ti=1,t=1N,L is uniformly distributed on the time and space domain, the orthogonal spatial basis function is φixi=1∞, the temporal basis function is ψiti=1∞, and the coefficient is αii=1∞. Then the output can be expressed as follows [25–27]:
(17)yx,t=∑i=1∞αiφixψit.The basis function is orthogonal; that is,(18)φix,φjx=∫Ωφixφjxdx=0,i≠j,1,i=j,ψit,ψjt=∫Tψitψjtdt=0,i≠j,1,i=j,where ·,· denotes the inner product.In practical application, only a few dominant bases which can reflect most system information are chosen. Then(19)ynx,t=∑i=1nαiφixψit,αi=φixψit,yx,t.To obtain the dominant bases, we can get the below minimization problem:(20)minφix,ψit1N1L∫T∫Ωyx,t−ynx,t2dxdt,s.t.φix,φix=1,φix∈L2Ω,i=1,…,n,ψit,ψit=1,ψit∈L2T,i=1,…,n.Bring in a Lagrangian multiplier and convert (20) into an unconstrained optimization problem:
(21-a)J=1N1L∫T∫Ωyx,t−∑i=1nαiφixψit2dxdt+∑i=1nλiφix,φix−1+∑i=1nηiψit,ψit−1=1N1L∫T∫Ωy2x,t−2yx,t∑i=1nαiφixψit+ΛΛTdxdt+∑i=1nλi∫Ωφi2xdx−1+∑i=1nηi∫Tψi2tdt−1,(21-b)Λ=α1α2⋮αnTφ1x0⋯00φ2x0⋮⋮0⋱00⋯0φnxψ1tψ2t⋮ψnt,(21-c)ΛΛT=∑i=1nαi2.Take the variation of (21-a) with respect to φix and ψit; then
(22)δJ=1N1L∫T∫Ω−2yx,t∑i=1nαiψitδφix+−2yx,t∑i=1nαiφixδψitdxdt+∑i=1n2λi∫Ωφixδφixdx+∑i=1n2ηi∫Tψitδψitdt.Equation (22) can be rewritten as
(23)δJ=∫Ω1N1L∫T−2yx,t∑i=1nαiψitdt+∑i=1n2λiφixδφixdx+∫T1N1L∫Ω−2yx,t∑i=1nαiφixdx+∑i=1n2ηiψitδψitdt.The necessary condition of the extreme value for this functional problem isδJ=0. Since φix and ψit can be an arbitrary function, the below equation can be developed:
(24)1N1L∫Tyx,t∑i=1nαiψitdt=∑i=1nλiφix,1N1L∫Ωyx,t∑i=1nαiφixdx=∑i=1nηiψit.For all sample points in the time domain and space domain,∫Tfx,tdt=∑t=1Lfx,t and ∫Ωfx,tdx=∑z=1Nfx,t can be drawn approximately. Here, z=1,2,…,N correspond to x1,x2,…,xN. Then (24) can be written as
(25-a)1N1L∑t=1Lyx,t∑i=1nαiψit=∑i=1nλiφix,(25-b)1N1L∑z=1Nyx,t∑i=1nαiφix=∑i=1nηiψit.In (25-a), αiψit can be regarded as the projection of yx,t on φix; that is, αiψit=yx,t,φix. Thus,
(26)1N1L∑t=1L∑i=1nyx,t∫Ωφiξyξ,tdξ=∑i=1nλiφix.LetRx,ξ=1/N1/L∑t=1Lyx,tyξ,t; then
(27)∑i=1n∫ΩRx,ξφiξdξ=∑i=1nλiφix.Considering the mutual independence fori=1,…,n, (27) is true if and only if the below condition is satisfied:
(28)∫ΩRx,ξφiξdξ=λiφix.Similarly, setRt,τ=1/N1/L∑z=1Nyx,tyx,τ, and (29) can be obtained according to (25-b).
(29)∫TRt,τψiτdτ=ηiψit.From the above, the necessary condition for the optimization problem in (20) is
(30)∫ΩRx,ξφiξdξ=λiφix,∫TRt,τψiτdτ=ηiψit,i=1,…,n.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).


#### 3.2.2. Snapshot Method

The snapshot method [28] can get the solution of (30) effectively. Assume that the spatial basis function φix and temporal basis function ψit are the linear combination of a series of snapshots:
(31)φix=∑t=1Lγityx,t,ψit=∑z=1Nωizyx,t.Substitute (31) into (30); then
(32)1N1L∫Ω∑t=1Lyx,tyξ,t∑k=1Lγikyξ,kdξ=λi∑t=1Lγityx,t,1N1L∫T∑z=1Nyx,tyx,τ∑z=1Nωizyz,τdτ=ηi∑z=1Nωizyx,t.Define the correlation function of two points as follows:(33)Ctk=1NL∫Ω∑k=1Lyξ,tyξ,kdξ,Cxz=1NL∫T∑z=1Nyx,τyz,τdτ.Then the optimization problem in (30) is converted into the below matrix eigenvalue problem:
(34)Ctkγi=λiγi,Cxzϖi=ηiϖi,where γi=γi1,…,γiLT and ϖi=ϖi1,…,ωiNT denote the ith eigenvectors.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 φix and temporal basis function ψit.


#### 3.2.3. Dimension Selection

In the above problem, suppose that the maximal number of the nonzero eigenvalue isK, which satisfies K≤minN,L. Sort these eigenvalues in a descending order, then λ1≥λ2≥⋯≥λK correspond to the spatial basis functions φ1x,φ2x,…,φKx, and η1≥η2≥⋯≥ηK correspond to the spatial basis functions ψ1t,ψ2t,…,ψKt.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
(35)Ei=μi∑j=1Kμj.Generally speaking, when the proportion is more than 99%, it can be concluded that these bases can reflect the most energy; the correspondingn 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·i=1n, the below equation can be obtained [29]:
(36)∑i=1nyx,t,ϕi·2=∑i=1nμi≥∑i=1nyx,t,θi·2.This shows that the KL method can get a lower order than other methods in modeling.


### 3.3. Biorthogonal Spatial-Temporal Model

For ASP flooding, the spatial basis functionφixi=1n and temporal basis function ψiti=1l can be obtained by decomposing the state sample information with the method in Section 3.2. And the input basis function ϕixi=1m can be obtained by decomposing the injection concentration of displacing agents cΘin combined with ψiti=1l. Then the expanded biorthogonal spatial-temporal model is presented as follows:The inputcΘin can be expressed as
(37)cΘinz,t=∑i=1mciϕizψit,where ci=cΘinz,t,ϕizψit denotes the coefficient.The unmodeled dynamics and errordx,t can be expressed as
(38)dnx,t=∑i=1ndiφixψit,where dit=dx,t,φixψit denotes the coefficient.Suppose that the blockhx,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 hx,z,τ is
(39)hn,lx,z,τ=∑i=1n∑j=1m∑k=1lαi,j,kφixϕjzψkτ,where αi,j,k∈ℝi=1,…,n,j=1,…,m,k=1,…,l denotes the constant coefficients corresponding to the basis functions φix, ϕjz, 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 thaty=fu,u∈D,andu=φx,x∈Q,u∈Z,D∩Z≠∅. If both y=fu and u=φx exist inverse functions and they are u=f−1y and x=φ−1u, respectively, the composite function y=fφx must exist an inverse function; the function is x=φ−1f−1y.Theorem 1.
For a biorthogonal spatial-temporal Wiener system in Figure3, suppose that ξix is the spatial basis function, ςit is the temporal basis function, ρi is the coefficients, and P·:ℝ⟶ℝ and P−1·:ℝ⟶ℝ are the inverse functions of each other. If υx,t=Nvx,t=∑i=1nρiξixςitPvx,t, for an arbitrary given set of spatial basis function ϖix and temporal basis function wit, there must be a set of coefficients μi which can make vx,t=N−1υx,t=∑i=1nμiϖixwitP−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 forN·.
Asυx,t=∑i=1nξixρitPvx,t, P· and P−1· are the inverse functions of each other, according to Lemma 1:
(40)vx,t=P−11∑i=1nρiξixςitυx,t.
Take the derivative of (40) with respect to yx,t; then
(41)dvx,tdυx,t=dP−1udu⋅dudυx,t=dP−1udu⋅1∑i=1nρiξixςit,where u=1/∑i=1nρiξixςitυx,t.
Forvx,t=N−1υx,t, take the derivative with respect to υx,t; then
(42)dvx,tdυx,t=∑i=1nμiϖixwit⋅dP−1udu⋅dudυx,t=∑i=1nμiϖixwit⋅dP−1udu,where u=υx,t.
To make sure thatvx,t=N−1υx,t and υx,t=Nvx,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ϖixwit satisfies the below relation:
(43)∑i=1nμiϖixwit=1∑i=1nρiξixςit.
From our experience, there are infinite pairs ofwit, ϖix and μi. Let χ=1/n∑i=1nρiξixςit; then
(44)∑i=1nμiϖixwit=1nχ=1n2χ+1n2χ+⋯+1n2χ⏟n.
Ifμiϖixwit≡1/n2χ, then for an arbitrary given set of wit and ϖix, there is the corresponding set of μi=1/n2χwitϖix, 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 wit, ϖix, and μi can be selected approximately under a given error. Among all the combinations, one group of wit, ϖix, and μi by the KL theory can be found.Assume thatN· is invertible; then vx,t can be expressed as (45) according to Theorem 1.
(45)vx,t=∑i=1nP−1υx,tBiφixψit,where Bi is the identification parameter.In addition,vx,t=hx,z,qcΘinz,t+dx,t. That is,
(46)vx,t=∑τ=0t∫Ω∑i=1n∑j=1m∑k=1lαi,j,kφixϕjzψkt−τ∑r=1mϕrzψrτcrdz+dx,t.Uniting (38), (45), and (46), then
(47)∑i=1nP−1υx,tBiφixψit=∑τ=0t∫Ω∑i=1n∑j=1m∑k=1lαi,j,kφixϕjzψkt−τ∑r=1mϕrzψrτcrdz+dx,t=∑i=1nφix∑j=1m∑k=1lαi,j,k∑r=1m∫Ωϕjzϕrzdz∑τ=0tψkt−τψrτcr+∑i=1ndiφixψit.Letϕj,r=∫Ωϕjzϕrzdz and Lk,rt=∑τ=0tψkt−τcr; then
(48)∑i=1nP−1υx,tBiφixψit=∑i=1nφix∑j=1m∑k=1lαi,j,k∑r=1mψrτϕj,rLk,rt+∑i=1ndiφixψit.Project (48) on the temporal basis function ψst and spatial basis function φhx with the Galerkin method. Considering the orthogonal property of temporal basis function, (48) can be inferred.
(49)∑i=1n∫ΩφhxφixdxP−1υx,tBi=∑i=1n∫Ωφhxφixdx∑j=1m∑k=1lαi,j,k∑r=1mϕj,rLk,rt+∑i=1n∫Ωφhxφixdxdi.Since the spatial basis function is orthogonal,n equations can be drawn from (49).
(50)P−1υx,tBi=∑j=1m∑k=1lαi,j,k∑r=1mϕj,rLk,rt+di,i=1,…,n.SetLj,kt=∑r=1mϕj,rLk,rt∈ℝ, B=B1,…,BnT∈ℝn, d=d1,…,dnT∈ℝn, and αj,k=α1,j,k,…,αn,j,kT∈ℝn; then
(51)P−1υx,tB=∑j=1m∑k=1lαj,kLj,kt+d.Without the loss of generality, takeα1,1=I. Equation (51) can be rewritten as
(52)d=B−ATL¯t,where L¯t=L1,2t,…,L1,lt,L2,1t,…,Lm,lt,−P−1υx,tT∈ℝml and A=α1,2,…,α1,l,α2,1,…,αm,l,B∈ℝn×ml.The only unknown parameter in (52) is A which can be obtained by a quadratic criterion on the prediction errors [31] as follows:
(53)A^=argminA1L∑t=1L−L1,1t−ATL¯t2.ParameterA can be identified by using the LSE method [32]:
(54)A^=1L∑t=1LL¯tL¯Tt−1−1L∑t=1LL¯tL1,1Tt.So far the value ofA has been determined, that is, α and B. Parameter ci 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.


### 3.4. Output Identification

To acquire the mathematical description between system stateυ and output fwt, the multivariable ARMA identification [31] is introduced which is identified by the recursive least squares (RLS) [33]. Suppose that there are ns states in ASP flooding. The ARMA model is
(55)fwt=U0q−1fwt+U1q−1υ1x,t+U2q−1υ2x,t+⋯+Unsq−1υnsx,t,where U0q−1=∑i=1N0U0,iq−i, U1q−1=∑i=1N1U1,iq−i, U2q−1=∑i=1N2U2,iq−i, and Unsq−1=∑i=1NnsUns,iq−i; q−1 is the backward shift operator which means backing up one sampling period; U0,U1,…,Uns are the matrix polynomial of q−1; N0 denote the time lag for output; and N1,…,Nns denote the time lag for all states.Convert the above problem into the linear regression model,(56)fwt=θ⋅F⌢x,t,where
(57)F⌢x,t=fwt−1T,…,fwt−N0TTυ1x,t−1T,…,υ1x,t−N1TT⋮υnsx,t−1T,…,υnsx,t−NnsTT,θ=U0,1,…,U0,N0,U1,1,…,U1,N1,U2,1,…,Uns,1,…,Uns,Nns.Identify parameterθ with RLS in (58).
(58)θ^t=θ^t−1+Ktfwt−F⌢x,tTθ^t−1,Kt=Pt−1F⌢x,tF⌢x,tTPt−1F⌢x,t+μ−1,Pt=1μI−KtF⌢x,tTPt−1,where Kt denotes the weight matrix, 0<μ<1 is the forgetting factor, and Pt 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.


## 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.Figure 3
Spatial-temporal Wiener system.The distributed parameter identification contains a linear dynamical block followed by a static nonlinear blockN·:ℝ⟶ℝ, the description is as follows.
(15)υx,t=Nvx,t=N∑τ=0t∫Ωhx,z,τcΘinz,t−τdz+dx,t,where hx,z,q1×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, vx,t denotes the system intermediate variable, cΘinz,t∈ℝ and υx,t∈ℝ stand for the input and state output of the system at time t separately, and dx,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 blockF·:ℝ⟶ℝ. Its aim is to build the relation between system states and target outputs. The formula description is shown in (16).
(16)f⌢wt=Fυx,t+d′x,t,where f⌢wt∈ℝ 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 identifyN and h according to the input-state data and F on the basis of state-output data.


## 3.2. 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 outputyxi,ti=1,t=1N,L is uniformly distributed on the time and space domain, the orthogonal spatial basis function is φixi=1∞, the temporal basis function is ψiti=1∞, and the coefficient is αii=1∞. Then the output can be expressed as follows [25–27]:
(17)yx,t=∑i=1∞αiφixψit.The basis function is orthogonal; that is,(18)φix,φjx=∫Ωφixφjxdx=0,i≠j,1,i=j,ψit,ψjt=∫Tψitψjtdt=0,i≠j,1,i=j,where ·,· denotes the inner product.In practical application, only a few dominant bases which can reflect most system information are chosen. Then(19)ynx,t=∑i=1nαiφixψit,αi=φixψit,yx,t.To obtain the dominant bases, we can get the below minimization problem:(20)minφix,ψit1N1L∫T∫Ωyx,t−ynx,t2dxdt,s.t.φix,φix=1,φix∈L2Ω,i=1,…,n,ψit,ψit=1,ψit∈L2T,i=1,…,n.Bring in a Lagrangian multiplier and convert (20) into an unconstrained optimization problem:
(21-a)J=1N1L∫T∫Ωyx,t−∑i=1nαiφixψit2dxdt+∑i=1nλiφix,φix−1+∑i=1nηiψit,ψit−1=1N1L∫T∫Ωy2x,t−2yx,t∑i=1nαiφixψit+ΛΛTdxdt+∑i=1nλi∫Ωφi2xdx−1+∑i=1nηi∫Tψi2tdt−1,(21-b)Λ=α1α2⋮αnTφ1x0⋯00φ2x0⋮⋮0⋱00⋯0φnxψ1tψ2t⋮ψnt,(21-c)ΛΛT=∑i=1nαi2.Take the variation of (21-a) with respect to φix and ψit; then
(22)δJ=1N1L∫T∫Ω−2yx,t∑i=1nαiψitδφix+−2yx,t∑i=1nαiφixδψitdxdt+∑i=1n2λi∫Ωφixδφixdx+∑i=1n2ηi∫Tψitδψitdt.Equation (22) can be rewritten as
(23)δJ=∫Ω1N1L∫T−2yx,t∑i=1nαiψitdt+∑i=1n2λiφixδφixdx+∫T1N1L∫Ω−2yx,t∑i=1nαiφixdx+∑i=1n2ηiψitδψitdt.The necessary condition of the extreme value for this functional problem isδJ=0. Since φix and ψit can be an arbitrary function, the below equation can be developed:
(24)1N1L∫Tyx,t∑i=1nαiψitdt=∑i=1nλiφix,1N1L∫Ωyx,t∑i=1nαiφixdx=∑i=1nηiψit.For all sample points in the time domain and space domain,∫Tfx,tdt=∑t=1Lfx,t and ∫Ωfx,tdx=∑z=1Nfx,t can be drawn approximately. Here, z=1,2,…,N correspond to x1,x2,…,xN. Then (24) can be written as
(25-a)1N1L∑t=1Lyx,t∑i=1nαiψit=∑i=1nλiφix,(25-b)1N1L∑z=1Nyx,t∑i=1nαiφix=∑i=1nηiψit.In (25-a), αiψit can be regarded as the projection of yx,t on φix; that is, αiψit=yx,t,φix. Thus,
(26)1N1L∑t=1L∑i=1nyx,t∫Ωφiξyξ,tdξ=∑i=1nλiφix.LetRx,ξ=1/N1/L∑t=1Lyx,tyξ,t; then
(27)∑i=1n∫ΩRx,ξφiξdξ=∑i=1nλiφix.Considering the mutual independence fori=1,…,n, (27) is true if and only if the below condition is satisfied:
(28)∫ΩRx,ξφiξdξ=λiφix.Similarly, setRt,τ=1/N1/L∑z=1Nyx,tyx,τ, and (29) can be obtained according to (25-b).
(29)∫TRt,τψiτdτ=ηiψit.From the above, the necessary condition for the optimization problem in (20) is
(30)∫ΩRx,ξφiξdξ=λiφix,∫TRt,τψiτdτ=ηiψit,i=1,…,n.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).


### 3.2.2. Snapshot Method

The snapshot method [28] can get the solution of (30) effectively. Assume that the spatial basis function φix and temporal basis function ψit are the linear combination of a series of snapshots:
(31)φix=∑t=1Lγityx,t,ψit=∑z=1Nωizyx,t.Substitute (31) into (30); then
(32)1N1L∫Ω∑t=1Lyx,tyξ,t∑k=1Lγikyξ,kdξ=λi∑t=1Lγityx,t,1N1L∫T∑z=1Nyx,tyx,τ∑z=1Nωizyz,τdτ=ηi∑z=1Nωizyx,t.Define the correlation function of two points as follows:(33)Ctk=1NL∫Ω∑k=1Lyξ,tyξ,kdξ,Cxz=1NL∫T∑z=1Nyx,τyz,τdτ.Then the optimization problem in (30) is converted into the below matrix eigenvalue problem:
(34)Ctkγi=λiγi,Cxzϖi=ηiϖi,where γi=γi1,…,γiLT and ϖi=ϖi1,…,ωiNT denote the ith eigenvectors.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 φix and temporal basis function ψit.


### 3.2.3. Dimension Selection

In the above problem, suppose that the maximal number of the nonzero eigenvalue isK, which satisfies K≤minN,L. Sort these eigenvalues in a descending order, then λ1≥λ2≥⋯≥λK correspond to the spatial basis functions φ1x,φ2x,…,φKx, and η1≥η2≥⋯≥ηK correspond to the spatial basis functions ψ1t,ψ2t,…,ψKt.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
(35)Ei=μi∑j=1Kμj.Generally speaking, when the proportion is more than 99%, it can be concluded that these bases can reflect the most energy; the correspondingn 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·i=1n, the below equation can be obtained [29]:
(36)∑i=1nyx,t,ϕi·2=∑i=1nμi≥∑i=1nyx,t,θi·2.This shows that the KL method can get a lower order than other methods in modeling.


## 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 outputyxi,ti=1,t=1N,L is uniformly distributed on the time and space domain, the orthogonal spatial basis function is φixi=1∞, the temporal basis function is ψiti=1∞, and the coefficient is αii=1∞. Then the output can be expressed as follows [25–27]:
(17)yx,t=∑i=1∞αiφixψit.The basis function is orthogonal; that is,(18)φix,φjx=∫Ωφixφjxdx=0,i≠j,1,i=j,ψit,ψjt=∫Tψitψjtdt=0,i≠j,1,i=j,where ·,· denotes the inner product.In practical application, only a few dominant bases which can reflect most system information are chosen. Then(19)ynx,t=∑i=1nαiφixψit,αi=φixψit,yx,t.To obtain the dominant bases, we can get the below minimization problem:(20)minφix,ψit1N1L∫T∫Ωyx,t−ynx,t2dxdt,s.t.φix,φix=1,φix∈L2Ω,i=1,…,n,ψit,ψit=1,ψit∈L2T,i=1,…,n.Bring in a Lagrangian multiplier and convert (20) into an unconstrained optimization problem:
(21-a)J=1N1L∫T∫Ωyx,t−∑i=1nαiφixψit2dxdt+∑i=1nλiφix,φix−1+∑i=1nηiψit,ψit−1=1N1L∫T∫Ωy2x,t−2yx,t∑i=1nαiφixψit+ΛΛTdxdt+∑i=1nλi∫Ωφi2xdx−1+∑i=1nηi∫Tψi2tdt−1,(21-b)Λ=α1α2⋮αnTφ1x0⋯00φ2x0⋮⋮0⋱00⋯0φnxψ1tψ2t⋮ψnt,(21-c)ΛΛT=∑i=1nαi2.Take the variation of (21-a) with respect to φix and ψit; then
(22)δJ=1N1L∫T∫Ω−2yx,t∑i=1nαiψitδφix+−2yx,t∑i=1nαiφixδψitdxdt+∑i=1n2λi∫Ωφixδφixdx+∑i=1n2ηi∫Tψitδψitdt.Equation (22) can be rewritten as
(23)δJ=∫Ω1N1L∫T−2yx,t∑i=1nαiψitdt+∑i=1n2λiφixδφixdx+∫T1N1L∫Ω−2yx,t∑i=1nαiφixdx+∑i=1n2ηiψitδψitdt.The necessary condition of the extreme value for this functional problem isδJ=0. Since φix and ψit can be an arbitrary function, the below equation can be developed:
(24)1N1L∫Tyx,t∑i=1nαiψitdt=∑i=1nλiφix,1N1L∫Ωyx,t∑i=1nαiφixdx=∑i=1nηiψit.For all sample points in the time domain and space domain,∫Tfx,tdt=∑t=1Lfx,t and ∫Ωfx,tdx=∑z=1Nfx,t can be drawn approximately. Here, z=1,2,…,N correspond to x1,x2,…,xN. Then (24) can be written as
(25-a)1N1L∑t=1Lyx,t∑i=1nαiψit=∑i=1nλiφix,(25-b)1N1L∑z=1Nyx,t∑i=1nαiφix=∑i=1nηiψit.In (25-a), αiψit can be regarded as the projection of yx,t on φix; that is, αiψit=yx,t,φix. Thus,
(26)1N1L∑t=1L∑i=1nyx,t∫Ωφiξyξ,tdξ=∑i=1nλiφix.LetRx,ξ=1/N1/L∑t=1Lyx,tyξ,t; then
(27)∑i=1n∫ΩRx,ξφiξdξ=∑i=1nλiφix.Considering the mutual independence fori=1,…,n, (27) is true if and only if the below condition is satisfied:
(28)∫ΩRx,ξφiξdξ=λiφix.Similarly, setRt,τ=1/N1/L∑z=1Nyx,tyx,τ, and (29) can be obtained according to (25-b).
(29)∫TRt,τψiτdτ=ηiψit.From the above, the necessary condition for the optimization problem in (20) is
(30)∫ΩRx,ξφiξdξ=λiφix,∫TRt,τψiτdτ=ηiψit,i=1,…,n.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).


## 3.2.2. Snapshot Method

The snapshot method [28] can get the solution of (30) effectively. Assume that the spatial basis function φix and temporal basis function ψit are the linear combination of a series of snapshots:
(31)φix=∑t=1Lγityx,t,ψit=∑z=1Nωizyx,t.Substitute (31) into (30); then
(32)1N1L∫Ω∑t=1Lyx,tyξ,t∑k=1Lγikyξ,kdξ=λi∑t=1Lγityx,t,1N1L∫T∑z=1Nyx,tyx,τ∑z=1Nωizyz,τdτ=ηi∑z=1Nωizyx,t.Define the correlation function of two points as follows:(33)Ctk=1NL∫Ω∑k=1Lyξ,tyξ,kdξ,Cxz=1NL∫T∑z=1Nyx,τyz,τdτ.Then the optimization problem in (30) is converted into the below matrix eigenvalue problem:
(34)Ctkγi=λiγi,Cxzϖi=ηiϖi,where γi=γi1,…,γiLT and ϖi=ϖi1,…,ωiNT denote the ith eigenvectors.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 φix and temporal basis function ψit.


## 3.2.3. Dimension Selection

In the above problem, suppose that the maximal number of the nonzero eigenvalue isK, which satisfies K≤minN,L. Sort these eigenvalues in a descending order, then λ1≥λ2≥⋯≥λK correspond to the spatial basis functions φ1x,φ2x,…,φKx, and η1≥η2≥⋯≥ηK correspond to the spatial basis functions ψ1t,ψ2t,…,ψKt.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
(35)Ei=μi∑j=1Kμj.Generally speaking, when the proportion is more than 99%, it can be concluded that these bases can reflect the most energy; the correspondingn 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·i=1n, the below equation can be obtained [29]:
(36)∑i=1nyx,t,ϕi·2=∑i=1nμi≥∑i=1nyx,t,θi·2.This shows that the KL method can get a lower order than other methods in modeling.


## 3.3. Biorthogonal Spatial-Temporal Model

For ASP flooding, the spatial basis functionφixi=1n and temporal basis function ψiti=1l can be obtained by decomposing the state sample information with the method in Section 3.2. And the input basis function ϕixi=1m can be obtained by decomposing the injection concentration of displacing agents cΘin combined with ψiti=1l. Then the expanded biorthogonal spatial-temporal model is presented as follows:The inputcΘin can be expressed as
(37)cΘinz,t=∑i=1mciϕizψit,where ci=cΘinz,t,ϕizψit denotes the coefficient.The unmodeled dynamics and errordx,t can be expressed as
(38)dnx,t=∑i=1ndiφixψit,where dit=dx,t,φixψit denotes the coefficient.Suppose that the blockhx,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 hx,z,τ is
(39)hn,lx,z,τ=∑i=1n∑j=1m∑k=1lαi,j,kφixϕjzψkτ,where αi,j,k∈ℝi=1,…,n,j=1,…,m,k=1,…,l denotes the constant coefficients corresponding to the basis functions φix, ϕjz, 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 thaty=fu,u∈D,andu=φx,x∈Q,u∈Z,D∩Z≠∅. If both y=fu and u=φx exist inverse functions and they are u=f−1y and x=φ−1u, respectively, the composite function y=fφx must exist an inverse function; the function is x=φ−1f−1y.Theorem 1.
For a biorthogonal spatial-temporal Wiener system in Figure3, suppose that ξix is the spatial basis function, ςit is the temporal basis function, ρi is the coefficients, and P·:ℝ⟶ℝ and P−1·:ℝ⟶ℝ are the inverse functions of each other. If υx,t=Nvx,t=∑i=1nρiξixςitPvx,t, for an arbitrary given set of spatial basis function ϖix and temporal basis function wit, there must be a set of coefficients μi which can make vx,t=N−1υx,t=∑i=1nμiϖixwitP−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 forN·.
Asυx,t=∑i=1nξixρitPvx,t, P· and P−1· are the inverse functions of each other, according to Lemma 1:
(40)vx,t=P−11∑i=1nρiξixςitυx,t.
Take the derivative of (40) with respect to yx,t; then
(41)dvx,tdυx,t=dP−1udu⋅dudυx,t=dP−1udu⋅1∑i=1nρiξixςit,where u=1/∑i=1nρiξixςitυx,t.
Forvx,t=N−1υx,t, take the derivative with respect to υx,t; then
(42)dvx,tdυx,t=∑i=1nμiϖixwit⋅dP−1udu⋅dudυx,t=∑i=1nμiϖixwit⋅dP−1udu,where u=υx,t.
To make sure thatvx,t=N−1υx,t and υx,t=Nvx,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ϖixwit satisfies the below relation:
(43)∑i=1nμiϖixwit=1∑i=1nρiξixςit.
From our experience, there are infinite pairs ofwit, ϖix and μi. Let χ=1/n∑i=1nρiξixςit; then
(44)∑i=1nμiϖixwit=1nχ=1n2χ+1n2χ+⋯+1n2χ⏟n.
Ifμiϖixwit≡1/n2χ, then for an arbitrary given set of wit and ϖix, there is the corresponding set of μi=1/n2χwitϖix, 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 wit, ϖix, and μi can be selected approximately under a given error. Among all the combinations, one group of wit, ϖix, and μi by the KL theory can be found.Assume thatN· is invertible; then vx,t can be expressed as (45) according to Theorem 1.
(45)vx,t=∑i=1nP−1υx,tBiφixψit,where Bi is the identification parameter.In addition,vx,t=hx,z,qcΘinz,t+dx,t. That is,
(46)vx,t=∑τ=0t∫Ω∑i=1n∑j=1m∑k=1lαi,j,kφixϕjzψkt−τ∑r=1mϕrzψrτcrdz+dx,t.Uniting (38), (45), and (46), then
(47)∑i=1nP−1υx,tBiφixψit=∑τ=0t∫Ω∑i=1n∑j=1m∑k=1lαi,j,kφixϕjzψkt−τ∑r=1mϕrzψrτcrdz+dx,t=∑i=1nφix∑j=1m∑k=1lαi,j,k∑r=1m∫Ωϕjzϕrzdz∑τ=0tψkt−τψrτcr+∑i=1ndiφixψit.Letϕj,r=∫Ωϕjzϕrzdz and Lk,rt=∑τ=0tψkt−τcr; then
(48)∑i=1nP−1υx,tBiφixψit=∑i=1nφix∑j=1m∑k=1lαi,j,k∑r=1mψrτϕj,rLk,rt+∑i=1ndiφixψit.Project (48) on the temporal basis function ψst and spatial basis function φhx with the Galerkin method. Considering the orthogonal property of temporal basis function, (48) can be inferred.
(49)∑i=1n∫ΩφhxφixdxP−1υx,tBi=∑i=1n∫Ωφhxφixdx∑j=1m∑k=1lαi,j,k∑r=1mϕj,rLk,rt+∑i=1n∫Ωφhxφixdxdi.Since the spatial basis function is orthogonal,n equations can be drawn from (49).
(50)P−1υx,tBi=∑j=1m∑k=1lαi,j,k∑r=1mϕj,rLk,rt+di,i=1,…,n.SetLj,kt=∑r=1mϕj,rLk,rt∈ℝ, B=B1,…,BnT∈ℝn, d=d1,…,dnT∈ℝn, and αj,k=α1,j,k,…,αn,j,kT∈ℝn; then
(51)P−1υx,tB=∑j=1m∑k=1lαj,kLj,kt+d.Without the loss of generality, takeα1,1=I. Equation (51) can be rewritten as
(52)d=B−ATL¯t,where L¯t=L1,2t,…,L1,lt,L2,1t,…,Lm,lt,−P−1υx,tT∈ℝml and A=α1,2,…,α1,l,α2,1,…,αm,l,B∈ℝn×ml.The only unknown parameter in (52) is A which can be obtained by a quadratic criterion on the prediction errors [31] as follows:
(53)A^=argminA1L∑t=1L−L1,1t−ATL¯t2.ParameterA can be identified by using the LSE method [32]:
(54)A^=1L∑t=1LL¯tL¯Tt−1−1L∑t=1LL¯tL1,1Tt.So far the value ofA has been determined, that is, α and B. Parameter ci 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.


## 3.4. Output Identification

To acquire the mathematical description between system stateυ and output fwt, the multivariable ARMA identification [31] is introduced which is identified by the recursive least squares (RLS) [33]. Suppose that there are ns states in ASP flooding. The ARMA model is
(55)fwt=U0q−1fwt+U1q−1υ1x,t+U2q−1υ2x,t+⋯+Unsq−1υnsx,t,where U0q−1=∑i=1N0U0,iq−i, U1q−1=∑i=1N1U1,iq−i, U2q−1=∑i=1N2U2,iq−i, and Unsq−1=∑i=1NnsUns,iq−i; q−1 is the backward shift operator which means backing up one sampling period; U0,U1,…,Uns are the matrix polynomial of q−1; N0 denote the time lag for output; and N1,…,Nns denote the time lag for all states.Convert the above problem into the linear regression model,(56)fwt=θ⋅F⌢x,t,where
(57)F⌢x,t=fwt−1T,…,fwt−N0TTυ1x,t−1T,…,υ1x,t−N1TT⋮υnsx,t−1T,…,υnsx,t−NnsTT,θ=U0,1,…,U0,N0,U1,1,…,U1,N1,U2,1,…,Uns,1,…,Uns,Nns.Identify parameterθ with RLS in (58).
(58)θ^t=θ^t−1+Ktfwt−F⌢x,tTθ^t−1,Kt=Pt−1F⌢x,tF⌢x,tTPt−1F⌢x,t+μ−1,Pt=1μI−KtF⌢x,tTPt−1,where Kt denotes the weight matrix, 0<μ<1 is the forgetting factor, and Pt 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.


## 4. Optimization Procedures

The biorthogonal spatial-temporal Wiener modeling method in Section3 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 fw of production wells is related to the water saturation Sw of corresponding wells directly, only the relation between Sw and fw 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:
(59-a)maxJNPV=∫0tf1+χ−t/ta∭ΩPoil1−fwqoutt−∑ΘPΘqincΘintdσ−Pcostdt,(59-b)Swx,y,z,t=Nvx,y,z,t=N∑τ=0t∫Ωhx,y,z,ξ,ζ,ς,τcΘinξ,ζ,ς,t−τdz+dx,y,z,t,(59-c)fwt=F^fwt−1,…,f^wt−N0,S^wx¯,y¯,z¯,t−1,…,S^wx¯,y¯,z¯,t−N1,where x¯,y¯,z¯ is the location of the production well, N0,N1 denote the time lag of f^w,S^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 intoP+1stages. The optimization variable is cΘink,k=1,…,P. Note that cΘinP+1≡0.(2)
In initialization, give the slug length and select the initial control feasible regionrΘin and the initial control cΘin0. The number of the initial discrete state grid is M. Select R values for every grid. Set contraction factor γ, maximal iteration lmax, and convergence accuracy ε. Start from generation l=1.(3)
Set the current control domainrΘl=rΘin.(4)
In the current feasible region, generateM−1 injection concentrations with a uniform strategy [35] in the region cΘin∗l−1k−rΘl−1k,cΘin∗l−1k+rΘl−1k in which cΘin∗l−1 is the optimal control strategy obtained from the previous generation. Specially, l=1 corresponds to the initial condition. Calculate the system as (59-b), and get M state trajectories. M state values Swk−1,k=1,…,P+1 at the beginning of each stage constitute M state grids of the corresponding stage.(5)
Start fromP+1th stage, and calculate the state and performance index as (59-a)–(59-c). Record the result.(6)
Start from thePth stage, and generate R injection concentration for every state grid Swk−2 as (60).
(60)cθinlP=cθin∗l−1P+ζrθl−1P,where ζ denotes a diagonal matrix in which all elements on diagonal are random numbers belonging to −1,1 and cθin∗l−1P 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:
(61)cΘinlP=0,cΘinlP<0,cΘmax,cΘinlP>cΘmax.Calculate the system state fromtP−1 to tP in (8). Then select the optimal injection concentration cΘin∗lP+1 which is closest to the value of current tP and calculate the system state from tP to tP+1. Compute the total performance index from tP−1 to tP+1. Compare and record the injection concentration cΘin∗lP with the maximal NPV for every grid at stage P.(7)
Shift the time step forward, and repeat step (6) until the first slug. Compute the optimal performance index of the whole time domain as (59-a), and save the corresponding injection concentration cΘin∗lk with the maximal NPV.(8)
Shrink the feasible domain of decision variables:(62)rΘl+1k=γrΘlk,k=1,2,…,P+1.Take the optimal injection concentration obtained from step (7) as the center of the feasible domain in the next iteration.(9)
Add the iteration timel=l+1. If Jnew−Jold≥ε, go to step (4), and continue computing and iterating. If Jnew−Jold<ε, 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.


## 5. 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 Figure4.Figure 4
The distribution diagram of well position.


#### 5.1.1. 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 × 106 m3. 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.Figure 5
The distribution diagram of initial permeability.Figure 6
The distribution diagram of initial pressure. The unit is MPa.Figure 7
The distribution diagram of initial water saturation.The water injection rate of each injection well is 83 m3/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].Table 1
The liquid volume quantity of production wells.
Wells
S1-1
S1-2
S1-3
S2-1
S2-2
S2-3
S3-1
S3-2
S3-3
m3/day
20.75
41.5
20.75
83
41.5
41.5
20.75
41.5
20.75


### 5.2. 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 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:(63)RMSE=∫∑ex,y,z,t2dxdydz∫dxdydz∑Δt,where ex,y,z,t=S^wx,y,z,t−Swx,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 the coefficients and calculate the averageRMSE which is shown in Table 2.Table 2
R
M
S
E for different numbers of basis function.
n
1
2
3
4
5
R
M
S
E
0.0751
0.0502
0.0296
0.0209
0.0197Note that the number of spatial basis function is equal to that of temporal basis function. This is decided by the matrix operations. According to Table2, 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 N0=5; the time lag of water saturation is N1=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:(64)cpin=2.7,1.4,1.1kg/m3,cOHin=3.7,2.8,1.3kg/m3,csin=2.9,1.3,0.6kg/m3.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 Figure8(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 in Figures 8 and 9.Figure 8
Comparison diagram of water saturation.
(a)
The injection mass concentration
(b)
The output of real system
(c)
The output of modelingFigure 9
Modeling error of water saturation.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.Figure 10
Comparison for the moisture content of nine production wells.From the above, the biorthogonal spatial-temporal Wiener modeling method can model for ASP flooding well. It has good modeling accuracy and generalization ability.


### 5.3. 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 Section4 to optimize the optimal injection strategy for alkali, surfactant, and polymer.


#### 5.3.1. Parameter Setting

P
=
3, Pp=6.85$/kg, Ps=8.31$/kg, Pa=2.86$/kg, the production cost per day is Pcost=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 is equal to 0.346$/L. r=0.85, cΘmax=4, R=5, M=3, ε=1×10−4, η=0.7, MOH=1300 t, Ms=1000 t, and Mp=700 t. The initial control strategy is cOHin=2.9 kg/m3, cpin=2.1 kg/m3, and csin=1.8 kg/m3.In order to illustrate the accuracy of the proposed method, the mechanism model is also solved with the IDP algorithm in Section4 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
(65)cOHin=3.5,2.59,1.8kg/m3,cpin=2.26,1.17,0.69kg/m3,csin=2.42,1.5,0.49kg/m3.Figure 11
Comparison of the optimal injection concentration.Figure 12
The moisture content and oil production of production wells for PM.Figure 13
The moisture content and oil production of production wells for MM.Figure 14
The average moisture content comparison of production wells.Figure 15
The average oil production comparison of 9 production wells.In Figure11, 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 Figures12 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 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 Table3.Table 3
Parameter comparison.
Parameters
Mechanism model (MM)
Proposed method (PM)
Consumption of polymer (t)
659.312
656.563
Consumption of surfactant (t)
701.184
702.777
Consumption of alkali (t)
1258.944
1257.350
Oil production (t)
143,371
143,060
NPV (million $)
47.05
46.97By 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 by311t, 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.


## 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 Figure4.Figure 4
The distribution diagram of well position.


### 5.1.1. 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 × 106 m3. 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.Figure 5
The distribution diagram of initial permeability.Figure 6
The distribution diagram of initial pressure. The unit is MPa.Figure 7
The distribution diagram of initial water saturation.The water injection rate of each injection well is 83 m3/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].Table 1
The liquid volume quantity of production wells.
Wells
S1-1
S1-2
S1-3
S2-1
S2-2
S2-3
S3-1
S3-2
S3-3
m3/day
20.75
41.5
20.75
83
41.5
41.5
20.75
41.5
20.75


## 5.1.1. 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 × 106 m3. 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.Figure 5
The distribution diagram of initial permeability.Figure 6
The distribution diagram of initial pressure. The unit is MPa.Figure 7
The distribution diagram of initial water saturation.The water injection rate of each injection well is 83 m3/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].Table 1
The liquid volume quantity of production wells.
Wells
S1-1
S1-2
S1-3
S2-1
S2-2
S2-3
S3-1
S3-2
S3-3
m3/day
20.75
41.5
20.75
83
41.5
41.5
20.75
41.5
20.75


## 5.2. 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 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:(63)RMSE=∫∑ex,y,z,t2dxdydz∫dxdydz∑Δt,where ex,y,z,t=S^wx,y,z,t−Swx,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 the coefficients and calculate the averageRMSE which is shown in Table 2.Table 2
R
M
S
E for different numbers of basis function.
n
1
2
3
4
5
R
M
S
E
0.0751
0.0502
0.0296
0.0209
0.0197Note that the number of spatial basis function is equal to that of temporal basis function. This is decided by the matrix operations. According to Table2, 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 N0=5; the time lag of water saturation is N1=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:(64)cpin=2.7,1.4,1.1kg/m3,cOHin=3.7,2.8,1.3kg/m3,csin=2.9,1.3,0.6kg/m3.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 Figure8(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 in Figures 8 and 9.Figure 8
Comparison diagram of water saturation.
(a)
The injection mass concentration
(b)
The output of real system
(c)
The output of modelingFigure 9
Modeling error of water saturation.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.Figure 10
Comparison for the moisture content of nine production wells.From the above, the biorthogonal spatial-temporal Wiener modeling method can model for ASP flooding well. It has good modeling accuracy and generalization ability.


## 5.3. 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 Section4 to optimize the optimal injection strategy for alkali, surfactant, and polymer.


### 5.3.1. Parameter Setting

P
=
3, Pp=6.85$/kg, Ps=8.31$/kg, Pa=2.86$/kg, the production cost per day is Pcost=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 is equal to 0.346$/L. r=0.85, cΘmax=4, R=5, M=3, ε=1×10−4, η=0.7, MOH=1300 t, Ms=1000 t, and Mp=700 t. The initial control strategy is cOHin=2.9 kg/m3, cpin=2.1 kg/m3, and csin=1.8 kg/m3.In order to illustrate the accuracy of the proposed method, the mechanism model is also solved with the IDP algorithm in Section4 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
(65)cOHin=3.5,2.59,1.8kg/m3,cpin=2.26,1.17,0.69kg/m3,csin=2.42,1.5,0.49kg/m3.Figure 11
Comparison of the optimal injection concentration.Figure 12
The moisture content and oil production of production wells for PM.Figure 13
The moisture content and oil production of production wells for MM.Figure 14
The average moisture content comparison of production wells.Figure 15
The average oil production comparison of 9 production wells.In Figure11, 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 Figures12 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 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 Table3.Table 3
Parameter comparison.
Parameters
Mechanism model (MM)
Proposed method (PM)
Consumption of polymer (t)
659.312
656.563
Consumption of surfactant (t)
701.184
702.777
Consumption of alkali (t)
1258.944
1257.350
Oil production (t)
143,371
143,060
NPV (million $)
47.05
46.97By 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 by311t, 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.


## 5.3.1. Parameter Setting

P
=
3, Pp=6.85$/kg, Ps=8.31$/kg, Pa=2.86$/kg, the production cost per day is Pcost=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 is equal to 0.346$/L. r=0.85, cΘmax=4, R=5, M=3, ε=1×10−4, η=0.7, MOH=1300 t, Ms=1000 t, and Mp=700 t. The initial control strategy is cOHin=2.9 kg/m3, cpin=2.1 kg/m3, and csin=1.8 kg/m3.In order to illustrate the accuracy of the proposed method, the mechanism model is also solved with the IDP algorithm in Section4 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
(65)cOHin=3.5,2.59,1.8kg/m3,cpin=2.26,1.17,0.69kg/m3,csin=2.42,1.5,0.49kg/m3.Figure 11
Comparison of the optimal injection concentration.Figure 12
The moisture content and oil production of production wells for PM.Figure 13
The moisture content and oil production of production wells for MM.Figure 14
The average moisture content comparison of production wells.Figure 15
The average oil production comparison of 9 production wells.In Figure11, 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 Figures12 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 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 Table3.Table 3
Parameter comparison.
Parameters
Mechanism model (MM)
Proposed method (PM)
Consumption of polymer (t)
659.312
656.563
Consumption of surfactant (t)
701.184
702.777
Consumption of alkali (t)
1258.944
1257.350
Oil production (t)
143,371
143,060
NPV (million $)
47.05
46.97By 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 by311t, 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.


## 6. 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 spatial-temporal 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 flooding which can avoid solving the complex mechanism equations. It is also suitable for solving other complex optimization problems in engineering.



---

*Source: 9248161-2018-10-17.xml*