Statistical inference of the mixed linear model with incorrect stochastic linear restrictions

1.   Introduction

In the later stages of oil field development, water injection for enhanced oil recovery becomes a critical extraction method ^[1,2,3]. The purpose of the oil field water injection pipeline network is to distribute water from injection stations to various injection wells according to production needs, meeting the flow rate and pressure requirements of different wells. Oil field water injection systems typically cover areas of tens of square kilometers, and their power consumption generally accounts for about 40% of the total electricity consumption of the oil field ^[4,5]. Therefore, establishing and solving an energy consumption optimization model for the water injection system can reduce electricity consumption while meeting production requirements.

In solving the energy consumption optimization model for the water injection system, the calculation of node pressure is closely related to parameters such as pipe roughness and diameter. Currently, the selection of these parameters is based on the values at the time of pipeline installation. However, since the oil field water injection network is a high-pressure pipeline system with relatively small diameters, and the water transported is treated oily wastewater, the corrosion of these pipes is more severe compared to other networks, and the pipelines have been in place for a long period. Therefore, the pipe roughness coefficient and diameter may have changed significantly from their values at the time of installation, leading to considerable errors when using the installation data for simulation and optimization ^[6,7]. Thus, it is necessary to conduct correction studies on the pipe roughness coefficient and diameter for the oil field water injection network.

Extensive research has been conducted on parameter estimation for oil field water injection networks as water distribution pipeline systems. Three primary methods have been proposed ^[8,9]:

(1) Trial and Error Method: This method requires multiple manual repetitions of judgment and adjustment, resulting in very slow convergence and no guarantee of achieving the desired results (see ^[10,11]).

(2) Explicit Calibration Method: This method involves solving a series of extended steady-state hydraulic equations. However, it requires the number of calibration parameters to match the number of observational (field measurement) data, which is difficult to achieve in practice ^[12].

(3) Implicit Calibration Method: This method establishes an implicit model based on optimization techniques. It primarily estimates calibration parameters by using optimization algorithms combined with hydraulic simulation models to minimize the difference between observed and simulated results. This method is currently the primary approach and has been widely studied ^[13,14] due to its effectiveness in handling complex hydraulic systems. The calibration variables for these models include parameters such as nodal demand and pipe roughness ^[15]. Typically, the objective function of the model is the error between measured and simulated pressures. Various optimization methods have been employed to solve the related optimization models; however, these algorithms cannot guarantee obtaining the global optimal solution ^{[16,17,18,19,20]}. Although optimization techniques using genetic algorithms (GA) for model calibration have been proposed to achieve the global optimal solution ^{[21,22,23,24,25]}, these methods also cannot ensure obtaining the global optimal solution.

The problem of correcting the pipe roughness coefficient in oil field water injection networks shares similarities with that in urban water supply networks, but there are also notable differences. In oil field water injection networks, the corrosive effects of oily wastewater have significantly altered the pipe roughness coefficient and diameter compared to their values at the time of installation. However, due to objective production constraints, it is challenging to obtain multi-condition data for oil field water injection networks. Currently, there is relatively limited research on correcting the pipe roughness coefficient in these networks.

Wang et al. ^[26] applied methods from urban water supply network roughness coefficient correction to inverse research on oil field water injection pipeline roughness coefficients. She proposed methods based on graph theory, sensitivity analysis, neural networks, and particle swarm optimization. However, graph theory and neural network methods can only handle small, ideal networks and exhibit low accuracy in practical applications.

Wang et al. ^[27] investigated the issue of inaccurate empirical roughness coefficient values in oil field water injection networks. These inaccuracies can severely impact the operational efficiency of the network. They used orthogonal experiments and sensitivity analysis to identify and analyze pipe combinations that require precise adjustments to the friction coefficient. An improved method was demonstrated, calculating these factors through node equations to ensure more accurate and efficient operation of the water injection system. However, these methods may have limitations and might not fully reflect real-world conditions, as other, more complex factors or interactions might not have been considered.

Ren et al. ^[28] proposed a mathematical model for correcting the pipe roughness coefficient in oil field water injection networks. They used particle swarm optimization and simulated annealing algorithms for iterative optimization of the multivariable, multiparameter roughness coefficient correction problem. However, despite establishing an optimization model, obtaining satisfactory results is challenging due to the existence of multiple global optimal solutions.

Neither traditional optimization algorithms nor intelligent optimization algorithms can guarantee finding the global optimal solution to the optimization problem. More importantly, due to the underdetermined nature of the node equations, the established optimization model has multiple global optimal solutions ^[29], making the correction results of such models inaccurate.

The contribution of this paper is the presentation of a mathematical model for correcting pipe roughness coefficient under a single operating condition, along with an efficient numerical method for solving this model. Additionally, the established model has a unique solution.

This paper first presents a mathematical model for correcting the pipe roughness coefficient under a single operating condition, with changes in pipe diameter being attributed to changes in roughness, thereby reducing the number of parameters for hydraulic model calibration without affecting the hydraulic calculations of the pipeline network. The mathematical model is solved using single-condition data, reducing the dependence on multi-condition data, which is common in traditional pipe roughness coefficient correction optimization algorithms. Additionally, the model considers the roughness coefficient values at the time of pipe installation and limits the range of roughness coefficients within the model, making it more realistic. Second, by using matrix singular value decomposition ^[30,31], the optimization model's solution is transformed into a positive definite quadratic programming problem. Since the solution to a positive definite quadratic programming problem exists and is unique, it is demonstrated that the solution to the mathematical model also exists and is unique. Finally, the interior-point method is used to solve the model, ultimately obtaining the global optimal solution of the optimization model.

This study uses several key symbols and parameters which are defined in Table 1.

2.   Materials and methods

2.1.   Hydraulic models

The oil field water injection system consists of injection stations, distribution rooms, injection wells, and the connecting pipeline network, forming a complex and extensive fluid network. Typically, the number of nodes in the pipeline network can reach thousands, requiring extensive computations to solve the nodal pressure equations. To reduce the dimensionality of the system equations while retaining the essential characteristics of the original system, simplification strategies are employed. The simplified water injection network consists of main injection lines, injection stations, and pipeline intersections, forming a looped network ^[32].

According to the principle of mass conservation, for any given node, the inflow to the node equals the outflow from the node, thus satisfying the node flow balance. The mathematical expression is as follows ^[12]:

where, the index $i$ refers to the node number in the pipeline network, where $i = 1, 2, \cdots, n$. The index $j$ represents the node that is part of the same pipeline segment as node $i$. $Q_i$ represents the flow at node $i$, with inflow being negative and outflow positive. $q_{ij}$ denotes the pipe flow from node $i$ to node $j$, with its sign determined by the pressure difference between two nodes. When the pressure at node $i$ is greater than that at node $j$, $q_{ij}$ is positive; when the pressure at node $i$ is less than that at node $j$, $q_{ij}$ is negative.

The pressure drop equation in pipeline hydraulic calculations represents the relationship between pipeline flow and head loss, which can be expressed by the exponential formula below ^[12]:

where, $H_i$ and $H_j$ are the pressures at the two nodes $i$ and $j$ of the pipeline; $s_{ij}$ is the coefficient term; and

varies depending on the formula used. This paper adopts the Hazen-Williams formula, which is widely used in pipeline network calculations, and its form is as follows ^[12]:

where, $h_{ij}$ represents the head loss in the pipeline between nodes $i$ and $j$, measured in meters (m); $l_{ij}$ is the length of the pipeline, measured in meters (m); $d_{ij}$ is the diameter of the pipeline, measured in meters (m); and $C_{ij}$ is the Hazen-Williams coefficient, which is the roughness coefficient to be calibrated in this paper; The numerical values 1.852 and 4.87 are constants from the Hazen-Williams equation, derived from experimental data.

2.2.   Pipe roughness coefficient correction

The direct problem of hydraulic calculation for pipeline networks is defined as follows: given the flow rates at each node, pipe lengths, pipe diameters, and a reference point pressure in the network, and assuming the roughness coefficients of each pipe are known. In this case, the pressures at each node can be determined by simultaneously solving the pressure drop equations and the node flow balance equations. Subsequently, these node pressures can be used to solve the energy consumption optimization model for the water injection system ^[33,34].

Currently, the pipe diameters and roughness coefficients employed in hydraulic calculations for oil field water injection networks primarily rely on data from the time of pipeline installation. However, over time, some pipes may experience corrosion, scaling, and other issues, leading to changes in their diameter and roughness coefficients. This can introduce significant errors when using the installation data for simulation calculations and optimization. Therefore, it is necessary to correct the pipe diameters and roughness coefficients.

Conventional methods typically involve solving an optimization problem using multi-condition data, and the established models theoretically lack a unique solution, making it difficult to obtain a global optimal solution ^[35,36]. Therefore, this paper explores the development of an optimization model for correcting pipe roughness coefficients using single-condition data. To reduce the number of correction parameters, in the Hazen-Williams formula (2.3), we retain the pipe segment diameter as the value from the time of installation, attributing its changes to the roughness coefficient. Thus, only the pipe roughness coefficient needs to be corrected, simplifying the calculation process without altering the results of the hydraulic simulation. Therefore, the correction problem addressed in this paper is as follows: under single-condition data, given the known pressures at each node, how to use the pressure drop equations and the Hazen-Williams formula to solve for the pipe roughness coefficients.

3.   Model for correcting pipe roughness coefficients and numerical solution

The approach adopted in this paper significantly diverges from the conventional methodologies in the field. This process is divided into two steps rather than directly solving for the pipe roughness coefficients. First, the flow rate for each pipe is determined under the current operating conditions. Then, given the known pressures at both ends of the pipe, the pipe roughness coefficient is determined using Eq (2.3).

3.1.   Determining the pipe flow optimization model

Assume the pipeline network consists of nodes, with pressure values denoted as $H_1, H_2, \cdots, H_n$. The network comprises $m$ pipelines, with roughness coefficients denoted as $C_1, C_2, \cdots, C_m$ for each pipeline and let

In a looped pipeline network, the number of nodes $n$ is less than the number of pipelines $m$. Since

the node continuity equations are denoted as:

There exists at least one equation in this system that can be linearly represented by the remaining equations. It is necessary to remove one redundant equation. Without loss of generality, by removing the last equation, the new continuity equations are ^[37]:

Suppose the two nodes at the ends of the $k$-th pipeline are numbered sequentially as $i$ and $j$, and let

It follows that

The system of continuity Eq (3.2), where $q_1, q_2, \cdots, q_m$ are the unknowns, consists of $m$ equations.

For the sake of clarity and convenience in presentation, we will express Eq (3.2) in matrix-vector form. Define

and $\mathbf{A}$ as the coefficient matrix. Obviously, the order of $\mathbf{A}$ is $(n - 1)\times m$. With this setup, the Eq (3.2) can be written in matrix-vector form as:

We denote the known initial values of the roughness coefficients at the time of pipeline installation as

Next, we will study how to utilize the system of Eq (3.3), and in combination with the initial roughness coefficient values $\mathbf{C}_0$ and the range of roughness values

to establish an optimization model for solving the pipeline flow rates.

Taking the initial roughness coefficient $\mathbf{C}_0$ as the pipeline roughness coefficient during the calculation, we can utilize the node pressure equations to determine the pressure at each node. Based on the calculated node pressures, we can use the pressure-drop Eq (2.3) to calculate the corresponding initial values $\mathbf{q}_0$ of pipeline flow rates.

Generally, pipeline roughness coefficients have a range of values $[c_1, c_2]$. Define

and

Regarding the pipeline roughness coefficients, there is a constraint

Thus, by using the method described above to determine $\mathbf{q}_0$, we can obtain the pipeline flow rates $\mathbf{q}_{\min}$ and $\mathbf{q}_{\max}$ corresponding to $\mathbf{C}_{\min}$ and $\mathbf{C}_{\max}$, ultimately deriving the constraint conditions for pipeline flow rates as follows

The pipeline flow model established in this paper aims to find $\mathbf{q}$ that satisfies the Eq (3.3), the constraint condition (3.4), and minimizes $\|\mathbf{q} - \mathbf{q}_0\|_2$. Therefore, the mathematical model for the pipeline roughness coefficient inversion is:

Since Eq (3.3) has a solution, and the rank of matrix $\mathbf{A}$ is less than the number of unknowns, Eq (3.3) has infinitely many solutions, and it has the same solution as $\min \|\mathbf{Q} - \mathbf{Aq}\|_2^2$. Therefore, the mathematical model of the problem can be transformed into: finding $\mathbf{q}$ such that:

Define

and

The final mathematical model is to find $\mathbf{q}'$ such that:

If we solve for $\mathbf{q}'$, then

3.2.   Numerical solution algorithm for pipeline flow model

The matrix $\mathbf{A}$ is decomposed using singular value decomposition as:

Then,

where $\mathbf{U}$ is an orthogonal matrix of order $n-1$, $\mathbf{V}$ is an orthogonal matrix of order $m$, and $\Sigma_r$ is a diagonal matrix of $r$ order. Therefore, we have:

where $\mathbf{c}_1$ is a vector composed of the first $r$ elements of vector $\mathbf{U}^T \mathbf{Q}'$, and $\mathbf{c}_2$ is a vector composed of the last $n-r-1$ elements of vector $\mathbf{U}^T \mathbf{Q}'$. Based on the above decomposition, it is evident that when

and $\mathbf{y}_2$ is chosen arbitrarily, the resulting

is guaranteed to be a solution for minimizing $\|\mathbf{Q}' - \mathbf{A} \mathbf{q}'\|_2^2$.

Further research is needed to explore how to find $\mathbf{y}_2$ such that the corresponding $\mathbf{q}'$ can satisfy as follows:

Let $\mathbf{W}$ denote the first $r$ columns of matrix $\mathbf{V}$, and let $\mathbf{M}$ denote the last $m-r$ columns. Then

thus,

Therefore, the constraint is as follows:

That is

In summary, let

then $\mathbf{y}_2$ satisfies as follows:

Then, the corresponding

is the solution to problem (3.7). Since Eq (3.15) is a positive definite quadratic programming problem with linear inequality constraints, it has a unique solution and the solving process is relatively simple. Since the objective function is strictly convex and the constraints are linear, the problem is a convex optimization problem, which theoretically guarantees the existence of a unique global optimal solution ^[38,39,40]. In this paper, we employ the well-established interior-point method, which is widely recognized for its good convergence and stability properties. For positive definite quadratic programming problems, where the objective function is strictly convex and the constraints are linear, the interior-point method ensures global convergence to the optimal solution. After obtaining the solution $\mathbf{y}_2$ of Eq (3.11), the solution to problem (3.7) is obtained as

3.3.   Calibration of pipeline roughness coefficient

After obtaining the solution $\mathbf{q}'$ for the pipeline flow optimization model, the actual pipeline flow

can be determined. When the pressures at all pipeline nodes are known, $s_{ij}$ is determined by the following formula:

Upon determining $s_{ij}$, the roughness coefficient $C_{ij}$ for each pipeline can be calculated based on the following equation:

Therefore, the calculation procedure for correcting the pipeline roughness coefficient under a single operating condition is as follows:

Step 1: Utilize the known data to determine the coefficient matrix $\mathbf{A}$.

Step 2: Utilize the known node pressures and initial roughness coefficient $\mathbf{C}_0$ to determine the flow rate $\mathbf{q}_0$.

Step 3: Utilize the known node pressures, pressure drop equations, and the range of roughness coefficients to determine $\mathbf{q}_{\text{min}}'$ and $\mathbf{q}_{\text{max}}'$.

Step 4: Obtain the solution $\mathbf{q}'$ to problem (3.7).

Step 5: Set

Step 6: Determine the value of $s_{ij}$ using Eq (3.17).

Step 7: Determine the roughness coefficient $C_{ij}$ for each pipeline using Eq (3.18).

If the pressures at some nodes are unknown, the initial roughness coefficient can be used to estimate the node pressures. The estimated pressures can then replace the unknown node pressures, making the pressures at all nodes known. Subsequently, the roughness coefficients for the pipelines can be determined using the method outlined above. It should be noted that if the pressures at some nodes are unknown, the estimation accuracy may decrease.

4.   Case study calculation

Case study: This is a simplified real-world water injection network consisting of 7 injection stations, 98 nodes, 131 pipelines, and 34 loops. Nodes 17, 34, 42, 48, 66, 79, and 83 represent the locations of the injection stations. A simplified diagram of the water injection network is shown in Figure 1.

For detailed pipeline parameters, node parameters, and initial roughness coefficients used in the network design, see Tables A.1–A.3 in the Appendix. Since the true roughness coefficients of the pipelines are unknown, we modified the initial roughness coefficients of some pipelines to simulate the actual conditions of the network. These modified coefficients were regarded as the true roughness coefficients for the purpose of this study, serving as a benchmark to evaluate the performance of our proposed calibration method.

Table 2 shows the details of the roughness coefficients: the first column represents the node numbers, the second column shows the initial pipe roughness coefficients, and the third column represents the modified roughness coefficients, which are considered the actual roughness coefficients used in the study. Hydraulic simulations were then conducted using these modified coefficients to calculate the node pressures. The obtained node pressures were used as known values, and the proposed method was applied to calibrate the pipe roughness coefficients. Finally, by comparing the calibrated roughness coefficients with the true roughness coefficients, the effectiveness and accuracy of the method were assessed.

According to Table 2, the roughness coefficients of 20 pipelines in the network were adjusted. Using the new roughness coefficients and the basic data of the network, we solved the nodal pressure equations to obtain the pressure values at each node.

With all node pressures known, we determined the roughness coefficient of each pipeline segment. Using the initial roughness coefficients, known node pressures, pipeline segment radii, lengths, and node flow data, the roughness coefficients of the pipeline segments were calculated using the correction method proposed in this paper, and the roughness coefficient correction results were obtained (see Table 3). The first column of Table 3 represents the node numbers, the second column shows the initial roughness coefficients, the third column represents the true roughness coefficients, and the fourth column shows the roughness coefficients calculated using the method proposed in this paper.

The calculation results indicate that if the actual roughness coefficient of a pipeline is equal to the initial roughness coefficient, the roughness coefficient obtained using the proposed method equals the actual roughness coefficient. Therefore, the correction results for the pipeline segments where the actual roughness coefficient differs from the initial roughness coefficient are presented.

After the calculation, the average error between the initial roughness coefficients and the actual roughness coefficients is 7.08%, while the average error between the roughness coefficients using the proposed method and the actual roughness coefficients is 2.18%, representing a reduction of 4.9%. When solving the hydraulic calculations and operational optimization problems, both of which require the use of roughness values, the roughness coefficients obtained by the proposed method yields better results compared to using the initial roughness coefficients.

5.   Discussion

(1) This paper establishes a constrained least squares mathematical model for calibrating the pipe roughness coefficients in oil field water injection networks under a single operating condition and proposes a global optimal solution method to address the rank-deficient least squares problem with constraints. Using this method, satisfactory results were obtained by simulating the calibration of the roughness coefficients for a large real-world pipeline network. The average error between the calibrated roughness coefficients and the actual roughness coefficients is 2.18%. The remaining error is due to the fact that the theoretical solution of the model established under a single operating condition does not necessarily guarantee the actual solution, which is a limitation of using single-condition data for roughness coefficient inversion.

If the proposed method is applied to multi-condition data, each additional condition would add another least squares constraint, significantly increasing the dimensionality of the mathematical model. This would lead to a substantial computational load when performing singular value decomposition. The key issue here is to design an effective block-diagonal matrix singular value decomposition algorithm. This is a crucial challenge for future research. If an efficient block-diagonal matrix singular value decomposition can be implemented, the proposed method could be extended to multi-condition data, potentially improving the accuracy of the pipe roughness coefficient further.

(2) This paper does not compare with other roughness coefficient correction methods, as other models typically require multi-condition data, which is not easily obtainable for oil field water injection systems. This paper primarily focuses on developing a mathematical model for calibrating pipe roughness under a single operating condition and exploring how to efficiently solve the model. Currently, mainstream methods are mainly focused on developing and solving calibration models for pipe roughness under multi-operating conditions. However, when handling single-condition data, existing methods theoretically lead to a multi-solution optimization problem, making them unsuitable for single-condition data.

(3) The method proposed in this paper requires all node pressures to be known. If some node pressures are unknown, the initial roughness coefficients can be used to estimate node pressures once, and the estimated pressures can replace the unknown values. Thus, all node pressures become known, although some pressures will be approximations rather than measured values. The method can then be used to correct the roughness coefficients, though the accuracy may decrease. Alternatively, fuzzy optimization techniques can be explored to handle situations where some node pressures are unknown.

6.   Conclusions

This paper establishes a mathematical model for correcting the pipe roughness coefficient in oil field water injection networks under a single operating condition. The solution to the mathematical model is unique, and by using matrix singular value decomposition, the original problem is transformed into a positive definite quadratic programming problem with linear inequality constraints, thereby obtaining the global optimal solution. In comparison with conventional optimization methods for determining the roughness coefficient, this method requires only one matrix singular value decomposition and solving a positive definite quadratic programming problem, thereby increasing computation efficiency. Simulation of roughness coefficient correction for a large real network shows that the average error between the corrected roughness coefficients and the actual roughness coefficients is 2.18%, while the average error between the initial roughness coefficients and the actual roughness coefficients is 7.08%, representing a reduction of 4.9%. This demonstrates the model's validity and the effectiveness of the calculation method.

Author contributions

Yuxue Wang: conceived the study, designed the methodology, performed the data analysis, and contributed to writing the manuscript; Songyu Bai: contributed to the study design, conducted the case study calculations, and assisted with the interpretation of results and manuscript revisions. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix