A Lagrange Relaxation Based Decomposition Algorithm for Large-Scale Offshore Oil Production Planning Optimization

Abstract

In this paper, a new Lagrange relaxation based decomposition algorithm for the integrated offshore oil production planning optimization is presented. In our previous study (Gao et al. Computers and Chemical Engineering, 2020, 133, 106674), a multiperiod mixed-integer nonlinear programming (MINLP) model considering both well operation and flow assurance simultaneously had been proposed. However, due to the large-scale nature of the problem, i.e., too many oil wells and long planning time cycle, the optimization problem makes it difficult to get a satisfactory solution in a reasonable time. As an effective method, Lagrange relaxation based decomposition algorithms can provide more compact bounds and thus result in a smaller duality gap. Specifically, Lagrange multiplier is introduced to relax coupling constraints of multi-batch units and thus some moderate scale sub-problems result. Moreover, dual problem is constructed for iteration. As a result, the original integrated large-scale model is decomposed into several single-batch subproblems and solved simultaneously by commercial solvers. Computational results show that the proposed method can reduce the solving time up to 43% or even more. Meanwhile, the planning results are close to those obtained by the original model. Moreover, the larger the problem size, the better the proposed LR algorithm is than the original model.

1. Introduction

Oil and gas resources are the blood of national development. Oil and gas production plays an extremely important role in promoting industrial development and social progress. In recent years, China′s dependence on foreign crude oil continues to increase, posing a serious threat to national energy security [1]. In order to alleviate the shortage of oil and gas resources in China, on the one hand, China has actively taken measures to save energy and reduce emissions and put forward the goals of green, low-carbon and economic development. On the other hand, the survey report of the United States Geological Survey (USGS) [2] indicates that global deep-water oil resources have great exploration and development potential. Exploration shows that China′s deep-water areas are rich in oil and gas resources [3]. The vast deep-water area of the South China Sea has the basic geological conditions for the formation of large and medium-sized oil and gas fields and good oil and gas exploration prospects and resource potential [4]. In the last 10 years, most of the newly discovered oil and gas fields in China have been located at sea, and the proportion of deep-sea oil and gas fields keeps increasing [5]. The rapid development of offshore engineering equipment has greatly enhanced the ability of deep-water oil and gas exploration and development. Policies have supported and guided the process of deep-water oil and gas exploration and development. However, persistently low oil prices have had a devastating impact on offshore oil development [6]. Reducing production costs is a priority in response to falling oil prices. Planning and scheduling optimization is an attractive alternative to reducing production costs.

Due to the high complexity and lack of connection between the above-water platform and underwater oil production, the research on modeling and planning optimization of the whole process hybrid system combining above-water and underwater is rare, or the specific research is only carried out for a certain oil production link, without considering all links together, such as injection and flooding, artificial lifting and crude oil storage. So, it is difficult to optimize and balance the offshore oil and gas production process, all resources and equipment. Moreover, the harsh underwater environment of deep-sea oil production will bring a series of pipeline transportation safety problems, so the flow safety guarantee is closely related to the whole production process, but there are few studies that consider both process optimization and flow safety guarantee. Based on this, in order to maximize the comprehensive economic benefits and ensure the smooth progress of the oil production process, it is necessary to carry out research on the flow assurance guarantee and plan optimization of offshore oil and gas production.

The early research fields of production planning mainly focus on batch production process and petroleum refining process in petrochemical enterprises [7,8,9,10,11]. Gupta et al. [12] set up a multi-period MINLP model for the development of offshore oil and gas, considering the three-phase mixing ratio of oil, gas and water, the installation and expansion of oil tankers, the connection between pipelines and FPSO and drilling and production speed. Vijay et al. [13] built a MINLP model aiming at the problem of offshore oil and gas infrastructure construction and production equipment sharing. Aseeri et al. [14] studied offshore oil and gas production optimization problem under the financial risk management strategy, considering the sequence and productivity of platform construction, drilling and oil well production and budget constraints. Ortiz-Gomez [15] studied different complex oil well production plans in the same reservoir with the fixed topology and described three mixed integer multiperiod optimization models to solve the oil production profiles and operation/shut of oil wells in each time period. At the level of production planning, Lu [16] comprehensively considered manufacturing resource planning (MRPII), customers′ orders and constraint theory to make production planning, so as to carry out effective production scheduling control. Saravanan [17] studied how to improve productivity in small-scale industries by implementing lean tools to reduce delivery times. At present, there are few studies on the planning optimization of oil and gas production.

In the previously published article, a comprehensive planning issue considering both oil well operation and flow assurance was given. In particular, a multi-period mixed integer nonlinear programming (MINLP) model was proposed to minimize the total operation cost, considering well production state, polymer flooding, energy consumption, platform inventory and flow assurance. By solving this integrated model, each well’s working state, flow rates and chemicals injection rates can be optimally determined [18]. Compared with individual optimization, integrated optimization can save more energy and greatly reduce costs. When confronted with the industrial cases, a large-scale model with plenty of variables and lots of constraints is obtained, finally resulting in solving difficulty. Lots of efforts and attempts have been done to overcome the solving difficulty. For now, the known solving strategies of MINLP problems are global optimization algorithm [19], branch and bound algorithm [20,21], external approximation algorithm [22], particle swarm optimization algorithm [23], decomposition algorithm [24,25,26] and so on. Despite the valuable progress in solving algorithms, the direct use of these algorithms in our MINLP model is not applicable because of the unacceptable solving time. To improve the solving efficiency, it is necessary to make optimization decisions for the real-world offshore oil and gas production case in a reasonable time.

After the Lagrange multiplier-based coupling constraints relaxation algorithm was first introduced in the 1970s, there have been many variations and successful applications in process system engineering community. Held [27] used this algorithm to solve the traveling salesman problem. Then Fisher et al. [28] applied Lagrange relaxation (LR) algorithms to solve the job shop scheduling problem. Although the optimality was theoretically guaranteed, it was proved not to be effective for practical problems. Until the 1990s, Luh et al. [29] applied LR algorithm to parallel machine scheduling. The mechanism of fast estimating rescheduling cost by solving Lagrange multipliers was discussed in detail. This was the first attempt in which LR algorithm was used in a practical scheduling issue. From then, studies have focused on improving the efficiency of the algorithm. Luh [30] and Yu [31] proposed an augmented LR algorithm to improve the lower bound quality and reduce the duality gap by introducing the penalty of the constraint destruction degree. Susara and Grossmann et al. [32] proposed a heuristic algorithm based on LR to solve the long-term design and planning of offshore hydrocarbon field infrastructures with complex economic objectives. They proposed an iterative format to speed up the solution of the problem, providing new ideas for subsequent research. Campanagara [33] presented a predictive control approach that obtains a high thermal comfort level optimizing the use of an HVAC (Heating, Ventilation and Air Conditioning) system by means of a cost function. Meanwhile, they used Lagrangian dual method to optimize the procedure, which allows the use of parallel programming paradigms in an easy way. At present, the LR algorithm has been applied to various industries, such as mining, high-speed railway, power industry, water conservancy, etc. To the best of our knowledge, there is little report about LR in oil and gas production optimization, especially in the offshore oil and gas domain. Motivated by the wide and successful application of LR, a batch decomposition method based on LR is proposed in this paper to accelerate the solving speed. The algorithm relaxes the multi-batch supply constraint, which is a coupling constraint in the original problem, to obtain a relaxation of the original problem. This relaxation problem is decomposed into several independent single-batch sub-problems to solve the difficult problem of long solution times for large-scale problems. The case study shows that the proposed algorithm can reduce the solution time to 43% (or even more) and the optimality gap is also reduced to 0.72%.

In our previous study, the research mainly focuses on the integrated planning optimization model considering the well operation and flow assurance. Following the previous studies, in order to cope with the difficult situation of solving the multi-cycle full process integrated MINLP model, this paper further investigates the batch decomposition strategy based on Lagrangian relaxation algorithm to improve the solving efficiency. The rest is organized as follows. First, Section 2 gives the problem statement and process description. Based on process analysis, Section 3 gives a detailed solution algorithm. Section 4 takes the actual production process as an example to verify the feasibility of the proposed Lagrange batch decomposition algorithm, and results directly solved by the solver are compared. Finally, a conclusion is drawn in Section 5.

2. Process Description and Problem Statement

2.1. Process Description

From the wells to the platform, the whole production process can generally be divided into three parts: the under-well reservoir process, the under-water production process and the over-water platform section (Figure 1).

Wells are widely distributed in oilfield. Typically, an oilfield contains a lot of reservoirs, each of which in turn contains many wells. Wells can be divided into different batches of oil wells by the wells’ geological properties and physical characteristics, as shown in Figure 2. Wells belonging to the same batch interconnect with each other through a complex comprehensive pipeline network to convey liquid to manifold (i.e., PLEM in Figure 2). Wells in the same reservoir are grouped into one batch and normally share a surface equipment, usually named floating production, storage and offloading unit (FPSO). The whole oil and gas production process include five parts: optimization of oil well start-up and shutdown operation, optimization of electric submersible pump energy consumption, optimization of FPSO storage, optimization of polymer injection and flooding and flow safety guarantee. The key and difficult points of system integration modeling are:

• Optimize the capacity consumption of electric submersible pump (ESP) under changing flow conditions;
• Use the flow guarantee mechanism to balance the optimal operation scheme of oil wells and ensure the flow safety;
• Allocate the injection strategy of each well with a given amount of polymer;
• Integrate well and platform operations to separate or store oil/gas delivered to the platform.

Considered in terms of practical industrial problems, the following issues need to be considered in planning modelling:

• To facilitate modeling, the entire offshore oil and gas production process is regarded as a continuous production process, and all production-related variables can be connected through time;
• The whole oilfield is divided into several blocks according to geographical location, product characteristics and other conditions, and the modeling is optimized according to the blocks;
• To ensure smooth production, start-up and shutdown operation of each underwater well shall be considered;
• In order to ensure safe production, considering the protective effect of flow assurance guarantee on the production process, the cost of single wax removal is considered.

Taking the process flow diagram of offshore oil and gas production process shown in Figure 2 as the research object, the system integration modeling is carried out in consideration of equipment optimization at the operational layer and flow safety assurance technology. Several assumptions are made in this study as follows:

• The production wells are separated and totally independent of each other. It is natural because each well has its own independent reservoir.
• During the middle and later periods of oilfield development, artificial lift technology and polymer flooding is indispensable;
• All the electric submersible pumps have the same working characteristic curve;
• Geological properties characterizing the well are available;
• In the absence of polymerization flooding, oil recovery rate remains the lowest;
• The location of easily blocked pipeline section is known;
• Ignore the pressure change in the pipe.

With the above assumptions, the model relies on the following given information:

• A planning horizon and planning period;
• Production tasks for each batch of oil wells along the planning horizon;
• Working load range of oil production wells;
• A set of storage bins, their minimum and maximum stock and initial inventories;
• The penalty of switching operations and stock out;
• A set of cost coefficient and model parameters.

The decision variables are:

• The production rate and operating state of each oil well in each time period;
• The detailed delivery quantity in each oil batch in each time period;
• The injection displacement volume of each well;
• Diesel fuel consumption within each planning period;
• The wax removal cycle of each well.

The integrated planning model defined as a multi-period MINLP has been developed in our previous report. To improve readability, the original model is detailed as follows.

Mathematically, the objective function is given as follows:

$$\min Z=Z_1+Z_2+Z_3+Z_4+Z_5+Z_6$$

(1)

The objective described in Equation (1) aims at minimizing the overall cost (Z), which includes the oil well open-close switching penalty(Z₁), energy consumption (Z₂), oil inventory (Z₃), and chemicals cost (Z₄), wax removal cost (Z₅) and the costs of stock out penalty (Z₆). Constraints related to Z₁ to Z₆ are shown in Equations (2)–(48). The specific meaning of each equation has been specified in our published articles ¹⁸.

$$Z_1={\displaystyle \sum }_i{\displaystyle \sum }_t{\alpha }_{i}w{f}_{i,t}$$

(2)

$$w{f}_{i,t}+w_{i,t}\ge w_{i,t+1}\hspace{1em}\forall i\in I,t\in T$$

(3)

$$w{f}_{i,t}+w_{i,t+1}\ge w_{i,t}\hspace{1em}\forall i\in I,t\in T$$

(4)

$$p_i^{end}=p_i^{in}-\frac{141.2x_{i,t}B\mu }{kh}\left(\frac{1}{2}\left[ln\frac{0.000246kt}{\Phi \mu c_i{r}_i^2}+0.80907\right]\right)\hspace{1em}\forall i\in I,t\in T$$

(5)

$$p_i^{end}=p_i^{in}-c_1{x}_{i,t}\left(lnDr+c_2\right)\hspace{1em}\forall i\in I,t\in T$$

(6)

$$X{P}_{i,t}={\pi }_i{x}_{i,t}{a}_{i0}\left(a_{i1}+\ln Dr\right)\hspace{1em}\forall i\in I,t\in T$$

(7)

$$p_{i,t}^{end}=p_{i,t}^{in}-X{P}_{i,t}\hspace{1em}\forall i\in I,t\in T$$

(8)

$$p_{i,t}^{in}-X{P}_{i,t}\ge p_i^{low}\hspace{1em}\forall i\in I,t\in T$$

(9)

$$S{P}_{i,t}=b_{i0}\left(b_{i1}+\ln Dr\right)\left(1-w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(10)

$$p_{i,t}^{end}-p_{i,t}^{in}+X{P}_{i,t}\ge -M\left(1-w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(11)

$$p_{i,t}^{end}-p_{i,t}^{in}+X{P}_{i,t}\le M\left(1-w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(12)

$$p_{i,t}^{in}-p_i^{low}-X{P}_{i,t}\ge -M\left(1-w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(13)

$$p_{i,t}^{end}-p_{i,t}^{in}-S{P}_{i,t}\ge -M\left(1-Y_{i,t}+w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(14)

$$p_{i,t}^{end}-p_{i,t}^{in}-S{P}_{i,t}\le M\left(1-Y_{i,t}+w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(15)

$$p_{i,t}^{in}-p_i^{up}+S{P}_{i,t}\le M\left(1-Y_{i,t}+w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(16)

$$p_{i,t}^{end}-p_i^{up}\ge -M\left(Y_{i,t}+w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(17)

$$p_{i,t}^{end}-p_i^{up}\le M\left(Y_{i,t}+w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(18)

$$p_{i,t}^{in}-p_i^{up}+S{P}_{i,t}\ge -M\left(Y_{i,t}+w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(19)

$$p_{i,t}^{in}=p_{i,t-1}^{end}\hspace{1em}\forall i\in I,t\in T$$

(20)

$$p_{i,1}^{in}=p_i^{initial}\hspace{1em}\forall i\in I$$

(21)

$$el{e}_{cost}={\displaystyle \sum }_i{\displaystyle \sum }_t{\beta }_i\left(x_{i,t},w_{i,t}\right)\hspace{1em}\forall i\in I,t\in T$$

(22)

$$w_{i,t}{x}_i^{min}\le x_{i,t}\le w_{i,t}{x}_i^{max}\hspace{1em}\forall i\in I,t\in T$$

(23)

$$Z_2=p{e}_1 · el{e}_{cost}$$

(24)

$$I_{k,t}=I_{k,t-1}+{\displaystyle \sum }_{i\in K}{x}_{i,t}-p{r}_{k,t}\hspace{1em}\forall k\in K,t\in T$$

(25)

$$I_{k,1}=I_k^{initial}\hspace{1em}\forall k\in K$$

(26)

$$I^{min}\le I_{k,t}\le I^{max}\hspace{1em}\forall k\in K,t\in T$$

(27)

$$Z_3={\displaystyle \sum }_k{\displaystyle \sum }_t\gamma ·I_{k,t}\hspace{1em}\forall k\in K,t\in T$$

(28)

$$\Delta E_{i,t}=\raisebox{1ex}{$w_{i,t}\left(x_{i,t}-x_i^{min}\right)$}\!\left/ \!\raisebox{-1ex}{$x_i^{min}$}\right.\hspace{1em}\forall i\in I,t\in T$$

(29)

$$\log P_{i,t}=A_i+B_i\Delta E_{i,t}\hspace{1em}\forall i\in I,t\in T$$

(30)

$${\mathrm{Z}}_4={\displaystyle \sum }_i{\displaystyle \sum }_t\delta ·P_{i,t}\hspace{1em}\forall i\in I,t\in T$$

(31)

$$Q_{in}-Q_{out}-Q_r=Q_{acc}$$

(32)

$$Q_{in}=\rho C_{p}vA{T}_L\Delta t$$

(33)

$$Q_{out}=\rho C_{p}vA{T}_{L+\Delta \mathrm{L}}\Delta t$$

(34)

$$Q_{acc}=\rho C_{p}A\Delta L\Delta Te$$

(35)

$$Q_r=\frac{2\pi r{k}_1\Delta L\Delta t\left(T{e}_{k,t}-T_{out}\right)}{R_t}\hspace{1em}\forall k\in K,t\in T$$

(36)

$$R_t=\frac{1}{h_{in}r}+\frac{1}{{\lambda }_{ins}}ln\frac{r+s+s_{tub}}{r+s_{tub}}$$

(37)

$$w_{i,t}T{e}_k^{\min }\le T{e}_{k,t}\le w_{i,t}T{e}_k^{\max }\hspace{1em}\forall k\in K,t\in T,\forall i\in I$$

(38)

$$Z_5=floor\left(\frac{TT}{T{l}_k}\right)·\tau$$

(39)

$$M{l}_i=2F_d{\displaystyle \sum }_k^K{x}_i\frac{e_k{}^2+D_k{e}_k{}^2}{D_k{}^2+2e_k{}^2+2D_k{e}_k}\hspace{1em}\forall i\in I$$

(40)

$$V{l}_k=\raisebox{1ex}{$M{l}_k$}\!\left/ \!\raisebox{-1ex}{$Gl$}\right.\hspace{1em}\forall k\in K$$

(41)

$$v_k=\frac{D_k-\sqrt{D_k{}^2-\frac{4V{l}_k}{\pi L_k}}}{2}\hspace{1em}\forall k\in K$$

(42)

$$\mathrm{T}{l}_k=\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2v_k$}\right.\hspace{1em}\forall k\in K$$

(43)

$$0<h\le h^{max}$$

(44)

$$Z_6=\sum \theta ·\left(d_t-p{r}_t\right)$$

(45)

$$d_t={\displaystyle \sum }_k{d}_{k,t}\hspace{1em}\forall t\in T$$

(46)

$$p{r}_t={\displaystyle \sum }_{k}p{r}_{k,t}\hspace{1em}\forall t\in T$$

(47)

$$p{r}_{k,t}\le d_{k,t}\hspace{1em}\forall t\in T,k\in K$$

(48)

2.2. Problem Statement

The MINLP model is theoretically difficult to solve especially for the industrial case. In addition, long-period planning time horizon and a large amount of parallel equipment will result in a large number of continuous and discrete variables, leading to the difficulty of low solving efficiency. In order to further improve the solving efficiency and meet the actual industrial production requirements, a batch decomposition algorithm based on Lagrange relaxation (LR) is proposed in this paper to solve the MINLP model. Problems that can use LR algorithm generally have the following characteristics: the objective function is differentiable with respect to the decision variables and constraints and there is a strong coupling between the constraints. The idea of LR [34] is to introduce Lagrange multipliers, which combine coupling constraints with multipliers as a form of relaxation into the original objective function as a penalty term. The original problem is transformed into a relaxation problem without the coupling constraint, thus decomposing the relaxation problem into sub-problems with multiple decision variables, fewer constraints and that is easier to solve. The batch decomposition algorithm based on Lagrange relaxation proposed in this paper decomposes the simultaneous optimization of multi-batch well packs into the optimization of single-batch well packs. Lagrange multipliers are introduced to relax the coupling constraints of multi-batch well packs, and dual problems are constructed to form cyclic iterations. Lagrange multiplier is introduced to relax the multi-batch coupling constraint to the objective function, the single-batch subproblem is obtained, and the solution is taken as the upper bound of the original problem. Secondly, the Lagrange dual problem is constructed, and the solution of the dual problem is taken as the lower bound of the original problem. The quality of the solution is judged by the duality gap between the upper and lower bounds, and the iterative algorithm is realized by continuously updating the multiplier by the step-by-step gradient method, which solves the large-scale problem caused by multi-batch coupling and simultaneous solution and greatly accelerates the solution speed.

2.3. LR Algorithm Implementation

The main steps of solving the planning problem are as follows: The multi-batch supply constraints as coupling constraints, i.e., Equation (48), are relaxed and embedded as a penalty term multiplied by Lagrange multipliers in the objective function, leading to a relaxation of the original problem. The relaxation problem can be decomposed into several single independent batches of sub-problems that are easier to solve. The Lagrange dual problem is represented by the LR method trying to maximize the dual objective function based on the duality theory. The solution of the dual problem is the lower bound of the original problem. Then, a feasible solution is constructed according to the solutions of the subproblems as the upper bound of the original problem. The goodness of the solution is judged by the duality gap of the upper and lower bounds. The sub gradient optimization method is used to update Lagrange multipliers. Until the iteration termination condition is met, the algorithm is terminated. The block diagram of LR algorithm is shown in Figure 3.

3. Multi-Well Batch Decomposition Algorithm

The wells in different blocks are divided into different well batches k. The production tasks of each well batch should meet the established annual production plan. Furthermore, each oil well batch relates to the others through underwater pipelines, with coupling constraints. Therefore, the multi-well batch coupling constraint is relaxed into the objective function of the original problem, and the relaxation problem is constructed to eliminate the coupling relationship between each well batch. At the same time, the relaxation problem is generally separable, which can be decomposed into sub-problems of each single well batch.

The objective function of the optimization problem is Equation (1), which means to find the minimum value of the total production cost.

3.1. Construction of Lagrange Relaxation LRP

Where the total cost is expressed as the superposition of costs of each link, as shown in Equation (49):

$$\begin{gathered} {\mathrm{Z}=Z}_1{+Z}_2{+Z}_3{+Z}_4{+Z}_5{+Z}_6 \\ ={\displaystyle \sum }_i{\displaystyle \sum }_t\underset{i}{\partial }w{f}_{i,t}+p{e}_1·el{e}_{cost}+{\displaystyle \sum }_k{\displaystyle \sum }_t\gamma I_{k,t}+{\displaystyle \sum }_i{\displaystyle \sum }_t\delta ·P_{i,t}+floor\left(\frac{TT}{T{l}_k}\right)·\tau +{\displaystyle \sum }_t\theta ·\left(d_t-p{r}_t\right) \\ ={\displaystyle \sum }_i{\displaystyle \sum }_t{\partial }_i\left[w_{i,t}\left(1-w_{i,t-1}\right)+w_{i,t-1}\left(1-w_{i,t}\right)\right]+p{e}_1·{\displaystyle \sum }_i{\displaystyle \sum }_t{\beta }_i\left(x_{i,t},w_{i,t}\right)+{\displaystyle \sum }_k{\displaystyle \sum }_t\gamma I_{k,t}+{\displaystyle \sum }_i{\displaystyle \sum }_t\delta ·P_{i,t} \\ +floor\left(\frac{TT}{T{l}_k}\right)·\tau +{\displaystyle \sum }_k{\displaystyle \sum }_t\theta ·\left(d_{k,t}-p{r}_{k,t}\right) \end{gathered}$$

(49)

The original optimization problem satisfies the constraint conditions from Equations (2)–(48).

A set of Lagrange multiplier ${\lambda }_{k,t}$ is introduced to relax the batch demand constraint (48) into the objective function to obtain the Lagrange relaxation problem (LRP), as shown in Equation (50):

$$Z_{LRP}=min\left\{Z\right.+\left.{\lambda }^T\left(p{r}_{k,t}-d_{k,t}\right)\right\}$$

(50)

where the batch demand constraint (49) can be converted into the form of Equation (51):

$$p{r}_{k,t}-d_{k,t}={\displaystyle \sum }_i{x}_{i,t}+I_{k,t-1}-I_{k,t}-d_{k,t}$$

(51)

The objective function of the relaxation problem can be translated into Equation (52):

$$\begin{array}{cc}{Z}_{LRP}\hfill & \\ =\min \{Z+{\lambda }^T\hfill & \left(p{r}_{k,t}-d_{k,t}\right)\}\hfill \\ =\mathrm{Z}+{\displaystyle \sum _k}{\displaystyle \sum _t}{\lambda }_{k,t}\hfill & \left({\displaystyle \sum _i}{x}_{i,t}+I_{k,t-1}-I_{k,t}-d_{k,t}\right)\hfill \\ ={\displaystyle \sum _i}{\displaystyle \sum _t}{\partial }_i[w_{i,t}\hfill & \left(1-w_{i,t-1}\right)+w_{i,t-1}\left(1-w_{i,t}\right)]+p{e}_1·{\displaystyle \sum _i}{\displaystyle \sum _t}{\beta }_i\left(x_{i,t},w_{i,t}\right)\hfill \\ & +{\displaystyle \sum _k}{\displaystyle \sum _t}\mathsf{\gamma }{I}_{k,t}+{\displaystyle \sum _i}{\displaystyle \sum _t}\delta ·P_{i,t}+floor\left(\frac{TT}{T{l}_k}\right)·\tau \hfill \\ & +{\displaystyle \sum _k}{\displaystyle \sum _t}\theta ·\left(d_{k,t}-p{r}_{k,t}\right)+{\displaystyle \sum _k}{\displaystyle \sum _t}{\lambda }_{k,t}({\displaystyle \sum _{i\in I}}{x}_{i,t}+I_{k,t-1}-I_{k,t}\hfill \\ & -d_{k,t})\hfill \end{array}$$

(52)

Similar terms in the relaxation problem are combined and converted into the form that can be decomposed into single-well batch sub-problems, as shown in Equation (53):

$$\begin{array}{cc}{Z}_{LRP}\hfill & \\ ={\displaystyle \sum _i}{\displaystyle \sum _t}{\partial }_i[w_{i,t}\hfill & \left(1-w_{i,t-1}\right)+w_{i,t-1}\left(1-w_{i,t}\right)+p{e}_1{\sigma }_i\left|x_{i,t+1}-x_{i,t}\right|]\hfill \\ & +{\displaystyle \sum _k}{\displaystyle \sum _t}({\displaystyle \sum _i}{\lambda }_{k,t}+p{e}_1{\beta }_i{w}_{i,t})x_{i,t}\hfill \\ & +{\displaystyle \sum _k}{\displaystyle \sum _t}\left[\left(\mathsf{\gamma }-{\lambda }_{k,t}\right)I_{k,t}+{\lambda }_{k,t}{I}_{k,t-1}\right]+{\displaystyle \sum _i}{\displaystyle \sum _t}\delta ·P_{i,t}\hfill \\ & +floor\left(\frac{TT}{T{l}_k}\right)·\tau +{\displaystyle \sum _k}{\displaystyle \sum _t}\left[\left(\theta -{\lambda }_{k,t}\right)d_{k,t}-\theta p{r}_{k,t}\right]\hfill \end{array}$$

(53)

The relaxation problem satisfies the constraints of Equations (2)–(47).

Each term in the above equation can be derived from the well batch k. Therefore, given a set of values of Lagrange multiplier ${\lambda }_{k,t}$, the Lagrange relaxation problem for multiple well batches can be decomposed into a single well batch problem for each well in each well batch, as shown in Equation (54) below:

$$\begin{gathered} \begin{array}{cc}min{L}_i\left(\lambda \right)\hfill & \\ L_i\left(\lambda \right)={\displaystyle \sum _t}\left\{\begin{array}{c}{\partial }_i\left[w_{i,t}\left(1-w_{i,t-1}\right)+w_{i,t-1}\left(1-w_{i,t}\right)+p{e}_1{\sigma }_i\left|x_{i,t+1}-x_{i,t}\right|\right]\\ +\left({\lambda }_{k,t}+p{e}_1{\beta }_i{w}_{i,t}\right)x_{i,t}+\delta ·P_{i,t}\\ +\left[\left(\theta -{\lambda }_{k,t}\right)d_{k,t}-\theta p{r}_{k,t}\right]\end{array}\right\}\hfill \end{array} \\ + floor\left(\frac{TT}{T{l}_i}\right)·\tau ,\forall i\in I \end{gathered}$$

(54)

The single well batch subproblem satisfies the constraint conditions of Equations (2)–(47), and the final production target only needs to complete the annual planned production task of the corresponding well batch. The feasible solution obtained by solving the subproblem is regarded as the upper bound of the optimal solution of the original problem.

3.2. Construct Lagrange Duality Problem

The Lagrange relaxation problem is constructed to obtain the lower bound of the optimal solution of the original problem. However, the quality of the approximate optimal solution obtained by using this lower bound can hardly be satisfied. So, the Lagrange dual problem needs to be constructed to improve the quality of the lower bound and get a tighter bound.

Lagrange relaxation problem (LRP) is considered, as shown in Equation (55):

$$\begin{array}{c}{Z}_{LRP}=\min L\left(\lambda \right)\hfill \\ ={\min }_{x,w}\left\{\begin{array}{c}{\displaystyle \sum _i}{\displaystyle \sum _t}{\partial }_i\left[w_{i,t}\left(1-w_{i,t-1}\right)+w_{i,t-1}\left(1-w_{i,t}\right)+p{e}_1{\sigma }_i\left|x_{i,t+1}-x_{i,t}\right|\right]\\ +{\displaystyle \sum _k}{\displaystyle \sum _t}\left({\displaystyle \sum _i}\left({\lambda }_{k,t}+p{e}_1{\beta }_i{w}_{i,t}\right)x_{i,t}\right)+{\displaystyle \sum _k}{\displaystyle \sum _t}\left[\left(\mathsf{\gamma }-{\lambda }_{k,t}\right)I_{k,t}+{\lambda }_{k,t}{I}_{k,t-1}\right]\\ +{\displaystyle \sum _i}{\displaystyle \sum _t}\delta ·P_{i,t}+floor\left(\frac{TT}{T{l}_k}\right)·\tau +{\displaystyle \sum _k}{\displaystyle \sum _t}\left[\left(\theta -{\lambda }_{k,t}\right)d_{k,t}-\theta p{r}_{k,t}\right]\end{array}\right\}\hfill \end{array}$$

(55)

Equation (55) is the function of decision variables $x_{i,t}$ and $w_{i,t}$, satisfying the constraint conditions Equations (2)–(47), with the purpose of obtaining the minimum cost of the objective function.

Then the dual problem can be expressed as the following Formula (56):

$$Z_{DP}=\max \left(minL\left(\lambda \right)\right)$$

(56)

$$s.t.\hspace{1em}\lambda \ge 0$$

(57)

The above equation is the function of Lagrange multiplier $\lambda$. The fixed decision variable obtained from the relaxation problem is taken as the input of the dual problem, and Equation (57) is taken as the constraint condition of the dual problem. The maximum value in the minimum problem is obtained to provide a tighter lower bound.

3.3. Algorithm Iteration

In order to make the feasible solution obtained by the decomposition algorithm approach to the optimal solution of the original problem step by step, the Lagrange multipliers given at the beginning need to be updated iteratively. So, the upper and lower bounds provided by the dual problem and the relaxation subproblem are closer and tighter until the desired requirements of the decision maker are met.

The updating principle of Lagrange multiplier adopts sub gradient descent algorithm, as shown in Equation (58):

$${\lambda }_{k,t}^{j+1}={\lambda }_{k,t}^j+{\epsilon }^j{g}_{k,t}\hspace{1em}\forall t\in T,k\in \mathrm{K}$$

(58)

The ${\epsilon }^j$ is the step size of the $j$ iteration, and $g_{k,t}$ is the sub gradient of coupling constraint in the Lagrange relaxation problem, which can be expressed as Equation (59):

$$g_{k,t}=d_{k,t}-p{r}_{k,t}=d_{k,t}+I_{k,t}-{{\displaystyle \sum }}_i{x}_{i,t}-I_{k,t-1}\hspace{1em}\forall t\in T,k\in \mathrm{K}$$

(59)

Step size ${\epsilon }^j$ can be defined as the following Equation (60) [35]:

$${\epsilon }^h=v\frac{L^{*}-L^j}{\Vert g_{k,t}^j{\Vert }_2}\hspace{1em}0\le v\le 2,\forall t\in T,k\in K$$

(60)

In the formula, $L^{*}$ represents the best objective function value found until the $h$ step of iteration. The size of duality gap is taken as the judgment criterion. $L^j$ is the dual function value obtained in the $j$ iteration. $v$ is the given initial step size, which is generally between 0 and 2. When the duality gap is less than a given requirement or reaches a certain CPU time (number of iterations), the algorithm terminates. The duality gap can be expressed as Equation (61):

$$\eta =\frac{L_{sub}-L_{DP}}{L_{DP}}\times 100\%$$

(61)

In Equation (61), $L_{sub}$ represents the value of the objective function corresponding to the feasible solution of the subproblem, and $L_{DP}$ represents the value of the objective function of the dual problem. All of the above symbols are detailed in Abbreviations.

4. Case Studies

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation and the experimental conclusions that can be drawn. In order to compare with the solution results of ALPHAECP and verify the advantages of this algorithm in solving large-scale problems, four cases with different model sizes are given. Case 1 includes 1 oil well batch, 4 production wells and 12 planning cycles. Case 2 includes 2 well batches, 8 production wells and 12 planning cycles. Case 3 includes 3 oil well batches, 12 production wells and 12 planning cycles, the same as the case of published papers [18]. Case 4 includes 3 well batches, 12 production wells and 24 planning cycles. The crude oil requirement per well batch per month is given separately.

Given four cases, the model size gets progressively larger, with Case 2 having one more well lot than Case 1, Case 3 having one more well lot than Case 2 and Case 4 having one more planning cycle than Case 3. Case 1 is a single oil well batch problem, which is not applicable to the proposed batch decomposition algorithm. It is directly solved by ALPHAECP solver, and the results of Case 1 are used for comparison of other cases. Other cases are solved by the proposed algorithm respectively, and the results of the unique variables are compared to verify the solution effect of the algorithm on different scale problems.

4.1. Results Presentation

As shown in

Table 1. Case 1.

Well Batch	Monthly Demand of Production
	1	2	3	4	5	6	7	8	9	10	11	12
1	12,600	15,000	15,000	16,200	9000	27,000	15,000	15,000	22,000	18,000	16,200	9000

,

Table 2. Case 2.

Well Batch	Monthly Demand of Production
	1	2	3	4	5	6	7	8	9	10	11	12
1	12,600	15,000	15,000	16,200	9000	27,000	15,000	15,000	22,000	18,000	16,200	9000
2	21,000	16,800	18,000	9000	11,400	15,000	9000	18,800	15,000	14,400	15,000	19,800

,

Table 3. Case 3.

Well Batch	Monthly Demand of Production
	1	2	3	4	5	6	7	8	9	10	11	12
1	12,600	15,000	15,000	16,200	9000	27,000	15,000	15,000	22,000	18,000	16,200	9000
2	21,000	16,800	18,000	9000	11,400	15,000	9000	18,800	15,000	14,400	15,000	19,800
3	19,200	16,200	9000	9000	10,200	9000	23,400	16,200	9000	24,000	13,200	21,000

and

Table 4. Case 4.

Well Batch	Monthly Demand of Production
	1	2	3	4	5	6	7	8	9	10	11	12
1	12,600	15,000	15,000	16,200	9000	27,000	15,000	15,000	22,000	18,000	16,200	9000
2	21,000	16,800	18,000	9000	11,400	15,000	9000	18,800	15,000	14,400	15,000	19,800
3	19,200	16,200	9000	9000	10,200	9000	23,400	16,200	9000	24,000	13,200	21,000
	13	14	15	16	17	18	19	20	21	22	23	24
1	21,000	16,800	18,000	9000	11,400	15,000	9000	18,800	15,000	14,400	15,000	19,800
2	19,200	16,200	9000	9000	10,200	9000	23,400	16,200	9000	24,000	13,200	21,000
3	12,600	15,000	15,000	16,200	9000	27,000	15,000	15,000	22,000	18,000	16,200	9000

, the monthly plan arrangement of each batch of oil wells is given in the table. There are four wells in each batch. The effectiveness of the algorithm is verified by four examples.

These cases are computed by GAMS win32 24.0.2, and subproblems are solved by the solver of ALPHAECP in an Intel core i5-7500 CPU, 3.41 GHz machine with 8 GB of RAM. The model size statistics directly solved by ALPHAECP solver are shown in

Table 5. Directly solved model size statistics.

Case	Formula for the Number	Number of Nonlinear Terms	Number of Discrete Variables	Number of Continuous Variables	Duality GAP (%)	CPU Run Time (S)
CASE 1	2219	180	768	1587	1	100
CASE 2	4438	394	1536	3174	1	1000
CASE 3	6656	591	2304	4760	1	7200
CASE 4	13,312	1167	4608	9520	5	14,400

. Furthermore, the comparison results between multi-well batch model solution and LR algorithm are shown in

Table 6. Comparison of algorithms.

Case		Cost	Relative Value of Difference (%)	Duality GAP (%)	CPU TIME (S)	Relative Value of Difference (%)
CASE 1	ALPHAECP	163,786,288		1	12.27
	LR
CASE 2	ALPHAECP	355,809,355	4.7	1	242.02	43.7
	LR	372,569,324		0.98	136.25
CASE 3	ALPHAECP	515,030,600	4.1	1	2515.13	48.9
	LR	535,930,049		0.96	1283.64
CASE 4	ALPHAECP	941,556,300	3.8	5	14,400	61.6
	LR	978,023,560		0.72	5531.25

.

4.1.1. Case 1: Single Oil Well Approval for 12 Months

The production schedule of Case 1 is shown in Table 1. Case 1 is a single-well batch problem, which can be solved directly by using the built-in ALPHAECP solver of GAMS. The solution gap is set at 1% and the CPU runtime is set at 100 s. The statistical results of Case 1 are shown in Table 5. The solution of single well batch is not applicable to the proposed batch decomposition algorithm, and the results are used to compare the results of the four cases to verify that large-scale problems are difficult to be solved directly. Thus, the validity of the algorithms in this chapter is highlighted by comparison.

4.1.2. Case 2: Two Oil Wells Were Approved for 12 Months

In Case 2, the size of the simulated reservoir can be divided into two different production blocks according to geographical conditions. Each block is the same batch of oil wells, and each block contains four production wells. Draw up the production plan for the next 12 months according to the market demand for different oil well batches. Table 2 shows the given monthly production schedule. The directly solved model size statistics are shown in Table 5. Then, the batch decomposition algorithm based on Lagrange relaxation proposed in this paper is used to decompose the two batches of the original problem into single batches to solve the sub-problems respectively. At the same time, LR dual problem is constructed to solve iteratively. Among them, the single-batch sub-problem and dual problem of Lagrange are solved directly by using the ALPHAECP solver built by GAMS. The comparison results of the two algorithms are shown in Table 6.

In Case 2, the total cost increased by 4.7% compared with the direct solution of LR algorithm. The solution time was significantly reduced by 43.7%.

Figure 4 shows the optimization solution results of Case 2. It can be seen from the figure that the output of each month plus the inventory of each month can meet the delivery requirements of that month. Therefore, there are no out-of-stock penalties occurring and the production planning requirements given in Case 2 can be met on time. So, applying the batch decomposition algorithm based on Lagrange relaxation to solve the planning optimization problem of the 12-month batch size of two oil wells, the market demand is satisfied, the process is not punished and the results are effective.

4.1.3. Case 3: Three Oil Wells Were Approved for 12 Months

In simulation case 3, the monthly production plan arrangement is shown in Table 3, the directly solved model size statistics are shown in Table 5 and the comparison of solution results of the two algorithms is shown in Table 6.

In Case 3, compared with the direct solution, the total cost of LR algorithm increased by 4.1%. The solution time was significantly reduced by 48.9%.

Figure 5 is the optimization results of Case 3. It can be seen from the figure that the output of each month plus the inventory of each month can meet the delivery requirements of that month. Therefore, there is no shortage penalty and the production plan requirements given in Case 3 can be completed on time. The results obtained by LR algorithm are effective and can meet the production requirements.

4.1.4. Case 4: Three Wells Group of 24 Months

In simulation case 4, the monthly production schedule arrangement is shown in Table 4, and the directly solved model size statistics are shown in Table 5, and the comparison of solution results of the two algorithms is shown in Table 6.

In Case 4, the total cost of solving LR algorithm increased by 3.8% compared with the direct solution. The solution time was significantly reduced by 61.6%.

Figure 6 shows the optimization solution results of Case 4, which is similar to Case 3. The result solved by LR algorithm can meet the production demand and be efficient.

4.2. Results Presentation

The three control case groups with different sizes are compared, and the relationship between solving gaps and computing time is shown in Figure 7. The blue line is the trend of solving gaps directly solved by the ALPHAECP solver for MINLP model at a specified time and a specified suspension gap. Furthermore, the red line is the trend of solving gaps solved by LR algorithm for MINLP model. Obviously, in three cases of different sizes, the GAP value of the model solved by LR algorithm converges faster and can reach the required optimal GAP value in a shorter time.

The results of the four cases are compared together to verify the proposed batch decomposition algorithm for solving problems of different sizes.

The final algorithm comparison results of the four cases are shown in Table 6 above. By comparing the solution results of each case in the table above, the following conclusions can be drawn. In terms of CPU running time, the Lagrange relaxation-based batch decomposition algorithm proposed in this chapter on is at least 40% better than solving directly using the solver. Whether the problem of larger model size is caused by the increase in the number of devices or the increase in the planning period, the algorithm can significantly improve the solving efficiency, with the average solution time increasing by 49%. Moreover, the larger the problem size to be solved is, the greater the improvement of the solving efficiency will be. However, the final cost obtained by this algorithm is not better than that obtained by ALPHAECP. The reason may be that the solution obtained by this algorithm is only approximately the optimal solution. By relaxing the constraint into the objective function, the final optimization cost may become higher. However, as the problem gets bigger, the solutions get closer and closer. The results obtained by the proposed LR algorithm may not be globally optimal. The focus is to improve solving efficiency and to obtain the solution in reasonable time for industrial-scale problem.

5. Conclusions

The LR algorithm is used to decompose the large-scale MINLP problem, which is solved optimally for multiple well batches simultaneously, into subproblems that are solved independently for a single well batch. The multi-well batch coupling constraint is relaxed into the objective function by introducing Lagrange multipliers to remove the multi-batch coupling from the constraint. The relaxed problem is still of the MINLP model and is still more difficult to solve than LP or MILP models. However, the problem uses the LR algorithm to increase the feasible domain range and reduce the time required to search for a solution. Moreover, the resulted small scale MINLP model is not approximated and its result is more precise. The original problem is decomposed into several single batch sub-problems and the feasible solution obtained by solving the sub-problems is taken as the upper bound of the original problem. The Lagrange dual problem is also constructed as a lower bound on the original problem. The pairwise gap that exists between the upper and lower boundaries is used as a criterion for judging the quality of the feasible solution. The Lagrange multiplier is updated using the sub gradient method, which iteratively updates the upper and lower bounds continuously, making them tighter and tighter until the decision maker′s expectations are met. Four cases of different sizes are solved using the algorithm in this paper and solved directly by the solver through the case of a given four different scale, using this algorithm and direct solver. By analyzing and comparing the two results, the algorithm proposed in this paper can effectively improve the solving efficiency for tackling large-scale problems, with an average saving of 49% in solving time. Furthermore, the greater the scale of the problem, the better the results and the more instructive it is for practical industrial production.