Optimization of load sharing for parallel compressors using a novel hybrid intelligent algorithm

Compressor stations, which usually consist of multiple compressors in parallel, are installed to power natural gas travel in pipelines. Compressor station optimization, which should be expressed as a mixed integer nonlinear programming (MINLP) problem, makes economic sense for the entire gas transmission system. However, it has often been simplified as a nonlinear programming (NLP) or mixed integer linear programming (MILP) problem in previous research. Most of existing solutions are based on discretization and a genetic algorithm (GA). This paper addresses the general MINLP problem for compressor station optimization without simplification; a novel hybrid intelligent algorithm is proposed to solve this problem. The proposed algorithm, DWOA, leverages advantages of the whale optimization algorithm (WOA) and differential evolution (DE). The proposed algorithm can balance exploration and exploitation to find the global optimal solution. An approach to handling constraints is also presented, where the original problem model is reformulated to be continuous by expanding the flow rate range of the compressor. A case study is performed to illustrate the performance of this approach. Results show that the continuous reformulated model is easier to solve, and DWOA produces a satisfactory solution that differs from theoretical results by only 1.61%. In addition, DWOA demonstrates better accuracy and stability than WOA, DE, and DE‐WOA, another hybrid algorithm. Therefore, this solution has the potential to promote comprehensive compressor station optimization.


| INTRODUCTION
Pipelines play an important role in oil and gas transportation. With the sharp demand for oil and gas consumption, pipelines have expanded throughout the world. 1,2 The concerns about global warming and shale gas revolution have also boosted the expansion of gas pipelines in recent years. 3,4 The number of long-distance oil and gas pipelines is approximately 3800 globally, of which the gas pipelines take a large proportion in terms of length. 5,6 For pipeline transportation of fluids (oil and gas), pressure gets lost due to friction and heat transfer. Pump or | 331 LI et aL.
compressor stations are typically set up along a pipeline to compensate the pressure loss. [7][8][9] The energy consumption of such stations is responsible for a considerable proportion of operating costs in the pipeline system. [10][11][12] This paper focuses on the optimal operation of the compressor station to reduce pipeline companies' budgets. In general, there are several compressors connected in parallel to increase the transport capacity in each station. The power required by the compressor to boost the gas is delivered by electrical drivers or gas turbines and changes with the gas flow rate and inlet conditions. 13,14 The air cooler installed at the compressor outlet is not considered in this work since it consumes much less energy than the compressor. 15 Compressor station optimization is often taken as a subproblem of natural gas pipeline network optimization. Generally, such studies are based on the assumption that all compressors within a station are of the same type and exhibit the same performance. [16][17][18][19][20] Alternatively, the station is treated as a whole, and the details of each compressor are ignored. [21][22][23][24][25][26][27] However, the performance of each compressor within a station varies, occasionally even for compressors of the same type. The optimization solution of a pipeline network is also impractical in the absence of a specific operational station scheme.
In light of these factors, several studies have recently considered compressor station optimization. Most have attempted to solve this problem via classical mathematical methods or by simplifying the model. Paparella et al 28 introduced an MINLP model where the switching costs of compressors were included in the objective function; specifically, the branch and bound method (BB) was employed in their study. Deng 29 considered two approaches based on discretization of the volume flow rate: mixed integer linear programming (MILP) and dynamic programming (DP). However, it is time-consuming to obtain a sufficiently accurate solution that depends on the discretization resolution. In addition, a simplified nonlinear programming (NLP) model without considering the on/off state of each compressor has been widely adopted. For example, Xeons et al 30 investigated the real-time optimization of compressors in parallel based on a fixed number of active compressors. Cortinovis et al 31 proposed a two-step framework to optimize load sharing among parallel compressors; in their framework, online performance tracking was the foundation of optimization, and the on/off state of each compressor was determined through performance tracking rather than a binary variable in the problem model. Later, Kumar and Cortinovis 32 extended this framework to compressors installed serially based on the assumption that the active machine is known.
Conversely, as an effective way to solve optimization problems, the meta-heuristic algorithm offers another option. However, this algorithm has generally remained limited to the genetic algorithm (GA). Nguyen et al 33 compared GA, MILP, and expert systems, ultimately suggesting that GA might be useful for large-scale problems but did not demonstrate superior performance to MILP in their cases. Based on a fixed number of active compressors, Hawryluk et al 34 investigated multi-objective optimization consisting of the minimum energy and maximum throughput of a compressor station using GA. Mahmoudimehr and Sanaye 35 compared GA with the exhaustive search method (ES) in cases with similar and dissimilar compressor units; findings revealed better performance of GA than ES in the latter case.
Based on a thorough literature review, this paper aims to present an accurate and reliable algorithm to solve the nonlinearity and nonconvex MINLP model for compressor station optimization. The main contributions of this paper are as follows: The remainder of this paper is organized as follows: Section 2 provides formulations associated with the general MINLP model for this problem. Section 3 introduces the model solution based on the proposed DWOA algorithm in light of constraint preprocessing. Section 4 outlines the results and discussions of a case study. Finally, Section 5 presents the main conclusions of this paper.

| Objective function
Power consumption generally emphasizes the optimization problem for a natural gas station. The objective function considered here aims to minimize the station power consumption P : In practice, there is often more than one compressor unit in a station, and units are usually arranged to form (1) minP a parallel topology (Figure 1). A compressor unit comprises three pieces of equipment: a driver, compressor, and cooler. The driver provides power for the compressor and can be considered an integral when analyzing power consumption. The power consumption of the cooler (typically an air cooler) is much smaller than that of the compressor; the cooler's power consumption is thus not considered in this paper. The station objective function can thus be rewritten as where P N,i is the power consumption of compressor i [MW], and n is the total number of compressors.

| Load balance
For compressors in parallel, all compressors are subject to the same suction conditions. The pressure ratio provided by these compressors should also be equal, 32 namely where 0 is the station compressor ratio, i is the compressor ratio for the ith compressor, and i is a binary variable representing the state of compressor i; "0" indicates the compressor is on, and "1" indicates the compressor is off.
In addition, the gas flow rate through all compressors should be equal to the total flow rate of the station, such that where Q 0 is the total flow rate [m 3 /s], and q i is the flow rate for compressor i [m 3 /s].
Note that it is common to set a cycle bypass for each compressor to avoid surge instabilities, as illustrated in Figure 1. However, surge instabilities are not considered in this paper, and the backflow is ignored in this load balance equation; these omissions make the problem formulations simpler and do not affect the universality of the solution algorithm. Moreover, only the electrical driver is considered here. This simplification helps draw attention to the proposed solution strategy, which can easily be extended to scenarios involving gas drivers.

| Compressor unit limitations
Performance governing equations for a typical natural gas compressor can be expressed as shown in Equations (5)- (11). 19,36 Gas compression through a compressor is a polytropic process. The polytropic head H can be calculated using Equation (5) once the suction temperature, suction pressure, and compression ratio are given: where H is the polytropic head [J/kg], Z s is the compression factor of natural gas at the suction side, R is the gas constant [J/ kg·K], T s is the gas temperature at the suction side [K], and k is the polytropic exponent, is the compression ratio.
By contrast, the polytropic head can be formulated as the function of the compressor rotation speed N and gas flow rate q: , q is the gas flow rate at the compressor inlet [m 3 /s] Combining the two equations of the polytropic head, the rotation speed can be determined based on the given suction conditions and compression ratio as well as the assumed flow rate.
Further, the corresponding polytropic efficiency can be calculated by the rotation speed and flow rate as follows: where b 4 , b 5 , b 6 are constants, is the polytropic efficiency.
Thus, the power consumption P N can be obtained using the following equation: where m is the mass flow rate of natural gas passing through the compressor [kg/s], P N is the power consumption of the compressor [MW].
The rotation speed of a compressor is limited to a certain interval as follows: where N min and N max are the minimum and maximum allowable rotation speed [rpm], respectively.
Besides the rotation speed limitation, the compressor capacity is also limited by two unstable situations (ie, surge and stagnation). The surge line is expressed as where a 1 , a 2 , a 3 are constants, Q surge is the minimum gas flow rate decided by surge [m 3 /s].
The stonewall line satisfies the following condition: where a 4 , a 5 , a 6 are constants, Q stonewell is the maximum gas flow rate decided by stonewall [m 3 /s].
In the H-q plan, four lines (ie, the minimum and maximum rotation speed lines, surge line, and stonewall line) define a closed region, denoted by the blue region in Figure 2. This region represents the working domain of a compressor. Working points with an equal compressor ratio will appear along an approximately horizontal line through the working domain, indicated by the dotted line in Figure 2. Accordingly, the flow rate range can be obtained based on the intersections of this horizontal line and boundary lines.

| Summary of optimization problem
Overall, the complete model for compressor station optimization can be expressed as where q min and q max are the minimum and maximum flow rate, respectively. Note that this model is based on the assumption that the suction conditions and compressor ratio are given and that all parameters in the compressor governing equations are available. Thus, the objective function and range of the flow rate for each compressor can be determined from the governing equations of the compressor.

| Model analysis and processing
This problem can be classified as a mixed integer nonlinear programming (MINLP) problem because the objective function is nonlinear and contains a binary variable. The problem is constrained by the load balance of the station and the compressor working domain under the given suction conditions and compressor ratio; that is, the candidate solutions obtained by random searching during the iteration process may be infeasible. This model therefore cannot be solved directly using conventional methods, in which case constraint handling is necessary. Moreover, the compressor's flow rate range depends on the compressor ratio required. The flow rate for each candidate solution can be enforced into the required interval during the solution process. Therefore, the primary focus in this paper is on the constraints of flow balance and the binary variable indicating the on/off state of the compressor.

| Fitness function
In this paper, flow balance is considered by building a new fitness function to replace the original objective function. According to the literature, 37 the fitness for feasible solutions is equal to the inverse of the objective function, whereas the fitness for infeasible solutions is equal to the inverse of FM plus the constraint violation value. The fitness function adopted herein is as follows: where FM is a fairly large constant. Therefore, feasible solutions are compared according to the value of the objective function of each, while infeasible solutions are compared according to their constraint violation values. The fitness value of an arbitrarily feasible solution is always larger than that of an infeasible one.

| Continuous reformulations
For a problem without a priori knowledge about potential solutions, the probabilities of each compressor being on and off are supposed to be theoretically equal, meaning that the binary variable δ has a 50% probability of being 0 and a 50% probability of being 1. However, the candidate volume flow rate for each compressor is uniformly distributed within the interval 0 ∪ q min , q max in the random search process; the likelihood of being 0 is quite small, which does not align with the original problem. Consequently, a single feasible solution may be unavailable when the initial population of the algorithm is insufficiently large, especially for a large-scale problem. To address this problem, a strategy based on the extension of the flow rate range is proposed instead of applying a binary variable. A symmetrical extension of the flow rate range is used in this paper, and the original flow rate range is replaced by the interval q min − 1 2 (q max − q min ), q max + 1 2 (q max − q min ) during the iteration process. Combined with the fitness function, the original problem can be rewritten as WOA was proposed by Mirjalili and Lewis, 38 documented as a powerful meta-heuristic algorithm inspired by the behavior of humpback whales. WOA includes an exploration phase and exploitation phase. Notably, there are two kinds of mechanisms to update the position in the exploitation phase: the shrinking encircling mechanism and spiral updating position. A random mechanism is adopted to search for prey during the exploration phase. The algorithm is mathematically formulated as follows.

A
Exploration phase (search for prey) In this phase, search agents (whales) update their positions randomly according to each other's position. This behavior is modeled as where � � → X is the position vector to be updated, is the position vector of a random whale chosen from the current population, t is the current iteration number, � � → A and � � � → C are coefficient vectors, | | is the absolute value operator, and .is the operator that implements element-by-element multiplication.

A Exploitation phase
In this phase, search agents tend to update their positions to the best search agent (ie, optimum candidate solution) in the current population with two alternative mechanisms.
• Shrinking encircling mechanism The following formulations are used to model the process of whales' shrinking encircling of prey: max Fitness A is greater than 1 in the exploitation phase. This is the condition necessary to switch between these two mechanisms.

• Spiral updating position
To mimic the helix-shaped movement of humpback whales when attacking prey, the following spiral equation is given: where b is a constant defining the shape of the logarithmic spiral, and l is a random number in [−1, 1].
Furthermore, switching schemes are designed among the three above mechanisms. Thus, the algorithm concludes as follows: where p is a random number.
Note that WOA places more emphasis on exploitation because the balance between exploration and exploitation depends on the linearly varying vector � � → A. Only one element of � � → A is effective for each dimension of an agent. Therefore, the following calculation is based on unidimensional analysis, and vector symbols are dismissed. Equation (19) can then be rewritten as Thus, the probability of performing Equation (15) can be calculated as follows: The likelihood of exploring the entire search space is fairly low. Although this circumstance can facilitate convergence and accelerate computation, it tends to promote trapping in local optima.

| Differential evolution (DE)
DE is a highly competitive algorithm in solving multimodal optimization problems. It includes three operators: mutation, crossover, and selection. Formulations for the basic differential evolution algorithm are as follows. 39

A Selection operator
The trial vector is not the terminal offspring individual unless its fitness value exceeds that of the original. The operator can be expressed as follows: For the sake of visualization, the process of implementing DE is depicted in Figure 3.

| Proposed DWOA
Based on the above analyses, there is a risk of prematurely converging to local optima when implementing WOA. Conversely, DE has an advantege in exploring the entire search space rather than exploiting the local region around the optimum. Therefore, a hybrid algorithm is proposed in this paper by combining WOA and DE, otherwise called DWOA. Exploration and exploitation each occur in the earlier stage of WOA, whereas only exploitation is considered in the later stage. Thus, DE is embedded into the second half of the evolutionary process of WOA. The specific operation seeks to generate the trial vector by conducting the crossover operator between the mutant vector from DE and the vector from WOA under the condition that |A| < 1 and p < .5. Then, the trial vector replaces the vector obtained by the shrinking encircling mechanism of WOA if its fitness is better. This step improves the ability of WOA to escape from local optima without much decline in computing speed because DE is performed, but not in every iteration. The proposed DWOA therefore includes the advantages of both algorithms and offers better performance to balance exploration and exploitation. The main DWOA procedure appears in Figure 4.
According to Figure 4, the computational complexity depends on three procedures: initialization, updating, and fitness evaluation.   N × max (D, f (D))). In addition, there are more temporary variables and parameters in DWOA since it is a combination of WOA and DE, resulting a slightly increase in the space complexity.

| Case description
A compressor station with six compressors is considered in this case study. The six compressors can be divided into four types: one A-type (Compressor #1), three B-type (Compressors #2-4), one C-type (Compressor #5), and one D-type (Compressor #6). The total volume flow rate is 14 m 3 /s under the suction conditions of 293.15 K and 3.3 MPa, approximately 37.2 million standard cubic meters per day (MMSCMD). The compressor ratio is taken as 1.5. Parameters for model compressors under this suction condition are listed in Table 1.

| Effectiveness of continuous reformulation
The exhaustive search method (ES) was used to solve the traditional MINLP model. The result of ES was taken as the standard to test the reliability of the proposed solution strategy because ES is certain to reach the global optimal solution as long as discretization is sufficiently precise. 35 Here, the step length of the volume flow rate in ES is 0.0001 m 3 /s. The original problem model and the continuous reformulated model were solved simultaneously using the basic WOA first. In both cases, the initial population size was 100, the number of iterations was 1000, and the calculation was repeated 50 times. Findings are shown in Table 2.
No feasible solutions were obtained from the 50 calculations for the original model, while the average minimum power consumption obtained from the continuous reformulated model was close to the standard value. These results indicate the effectiveness of the continuous reformulation of the original model. Accordingly, the following calculations are based on the continuous reformulated model.

| Performance of different algorithms
To highlight the superiority of DWOA, comparisons with basic DE, WOA, and DE-WOA from the literature 40 were performed. Calculations were repeated 50 times for each algorithm with an initial population of 100 and 1000 iterations per calculation. The parameter settings of optimization algorithms are listed in Table 3, and the results are shown in Table 4 and Figure 5. Four error metrics in Figure 5 are calculated as follows: where V represents the real value and V′ denotes the predicted value.
In addition, the computation speed and convergence characteristics were also compared, as shown in Table 5 and Figure 6, respectively.
According to statistics from the 50 calculations, the average value, maximum value, and minimum value of minimum power consumption obtained by DWOA were closest to the standard value. Figure 5 also indicates that the proposed DWOA demonstrated better accuracy and stability given lower values of its four metrics. Although DE and WOA took less time than DWOA, the accuracy of DE was lower and WOA showed a tendency of premature convergence in Figure 6, which might lead to a local optimal solution. The combination of DE and WOA in DWOA was more efficient than in DE-WOA, as evidenced by the lower power consumption and reduced computation time of DWOA. Overall, DWOA performed high competitively in finding the global optimum with reasonable time consumption. The superiority of DWOA is attributed to the combination mode of DE and WOA. Despite the good performance of DE and WOA in many applications, [41][42][43][44][45][46] there are some intrinsic drawbacks of both algorithms. The low probability of exploration calculated as Equation (27) suggests a rapid diversity loss of the WOA population, leading to prematurely converging. 47 DE suffers from the problem of stagnation, where the population does not converge to any point. 48 On the positive side, however, it benefits the population diversity. Therefore, the two algorithms are complementary to each other. In the earlier stage of DWOA, WOA operator is employed. The global exploration and local exploitation are performed randomly, showing a similar rapid convergence characteristic to WOA. In the later iteration process, DE is introduced into the phase of shrinking encircling of WOA. The new individual generated from three randomly selected parents using Equation (28) allows searching in a direction that deviates from the current optimal solution. In this way, DE contributes to jump out of the local optimum while WOA emphasizes on the local exploitation to improve the accuracy of the solution. A good balance of exploration and exploitation is achieved.

| Optimal operation scheme
The optimum values of decision variables obtained by DWOA and ES are listed in Table 6.
Compressors #2-4 are of the same type; therefore, the optimal unit commitment from DWOA is essentially the same as that from the exhaustive method. Specifically, four compressors (one A-type, two B-type, one C-type, and one D-type) were turned on in this case.

| CONCLUSIONS
This paper describes a model and solution to the optimization operation of a compressor station with several compressors in parallel. Specifically, the general MINLP model for a natural gas compressor station was reformulated as a continuous model by expanding the compressor's working domain. The superiority of the continuous reformulated model was demonstrated through a case study. Further, an improved version of WOA was proposed to solve this model. The enhanced algorithm, named DWOA, integrates the advantage of DE exploration to avoid prematurity due to excessive focus on exploitation in WOA. By comparing the results of DWOA with those of the original WOA and DE as well as another hybrid version (DE-WOA), the accuracy and stability of DWOA were verified. In the case presented in this paper, the optimal power consumption obtained by DWOA showed an error within 2% compared to the standard value. The optimal unit commitment was essentially identical to the standard.