# A Numerical Computation Approach for the Optimal Control of ASP Flooding Based on Adaptive Strategies

Authors: Shurong Li; Yulei Ge
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A numerical computation approach based on constraint aggregation and pseudospectral method is proposed to solve the optimal control of alkali/surfactant/polymer (ASP) flooding. At first, all path constraints are aggregated into one terminal condition by applying a Kreisselmeier-Steinhauser (KS) function. After being transformed into a multistage problem by control vector parameter, a normalized time variable is introduced to convert the original problem into a fixed final time optimal control problem. Then the problem is discretized to nonlinear programming by using Legendre-Gauss pseudospectral method, whose numerical solutions can be obtained by sequential quadratic programming (SQP) method through solving the KKT optimality conditions. Additionally, two adaptive strategies are applied to improve the procedure:(1) the adaptive constraint aggregation is used to regulate the parameter ρ in KS function and (2) the adaptive Legendre-Gauss (LG) method is used to adjust the number of subinterval divisions and LG points. Finally, the optimal control of ASP flooding is solved by the proposed method. Simulation results show the feasibility and effectiveness of the proposed method.


---

## Body

## 1. Introduction

With the exploitation of oil, most oil fields of China have stepped into the high water cut period [1]. The production of oil cannot satisfy our demands for daily life and economic development. A series of tertiary oil recovery technologies such as chemical flooding [2], microorganism flooding [3], and carbon dioxide flooding [4] are put to use. The newly emerging alkali/surfactant/polymer (ASP) flooding which is an important tertiary oil recovery technology can enhance oil production obviously. The basic idea is to utilize three displacing agents (alkali, surfactant, and polymer), whose synergistic effects are important to enhanced oil recovery, to change the physicochemical property [5]. Since the price of displacing agents is high, how to determine the injection strategy to maximize the profit as much as possible is always a challenging problem.The essence of determining the injection strategy for ASP flooding is an optimal control problem. In current industrial application, the index comparison method is usually adopted to select the injection strategy, in which the best one is chosen among many given feasible strategies by simulating on the numerical simulation software according to a defined index [6, 7]. This method is simple and easy to be operated, but it is too dependent on the experience of manipulators. The optimization result is not the optimum. The optimal control technology, which is first used in oil exploitation in [8], works on searching for the optimal control strategy from all feasible solutions with considering the process dynamic characteristics, and it can realize the simultaneous optimization of all variables. In addition, the whole theory is proved by mathematical method; it is scientific and reasonable. Therefore, it is suitable for the control of ASP flooding.Some scholars have studied the optimal control for oil exploitation. Jansen et al. [9] studied the problems of optimal control and nonlinear model predictive control for water flooding, in which the control variables are bottom hole pressure and flooding rate, and the index is the maximum of net present value (NPV). Lei et al. [10] presented mixed-integer iterative dynamic programming (IDP) to optimize the polymer flooding and got good result. Ramirez et al. [11, 12] optimized the injection strategy for enhanced oil recovery of surfactant flooding with optimal control theory, in which the necessary condition of optimal control was deduced on the basis of maximum principle before being solved by gradient method. Furthermore, this method was applied to carbon dioxide flooding, nitrogen flooding, and binary system flooding [13]. Zerpa et al. [14] used field scale numerical simulation and multiple surrogates to optimize alkaline–surfactant–polymer flooding process based on UTCHEM. Ge et al. [15] developed an approximate dynamic programming method to solve the optimization of ASP flooding, in which an Actor-Critic algorithm is introduced to search the optimal injection strategy. Most of the researches are about water flooding and single chemical flooding. The optimal control for ASP flooding needs to be further studied urgently.With regard to the research on optimal control methods, although extensive research on optimal control has been conducted, large-scale nonlinear dynamic models processing and efficient solution methods for optimal control problems are still two main barriers for the widespread industrial application of optimal control [16, 17]. Numerical methods for optimal control can be divided into two classes, which are direct methods and indirect ones. Direct methods, which include the control vector parameterization (CVP) [18], the direct multiple shooting [19], and the full discretization [20], are closely related to the discretization of the original control problem and the application of nonlinear programming technique. These methods are very popular and well suited for many practical problems [21]. What differs in these methods is how to select the kind of variables to be discretized and how to approximate the system equations. However, direct methods may bring about unsatisfactory results especially when the control variables are discontinuous, such as with switching points and singular arcs [22]. For these reasons, the reasonable discretization methods for variables and approximation methods for system equations are quite important.As to the discretization, if the discretization grid is too coarse, this will lead to some local controls being not in the right place and the accuracy cannot meet the prespecified requirements. What is more, if the discretization grid is too fine, the computational cost will be very high and the robustness will be poor. Generally speaking, the discretization grid is determined by visually examining the obtained result which is balanced between efficiency and precision. For example, an adaptive CVP method was proposed in [23], in which the problem was discretized adaptively over time spans. The discretization was sequentially refined based on a wavelet-based analysis of the optimal solution which was obtained in the previous optimization step. Thus, the efficiency is enhanced with less loss of accuracy.As to the approximation methods for system equations, pseudospectral method, which is a major kind of direct method, has been widely used in recent years [24]. In this method, states and controls are approximated by a series of orthogonal basis functions, such as Lagrange, Laguerre, and Legendre polynomials. Some famous pseudospectral methods have been developed in the past decades, for example, the Gauss pseudospectral method (GPM) [25], the Lobatto pseudospectral method (LPM) [26], and the Radau pseudospectral method (RPM) [27]. The most obvious difference between them is the selection of interpolation point. For smooth optimal control problems, these methods can obtain good accuracy and performance compared with traditional direct methods. But when the control is discontinuous, for example, the three-slug injection strategy of ASP flooding, there may be some problems especially in the discontinuous points. So, many adaptive methods are proposed to cope with the interpolation points and subintervals. Darby et al. [28] developed an adaptive method, in which the difference between a Legendre-Gauss-Radau approximation to the integral of the dynamics and an interpolated value of state was used to estimate the approximation error. According to this error, the interpolation points were selected under a given accuracy requirement.Constraint handling is an important step for a control problem before being optimized. In general, the constraints are disposed by penalty function method, in which a penalty factor is introduced to convert the original problem into an unconstrained optimization problem [29]. But there is no definite theory to determine the penalty factor, which is usually obtained by experience. If the value of the factor is not chosen scientifically, this will lead to adverse impact on the optimization. In [30], the constraint aggregation method, which was expressed as the Kreisselmeier-Steinhauser (KS) function, was first proposed to transform the path constraints into a terminal condition. To regulate the parameter ρ, an adaptive constraint aggregation method was presented by Zhang et al. [31], in which the value of ρ was changed according to the partial differential automatically.Many researches about the optimization algorithms have been studied. Fonseca proposed a stochastic simplex approximate gradient (StoSAG) for optimization under uncertainty, in which the gradient is estimated approximately by an ensemble of randomly chosen control vectors, known as Ensemble Optimization (EnOpt) in the oil and gas reservoir simulation community [32]. Chen and Reynolds studied the optimal control of inflow control valves and well operating conditions for the water-alternating-gas injection process [33]. Shirangi et al. developed a new methodology for the joint optimization of economic project life and time-varying well controls [34]. Sequential quadratic programming (SQP) is an effective method to solve constrained optimization problems, which has good global convergence and more than one order of local convergence and does not need to construct the penalty function factor when dealing with constraints [35]. It is suitable for the optimal control of ASP flooding. In [36], an approximate feasible direction method was put forward in which the constraint aggregation was adopted to transform all path constraints into a generalized constraint. Then SQP was applied to solve this optimal control problem. Bernardo et al. [37] presented a slug optimization method for water flooding in which the Kriging interpolation was utilized to build the model between variables and index before being optimized by SQP.To solve the optimal control for alkali/surfactant/polymer (ASP) flooding which has discontinuous control variables, a methodology based on constraint aggregation and pseudospectral method with adaptive strategies is proposed in this paper. The rest of this article is outlined as follows. In Section2, the specific model of ASP flooding is given and the optimal control for ASP flooding is converted into a more general problem in mathematical form. Section 3 introduces the details of the proposed methodology in this paper, such as processing the path constraints with constraints aggregation method, converting original problem into a multistage problem with CVP, introducing the new time variable, and discretizing the control problem with LG pseudospectral method. The KKT optimality condition of this optimal control problem is provided in Section 3.4. In Section 4, the adaptive strategies for the constraints aggregation and LG pseudospectral method are presented to regulate the control method. Then SQP is introduced to solve the problem. Section 5 applies the method proposed in this paper to solve the optimal control problem for ASP flooding and Section 6 contains some discussion and concluding remarks. Finally, some essential information for ASP flooding is given in the Appendix.


## 2. Problem Formulations

As to the optimal control problem for ASP flooding, the control variables are the injection concentrations of three displacing agents (alkali, surfactant, and polymer) at all injection wells; the states cover pressure, grid concentration, and water saturation; and the performance index is the maximal NPV. This is a complex distributed parameter control problem aiming at getting the optimal injection strategy to fulfill the maximal profit. Since the three-dimensional model for ASP flooding, which includes a series of divergences and cross terms, is too complex, we only consider the one-dimensional model in this paper. In this section, we will obtain the optimal injection strategy with the method proposed in this paper.


### 2.1. Optimal Control Model Description for One-Dimensional ASP Flooding

In view of [5, 6, 38–40], we can make the following descriptions.Considering a long tube core with diameterd and length L, inject the displacing agents liquid with flow Q from one side at a constant speed; the core porosity is ϕ and the residual oil saturation is Sor. On the basis of the model assumptions in [39], we add the following assumptions:(a)
All adsorption processes satisfy Langmuir isothermal adsorption equation.(b)
The displacing agents exist in water phase, the adsorption satisfies the generalized FICK law, and the balance is established momentarily.On the basis of the oil/water seepage continuity equations and adsorption diffusion equations of displacing agents, combining mass balance conditions, we can obtain the following model with interaction of alkali, surfactant, and polymer fully considered.The seepage continuity equation is(1)∂Sw∂t=-QAϕ∂fw∂z.The adsorption diffusion equation of surfactant is(2)ϕSw∂Cs∂t=-vw∂Cs∂z+∂∂zDsϕ∂Cs∂z-ρr∂Γs∂t.The adsorption diffusion equation of polymer is(3)ϕSw∂Cp∂t=-vw∂Cp∂z+∂∂zDpϕ∂Cp∂z-ρr∂Γp∂t.The adsorption diffusion equation of alkali is(4)ϕSw∂Ca∂t=-vw∂Ca∂z+∂∂zDaϕ∂Ca∂z-ρr∂Γa∂t-Ra,where a denotes the alkali, s denotes the surfactant, p denotes the polymer, Sw is the water saturation, A is the core cross section area, fw is the moisture content, vw denotes the seepage speed of water phase, and Ca,Cs,Cp and Da,Ds,Dp denote the concentration and diffusion coefficient of alkali, surfactant, and polymer, respectively. ρr denotes the core density, Γa,Γs,Γp are the adsorbing capacity of core for different displacing agents, and Ra is the alkali consumption.The initial conditions are(5a)Swz,tt=0=1-Sor,(5b)CΘt=0=0,Θ=a,s,p.The boundary conditions are(6a)fwz=0=1.0,(6b)Swz=0=1-Sor,(6c)∂CΘ∂zz=0=vwDΘϕCΘ-uΘ,∂CΘ∂zz=L=0,where uΘ denotes the injection concentration of displacing agents, which is the control variable.In application, the slug injection strategy is usually adopted. Suppose that there areP slugs,(7)uΘt=uΘi,ti-1≤t≤ti,i=1,2,…,P,0,tP≤t≤tf,where ti denotes the time node, the length of every slug is Ti=ti-ti-1, and uΘi denotes the injection concentration of displacing agent Θ in slug i.Furthermore, the dosage limit of displacing agents is(8)∑i=1PuΘi·Ti≤MΘp,where MΘp denotes the maximum usage of displacing agents Θ.The injection concentration and slug size limitations are(9)0≤uΘi≤uΘmax,∑i=1PTi=tP.The other physicochemical algebraic equations can be found in the Appendix.The maximum net present value (NPV) is chosen as the index; the specific description is(10)max⁡JNPV=∫0tf1+χ-tξo1-fwqout-ξΘqinuΘdt,where χ denotes the discount rate, qin and qout denote the volume flow rate of injection and production, and ξΘ denotes the price of three displacing agents.


### 2.2. Model Transformation

To make the problem more concise, we usex∈Rnw (CΘ,Sw,ΓΘ) to denote all states and u∈Rnu (the injection concentration of all displacing agents uΘ) to denote all control variables. Since the ASP flooding model contains many partial differentials with respect to time variable and spatial variable, it is difficult to solve with conventional methods. For simplicity, the finite difference method is used to discretize the spatial grid [6]. Thus there is only the differential with respect to the time variable. Then the optimization for ASP flooding can be reformulated as in the following continuous Bolza problem [41]:(11a)maxut,tfJutf,tf=θxtf,tf+∫t0tfϕxt,ut,tdt,(11b)s.t.x˙t=fxt,ut,t,xt0=x0,(11c)hxt,ut,t≤0,(11d)Nxtf,tf≤0,(11e)uLB≤ut≤uUB,where xt∈Rnx denotes the state variables including water saturation, pressure, and grid concentration. The initial state is x0. ut∈Rnu is the control variable including the injection concentration of displacing agents (alkali, surfactant, and polymer). The system process is described by function f∙∈Rnf. The path constraints and terminal constraints are expressed as h∙∈Rnh and N∙∈RnN, respectively. The whole process is optimized on the time interval t∈t0,tf. The interval of control is uLB,uUB associated with the lower bound and upper bound, which can be reformulated to two path constraints, -ut+uLB≤0 and ut-uUB≤0. The equality constraints can be converted into the inequality constraints by introducing the relaxing factors. For example, an equation gx=0 can be expressed as gx-ς≤0, in which ς denotes the relaxing factor.Then we take a series of measures to process this problem. The details will be presented in the following contents.


## 2.1. Optimal Control Model Description for One-Dimensional ASP Flooding

In view of [5, 6, 38–40], we can make the following descriptions.Considering a long tube core with diameterd and length L, inject the displacing agents liquid with flow Q from one side at a constant speed; the core porosity is ϕ and the residual oil saturation is Sor. On the basis of the model assumptions in [39], we add the following assumptions:(a)
All adsorption processes satisfy Langmuir isothermal adsorption equation.(b)
The displacing agents exist in water phase, the adsorption satisfies the generalized FICK law, and the balance is established momentarily.On the basis of the oil/water seepage continuity equations and adsorption diffusion equations of displacing agents, combining mass balance conditions, we can obtain the following model with interaction of alkali, surfactant, and polymer fully considered.The seepage continuity equation is(1)∂Sw∂t=-QAϕ∂fw∂z.The adsorption diffusion equation of surfactant is(2)ϕSw∂Cs∂t=-vw∂Cs∂z+∂∂zDsϕ∂Cs∂z-ρr∂Γs∂t.The adsorption diffusion equation of polymer is(3)ϕSw∂Cp∂t=-vw∂Cp∂z+∂∂zDpϕ∂Cp∂z-ρr∂Γp∂t.The adsorption diffusion equation of alkali is(4)ϕSw∂Ca∂t=-vw∂Ca∂z+∂∂zDaϕ∂Ca∂z-ρr∂Γa∂t-Ra,where a denotes the alkali, s denotes the surfactant, p denotes the polymer, Sw is the water saturation, A is the core cross section area, fw is the moisture content, vw denotes the seepage speed of water phase, and Ca,Cs,Cp and Da,Ds,Dp denote the concentration and diffusion coefficient of alkali, surfactant, and polymer, respectively. ρr denotes the core density, Γa,Γs,Γp are the adsorbing capacity of core for different displacing agents, and Ra is the alkali consumption.The initial conditions are(5a)Swz,tt=0=1-Sor,(5b)CΘt=0=0,Θ=a,s,p.The boundary conditions are(6a)fwz=0=1.0,(6b)Swz=0=1-Sor,(6c)∂CΘ∂zz=0=vwDΘϕCΘ-uΘ,∂CΘ∂zz=L=0,where uΘ denotes the injection concentration of displacing agents, which is the control variable.In application, the slug injection strategy is usually adopted. Suppose that there areP slugs,(7)uΘt=uΘi,ti-1≤t≤ti,i=1,2,…,P,0,tP≤t≤tf,where ti denotes the time node, the length of every slug is Ti=ti-ti-1, and uΘi denotes the injection concentration of displacing agent Θ in slug i.Furthermore, the dosage limit of displacing agents is(8)∑i=1PuΘi·Ti≤MΘp,where MΘp denotes the maximum usage of displacing agents Θ.The injection concentration and slug size limitations are(9)0≤uΘi≤uΘmax,∑i=1PTi=tP.The other physicochemical algebraic equations can be found in the Appendix.The maximum net present value (NPV) is chosen as the index; the specific description is(10)max⁡JNPV=∫0tf1+χ-tξo1-fwqout-ξΘqinuΘdt,where χ denotes the discount rate, qin and qout denote the volume flow rate of injection and production, and ξΘ denotes the price of three displacing agents.


## 2.2. Model Transformation

To make the problem more concise, we usex∈Rnw (CΘ,Sw,ΓΘ) to denote all states and u∈Rnu (the injection concentration of all displacing agents uΘ) to denote all control variables. Since the ASP flooding model contains many partial differentials with respect to time variable and spatial variable, it is difficult to solve with conventional methods. For simplicity, the finite difference method is used to discretize the spatial grid [6]. Thus there is only the differential with respect to the time variable. Then the optimization for ASP flooding can be reformulated as in the following continuous Bolza problem [41]:(11a)maxut,tfJutf,tf=θxtf,tf+∫t0tfϕxt,ut,tdt,(11b)s.t.x˙t=fxt,ut,t,xt0=x0,(11c)hxt,ut,t≤0,(11d)Nxtf,tf≤0,(11e)uLB≤ut≤uUB,where xt∈Rnx denotes the state variables including water saturation, pressure, and grid concentration. The initial state is x0. ut∈Rnu is the control variable including the injection concentration of displacing agents (alkali, surfactant, and polymer). The system process is described by function f∙∈Rnf. The path constraints and terminal constraints are expressed as h∙∈Rnh and N∙∈RnN, respectively. The whole process is optimized on the time interval t∈t0,tf. The interval of control is uLB,uUB associated with the lower bound and upper bound, which can be reformulated to two path constraints, -ut+uLB≤0 and ut-uUB≤0. The equality constraints can be converted into the inequality constraints by introducing the relaxing factors. For example, an equation gx=0 can be expressed as gx-ς≤0, in which ς denotes the relaxing factor.Then we take a series of measures to process this problem. The details will be presented in the following contents.


## 3. The Numerical Computation Approach Based on Constraint Aggregation and Pseudospectral Discretization

In this section, all path constraints are transformed into one terminal constraint by the constraints aggregation method. Furthermore, CVP is introduced to convert the original problem into a multistage problem (MSP) to cope with the situation of discontinuous control variables after the new time variable is brought in. Finally, to solve the optimal control problem, LG pseudospectral method is adopted to transform the MSP into a general discrete nonlinear programming problem (NLP) [42].


### 3.1. Constraint Aggregation

The main idea of constraints aggregation is the KS function, which was developed by Kreisselmeier and Steinhauser in 1976 [30]. The mathematical form is(12)KShx,ρ=1ρln⁡∑j=1Nceρhjx,where Nc denotes the number of path constraints hjx. ρ is an approximation parameter.KS function can generate one envelope surface which isC1-consecutive. This can estimate the maximal function of hx conservatively. The properties are as follows [43].(a)
For an arbitraryρ>0, there exist the following formulas:(13)max⁡hx≤KShx,ρ≤max⁡hx+1ρln⁡Nc,limρ→∞⁡KShx,ρ=KShx,ρ.(b)
Ifρ2>ρ1>0, then KShx,ρ1≥KShx,ρ2.(c)
KShx,ρ is a convex function, if and only if it is convex for all constraints hx.From the above properties, KS function is the underestimation for feasible region. The bigger the value ofρ is, the more accurate the estimation for the maximal constraint is. This can be obviously shown in Figure 1. What is more, if the original problem is convex, the problem after being processed by KS function is convex too.Figure 1
Estimate ability of KS function formax⁡g1x=5/log⁡x-0.2x-3,g2x=x2/40+x/5-1.The path constraints in (11a), (11b), (11c), (11d), and (11e) can be transformed into a terminal constraint using KS function; the detailed process is as follows [44].Considering the point constraints on the path constraints in (12), ti∈t0,tf, then(14)KShx,ρ=1ρln⁡∑j=1Nc∑i=1eρhjxti.Transform it into the continuous form by introducing the impulse function,(15)KShx,ρ=1ρln⁡∫t0tf∑j=1Nc∑i=1eρhjxtiδt-tidt,where δt-ti is the impulse function.Define a new statep,(16)p˙=fp=∑j=1Nc∑i=1eρhjxtiδt-ti,pt0=0.Substitute (16) into (15); then (17)KShx,ρ=1ρln⁡ptf-pt0.Sincehxt,ut,t≤0, the constraints processed meet the following condition:(18)KShx,ρ≤0.In consequence, the path constraints can be transformed into the following terminal time constraint:(19)pxtf,utf,tf-1≤0.Moreover, a new state variableq is introduced to the index in (11a), which satisfies(20)q˙=fq=ϕxt,ut,t,qt0=0.Then(21)qxtf,utf,tf=∫t0tfϕxt,ut,tdt.Letω=x,p,qT, F=f,fp,fqT; the original problem can be transformed into the following form:(22a)maxut,tfJut,tf=θxtf,tf+qxtf,utf,tf,(22b)s.t.ω˙t=Fxt,ut,t,(22c)pxtf,utf,tf-1≤0,(22d)Nxtf,tf≤0,(22e)xt0=x0,pt0=qt0=0.


### 3.2. Multistage Problem Formulation

To cope with the discontinuity points of control variables preferably, CVP is adopted to convert the original problem into a MSP [23]. The main idea of this process is to approximate the control with piecewise function according to the slug injection strategy. Discretize the whole time domain t0,tf into K uneven time periods. The specific time control nodes are(23)t0<t1<t2<⋯<tK=tf.Approximate the continuous control variable with piecewise function in every period; that is,(24)ut=u^k,t∈tk-1,tk,k=1,2,…,K.Furthermore, we deal with the time variable by defining a new time variableτ:(25)τ=t-tk-1tk-tk-1,k=1,2,…,K.Then every original periodt∈tk-1,tk is transformed into τ∈0,1.Take derivation of (25) with respect to τ; then,(26)dτ=1tk-tk-1dt,1≤k≤K.Definexkτ,ukτ as the state and control on kth period, respectively. Then, (22a), (22b), (22c), (22d), and (22e) can be reformulated as(27a)maxukτ,tkJukτ,tk=θxK1,1+qxK1,uK1,1,(27b)s.t.ω˙kτ=tk-tk-1Fxkτ,ukτ,τ,(27c)pxK1,uK1,1-1≤0,(27d)NxK1,1≤0,(27e)x10=x0,p10=q10=0,(27f)xk1=xk+10,1≤k≤K-1,where ω=x,p,qT and F=f,fp,fqT. Equation (27f) is the linkage constraint which is used to keep the continuity of state variables at the interface of all subintervals. The MSP in (27a), (27b), (27c), (27d), (27e), and (27f) consists of K subproblems corresponding to K subintervals. The final control is obtained by solving K optimal control problems.


### 3.3. Legendre-Gauss Pseudospectral Discretization

For being solved by well-developed optimization algorithms, Legendre-Gauss (LG) pseudospectral discretization is applied to transform the original optimal control problem into a NLP on every subinterval [26, 28, 42].ConstructNk+1-order Lagrange interpolation polynomials with Nk LG points to approximate the states and controls in the subinterval k∈1,2,…,K; then,(28a)ωkτ≈Wkτ=∑j=0NkWjkLjkτ.(28b)Ljkτ=∏l=0,l≠jNkτ-τlkτjk-τlk,j=0,1,2,…,Nk.(29a)ukτ≈Ukτ=∑i=0NkUikL^ikτ,(29b)L^ikτ=∏l=0,l≠iNkτ-τlkτik-τlk,i=1,2,…,Nk,where τ0k,τ1k,…,τNk-1k denote the LG interpolation points that are defined on every subinterval τk∈tk-1,tk. Wk(τ),Uk(τ) are the approximations of state and control on corresponding subinterval. Note that the interpolation points of state are one more than that of control.Take differential of state with respect toτ; then,(30)W˙kτ=∑j=0NkWjkL˙jkτ.The differential of Lagrange interpolation polynomials at every LG point can be expressed by a differential calculation matrixD∈RNk×Nk+1 which can be described as(31)Dsjk=L˙jkτ=∑c=0Nk∏l=0,l≠j,cNkτs-τlk∏l=0,l≠jNkτjk-τlk,s=0,1,2,…,Nk-1.Then the dynamic constraints in (27b) can be transformed into the following algebraic constraints:(32)∑j=0NkDsjkWjk-tk-tk-1FXskτ,Uskτ,τ=0.Then the optimal control problem can be reformulated into a new NLP:(33a)maxUkτ,tkJUkτ,tk=θXNkK1+qXNkK1,UNkK1,(33b)s.t.∑j=0NkDsjkWjk-tk-tk-1FXskτ,Uskτ,τ=0,s=0,1,2,…,Nk-1,(33c)pXNkK1,UNkK1-1≤0,(33d)NXNkK1≤0,(33e)X010=x0,p010=q010=0,(33f)XNkk1=X0k+10.Since the same variable is calculated two times atXNkk1 and X0k+10, respectively, (33b) can be reformulated as follows:(34a)∑j=0Nk-1DsjkWjk+DsNkkW0k+1-tk-tk-1FXskτ,Uskτ,τ=0,1≤k≤K-1,(34b)∑j=0Nk-1DsjKWjK+DsNkKWNkK-tK-tK-1FXsKτ,UsKτ,τ=0.After being discretized by LG pseudospectral method, the final NLP is obtained which consists of equations (33a), (33c)~(33f), (34a), and (34b). The solution of this NLP is the approximation solution of original optimal control problem in (11a), (11b), (11c), (11d), and (11e). This is easy to be solved by general optimization algorithms.


### 3.4. The KKT Optimality Conditions for Nonlinear Programming

Introducing the Lagrange multiplier vectors, the following augmented performance index can be drawn:(35)J∗Ukτ,tk=θXNkK1+qXNkK1,UNkK1-υT∙NXNkK1-α∙pXNkK1,UNkK1-1-∑k=1K-1∑s=0Nk-1νsk∙∑j=0Nk-1DsjkWjk+DsNkkW0k+1-tk-tk-1FXskτ,Uskτ,τ-∑s=0Nk-1νsK∙∑j=0Nk-1DsjKWjK+DsNkKWNkK-tK-tK-1FXsKτ,UsKτ,τ,where υ∈RN, α∈Rn, and ν∈Rnw denote the Lagrange multipliers associated with discretized original terminal constraints of (33d), discretized aggregated terminal constraints of (33c), and discretized process constraints of (34a) and (34b) in corresponding subintervals, respectively.Take partial derivative of (35) with respect to Xkτ,XNkK1,Ukτ, UNkK1,υ,tf,vk,vK,α, and equating them to zero, then the KKT optimal condition is concluded as (36a), (36b), (36c), (36d), (36e), (36f), and (36g).(36a)∂q∂XNkK1-υT∙∂N∂XNkK1-α∙∂p∂XNkK1+νNk-1KDNk-1NkK+∂θ∂XNkK1=0,(36b)∂q∂UNkK1-υT∙∂N∂UNkK1-α∙∂p∂UNkK1=0,(36c)∑j=0Nk-1Dsjk=tk-tk-1∂Fk∂Xskτ,1≤k≤K-1,0≤s≤Nk-1,(36d)∑j=0Nk-1DsjK=tK-tK-1∂FK∂XsKτ,0≤s≤Nk-1,(36e)∂Fk∂Uskτ=0,1≤k≤K,0≤s≤Nk,(36f)υT∙NXNkK1=0,(36g)α∙pXNkK1,UNkK1-1=0,where Fk≡F(Xskτ,Uskτ,τ), 1≤k≤K.


## 3.1. Constraint Aggregation

The main idea of constraints aggregation is the KS function, which was developed by Kreisselmeier and Steinhauser in 1976 [30]. The mathematical form is(12)KShx,ρ=1ρln⁡∑j=1Nceρhjx,where Nc denotes the number of path constraints hjx. ρ is an approximation parameter.KS function can generate one envelope surface which isC1-consecutive. This can estimate the maximal function of hx conservatively. The properties are as follows [43].(a)
For an arbitraryρ>0, there exist the following formulas:(13)max⁡hx≤KShx,ρ≤max⁡hx+1ρln⁡Nc,limρ→∞⁡KShx,ρ=KShx,ρ.(b)
Ifρ2>ρ1>0, then KShx,ρ1≥KShx,ρ2.(c)
KShx,ρ is a convex function, if and only if it is convex for all constraints hx.From the above properties, KS function is the underestimation for feasible region. The bigger the value ofρ is, the more accurate the estimation for the maximal constraint is. This can be obviously shown in Figure 1. What is more, if the original problem is convex, the problem after being processed by KS function is convex too.Figure 1
Estimate ability of KS function formax⁡g1x=5/log⁡x-0.2x-3,g2x=x2/40+x/5-1.The path constraints in (11a), (11b), (11c), (11d), and (11e) can be transformed into a terminal constraint using KS function; the detailed process is as follows [44].Considering the point constraints on the path constraints in (12), ti∈t0,tf, then(14)KShx,ρ=1ρln⁡∑j=1Nc∑i=1eρhjxti.Transform it into the continuous form by introducing the impulse function,(15)KShx,ρ=1ρln⁡∫t0tf∑j=1Nc∑i=1eρhjxtiδt-tidt,where δt-ti is the impulse function.Define a new statep,(16)p˙=fp=∑j=1Nc∑i=1eρhjxtiδt-ti,pt0=0.Substitute (16) into (15); then (17)KShx,ρ=1ρln⁡ptf-pt0.Sincehxt,ut,t≤0, the constraints processed meet the following condition:(18)KShx,ρ≤0.In consequence, the path constraints can be transformed into the following terminal time constraint:(19)pxtf,utf,tf-1≤0.Moreover, a new state variableq is introduced to the index in (11a), which satisfies(20)q˙=fq=ϕxt,ut,t,qt0=0.Then(21)qxtf,utf,tf=∫t0tfϕxt,ut,tdt.Letω=x,p,qT, F=f,fp,fqT; the original problem can be transformed into the following form:(22a)maxut,tfJut,tf=θxtf,tf+qxtf,utf,tf,(22b)s.t.ω˙t=Fxt,ut,t,(22c)pxtf,utf,tf-1≤0,(22d)Nxtf,tf≤0,(22e)xt0=x0,pt0=qt0=0.


## 3.2. Multistage Problem Formulation

To cope with the discontinuity points of control variables preferably, CVP is adopted to convert the original problem into a MSP [23]. The main idea of this process is to approximate the control with piecewise function according to the slug injection strategy. Discretize the whole time domain t0,tf into K uneven time periods. The specific time control nodes are(23)t0<t1<t2<⋯<tK=tf.Approximate the continuous control variable with piecewise function in every period; that is,(24)ut=u^k,t∈tk-1,tk,k=1,2,…,K.Furthermore, we deal with the time variable by defining a new time variableτ:(25)τ=t-tk-1tk-tk-1,k=1,2,…,K.Then every original periodt∈tk-1,tk is transformed into τ∈0,1.Take derivation of (25) with respect to τ; then,(26)dτ=1tk-tk-1dt,1≤k≤K.Definexkτ,ukτ as the state and control on kth period, respectively. Then, (22a), (22b), (22c), (22d), and (22e) can be reformulated as(27a)maxukτ,tkJukτ,tk=θxK1,1+qxK1,uK1,1,(27b)s.t.ω˙kτ=tk-tk-1Fxkτ,ukτ,τ,(27c)pxK1,uK1,1-1≤0,(27d)NxK1,1≤0,(27e)x10=x0,p10=q10=0,(27f)xk1=xk+10,1≤k≤K-1,where ω=x,p,qT and F=f,fp,fqT. Equation (27f) is the linkage constraint which is used to keep the continuity of state variables at the interface of all subintervals. The MSP in (27a), (27b), (27c), (27d), (27e), and (27f) consists of K subproblems corresponding to K subintervals. The final control is obtained by solving K optimal control problems.


## 3.3. Legendre-Gauss Pseudospectral Discretization

For being solved by well-developed optimization algorithms, Legendre-Gauss (LG) pseudospectral discretization is applied to transform the original optimal control problem into a NLP on every subinterval [26, 28, 42].ConstructNk+1-order Lagrange interpolation polynomials with Nk LG points to approximate the states and controls in the subinterval k∈1,2,…,K; then,(28a)ωkτ≈Wkτ=∑j=0NkWjkLjkτ.(28b)Ljkτ=∏l=0,l≠jNkτ-τlkτjk-τlk,j=0,1,2,…,Nk.(29a)ukτ≈Ukτ=∑i=0NkUikL^ikτ,(29b)L^ikτ=∏l=0,l≠iNkτ-τlkτik-τlk,i=1,2,…,Nk,where τ0k,τ1k,…,τNk-1k denote the LG interpolation points that are defined on every subinterval τk∈tk-1,tk. Wk(τ),Uk(τ) are the approximations of state and control on corresponding subinterval. Note that the interpolation points of state are one more than that of control.Take differential of state with respect toτ; then,(30)W˙kτ=∑j=0NkWjkL˙jkτ.The differential of Lagrange interpolation polynomials at every LG point can be expressed by a differential calculation matrixD∈RNk×Nk+1 which can be described as(31)Dsjk=L˙jkτ=∑c=0Nk∏l=0,l≠j,cNkτs-τlk∏l=0,l≠jNkτjk-τlk,s=0,1,2,…,Nk-1.Then the dynamic constraints in (27b) can be transformed into the following algebraic constraints:(32)∑j=0NkDsjkWjk-tk-tk-1FXskτ,Uskτ,τ=0.Then the optimal control problem can be reformulated into a new NLP:(33a)maxUkτ,tkJUkτ,tk=θXNkK1+qXNkK1,UNkK1,(33b)s.t.∑j=0NkDsjkWjk-tk-tk-1FXskτ,Uskτ,τ=0,s=0,1,2,…,Nk-1,(33c)pXNkK1,UNkK1-1≤0,(33d)NXNkK1≤0,(33e)X010=x0,p010=q010=0,(33f)XNkk1=X0k+10.Since the same variable is calculated two times atXNkk1 and X0k+10, respectively, (33b) can be reformulated as follows:(34a)∑j=0Nk-1DsjkWjk+DsNkkW0k+1-tk-tk-1FXskτ,Uskτ,τ=0,1≤k≤K-1,(34b)∑j=0Nk-1DsjKWjK+DsNkKWNkK-tK-tK-1FXsKτ,UsKτ,τ=0.After being discretized by LG pseudospectral method, the final NLP is obtained which consists of equations (33a), (33c)~(33f), (34a), and (34b). The solution of this NLP is the approximation solution of original optimal control problem in (11a), (11b), (11c), (11d), and (11e). This is easy to be solved by general optimization algorithms.


## 3.4. The KKT Optimality Conditions for Nonlinear Programming

Introducing the Lagrange multiplier vectors, the following augmented performance index can be drawn:(35)J∗Ukτ,tk=θXNkK1+qXNkK1,UNkK1-υT∙NXNkK1-α∙pXNkK1,UNkK1-1-∑k=1K-1∑s=0Nk-1νsk∙∑j=0Nk-1DsjkWjk+DsNkkW0k+1-tk-tk-1FXskτ,Uskτ,τ-∑s=0Nk-1νsK∙∑j=0Nk-1DsjKWjK+DsNkKWNkK-tK-tK-1FXsKτ,UsKτ,τ,where υ∈RN, α∈Rn, and ν∈Rnw denote the Lagrange multipliers associated with discretized original terminal constraints of (33d), discretized aggregated terminal constraints of (33c), and discretized process constraints of (34a) and (34b) in corresponding subintervals, respectively.Take partial derivative of (35) with respect to Xkτ,XNkK1,Ukτ, UNkK1,υ,tf,vk,vK,α, and equating them to zero, then the KKT optimal condition is concluded as (36a), (36b), (36c), (36d), (36e), (36f), and (36g).(36a)∂q∂XNkK1-υT∙∂N∂XNkK1-α∙∂p∂XNkK1+νNk-1KDNk-1NkK+∂θ∂XNkK1=0,(36b)∂q∂UNkK1-υT∙∂N∂UNkK1-α∙∂p∂UNkK1=0,(36c)∑j=0Nk-1Dsjk=tk-tk-1∂Fk∂Xskτ,1≤k≤K-1,0≤s≤Nk-1,(36d)∑j=0Nk-1DsjK=tK-tK-1∂FK∂XsKτ,0≤s≤Nk-1,(36e)∂Fk∂Uskτ=0,1≤k≤K,0≤s≤Nk,(36f)υT∙NXNkK1=0,(36g)α∙pXNkK1,UNkK1-1=0,where Fk≡F(Xskτ,Uskτ,τ), 1≤k≤K.


## 4. The Adaptive Strategies for Solving the Control Problem

Although the original optimal problem has been transformed into a general NLP, there are still many problems that need to be coped with, such as how to adjust the parameterρ in KS function, how to determine the subintervals, and how to regulate the degree of approximating in LG pseudospectral discretization. In view of these situations, two adaptive strategies are applied in constraint aggregation and LG pseudospectral discretization, respectively.


### 4.1. Adaptive Constraint Aggregation Strategy

Whenρ tends to infinity, the value of KS function is the value of maximal constraint in theory. But if ρ is too large, there will be angular points at intersection points of multiple constraints of original problems. The continuity of first-order derivative will be decreased. After being aggregated, the second-order derivative will become morbid at angular points which will lead to slow convergence or even no convergence. This can be well explained in Figure 2. To keep good smoothness of KS function, the value of ρ should be decreased adequately. But this may lead to the feasible region of KS function being less than the real feasible region and the actual optimal control being unobtainable when the actual optimal state is on the boundary. The adaptive strategy in [44] can cope with this well. In this method, the calculation can be carried out with a big value of ρ in the feasible region or near a single constraint, and ρ can be decreased adaptively according to the first-order derivative at intersection points of multiple constraints.Figure 2
The curvature of KS function corresponding tomax⁡g1x=5/log⁡x-0.2x-3,g2x=x2/40+x/5-1.Take partial differential of (12) with respect to ρ; then,(37a)∂KS∂ρ=1ρhmaxx-KS-Mhx,(37b)Mhx=∑j=1Nceρhjxhmaxx-hjx∑j=1Nceρhjx.Here, hmaxx is the maximal constraint among all path constraints.From property (a) of KS function in Section3.1, if all constraints satisfy hjx=hmaxx, then(38)KS=hmaxx+1ρln⁡m.Equation (37a) can be written as(39)KS1′=1ρ1ρln⁡m-Mhx.If the current state is near the intersection point ofk constraints, (37a) can be similarly written as(40)KSu′=-1ρMhx-1ρ2ln⁡k.Likewise, if the current state approximates only one constraint, (37a) can be transformed into another form:(41)KSd′=-1ρMhx.Defineκx,ρ=KSu′-KSd′/KS1′-KSd′, κx,ρ∈0,1. It reflects the relation between the derivative of KS function and the number of constraints. When in use, we can give one ρ∗ as the fiducial value and calculate κ as the distance parameter between the current state and the intersection point of multiple constraints at each iteration. Then the following adaptive constraints aggregation strategy can be defined:(42)ρ=ρ_+Δρ1-κx,ρ∗,where ρ_ is the lower bound of ρ and Δρ>0 is the regulating range. The choice of ρ∗ has to ensure that enough regulating range can be obtained near the intersection point of multiple constraints.


### 4.2. Adaptive LG Discretization Strategy

In LG pseudospectral method, enough interpolation points have to be chosen to meet given accuracy, but this will lead to the burdensome calculation. To regulate the subintervals reasonably, a mesh refinement algorithm that refers to [28, 42] is adopted in this section.For subintervalk, select Nk LG points (τ⌢1k,τ⌢2k,…,τ⌢Bkk). Let Bk=Nk+1, τ⌢1k=τ1k=tk-1, and τ⌢Bkk=tk. According to (27b) and (28a), if the differential of state W⌢k accords with the dynamics at every LG point τ⌢ik, we have(43)W⌢kτik=Wkτk-1+tk-tk-1∑l=1BkI⌢jlkFXkτ⌢lk,Ukτ⌢lk,τ⌢lk,where I⌢jlk,j,l=1,2,…,Bk, denotes the Bk×Bk LG integration matrix which is similar to (30).Define the absolute error and relative error as follows:(44)Eikτ⌢lk=W⌢ikτ⌢lk-Wikτ⌢lk.(45)eikτ⌢lk=Eikτ⌢lk1+maxj∈1,2,…,Bk⁡Wikτ⌢jk,where l=1,2,…,Bk; i=1,2,…,nx+2.Then the maximal relative error can be written as(46)emaxk=maxi∈1,2,…,nx+2,l∈1,2,…,Bk⁡eikτ⌢lk.Assume that the given accuracy tolerance isε. For subinterval k, if emaxk≤ε, we can conclude that the current LG points are reasonable. If emaxk>ε, we can adopt the following equation to update the LG points [28]:(47)Nk′=Nk+ceillogNk⁡emaxkε.Here ceil is the operator that rounds to the next highest integer.Assume thatNmin,Nmax are the lower bound and upper bound of the number of LG points in subinterval k. If the regulated number Nk′ satisfies Nk′≤Nmax, then Nk′ is adopted to update the number of LG points. On the other hand, if Nk′>Nmax which means the LG points are too many, then the current subintervals should be further refined. We use K′ to update the number of subintervals [42]:(48)K′=maxceilNk′Nmin,2.From the above, the adaptive LG discretization strategy has been achieved. It balances the accuracy and computation burden by updating the number of LG interpolation points and refining the subintervals.


### 4.3. Control Method

In this section, the SQP [35, 36] is adopted to solve the NLP that is described as (33a), (33b), (33c), (33d), (33e), (33f), (34a), and (34b). The specific procedures are as follows:(1)
Initialize the parametersρ∗,ρ_,Δρ,Nk,Nmax,Nmin, the number of subintervals K, the initial state W010=0, initial control U010=0, the accuracy tolerance εJ,εs, and the current iteration r=1.(2)
Transcribe to NLP of (33a), (33b), (33c), (33d), (33e), (33f), (34a), and (34b) with current subintervals.(3)
Calculate the value ofρ in step r according to (42) on the basis of the partial differentials of KS function.(4)
Calculate the states and performance indexJk of each subinterval at each LG point in step r according to (35) and (43).(5)
Estimate the erroremaxk for each subinterval.(6)
Assess the number of LG points: ifemaxk≤εs, go to step (7); if not, update the number of LG points as (47).(7)
Assess the refinement of subintervals: ifNk′≤Nmax, go to step (8); if not, update the subintervals as (48).(8)
Determine the search direction and search step length, solve the problem in (35) and (36a), (36b), (36c), (36d), (36e), (36f), and (36g), and calculate the optimal control Uk+1 and index Jk+1.(9)
IfJk+1-Jk<εJ, finish the computing process and output the optimal control Uk+1 and index Jk+1. If not, go to step (2) and let k=1; carry on the procedure until Jk+1-Jk satisfies the given accuracy.


## 4.1. Adaptive Constraint Aggregation Strategy

Whenρ tends to infinity, the value of KS function is the value of maximal constraint in theory. But if ρ is too large, there will be angular points at intersection points of multiple constraints of original problems. The continuity of first-order derivative will be decreased. After being aggregated, the second-order derivative will become morbid at angular points which will lead to slow convergence or even no convergence. This can be well explained in Figure 2. To keep good smoothness of KS function, the value of ρ should be decreased adequately. But this may lead to the feasible region of KS function being less than the real feasible region and the actual optimal control being unobtainable when the actual optimal state is on the boundary. The adaptive strategy in [44] can cope with this well. In this method, the calculation can be carried out with a big value of ρ in the feasible region or near a single constraint, and ρ can be decreased adaptively according to the first-order derivative at intersection points of multiple constraints.Figure 2
The curvature of KS function corresponding tomax⁡g1x=5/log⁡x-0.2x-3,g2x=x2/40+x/5-1.Take partial differential of (12) with respect to ρ; then,(37a)∂KS∂ρ=1ρhmaxx-KS-Mhx,(37b)Mhx=∑j=1Nceρhjxhmaxx-hjx∑j=1Nceρhjx.Here, hmaxx is the maximal constraint among all path constraints.From property (a) of KS function in Section3.1, if all constraints satisfy hjx=hmaxx, then(38)KS=hmaxx+1ρln⁡m.Equation (37a) can be written as(39)KS1′=1ρ1ρln⁡m-Mhx.If the current state is near the intersection point ofk constraints, (37a) can be similarly written as(40)KSu′=-1ρMhx-1ρ2ln⁡k.Likewise, if the current state approximates only one constraint, (37a) can be transformed into another form:(41)KSd′=-1ρMhx.Defineκx,ρ=KSu′-KSd′/KS1′-KSd′, κx,ρ∈0,1. It reflects the relation between the derivative of KS function and the number of constraints. When in use, we can give one ρ∗ as the fiducial value and calculate κ as the distance parameter between the current state and the intersection point of multiple constraints at each iteration. Then the following adaptive constraints aggregation strategy can be defined:(42)ρ=ρ_+Δρ1-κx,ρ∗,where ρ_ is the lower bound of ρ and Δρ>0 is the regulating range. The choice of ρ∗ has to ensure that enough regulating range can be obtained near the intersection point of multiple constraints.


## 4.2. Adaptive LG Discretization Strategy

In LG pseudospectral method, enough interpolation points have to be chosen to meet given accuracy, but this will lead to the burdensome calculation. To regulate the subintervals reasonably, a mesh refinement algorithm that refers to [28, 42] is adopted in this section.For subintervalk, select Nk LG points (τ⌢1k,τ⌢2k,…,τ⌢Bkk). Let Bk=Nk+1, τ⌢1k=τ1k=tk-1, and τ⌢Bkk=tk. According to (27b) and (28a), if the differential of state W⌢k accords with the dynamics at every LG point τ⌢ik, we have(43)W⌢kτik=Wkτk-1+tk-tk-1∑l=1BkI⌢jlkFXkτ⌢lk,Ukτ⌢lk,τ⌢lk,where I⌢jlk,j,l=1,2,…,Bk, denotes the Bk×Bk LG integration matrix which is similar to (30).Define the absolute error and relative error as follows:(44)Eikτ⌢lk=W⌢ikτ⌢lk-Wikτ⌢lk.(45)eikτ⌢lk=Eikτ⌢lk1+maxj∈1,2,…,Bk⁡Wikτ⌢jk,where l=1,2,…,Bk; i=1,2,…,nx+2.Then the maximal relative error can be written as(46)emaxk=maxi∈1,2,…,nx+2,l∈1,2,…,Bk⁡eikτ⌢lk.Assume that the given accuracy tolerance isε. For subinterval k, if emaxk≤ε, we can conclude that the current LG points are reasonable. If emaxk>ε, we can adopt the following equation to update the LG points [28]:(47)Nk′=Nk+ceillogNk⁡emaxkε.Here ceil is the operator that rounds to the next highest integer.Assume thatNmin,Nmax are the lower bound and upper bound of the number of LG points in subinterval k. If the regulated number Nk′ satisfies Nk′≤Nmax, then Nk′ is adopted to update the number of LG points. On the other hand, if Nk′>Nmax which means the LG points are too many, then the current subintervals should be further refined. We use K′ to update the number of subintervals [42]:(48)K′=maxceilNk′Nmin,2.From the above, the adaptive LG discretization strategy has been achieved. It balances the accuracy and computation burden by updating the number of LG interpolation points and refining the subintervals.


## 4.3. Control Method

In this section, the SQP [35, 36] is adopted to solve the NLP that is described as (33a), (33b), (33c), (33d), (33e), (33f), (34a), and (34b). The specific procedures are as follows:(1)
Initialize the parametersρ∗,ρ_,Δρ,Nk,Nmax,Nmin, the number of subintervals K, the initial state W010=0, initial control U010=0, the accuracy tolerance εJ,εs, and the current iteration r=1.(2)
Transcribe to NLP of (33a), (33b), (33c), (33d), (33e), (33f), (34a), and (34b) with current subintervals.(3)
Calculate the value ofρ in step r according to (42) on the basis of the partial differentials of KS function.(4)
Calculate the states and performance indexJk of each subinterval at each LG point in step r according to (35) and (43).(5)
Estimate the erroremaxk for each subinterval.(6)
Assess the number of LG points: ifemaxk≤εs, go to step (7); if not, update the number of LG points as (47).(7)
Assess the refinement of subintervals: ifNk′≤Nmax, go to step (8); if not, update the subintervals as (48).(8)
Determine the search direction and search step length, solve the problem in (35) and (36a), (36b), (36c), (36d), (36e), (36f), and (36g), and calculate the optimal control Uk+1 and index Jk+1.(9)
IfJk+1-Jk<εJ, finish the computing process and output the optimal control Uk+1 and index Jk+1. If not, go to step (2) and let k=1; carry on the procedure until Jk+1-Jk satisfies the given accuracy.


## 5. Numerical Solution of Optimal Control for ASP Flooding

For a real ASP flooding optimization problem in Section2.1, carry out a series of operations—constraints aggregation, multistage problem formulation, Legendre-Gauss pseudospectral discretization, and adaptive strategies—and get the KKT optimality conditions. Then SQP is applied to solve the KKT conditions as in Section 4.3 to get the optimal injection strategy. After dimensionless processing, the optimal injection strategy for ASP flooding is solved by the proposed method in this paper. To explain the result better, the experiment optimization method and CVP method are adopted to solve the same problem and to be compared with the proposed method. All the simulations are carried out on software Matlab R2011b.In this section, three-slug injection is used in simulation. There are twelve optimization variables which areuΘ1,uΘ2,uΘ3,TΘ. Set tf=2, tP=0.5, uamax=2, usmax=0.5, upmax=1.5, Mp=0.6, Ma=0.6, and Ms=0.2. Give lots of different control inputs, and compute the performance index. Through the comparison of experiment results many times, we select the best injection concentration and slug size with the maximum profit. We call these results obtained with this method as experiment optimized results (EOS). The injection strategy is ua0=(1.2,1,0.8), us0=(0.5,0.4,0.2), and up0=(1,0.6,0.4), the slug size is T0=(0.2,0.2,0.1), and the performance index is J0=0.2538.Solve this optimization problem with CVP method and the method proposed in this paper, respectively. The parameters are set as follows:(49)ρ_=200,Δρ=100,ρ∗=50,Nmin=4,Nmax=25,εs=10-3,εJ=10-4,K=3.Some essential data for solving the optimization model is shown in Table1.Table 1
Partial data for solving ASP flooding model.
Sign Value Sign Value Sign Value M 8.3 N 8.33 L 1 Sor0 0.35 S o r 0.02 ϕ 0.4728 a p 1 7 a p 2 4 a p 3 40 μ w 0 0.5 μ w m a x 18 μ o 20 v w 0.0045 ρ r 2.0 ρs 1.0 Note. Table 1 is reproduced from Ge et al. (2016), (under the Creative Commons Attribution License/public domain). The specific variable declaration and unit can be found in reference [6].After being optimized, the optimal injection strategy with CVP isua∗=(1.862,0.982,0.826), us∗=(0.461,0.350,0.403), up∗=(1.415,1.171,1.024), and T0=(0.22,0.21,0.07), and the dimensionless profit is J∗=0.26385. The optimal injection strategy with the proposed method is ua∗=(1.869,0.981,0.829), us∗=(0.463,0.35,0.402), up∗=(1.418,1.168,1.024), and T0=(0.23,0.19,0.08), and the dimensionless profit is J∗=0.26467. The results are shown in Figures 3–6.Figure 3
Comparison results of injection concentration for alkali.Figure 4
Comparison results of injection concentration for polymer.Figure 5
Comparison results of injection concentration for surfactant.Figure 6
Comparison results of performance index.In these figures, the blue chain dotted line, the red imaginary line, and the black dotted line denote the results for EOS, CVP, and proposed method, respectively. Particularly, to distinguish the results in Figure6, the black line with circle is used to denote the result of proposed method. According to these figures, the result of the proposed method in this paper is very close to that of CVP method, which are both better than EOS method. Regardless of the optimized injection concentration for three displacing agents or the optimized performance index, the values are nearly the same. The performance index optimized by the proposed method is increased by 0.01087 compared with that of EOS. This can fully prove the effectiveness of the proposed method. In addition, the switch time points and the running time are displayed in Table 2. From Figures 3–5, we can find that the LG points are very intensive at the discontinuous points. In view of the switch time points in Table 2, the disposal of the proposed method is much better than CVP method and EOS. The index is further enhanced by 0.00082. This is because of the introduction of adaptive strategies and multistage problem formulation, which can adjust the discontinuous point more precisely. What is more, the running time of the proposed method is less than the other two methods. What is noteworthy is that because of the introduction of approximation, the results may have tiny error compared with the theoretical value.Table 2
Comparison of three methods (proposed method, CVP, and EOS).
Switch time points T 1 T 2 T 3 Performance index CPU times(s) Proposed method 0.23 0.42 0.5 0.26467 89.57 CVP method 0.22 0.43 0.5 0.26385 126.36 EOS method 0.2 0.4 0.5 0.2538 Note. For the need of production, the switch time points for alkali, surfactant and polymer have to be kept the same.To discuss the performance of the proposed method further, ASP flooding problem is optimized by different accuracy tolerances. The corresponding results are shown in Table3.Table 3
Comparison of the proposed method on accuracy and speed under different accuracy tolerancesε.
ε Number of LG points Number of iterations e m a x k Performance index CPU times (s) 1 0 - 3 20 15 4.78×10-4 0.26467 89.57 1 0 - 4 23 16 5.19 × 1 0 - 5 0.26467 99.32 1 0 - 5 29 19 2.97 × 1 0 - 6 0.26468 127.11 1 0 - 6 38 23 7.08 × 1 0 - 7 0.26468 184.59 Note. Keep Nmin=4 for all calculations.From Table3, the LG points, iterations, and running times increase evidently with the accuracy tolerance refining. And the maximum of error in every step increases at the same time. This is because with the increasing of accuracy more LG points are needed to be introduced to realize the approximation of states and controls and the disposal of discontinuous points, and more iterations are needed to search for optima. But the performance index is almost invariant. So we can conclude that the result is the optimal one. Though a better accuracy is desired, this will lead to an increase in the computation burden. So, there is a need a balance accuracy and running time in application.


## 6. Conclusions

In this paper, a method based on constraint aggregation and LG pseudospectral method with adaptive strategies has been proposed to solve the optimal control for ASP flooding whose control variables are discontinuous. Firstly, we extended this problem to a class of optimal control problems. Then we introduced the constraint aggregation to convert all path constraints into one terminal condition and discretized the system equations with LG pseudospectral method. For simplicity, we also developed a new time variable and transformed the original problem into a multistage problem with CVP. Furthermore, we used two adaptive strategies to decrease the approximation error with regulatingρ of KS function, the LG points, and the refinement of subintervals adaptively. Finally, the optimal control problem for ASP flooding was solved by the proposed method. To demonstrate the performance better, we also solved this problem with CVP and experiment optimized result method. Numerical simulations showed that the proposed method can complete this well with less time and no loss of accuracy. The multiobjective optimization problem (the maximal oil production, minimal cost, etc.) and uncertainty problem (the uncertainty of geological parameters, oil price, external disturbance, etc.) for ASP flooding are very important to enhance the oil recovery and make the optimal injection strategy. These are the next research direction in the future.



---

*Source: 1040476-2018-05-31.xml*