Abstract
As an important aspect of oil and gas field construction, oil and gas field surface engineering often affects the efficiency and safety of oil and gas production, and it requires a large investment. In the past, much work has been done on layout optimization for gathering pipeline networks. However, few studies have considered multiway valves in the optimization of pipeline networks. Compared to traditional metering processes, a process using multiway valves can reduce construction and operation costs and enable automation of the well-selection operation. In this paper, an MINLP model is established in which the number of multiway valves and their numbers of channels are considered special constraints and the number of multiway valves and the associations between wells and multiway valves are treated as optimization variables. A specific heuristic algorithm for solving this problem is also proposed in this paper. We consider the coordinates of real-world wells and of randomly generated well locations as different examples to analyze the performance of this algorithm. These examples demonstrate that the algorithm can initially adjust the network through single-step iteration, while double-step iteration is efficient when all channels of a multiway valve have been associated with pipelines, and multistep iteration can help the objective function escape from local optima. Finally, numerical analysis results prove that the proposed algorithm can be used to efficiently solve the problem of interest and exhibits stable convergence.
1. Introduction
In recent years, as the production scale of oil and gas fields has expanded, the construction and labor costs incurred for the measurement of oil and gas production per well have also increased. In addition, the practical implementation of conventional metering methods is subject to certain limitations, especially in areas such as deserts and deep-sea regions. There are three main metering methods for oil and gas fields: single-well metering, metering by manual rotation, and automatic rotation metering. The single-well metering process requires the installation of metering equipment for each well, as shown in Figure 1(a). The initial investment is relatively expensive, especially when the well production rate is below the initial estimate. In the manual rotation measurement process, as shown in Figure 1(b), each valve that is to be opened or closed needs to be operated manually. Although it is not necessary to install metering devices at each well site, metering a well will require the manual operation of more than a dozen valves, and after metering, these valves will need to be operated again to discharge the liquid in the metering separator, thus requiring more operation time and incurring a higher cost. In addition, the metering device occupies a large area, and such devices are not suitable for offshore platforms. Regarding the third method, retrofitting an actuator onto the original manual valve of an existing manifold will greatly increase the construction cost. Therefore, it is necessary to find a device that can enable automatic well-selection metering while incurring lower construction and operation costs.
(a)
(b)
(c)
At present, many oil and gas field production enterprises use multiway valves (MVs) to achieve automatic rotation metering. The MV-based metering process can reduce the amount of manual operation required in remote areas, which is greatly beneficial for reducing operation costs. The metering process with an MV is shown in Figure 1(c). The blue dashed arrow indicates the flow that needs to be metered at present, and the black arrows indicate the flows from production wells. When switching to metering another well, the control system requires only remote operation, thus leading to considerable savings in operation cost. Compared to the traditional methods, the MV-based metering process reduces the numbers of valves, gauges, and actuators that need to be installed on the pipelines. Each device also occupies a smaller space, making these devices more suitable for a wider variety of oil and gas fields, especially for coalbed methane mining, shale gas, and other unconventional oil and gas resources. From here on, when the term “valve” is used in this paper, it refers to an MV.
The number of channels of an MV is a characteristic that is defined by the manufacturer and determines how many pipelines can access the MV. If an MV has 16 channels, this means that at most 16 wells can be connected to the MV. Manufacturers make MVs with various channel counts, such as 12, 14, and 16. Depending on the conditions in the field or other requirements such as consistency principle, the designer needs to select MVs with an appropriate number of channels to meet the requirements of the project. Therefore, during the design process, the number of wells connected to an MV is constrained by the number of channels of the MV, which depends on the designer’s selection. Depending on the constraint relaxation degree, there are three possible ways to specify MV-related constraints:(1)The number of MVs is unrestricted(2)Allowable ranges are given for the number of MVs and their channel counts(3)The channel counts of the MVs and the number of MVs are fixed in advance
For a given target oil and gas reservoir, the traditional solution is to invite experts to design feasible solutions based on real scenarios, but this approach does not necessarily guarantee the most economical cost. By contrast, optimization methods can help a designer to obtain the best solution. However, in the past, researchers have rarely considered constraints related to the channel counts of the MVs and the number of MVs, although such constraints are important to ensure that the optimal solution satisfies the required service conditions. In this paper, the problem of interest is formulated as a special P-median problem, and the constraints mentioned above are integrated into the mathematical model. By treating the associations between wells and MVs and the number of MVs as variables and the minimum cost as the objective function, a general mixed integer nonlinear programming (MINLP) model is established. Based on the connection and foresight matrices of the production wells and MVs, a foresight matrix search algorithm (FMSA) is proposed in this paper to solve the MINLP problem. In this algorithm, the nonlinear part of the problem is treated as a subproblem, and the optimization results are obtained through iteration. By examples, the convergence of the algorithm can be ensured, as can the quality and efficiency of the solution. In addition, it is proven that, with this algorithm, the solution process can be completed in a reasonable time even if the problem is large in scale.
This paper is divided into six sections: In Section 1, we have introduced the problem motivating this paper and what we have done to address it. In Section 2, we review the related research contributions that have been reported in the past. In Section 3, the detailed mathematical model is presented, and the special constraints are formulated. Subsequently, the details of solving the mathematical model using the proposed algorithm are introduced in Section 4. The results are analyzed in Section 5. Finally, conclusions and discussions are presented in Section 6. The data used in this paper can be found in Tables 1–5.
2. Literature Review
The complexity and large scale of the corresponding mathematical models make MINLP problems a popular and difficult topic to address for scholars around the world. The methods for solving MINLP problems can be generally divided into deterministic methods and heuristic methods. Methods of the former type can yield globally optimal solutions for small-scale problems but are often inefficient for large-scale problems, whereas the quality of the solutions obtained with the heuristic methods that have been proposed to date usually cannot be guaranteed, and these methods cannot be used to solve large-scale problems within a limited time.
The problem of interest here is different from the well group partitioning (WGP) problem, which usually considers factors such as the gathering radius and pipeline diameter. The purpose of the conventional WGP model is to determine the optimal mapping between wells and stations. In other words, subject to well constraints and other constraints, the wells are divided among different gathering stations to save construction costs. This problem is a kind of location-allocation problem. In a discrete location-allocation problem, P facilities are selected from a set of available candidate sites such that a given objective function is minimized. This problem is specifically referred to as the P-median problem. However, there are some differences between the model established in this paper and the general P-median model. In the problem considered here, the alternative facility (MV) locations are the same as the well locations, meaning that , whereas in the conventional P-median problem, is often less than n. Moreover, the optimization target of the conventional P-median problem is usually to minimize the total distance from facilities to customers, while in this paper, the target is to minimize the cost of all equipment, including pipelines and devices, and the model is constrained by the special constraints of MVs. In addition, special constraints on the numbers of channels and facilities (MVs) are imposed in the model considered here. All of these factors increase the scale of the optimization problem, making the problem more complicated and difficult to solve. In [1, 2], the authors have proven that the P-median problem is an NP-hard problem when the number of facilities is uncertain. That is, the MINLP problem solved in this paper is a special case of the location-allocation problem. In this problem, the number of MVs m and the associations between wells and valves are treated as optimization variables. The associations between pipelines and various facilities are precisely the elements of interest in oil and gas field surface engineering, and reasonable optimization of the pipeline network structures between facilities can effectively reduce the engineering and construction costs. There are many possible pipeline network structures, including star, tree, and ring structures. With the development of methods for exploiting unconventional oil and gas resources, tree-tree network structures [3, 4] have been used in marginal gas reservoirs for resources such as coalbed gas and shale gas. Such a structure can considerably reduce costs by virtue of the connection configuration between the wells and co-construction with existing wells and stations. Similarly, research on star-tree [5] and star-star [6] pipeline networks has become quite mature.
A pipeline network also includes facilities such as manifolds, platforms, and compressors. Many studies have investigated oil and gas network models that contain various types of facilities. Rodrigues et al. [7] proposed a 0-1 programming model to determine the connections between FPSO units, manifolds, and facilities as well as the sizes and locations of platforms. In [8, 9], the authors proposed a mixed integer linear programming (MILP) model to optimize the connections between wells and platforms while considering not only the installation time but also the linear drop of the reservoir pressure. Furthermore, Ramos Rosa et al. [10] considered the reservoir dynamics and secondary development to establish an MILP model seeking the maximum net present value (NPV) while simultaneously optimizing the manifold layout, routing, and pipeline diameter. Most of the oil and gas field layout optimization problems mentioned above can be transformed into mixed integer programming (MIP) problems, i.e., programming problems in which some of the independent variables take integer values while others take continuous values. MIP first became an independent branch of mathematical problem-solving when Gomory [11] proposed the cutting-plane method. Generally, MINLP problems can be solved using two types of methods: deterministic algorithms and heuristic algorithms. A deterministic method of solving an MINLP problem is to simplify it to an MILP problem, which can be effectively solved by means of dual simplex or sequential quadratic programming. However, the main difficulty is that appropriate linear approximations must be adopted for nonlinear functions. To this end, Möller [12], Martin et al. [13], Tomasgard [14], and Nørstebø [15] et al. used piecewise approximation to deal with the nonlinear functions in the model. Mikolajková et al. [16] proposed a linearization method for solving MINLP problems. Deterministic methods can guarantee a globally optimal solution when solving an MINLP problem, but they usually require considerable resources to solve large-scale problems.
Because of the complexity of MINLP problems, metaheuristic algorithms, including genetic algorithms (GAs) [17] and particle swarm optimization (PSO) [18], are also used to solve problems related to oil and gas pipeline networks. However, the computational burden of heuristic algorithms is usually heavy, making them an inefficient means of solving large-scale problems. In addition, either heuristic methods converge slowly or the quality of their results cannot be guaranteed. Because of the special characteristics of MVs, no related research on optimization problems involving MVs has been reported in recent years; therefore, to address this lack, a corresponding MINLP model and an algorithm for solving it are proposed in this paper.
3. Problem Definition and Formulation
3.1. Connection Matrix and Assumptions
In general, optimal solutions to various problems can be found through various optimization methods, but the solution found may not necessarily be suitable for the actual conditions and engineering design requirements of a particular situation. For a given oil and gas field, especially an unconventional reservoir with scattered wells, it is necessary to reserve some MV channels for future connections to additional production wells; therefore, either the channel counts of the MVs and the number of MVs need to be fixed in advance or ranges should be defined for these parameters to satisfy the need to reserve valves or channels for future use. If experts define several possible schemes in accordance with the given requirements, the most appropriate solution can be obtained by comparing the results of these schemes, but the quality of a solution chosen in this way cannot be guaranteed; thus, it is necessary to formulate and solve the corresponding optimization problem with constraints on the number of MVs and their channel counts.
Suppose that there are n wells that need to be divided into m sets and connected to the corresponding MVs. The MVs will be placed at well sites because no additional land acquisition is allowed; therefore, all well sites are possible alternative locations for MVs, and the number of MVs m is unknown. Let be the connection matrix, which reflects the connection relationships between the wells and the MVs. If , then well i is in well set j, which means that well i is connected to MV j. Thus, the remaining elements in the ith column of the matrix are all equal to 0. Here, a matrix , with dimensions of is introduced to indicate whether an MV is placed at the ith well site. If , this means that an MV is placed at the ith well site. If , this indicates that there is no MV placed at the ith well site. By incorporating into the matrix , we obtain the matrix , as shown in Figure 2. If in the matrix , then the corresponding MV can be identified by finding the row j that contains the element in the ith column of (below the first row) that has a value of 1; that is, MV j is placed at the ith well site. Thus, and the other wells connected to valve j can be found at the same time. Now, considering the characteristics of MVs, the following assumptions are adopted before solving the mathematical model in this paper:(1)Generally, the number of channels of an MV ranges from 6 to 16. According to the designer’s requirements, the corresponding MV specification will be given as either a constant or a range before optimization. Thus, the special constraints on the MVs will be specified beforehand.(2)All wells need to connect directly to an MV because this is the only way in which wells can be selected for metering; there are no pipeline connections between wells.(3)The possible locations of the MVs are preselected by the engineering company. One of the strategies that is usually adopted is to place the MVs at well sites to reduce land acquisition, which means that all well sites are potential alternative locations for MVs.(4)In this paper, we use the Euclidean distance to quantify the length of the pipeline between a well and an MV. Additionally, terrain and obstacles are not considered.(5)The costs of valves and instruments other than the MVs are not considered. These costs are neglected in the calculation since the use of MVs simplifies the metering process.(6)We assume that the same material and diameter are used for all pipelines. Therefore, factors such as pressure and flow rate are not considered.
3.2. MINLP Model
3.2.1. Objective Function
Given the above assumptions, the objective function can be simply expressed as follows:where represents the cost of all pipelines, which is equal to the sum of the costs of the pipelines connected to each MV:where represents the total cost of all pipelines connected to MV j, which can be expressed as follows:where “” denotes the matrix dot product operation.
The total cost of the pipelines connected to valve j is the dot product of , , and , where is the pipeline cost coefficient matrix, is the jth row of the connection matrix, and is a matrix consisting of the pipeline lengths from the wells to the valve. represents the coordinates of the ith well, . Specifically, , means that the ith well is connected to the jth MV. Then, represents the pipeline length from well i to the corresponding MV under the assumption that this well is connected to valve j; this pipeline length can be expressed as
represents the coordinates of the jth MV, which are unknown and are determined by solving a subproblem. The detailed solution method will be described in Section 4. Then, the cost of all pipelines can be written as follows:
Considering the costs of the MVs themselves, the total cost of all MVs can be expressed aswhere is a matrix with dimensions of that is used as a compact representation of the connection matrix, in which is a binary variable. If all elements in the jth row of the connection matrix are 0, then ; otherwise, . is the matrix of the MV cost coefficients. is the dot product of and . By summing the costs of the pipelines and valves, the objective function is obtained as follows:
3.2.2. Constraints
The objective function has been established, as shown in equation (7). Now, we will describe the constraints of the model. In this paper, the variable m is constrained by and , as follows:where is the upper bound on the number of MVs, while is the minimum value that m can take; both and are integers. As the range of m shrinks, the number of MVs will finally be specified; that is, . Only one nonzero element is allowed to exist in each column of the connection matrix because each well can be connected to only one MV. This constraint can be expressed as follows:
Given that row in the connection matrix represents one MV, the sum of the elements in each row of this matrix cannot be greater than the channel count of the corresponding MV. Since the channel counts of the MVs selected by the designer might differ, is used to represent the maximum MV channel count provided by the manufacturer, and is used to represent the specified channel count of the jth MV:
The following constraints should be satisfied:
If is the channel count of MV j as chosen or required by the designer and is a binary variable that represents whether a specified channel count for the jth valve is actually used to constrain the mathematical model, according to the constraint relaxation degree, then the maximum number of wells that can actually be connected to MV j is defined as
By combining this equation with equation (11), we obtain
Usually, it is not permitted to connect zero wells to a given MV. Therefore, taking 1 as the default lower bound for each row of the matrix, we have
Sometimes, it is unreasonable to connect only a few wells to an MV. In this case, we use to indicate a specified lower bound on the number of wells to which an MV is allowed to be connected. The corresponding constraint can be expressed as
The minimum number of wells that can be connected to MV j is determined by . is a binary variable that controls the value of . If , a specific lower bound is not required and the default is used, whereas if the value of this variable is 0, a specific lower bound is applied. Thus, using the same formulation used for the upper bound, we obtain the following constraint:
By combining equations (16)∼(18), we can write the constraint related to the lower bound as follows:
The following relationships should be satisfied:
Equations (20) and (21) are used to ensure that the lower bound will not exceed the upper bound. Equations (22) and (23) ensure that both the required lower bound and the bound that is actually applied are at least as large as the default lower bound. Finally, equation (24) stipulates that the real lower bound should not exceed the required lower bound. Since the matrix represents the MVs, the sum of the elements of the matrix should satisfy
Because only one MV is allowed in each well group, when we subtract from any row of the connection matrix L and then add the absolute value of the results together, we will have for each row of . Then, the constraint that only one MV is allowed in each well group can be written as
According to the above equations, the mathematical model can be expressed as
The decision variables used in the mathematical model represent the associations between the wells and the MVs, and the number of valves and their channel counts are treated as constraints. The binary variables and determine the range of the channel counts of the MVs, while and determine the range of the number of MVs.
4. The Proposed Algorithm for Solving the MINLP Problem
The underlying principle of methods such as steepest gradient descent [19] is to find the direction in which the objective function decreases the fastest at the present point, take a step in that direction, and then repeat this process, thus iteratively reducing the value of the objective function. We use the same idea to search for the optimal solution in a discrete space. In this paper, the mathematical model is described in the matrix form; thus, the step size and the number of steps must be expressed in terms of a measure of the distance between two matrices. If we treat the whole matrix as a point, then we can move one nonzero element in the matrix to obtain a new matrix, and this new matrix can be treated as a point adjacent to the original point; that is, the new matrix is one step away from the original matrix. All the points obtained by changing one individual element are the points closest to the current point. As an example, let us consider a connection matrix that is equivalent to the 5 × 5 unit matrix; then, we can obtain a foresight matrix corresponding to this connection matrix, as shown in Figure 3. The square columns represent the objective function values, where each row of the matrix represents a well and each column represents an MV. Then, the values in the jth column of the matrix represent the values of the objective function after the corresponding nonzero elements are moved to the jth valve; that is, the square column in the ith row and jth column represents the objective function value after well i is connected to MV j.
If the initial connection matrix is a unit matrix, then initially, the values on the diagonal of the foresight matrix are equal to the objective function value of the original matrix because the values on the diagonal of the foresight matrix, represented by square columns with stripes, are obtained only when no element is moved. For a case in which the constraints are not satisfied after a single nonzero element is moved to another position, the corresponding objective function value is represented by a dark gray square column; all such objective function values are replaced by a sufficiently large penalty value M to ensure that the algorithm will not select such a point. White columns are used to represent the feasible domain, meaning that these positions in the matrix represent points that satisfy the constraints of the mathematical model. The minimum value in the feasible domain, denoted by , is found as the result of the current iteration, as shown in Figure 3. As shown in the figure, the minimum objective function value is found in the 4th row and 1st column; therefore, we set the element in the 4th row and 1st column to 1 and set the remaining elements in the 4th row to 0. In this way, we can find the path that causes the objective function value to drop the fastest from the current point, thus completing one single-step iteration. When all values in the feasible domain are greater than or equal to the current value, the current point is considered to be a locally optimal solution.
We can always construct a foresight matrix such as that shown in Figure 3 after every step and then find the point that minimizes the objective function for the next step until no value less than the current value of the objective function is found. To obtain a better solution once the current point has fallen into a local optimum, we can consider increasing the step size, thus enabling the algorithm to escape from the local optimum and find a better feasible solution after further iterations.
4.1. Solving the Subproblem
All wells connected to MV j are assigned to the same well set. The MV will be placed at one of the wells in this set such that the total weighted length from the wells to the MV is as short as possible, as expressed in equation (23). Thus, a subproblem must be solved to find the placement of valve j that results in the shortest pipeline length in well group j:
This subproblem can be solved using a simple algorithm. The coordinates of the wells connected to valve j are represented by . Because the location of MV j is unknown, is assumed to represent the valve coordinates that minimize the sum of the pipeline lengths connected to valve j. The valve will be placed at each well site in turn; then, we can compare the results to find the optimal placement position for the valve.
4.2. Details of the Algorithm
The algorithm for solving the MINLP problem obtains the solution through iteration. The iterative process mainly comprises 4 steps; this process is depicted in Figure 4 and described in detail below.
4.2.1. Step 1: Initialization
If the constraint on the number of MVs m is relaxed, m can range from 1 to n. In the extreme case, a unit matrix is taken as the initial feasible solution, meaning that each MV is initially assumed to be connected to only one well. Generally, if , we can initially set , and the wells will be sequentially assigned to the n MVs in accordance with the constraint conditions.
4.2.2. Step 2: Iteration
The initial costs can be obtained by substituting the initial distribution matrix into the objective function. Then, a blank matrix with dimensions of is initialized as the foresight matrix.
Starting from the first column, each nonzero element 1 in the connection matrix is moved up or down by one row; each movement will change the distribution, and the objective function value will change as the connection matrix changes. The corresponding objective function value is calculated for every movement, and the new objective function value is filled into the corresponding position of the foresight matrix.
In the example shown in Figure 5, element 1 in the first row and first column is moved to the second row and first column to obtain a new connection relationship. After substituting this new connection relationship into the objective function, the value corresponding to the new connection relationship is obtained. is then filled into the second row and first column of the foresight matrix, corresponding to the new position of the nonzero element in the connection matrix. If the constraints are not satisfied by this new connection relationship, a sufficiently large penalty value M is placed in the corresponding position instead. Once the foresight matrix has been completely filled, the path of fastest descent can be determined. This process is illustrated in Figure 5.
According to the constraints, the nonzero elements in the connection matrix can move only up and down; they cannot move left or right. Since each column of the matrix represents a well, when a nonzero element is moved from the jth row to the kth row in column i, this represents that this well is disconnected from MV j and connected to MV k, and the value of the element in the ith row and jth column returns to 0.
In this way, each 1 value in the connection matrix is moved from the top row to the bottom. Each time a 1 value is moved to a new position, a new objective function value is generated, filling in the corresponding position in the foresight matrix. When the new matrix represents a connection relationship that does not satisfy the constraints, a sufficiently large penalty value M is filled into the foresight matrix instead. Once the foresight matrix has been completed, we then search for the minimum value in this matrix that is smaller than the current value. For example, if is the smallest value in the foresight matrix, then for the next iteration, 1 in the ith column of the connection matrix is moved to the jth row, and the value at the original position is set to 0; thus, one single-step iteration is completed. By repeating this process, the current value of the objective function can be iteratively reduced, and new connection relationships can be formed. Figure 5 illustrates this iterative process for a 5 × 5 foresight matrix.
If all nonzero elements in a given row of the connection matrix leave that row by moving up or down, this means that the corresponding MV does not exist; in this case, the empty row should be removed from the connection matrix, decreasing the feasible domain. In this way, the number of MVs can be found. Specifically, if the elements of the jth row are all 0, then the jth row should be removed from the connection matrix, as shown in Figure 6.
4.2.3. Step 3: Increasing the Step Size
If a lower cost cannot be obtained by moving a single element in the connection matrix up or down, this means that the solution has become stuck in a local optimum. Therefore, the step size must be gradually increased to allow the algorithm to escape from this local optimum. For example, suppose that the matrix L is initially set as the unit matrix and that the initial step size is 1. When the objective function value cannot continue to decrease after a certain number of iterations, then increasing the step size can help to find a new connection relationship that will allow the objective function to continue to decrease. More specifically, suppose that the ith row of the matrix L contains two nonzero elements, meaning that MV i is connected to two wells, while the jth row contains three nonzero elements. In this case, both nonzero elements in the ith row are added to the jth row, setting all elements in the ith row to zero, and the new total cost under the new connection relationship is obtained. This cost value is filled into the ith row and jth column of the foresight matrix. Then, the matrix is restored to its previous state, and the algorithm proceeds to fill in the other elements in the foresight matrix in the same way. Each element in the foresight matrix is an objective function value, and these values are used as the basis for iteration.
In the foresight matrix for multistep iteration, the values on the diagonal should be equivalent to the objective function value when no wells are moved because for each element on the diagonal, the row and column correspond to the same well; therefore, the penalty value M is filled into the positions on the diagonal of the matrix. Moreover, the matrix elements below the diagonal to the left are symmetric to the matrix elements above the diagonal to the right; therefore, all matrix elements in either the lower left or the upper right can be directly set equal to M to avoid redundant calculations. As before, if the constraints are not satisfied after the elements of the connection matrix are moved, the penalty value M is also directly filled into the corresponding position. After the foresight matrix has been fully filled in, the smallest value in this matrix can be found. If this minimum value is smaller than the cost obtained in the previous iteration, then we can find the row and column corresponding to this minimum value, add the elements of the ith row to the jth row, and remove the ith row, whose elements are now all 0. At this time, the number of MVs m is set to m minus 1.
If the remaining channels of an MV do not allow connection to any further wells, the wells will be sequentially assigned to other MVs. This iterative process will repeat until the objective function value can no longer be reduced by such multistep iterations. Then, the algorithm will either return to Step 1 or proceed to Step 4.
4.2.4. Step 4: Adjustment
If there is no MV whose channels are fully occupied or for which the number of wells connected to that valve has reached the lower bound , then the optimal solution can be obtained through single-step or multistep iteration. However, according to the constraints, if an MV is fully occupied, then it cannot be connected to any new wells, whereas if a valve has reached its lower bound , then no more wells can be moved out of the corresponding row. In this case, the iterative process can be continued as follows:(1)The step size can be increased to 2 to test whether the objective function value will continue to decrease. First, we look for the fully occupied MV. Suppose that wells are connected to this MV; then, the remaining wells will each be moved to this MV in turn, replacing one of the existing wells, which is moved to another MV that meets the constraints. By repeating this process, the new values of the objective function are obtained and written into the foresight matrix based on the corresponding well numbers i and j. Finally, the path of steepest descent can be determined in the same way as before.(2)For an MV j that has been fully connected, we look for an MV t with more channels than MV j, where denotes the number of channels of MV j. Then, , and the total number of wells connected to each valve satisfies . Since , must be satisfied to allow the jth and tth valves to be exchanged; that is, the following conditions should be satisfied simultaneously:
If these conditions are satisfied, then the channel counts of the jth and tth MVs are exchanged, and we continue iterating steps 2 and 3 mentioned above. However, if i and t do not satisfy , , and , then the penalty value M is written into the corresponding position of the foresight matrix. The process above is repeated; if, after the Nth iteration, the objective function value cannot continue to decrease, meaning that , then the current objective function value corresponds to the optimal solution.
5. Analysis of Examples
5.1. Verification of the Validity of the Algorithm
First, a simple example is presented to demonstrate that the proposed algorithm can effectively solve the problem of interest. The associations between the wells and the MVs and the number of MVs are taken as variables, and their values are found by the algorithm via iteration; in particular, during the iterative process, multistep iteration can help the objective function value to escape from a local optimum. Let us consider a simple example for which we can make a reasonable prediction of the expected solution. Then, we can apply the proposed algorithm to solve the example and compare its solution with our expected result to prove that the algorithm can effectively solve the presented MINLP model.
For all examples in this paper, we consider steel pipes with an inner diameter of 53 mm and a wall thickness of 3.5 mm as the pipes linking the wells and valves; thus, the cost coefficient is a constant. Moreover, the cost of an MV is a known constant once the vender is decided. In this paper, the cost of an MV is treated as an average value that does not vary with the number of channels. Here, we consider an example of a field with 26 wells (Table 1; W10∼W35), for which the channel counts of the MVs and the number of MVs are not limited; therefore, the number of channels of each MV is uniformly set to 16, which is the maximum number of channels per MV that can be provided by general manufacturers.
We first perform single-step iteration, during which the number of MVs will usually stay at a relatively large value, as shown in Figure 7. When the identity matrix is used as the initial matrix, the number of MVs cannot be decreased below 8 via this process. However, increasing the step size can lead to further reduction in the value of the objective function, demonstrating that the reason that the objective function cannot be further reduced via single-step iteration is that the step size of each advance is too small. Thus, the objective function becomes stuck in a local optimum, and it is necessary to increase the step size to escape from that local optimum to allow the objective function to continue to decrease. The objective function decreases almost linearly during the initial 18 steps because the main factor affecting the decrease in the objective function value during the initial iterations is the number of MVs rather than the pipeline cost. However, the value of the objective function no longer decreases once it drops into a local optimum. After the 18th iteration, the step size must be increased to continue to reduce the objective function value. The single-step iteration process converges prematurely after the 18th iteration; thus, the findings demonstrate that the strategy of increasing the step size can efficiently help the objective function to escape from a local optimum.
The result of optimization with an increase in step size is shown in Figure 8. The hollow square symbols represent the positions of the MVs, the triangular symbols represent the wells connected to MV 1, and the circular symbols represent the wells connected to MV 2. For the 26 wells considered in this calculation, the algorithm completely divides them into two groups, and the solution satisfies all the constraints, which is obviously consistent with our expectations. This example proves that the proposed algorithm can be applied to solve the problem of interest and obtain reasonable optimization results.
5.2. Convergence of the Algorithm
Usually, we can judge whether an algorithm is valid according to whether the value of the objective function converges during the iterative process. If the algorithm is divergent, no matter how many times it is iterated, the value of the objective function will always fluctuate and will not approach a constant. In contrast, convergence of the algorithm means that the objective function value will tend toward stability after some numbers of iterations. In summary, convergence is a necessary condition for the algorithm to be valid.
Next, various examples are presented to analyze the convergence of the algorithm. First, we consider a field with 35 wells (Table 2; W1∼W35) as an example. MVs with four different channel counts are specified: . Thus, the number of channels of each MV is fixed in advance, the limitations on the numbers of pipelines that are allowed to be connected to different valves are different, and the number of MVs is also given, i.e., m = 4.
Figure 9 shows the iterative process for these 35 wells, which includes single-step iteration, double-step iteration, and redistribution of different kinds of valves. The iterative process still converges, and the optimal solution is obtained after the 21st iteration. After 14 iterations, the redistribution of different valve types is necessary to allow the objective function value to continue to decrease. Finally, a difference of 0 between two consecutive steps indicates that the final solution has been obtained.
As seen from Figure 9, the objective function value gradually tends toward stability and converges as iteration continues. The result is shown in Figure 10. Red circular symbols, blue triangular symbols, purple diamond symbols, and yellow triangular symbols are used to represent the individual wells connected to MVs 1, 2, 3, and 4, respectively, and square symbols are used to represent the MVs. The initial channel counts of the MVs corresponding to these different well groups are 14, 12, 10, and 8, respectively, while the final channel counts are 12, 10, 14, and 8, in sequence. In the case of tight constraints, the minimum value of the objective function can be obtained by using the proposed optimization algorithm.
In addition to the example above, an example with 40 wells (Table 1; W69∼108) and various channel counts is considered to prove the convergence of the algorithm; the results are shown in Figures 11(a) and 11(b). It can be seen from the figures that, during the initial iterations, the objective function value decreases at a relatively fast rate. Then, the rate of decrease slowly declines as iteration proceeds. The algorithm eventually converges when .
(a)
(b)
Again, the value of the objective function gradually converges and tends toward stability. In Figure 11(a), the iterative processes with a total of 46 and 48 channels are identical, indicating that once the number of channels exceeds a certain threshold, a further increase in the number of channels has no effect on the iterative process or the result. Additionally, although the iterative process for an initial feasible solution with 40 channels is different from the other 4 cases, the final result is not significantly different from the others, indicating that the initial solution does not have a great influence on the result.
The iterative processes for total channel counts of 46 and 48 are also fundamentally similar. Once the total number of channels exceeds 46, the iterative processes are basically the same. The results show that the algorithm presented above can stably converge. Therefore, these examples demonstrate that even under different constraints and in different scenarios, the algorithm can drive the objective function value toward stable convergence.
5.3. Parameter Analysis
In general, the designer of a pipeline system should consider the possibility of reserving channels for new wells in the future along with other factors when specifying the channel counts of the MVs and the number of MVs. However, different combinations of MVs with different numbers of channels might affect the results; therefore, we will analyze the effects of different channel counts and different numbers of MVs on the results of our algorithm. We consider a field with 40 wells (Table 1; W69∼108) as our example here, varying the constraints on the numbers of channels and valves to obtain different results. The results are shown in Table 6, from which we can see the influence of these different constraint parameters on the cost of the pipeline system.
Six different numbers of MVs are specified: 3, 4, 5, 6, 7, and 8. The table shows that the requirement of the use of too few MVs will lead to additional costs when the total number of channels is the same. The cost also rises if more than 4 MVs are used, and the smaller the total number of channels is, the faster the increase of cost is. When the total number of channels is fewer than or equal to 42, the cost rises at a significantly higher rate than in the other cases. Therefore, the total number of channels should be set to greater than 42 to make it easier to achieve lower costs. However, increasing the number of channels does not significantly reduce the overall cost once the total number of channels exceeds 48. The number of MVs has a greater impact on the cost. These findings suggest that the optimization of the number of MVs is the dominant factor affecting the value of the objective function, while the number of channels is the secondary influencing factor. These observations can provide a reasonable basis for practical design.
Next, a set of well distribution data from the field and a set of data randomly generated by MATLAB are used as examples to analyze how the well distribution density affects the optimization result. In the first of these examples, the coordinates of 40 wells (Table 1; W69∼108) that are densely distributed in an actual production field are considered for optimization. The number of MVs m is treated as an optimization variable, and each MV is specified to have 16 channels; in the optimization result obtained, four MVs are used. It can be seen from Figure 12(a) that the optimization result is reasonable.
(a)
(b)
In the second example, the positions of the 40 wells are randomly generated using the rand function (Table 3; W1∼W40), and the number of MVs is constrained to a range. The upper limit on the number of channels is uniformly set to 16 for each MV, and the problem is solved using the proposed algorithm. The well connections of the six MVs are shown in Figure 12(b). These optimization results are also reasonable.
5.4. Stability of the Algorithm
In this section, scatter plots are presented in Figure 13 to show the well and MV distribution results calculated using the proposed algorithm under the condition that the maximum number of channels per MV is constrained to 16 and the number of MVs is constrained to a certain range. The results prove that the algorithm can stably solve the MINLP problem and obtain reasonable results. Figures 13(a) and 13(b) show the result generated using the data based on W140∼W180 (Table 1), Figures 13(c) and 13(d) show the result generated using the data based on W67∼W128 (Table 1), Figure 13(e) shows the result generated using the data based on W1∼W80 (Table 1), and Figure 13(f) shows the result generated using the data based on W1∼W80 (Table 4).
(a)
(b)
(c)
(d)
(e)
(f)
5.5. Algorithm Performance
The scale of the computations performed in the proposed algorithm mainly depends on the size of the foresight matrix and the number of iterations. Under the assumption that all solutions are feasible, the maximum time complexity of the algorithm is , where m and n denote the numbers of rows and columns, respectively, of the matrix and N is the number of iterations. If the coordinate data are given in order, the time consumption is smaller, and if they are random, the time consumption is larger, as seen by comparing the examples with 191 wells in Table 7. The random data for 191 wells are shown in Table 5. For the examples with 26 and 41 wells, for the first example of each type listed in Table 7, the identity matrix is considered the initial feasible solution, and the computational scale is then reduced through multistep iteration. By comparison with the second example in each case, it is confirmed that the multistep iteration strategy can efficiently reduce the scale of the computations performed, thus saving time. However, regardless of the type of iteration used, the same optimal solutions can be found. In the second listed examples with 35, 60, and 80 wells, there is a specified upper bound on the number of channels for each MV, and the number of MVs is constant. The computation times are increased in these cases mainly because the single-step iteration process requires more iterations to approach the optimal solution. The solution times shown in Table 7 essentially meet our expectations with regard to the algorithm complexity. Thus, it is proven that the proposed algorithm can be used to solve various problem instances within a short time.
Finally, a large-scale problem with 192 wells (Table 1; W1∼W192) is considered, and the result is shown in Figure 14. Here, the channel count for each valve is constrained to 16, and the number of valves is allowed to vary from 15 to 39. It can be seen that the algorithm can still successfully separate different wells into different blocks in the face of a large number of wells with multiple valves. The results show that the algorithm can still effectively solve such a large-scale problem and yield reliable results. All of the above solution processes were implemented in MATLAB on an ASUS GX501VSK laptop (Intel Core i7-7700HQ Processor, 16 GB of DDR4 RAM).
6. Conclusion
This paper proposes the replacement of traditional metering processes with a process based on MVs, thus simplifying and automating the necessary operations and reducing the investment in redundant valve instrumentation and other equipment, the amount of space occupied by this equipment, and the cost of the station. A generalized MINLP model is presented, and a corresponding heuristic algorithm for MV optimization design is proposed for solving the presented MINLP problem. The numbers of channels of the MVs and the number of MVs are treated as special constraints in the MINLP problem.
The coordinates of actual production well sites and randomly generated data are used as examples for calculation and analysis. From the analysis of these examples, it can be seen that the behavior of the objective function value in the proposed algorithm is strictly convergent. A large collection of optimization results demonstrates that the algorithm shows stable convergence and good reliability. The analysis results show that the proposed heuristic algorithm can quickly and effectively solve the presented MINLP model for the design of a pipeline system with MVs. In this paper, we do not consider the optimization of the MV process itself in an integrated gathering system; this problem will be addressed in the future.
Nomenclature
| : | Binary variable. If , then a specific number of channels are given. By default, if , the number of channels of the jth MV is given by the manufacturer |
| : | Binary variable. If , then a specific lower bound on the number of wells to be connected to the jth MV is given, whereas a value of 0 indicates that the default lower bound is used |
| : | Number of channels of a fully connected MV |
| : | Number of channels of MV j |
| : | Matrix used as a compact representation of the connection matrix |
| : | Binary variable used to indicate MVs without any connected channels. If all elements in the jth row of are 0, then ; otherwise, |
| : | Total cost |
| : | Total cost of pipelines |
| : | Total cost of MVs |
| : | Total cost of the pipelines connected to the jth MV |
| : | Total length of pipelines |
| : | Length of the pipelines connected to MV j |
| : | Set of the total pipeline lengths obtained by placing MV j at all possible different well sites |
| : | Matrix representing the connection relationship between the wells and the valves |
| : | Element in the ith row and jth column of . indicates that the ith well is connected to the jth valve, and indicates that there is no connection between the ith well and the jth valve |
| : | Matrix with dimensions of used to express where the valves are placed. indicates that an MV is placed at the ith well site, and a value of 0 indicates no MV is placed at the corresponding well site |
| : | Connection matrix incorporating the information about the placement of the MVs |
| : | Penalty value |
| m: | Number of MVs |
| n: | Number of wells |
| : | Number of channels ultimately selected for MV j |
| : | Lower bound actually applied on the number of wells connected to MV j |
| : | Specified number of channels for MV j |
| : | Specified lower bound on the number of wells to be connected to valve j |
| : | Default lower bound on the number of wells to be connected to valve j |
| : | Maximum number of channels per MV that a manufacturer can provide and also the maximum value that can take |
| : | Set of well coordinates |
| : | Indicator that the ith well is connected to the jth MV |
| : | Set of MV coordinates |
| : | Objective function value obtained by moving the nonzero element in the jth column to the ith row |
| : | Foresight matrix |
| : | Minimum objective function value after N iterations |
| : | Lower bound on the number of MVs |
| : | Upper bound on the number of MVs |
| : | Cost coefficient matrix for pipelines |
| : | Cost coefficient matrix for MVs. |
Data Availability
The coordinate data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (51704253).