Two-stage robust power cost minimization in a natural gas compressor station

Yize Meng,Ruoran Chen,Keren Zhang,Tianhu Deng

Department of Industrial Engineering,Tsinghua University,Beijing,100084,China

Keywords:Natural gas Single station power minimization Nonconvex robust optimization C&CG algorithm

ABSTRACT Optimal operation of a compressor station is important since it accounts for 25% to 50% of a company's total operating budget.In short-term management of a compressor station,handling demand uncertainty is important yet challenging.Previous studies either require precise information about the distribution of uncertain parameters or greatly simplify the compressor model.We build a two-stage robust optimization framework of power cost minimization in a natural gas compressor station with nonidentical compressors.In the first stage,decision variables are the ON/OFF state of each compressor and discharge pressure.The worst-case cost of the second stage is incorporated in the first stage.Firststage decision variables feasibility is discussed and proper feasibility cuts are also proposed for the first stage.We employ a piece-wise approximation and propose accelerate methods.Our numerical results highlight two advantages of robust approach when managing uncertainty in practical settings:(1) the feasibility of first-stage decision can be increased by up to 45%,and (2) the worst-case cost can be reduced by up to 25% compared with stochastic programming models.Furthermore,our numerical experiments show that the designed accelerate algorithm has time improvements of 1518.9%on average(3785.9% at maximum).

1.Introduction

It is predicted that natural gas production of China will maintain a rapid growth(Liang et al.,2019)with peak gas of 323 billion cubic meters a year coming in 2036 (Li et al.,2016).Oil and gas tend to still play an important role in the coming three decades(Pan et al.,2020;Zhao et al.,2019;Huang et al.,2019).With a 4.6%increase in consumption in 2018,it accounted for nearly half of the increase in global energy demand according to the report of International Energy Agency (2019).Natural gas pipeline transmission is one important way to transmit natural gas.As natural gas travels along a pipeline over long distances,its pressure is reduced by the friction between itself and the pipeline wall.To compensate for the pressure loss,natural gas pressure is increased by compressing in a compressor station.

The optimal operation of a compressor station is important,since its operating cost accounts for 25%to 50%of a company's total operating budget(Luongo et al.,1989).We refer the reader to Ríos-Mercado and Borraz-S′anchez (2015) for a state-of-the-art survey on natural gas fuel cost minimization problems.In particular,the authors mentioned that one of the research challenges is to design stochastic models to handle cases when the variation of demand and supply is significant.

We consider operation management of a single compressor station facing short-term demand uncertainties.Many short-term factors may cause uncertainty in the total flow passing through the compressor station.For example,on the customer demand side,natural gas consumption can be largely influenced by weather,fuel switching,price and market variability.Thus,customer demand often exhibits large daily variations (Hellemoet al.2012) that are difficult to predict accurately.Other factors include unpredictable events due to malfunctions and failures and shortfall in upstream and downstream capacity (Ríos-Mercado and Borraz-S′anchez,2015).

To cope with the uncertainty issue in natural gas networks,in most of the current studies,assumptions are made about the probability characteristics.Methods such as stochastic programming (Padberg and Haubrich,2008;Ding et al.,2017) and chanceconstrained programming (Wintergerst,2017;Behrooz,2016) are utilized.In stochastic programming,the uncertainty is assumed to have a probabilistic description.Chance-constrained programming allows constraints to be violated with a specific probability.Both streams of work require heavily on the characteristics of the uncertain parameters.A third method is to characterize the load uncertainties in natural gas networks with multiple scenarios(Behrooz and Boozarjomehry,2017).

We adopt a two-stage robust optimization method for the following two reasons.First,uncertain parameters are assumed to belong to one uncertainty set.Therefore,it is more applicable in practical scenarios where the full distributional information of uncertain parameters is unknown or very difficult to determine accurately.Second,two-stage robust optimization overcomes the shortcomings of over-conservatism in single-stage robust optimization.This method has been widely applied to electricity power system,refinery,chemical,but has few applications in the natural gas field.

The model has two stages.In the first stage,the discharge pressure and ON/OFF state of each compressor are decided.Then the uncertainty is revealed.The uncertainty mainly refers to the flow rate passing through the compressor station.Also,suction pressure of the compressor station is uncertain,since it depends on the flow rate realization and the parameters of the previous station.Based on the realized uncertainty,the second stage decides the rotational speed of each compressor.The objective is to minimize the compressor station's worst-case power cost in the second stage.

To solve the problem,column-and-constraint generation(C&CG) algorithm (Zeng and Zhao,2013) is utilized,which is proved to be efficient for two-stage robust optimization problem.It composes of a master problem and a sub-problem.We combine it with a gradient-ascent searching scheme in solving the bi-level sub-problem.

We contribute the literature in four ways:

1.We formulate a two-stage robust power cost minimization model in a compressor station.To the best of our knowledge,we provide one of the first studies on robust optimization with detailed natural gas compressor model formulation.Natural gas managers strongly suggest that when the solution directs the operation,the model should be as accurate to the real case as possible (Deng,2015).

2.The convexity/concavity of the constraints are examined.We provide conditions under which the first-stage decisions are feasible.The condition is then formulated as a feasibility cut to reduce search space.For the sub-problem,a gradient-ascent algorithm combined with mixed-integer linear programming models is proposed.As for the master problem,an accelerate algorithm for original piecewise linear approximation with same optimal value is provided,which is also applicable in twostage stochastic programming model.

3.We conduct a set of numerical experiments to verify the benefits of the robust model compared to the deterministic one.The results show that feasible possibility of first-stage decision variables is significantly increased and worst-case cost is reduced.However,extra power cost is the price to gain robustness and the price increases with uncertainty set size.

4.The proposed robust approach and analysis can offer insights on the trade-off between robustness and operational cost.According to the feedback of natural gas managers,adjustment of key parameters such as rotational speed of each compressor as flow rate changes is reasonable and necessary.Other less critical parameters should avoid frequent adjustments.Current adjustment strategy is mainly qualitative and lacks quantitative analysis.The proposed model can help them find better strategy through quantitative analysis.

In Section 2,we review the relevant literature.Section 3 describes the deterministic and robust optimization model.In Section 4,we discuss solution approaches for the robust optimization problem.Numerical experiments are studied in Section 5 to evaluate the algorithm and compare the robust optimization results with the deterministic ones.In Section 6,we summarize our conclusions.

2.Literature review

The deterministic model has nonconvex objective function and feasible region,making this problem difficult to solve.Methods to solve the deterministic fuel cost minimization problem include dynamic programming,gradient search,geometric programming,particle swarm optimization algorithm (Zheng and Wu,2012) and linear approximation.Among these methods,particle swarm optimization algorithm and gradient search can not guarantee optimality.Dynamic programming and piecewise linear approximation of MILP guarantee optimality but often demand computing resource.Ríos-Mercado and Borraz-S′anchez (2015) gave an overview of related methods as well as other relevant optimization problems arising in natural gas networks.

Approximation with a piecewise linear function is widely researched in deterministic formulation.De Wolf and Smeers(2000) approximated the nonlinear relations between flow rate and pressure and used an extension of the simplex method to solve the problem.Geißler et al.(2015) incorporated heat power constraints,natural gas-mixing constraints,and detailed physical compressor working domain constraints.They guaranteed convergence to a partially optimal value.Martin et al.(2006)approximated the model without binary variables,and the subpolyhedra associated with the model have a computationally tractable number of vertices.

Another stream in deterministic fuel cost minimization problem utilizes dynamic programming (DP) methods.Since it guarantees optimality and meets the requirement of high solution accuracy from practitioners.Wong and Larson(1968)first applied DP models to pipeline networks.Borraz-S′anchez and Haugland(2011)enabled a significant speed-up for the DP formulation as they adopted an adaptive discretization scheme.Deng et al.(2019) approximated the solution using state dimension reduction and neighborhood searching and significantly accelerated the calculation of the DP model.Zhang et al.(2020) considered the transient-state natural gas transmission problem in which adjacent periods' natural gas flow had a relationship.They utilized an approximate dynamic programming(ADP)method to solve the problem.In this paper,we mainly focus on the nonlinear programming formulation for steady-state networks and related methods.

Considering uncertainty in natural gas operation is challenging since the problem is even more complicated than deterministic model,which is already quite difficult.Gotzes et al.(2016) took a step in studying the probability of a feasible load in a passive network with random demands.Feasibility here means that the existing flow meets the loads while the pressure in the pipelines is within given bounds.In addition,there are current three growing streams of optimization under uncertainty:stochastic programming,chance-constrained optimization and robust optimization.The first two methods require assumption about probability information of uncertainty distribution and the probability may significantly affect the outcome.These methods put high demands on the accuracy of probability prediction.

At the optimal solution of chance-constrained programming,constraint violation is allowed to a certain degree.Wintergerst(2017) applied a chance-constrained optimization methodology to a natural gas network.He considered uncertain boundary flows and nonconstant compressibility factors of the quasi-linear isothermal Euler equation.He showed that medium-sized network optimization can be solved in reasonable time.Behrooz and Boozarjomehry considered the demand uncertainty and utilized a chance-constrainted programming formulation to decide compressor stations discharge pressure in steady state (Behrooz,2016) and transient state (Behrooz and Boozarjomehry,2017).They characterized the user load with a Gaussian random variable and assumed the mean and variance were known.

In stochastic programming models,the uncertainties in demand and supply are represented by several scenarios,and the forecast quality influences the accuracy of optimal solutions.Padberg and Haubrich (2008) optimized natural gas portfolios of commercial enterprises under demand and price uncertainties.They generated a tree consisting of all possible scenarios with defined probabilities.Zavala (2014) utilized stochastic optimal control model with scenario based uncertainty formulation.Ding et al.(2017) proposed a multi-stage stochastic programming model for the expansion coplanning of natural gas and power networks,and the uncertainties in demand were considered.

Compared with stochastic optimization and chance-constrained programming,robust optimization dose not require probability distribution assumption,which is more applicable in practical scenarios.Gabrel et al.(2014) and Gorissen et al.(2015) gave comprehensive overview of recent developments in robust optimization.Two-stage robust optimization allows some decision variables to be adjusted after the uncertainty realizes.C&CG method is an efficient algorithm firstly proposed by Zeng and Zhao(2013)and is widely used in two-stage robust optimization in other fields.

It is worth pointing out that there are only a few researches relevant to robust optimization application in natural gas operation field.In Aßmann et al.(2018),they managed to answer whether,for any possible realization of the uncertainty in demand,there would be a feasible adjustable second-stage decision that can be satisfied by the network.They also formulated a two-stage robust optimization model and exploited the decomposable problem structure and gained large speedups (Aßmann et al.,2019).However,either the cost function was simplified as linear or the working domain constraints were not considered.Most details of the original problem were omitted.

Different with Aßmann et al.(2019),we take into account detailed compressor model in a compressor station with robust optimization.The work is most closely related to the study by Deng(2015),who compared dynamic programming and MILP in optimization of single compressor station.Mahmoudimehr and Sanaye(2014) also considered deterministic optimization in a single station.We classify the most relevant literature and our model by topology (single compressor or network),compressor model(detailed or simplified,where to consider on-off status of compressors as decision variables),uncertainty (deterministic or uncertain) and solution approaches in Table 1.

Table 1 Summary of literature review.

Table 2 Notation and symbols.

In power system research field,which is similar to natural gas network,accommodating the uncertainty in demand is widely studied.Jiang et al.(2014) studied two-stage robust unit commitment problem and solved it by Benders' decomposition based on exact and heuristic approaches.Dashti et al.(2016) studied a twostage robust generation scheduling for a hydrothermal power system with water inflow uncertainty.Ding et al.(2015)addressed the uncertainties of wind power output in active distribution networks using two-stage robust optimization.Methods to handle uncertainties in medical field are studied in Liu et al.(2019),Yang et al.(2019).

3.Model formulation

The two-stage working process under uncertainty in practice is as follows.Firstly the ON/OFF state of each compressor and discharge pressure of compressor station are decided,which are the first-stage decisions.After the uncertainties of volumetric flow rate and suction pressure realize,the rotational speed of each compressor is adjusted to meet the set point of discharge pressure decided in the first stage.Similar process is also introduced in Behrooz (2016).

In deterministic formulation,fluctuation is not considered.The problem is to find the optimal flow rate allocation and compressor rotational speed in a station with given suction pressure and total volumetric flow rate.The objective is to minimize the power cost.We start with the deterministic model and the following assumptions:

(A1) Only one compressor station is considered.It consists ofNnonidentical parallel compressors with the same suction temperature (Ts),suction pressure (ps),and discharge pressure (pd).

(A2) The total natural gas inflow of the station equals its total outflow.The compressors are not powered by the natural gas passing through them.

We use superscriptnto indicate a compressor,where 1 ≤n≤N.Each compressor has specific physical parameters and operationspecific values.These settings determine the mass flow rate range of each compressor.The aim is to minimize the total power cost with an optimal combination of discharge pressure,the working status and rotational speed of each compressor.All the notation and symbols are listed in Table 2.

3.1.Deterministic model

In this part,we introduce the power cost minimization problem(PCMP).The objective is the sum of power cost consumed by working compressorn,which depends on its mass flow ratefn,isentropic headHnand isentropic efficiency ηn:

The isentropic head within the compressor station are the same since compressors are installed in parallel.The isentropic head is determined by pressure and temperature as follows:

whereRis the natural gas constant andm=（σ 1）/ σ,σ is the isentropic exponent,Zsis the compression factor at the suction side of compressor station,both are related to the suction pressurepsin complex forms (Mari′c et al.,2005;Mohring et al.,2004;Mahmoudimehr and Sanaye,2014).

LetQnbe the actual volumetric flow rate andSnbe the compressor rotational speed of compressorn.These two variables need to satisfy:

The isentropic efficiency ηnhas a quadratic representation in terms ofQnandSn:

The relationship betweenQnandfnis

By combining equations(4)-(6),the isentropic efficiency ηncan be expressed as a function of rotational speedSn，discharge pressurepdand suction pressureps,therefore,we can express the consumption power function in equation(1)asPn（Sn，ps，pd）,sinceQnis intermediate variable,we can also usePn（Qn，ps，pd）.The total power consumption of the station is the sum of all the compressor power consumption levels.

The working domain of compressornis given by the following four constraints.Firstly the constraints on rotational speedSnare:

In addition,we have

Any violation of these constraints can lead to impractical solutions.Therefore,it is necessary to ensure the feasibility as much as possible.According to our conversations with natural gas managers,to ensure the feasibility of the operation point is the most important goal for the operators.

Minimization of the power cost under these constraints with deterministic total volumetric flow rateQ=and suction pressureps=is denoted as PCMP,which is also referred to as deterministic model in the following text:whereYnequals to 1 if compressornis on andMis a big positive number.The optimization problem is transformed into a robust optimization problem in section 3.2.The problem is nonlinear and nonconvex.

3.2.Uncertainty sets and robust approach

In this section,we discuss how to optimize the power cost under uncertainty.The natural gas is supplied from upstream storage facilities,production sites or natural gas import terminals and delivered to customers.Natural gas may be extracted out of or injected into the network in an intermediate node.Uncertain customer loads imposed on the system can be represented by the change of flow rate in these nodes.

Fig.1 is a part of typical natural gas pipeline transmission network in both China (Deng et al.,2019) and Belgium (Bakhouya and De Wolf,2007) with at least one intermediate node placed before the selected compressor station.This frame is applicable in any kind of network topology including gun-barrel,tree-shaped and cyclic.In node N2,a customer demand results in volumetric flow rateQextractextracted from the network,which is uncertain as well as volumetric flow rate into the networkQinandQextract.We assume thatQinandQextractvary independently.Then since the suction pressurepscan be calculated from pipeline propertiesW12，W23，QinandQ=QextractQin,we can separate the variation range ofQandps.

We consider two main uncertainty sources:demand fluctuation and suction pressure fluctuation.The pressure loss in the pipeline can be affected by a number of physical properties and gas mixture quality.Many of them can be affected by uncertainty or hard to measure (Aßmann et al.,2019).Assuming the varying range of demand is=[Ql，Qu,the varying range of suction pressure is,similar with Aßmann et al.(2019),we use an uncertainty set of the form.The minimization of the worst-case scenario cost under MIcan be formulated as the following two-stage model:

Fig.1.A typical natural gas pipeline transmission network part with intermediate demand node.

The method to construct MIin practice is as follow.The compressor station operators can predict the volumetric flow rate in the short-term future.Based on the average historical prediction accuracy,a range of possible realization ofQcan be calculated.To obtain the range of the suction pressure,some analysis from historical records of the value are necessary.We can setas the maximal value ofpsin the most recent time window.The calculation ofis similar.In this way,set MIis constructed.

4.Solution approaches for the robust optimization problem

In this section,we introduce solution approaches for problem(11).We first examine the constraint characteristic and find out that when parameters satisfy some conditions,the constraints are all convex or concave.In Section 4.2,we present condition under whichrelatively complete recourseassumption holds to ensure that the first-stage decisions are feasible in all possible second stage uncertainty realizations.In Section 4.3,we give the C&CG algorithm(Zeng and Zhao,2013)and several techniques to solve the problem.Finally,in Section 4.4 we discuss the reduced formulation under convexity assumption of the objective function.

4.1.Constraint characteristic

In this section,we assume that all the first-stage decision variables and uncertainties are given as in the fourth level of (12),which means that we only focus on the second stage optimization.We show that the convexity or concavity of constraints (7)-(9)holds in this stage under certain conditions.

The property can be revealed by equivalent transform of the model,after we transform decision variableSntofn.The objective in model (11) is equal to the following optimization objective:

In the fourth level of problem(12),the first-stage decisions and the suction pressure,total volumetric flow rateQare fixed.Thus to decideSnis equivalent to decideQn.Since

wherek=is given in this innermost problem,then the four constraints mentioned above are only functions of mass flow ratefnby combining equations (2),(3) and (6):

In this section,we examine the formulation of the working domain of each compressor and show that under some parameter conditions the constraints are convex/concave.These conditions are satisfied as shown in several data settings(Deng,2015;Mahmoudimehr and Sanaye,2014).Thefindings in this section enable us to calculate the feasible range of Qn with gradient-based algorithm,which can simplify the calculation to some extent.Otherwise,two quartic equations need to be solved(Deng et al.,2019).For example,when calculating the minimal volumetricflFig.2,if H3≤Hn≤H1,a system of equations needs to be solved,which is equivalent to the following unary quartic equation with Sn being independent variable:

Fig.2.The working domain of a working compressor under uncertainty

Even though there are formulas forfinding roots of this type of equation,the formulas are quite complicated.Numerical methods such as Newton-Raphson method are utilized to quickly obtain solutions when analytical solutions are not necessary,which is common in practical scenarios.Thefindings can replace the Newton’s method or combine with it to overcome shortcomings of Newton’s method such as initial value dependence,solution incompleteness,etc.In the following text,we denote the minimal and maximal volumetricflow rate at given suction and discharge pressure as

4.2.First-stage decision feasibility

The first-stage decision variables should ensure that for any possible realization of the uncertainty,there is always a feasible second-stage solution.In this section,we manage to answer under what condition the relatively complete recourse assumption holds,which is for any given first-stage decision variablesand any realization of the uncertaintythe innermost problem(14) is feasible.If not,a feasibility cut is required.

Fig.3.The isentropic head and effects of suction and discharge pressure.

For a compressor with first-stage decisionYnbeing 1,the working domain of the compressor is shown in Fig.2.Withgiven,the isentropic head is given,thus the volumetric flow rate range can be easily calculated.For any given value ofH,algorithm to calculate minimal and maximal value of volumetric flow rate can be seen at supplementary materials of Deng et al.,(2019).

We propose a method which builds on top of an assumption about the monotonicity of the isentropic head.The calculation ofHin equation (3) takes a complex form.It is because that both compression factorZsand isentropic exponent σ are related to suction and discharge pressure in complicated form.There are articles focusing on calculating the isentropic head (Mari′c et al.,2005).

There is no analytical form condition against which we can derive monotonicity.However,we test numerically and find that the monotonicity holds under a wide range of pressure values.For a given suction temperature of 293.15 K,we test seven settings of discharge pressure ranging from 10.6 MPa to 11.8 MPa,suction pressure ranging from 8.0 MPa to 9.2 MPa.In each curve,suction pressure is discretized at 0.01 MPa.Fig.3 suggests that isentropic head is monotonous in pressure.We introduce the following assumption.

Assumption 1.For a given suction temperature,the isentropic head H is decreasing with suction pressurepsand increasing with discharge pressurepd.

Assumption 2.With fixed suction pressure,the minimal volumetric flow rateis quasi-convex with respect topd,and the maximal volumetric flow rateis quasi-concave with respect topd.

Assumption 2 seems to be intuitive in Fig.2,what's more,it has been numerically implied by massive literature (Wu et al.,2000;Deng,2015;Deng et al.,2019;Behrooz,2016;Mahmoudimehr and Sanaye,2014).Based on this assumption,we discuss the working domain under uncertainty.We have the following proposition about the minimum span ofQn.

where the inequality is due to the definition of quasi-convex.Then we can have the following proposition:

Proposition 3.For a given first-stage decisioncalculate the range of volumetric flow rate of each compressor as inequation(15)and denote the uncertainty set of total volumetric flow rate[Ql，Quof compressor station as.If

4.3.C & CG algorithm and accelerate method

The column and constraint generation (C&CG) algorithm decomposes the two-stage robust optimization problem into a master problem (MP) and a subproblem (SP).Two-stage robust optimization problems are usually hard to solve.Ben-Tal et al.(2004)studied two-stage robust linear programs and showed that the problems in general were NP-hard.

For a given first-stage pair of variablethe SP is solved to identify the worst-case scenario and load conditions:

Algorithm 1.Gradient-ascent Algorithm for SP

At each point,the functiong（Q，ps）can be calculated by mixedinteger linear programming (MILP),dynamic programming (Deng,2015),genetic algorithm (GA) (Mahmoudimehr and Sanaye,2014),etc.In this paper,we choose the method of MILP to solveg（Q）since it is not sensitive to the discretization length and is much faster than the dynamic programming and GA in smaller intervals(Deng,2015).The overall algorithm for solving the SP can be seen in Algorithm 1.

The C&CG algorithm procedure utilizes a master-subproblem framework.The SP returns either a finite optimal solution or+∞with some uncertainty realizations if the first-stage decision variables are infeasible.The MP considers only several selected uncertainty realizations with index l for each realization.The original formulation of MP is:

which is nonlinear.

We utilize piece-wise linear approximation to reformulate the problem into a mix-integer linear problem (MILP).The benefit of this approximation is that the MP can be as accurate as possible,since the compression factorZs，land m are related to the decision variables based on simulation functions which are too complicated to be expressed in explicit form.We discretize the discharge pressurepdin the first stage and defineIpdnew binary decision variablesXi,withXi=1 corresponding to:

where[Ji，l，n=｛1，2，...Ji，l，n｝ for notational consistency.

Through this modeling method,the discharge pressure can also be expressed with constraints:

Thus only when compressor n is working,the actual discharge pressure has relationship withZi，n，j，l.Similarly,the actual volumetric flow rate expressed withZi，n，j，lshould be within the working domain only whenXiandYnequal to 1.

We first list the constraints of the approximation model of MP:

then the whole approximation model of MP is as follow:

(MP-MILP)

The acceleration is based on four thoughts.First,with given first-stage decisions and suction pressure,the feasible volumetric flow range of each compressor on can be calculated.Thus similar with the thoughts proposed in Proposition 3,the possible range of total volumetric flow range is compared with the uncertainty realization.Infeasible first-stage decision can be cut.Second,if we only consider a realization of uncertainty in problem,we can get a lower bound ofIf the lower bound is larger than the current optimal value,there is no need to consider calculatingThird,similarly,relaxing binary variablesZn，j，linto continuous variables will also return a lower bound.In that case,problemturns to a linear programming problem,which is fast to solve.Last,according to Assumption 1 and 2,as discharge pressure increases,after a certain point,the feasible volumetric flow rate range of a compressor starts to shrink,which means that it’s likely that the optimal power cost increases as well.Thus under givens,if after a pointwe find the power cost starts to increas e,all the discharge pressure larger thancan be added into the relaxed formulation to get a lower bound of all left discharge pressure search space.

Algorithm 2.

Accelerated Downward Algorithm for MP

What's more,we have the following proposition:

Proposition 4.TheAlgorithm 2 has the same optimal value as MPMILP.

4.4.Reduced formulation under convexity assumption

Thefirst inequality is based on the definition of g（Q3）,and the second inequality is a direct result ofAssumption 3.Thefinal equality is based on the definition of g（Q1）,g（Q2）.

With the conclusion drawn inProposition 5,the maximal value of a convex function can be found by comparing the function value in minimal and maximal Q.To ensure feasibility,the uncertain realization and should be considered.According toWu et al.(2000),it is typical that the single-compressor power decreases with the suction pressure ps.If this condition holds,to solve the two-stage robust optimization model,we only need to solve one time of the master problem in C&CG Algorithm with four specific scenarios.

In general cases where the monotonicity about psandAssumption 3do not hold,we can not ensure the objective function is convex and reduced formulation in this section is no longer applicable.For instance,if the parameter settings are the same withWu et al.(2000),the C&CG Algorithm frame mentioned inSection 4.3is necessary.In other cases where the convexity of the power cost function is hard to verify,the C&CG Algorithm is also necessary.

5.Numerical results

In this section,we consider six compressors of four different types in a station.The physical parameters are the same as in Deng(2015).We first show the numerical setups and verify the Assumption 3.Then we conduct experiments by the following steps:(1)computation results of the proposed MP-MILP model for the two-stage robust optimization problem under different discretization scheme and uncertainty sets;(2) simulation results comparison between the deterministic power minimization model and the robust model under assumed distribution of uncertainties;(3) cost with varying uncertainty sizes of suction pressure and volumetric flow rate.The codes were developed using C++.We use commercial software Gurobi 9.1.1 for the models.The experiments were performed on an Intel Core i90-9900K 3.60 GHz CPU.

5.1.Input parameters and power cost convexity

The input parameters are listed in Table 3.The compressor station includes one Type A,three Type B,one Type C,and one Type D compressors.The suction temperature is 293.15 K.We examine the convexity of the power costPn（Qn，ps，pd）with respect toQn.The suction pressure varies from 6.75 MPa to 9.14 MPa,and discharge pressure varies from 6.80 MPa to 11.95 MPa,both are discretized at 0.1 MPa.It turns out that for Type A,Type B and Type C compressors,the convexity hold.For Type D compressor,except in certain condition,power costPn（Qn，ps，pd）is nonconvex.Actually,for Type D compressor,in most cases the power cost is nonconvex.We show in Fig.4 the selected plots convex power cost for all the four types,respectively.In Fig.5 we show the selected plots of nonconvex power cost for Type D compressor with corresponding changing of gradient.In both figures the unit of thevolumetric flow rate is under normal condition,that is,the mass flow rate over the density under standard conditions.

Table 3 Physical parameters of the four types of compressors.

The results imply that Assumption 3 does not hold.Thus we need C&CG algorithm to solve the problem.However,even though the power cost is nonconvex,we find thatfor any compressorn,which means that Corollary 1 holds.Thus,with given pressure values,we can use gradient-based algorithm to calculate the range of flow rate.

5.2.Computation time and optimality of MP-MILP and Algorithm 2

In this part,we compare the MP-MILP with Algorithm 2 by two assessment metrics,optimality gaps (Gap) and computation time improvement (Imp).Supposeandare the optimal values of MP-MILP and Algorithm 2.The optimal gap ofwith respect tois defined asThe improvement of computation time of Algorithm 2 with respect to MP-MILP is calculated bywheret0andt1are the computation time of MP-MILP and Algorithm 2,respectively.

We test the performance of MP-MILP and Algorithm 2 under a varying discretizations,uncertainty sets and number of compressors.To better understand the performance of Algorithm 2,we delete the ninth to fifteenth lines in Algorithm 2,only whenwe need to continue from sixteenth line.We denote it as algorithm 3.The performance of algorithm 3 is recorded so as to test the effects of relaxation related toon the Algorithm 2.

The suction pressure ranges from 8.2 MPa to 9.0 MPa and the volumetric flow rate ranges from 32.0m3/s to 36.0m3/s.The search scope of discharge pressure is 9.0 MPa-12.0 MPa.We set the discretization interval of suction pressure as 0.05 MPa,the discretization interval of volumetric flow rate as 0.01m3/s unless otherwise specified.Default uncertainty set has six scenarios with suction pressure being ｛8.2，9.0｝ MPa and volumetric flow rate being ｛32.0，34.0，36.0｝m3/s unless otherwise specified.

The results are shown in Tables 4-6.The optimal value of MPMILP,Algorithm 2 and algorithm 3 are the same in all numerical instances.Compared with MP-MILP,algorithm 3 already has significant computational time improvements with a minimum speed improvement of 224.7%.Algorithm 2 further promote the speed to a certain extent with time improvements of 1518.9% on average(3785.9%at maximum).The results also imply that among the four accelerate techniques,the first three techniques excluding cutting related toare crucial to outperform MP-MILP.We utilize Algorithm 2 in the subsequent calculation.

5.3.Computation results of C&CG algorithm

In this part,we solve the whole two-stage robust model with C&CG algorithm and show the results.The gradient-ascent algorithm for SP terminates when the gap between upper bound and lower bound is smaller than 0.001 kW.The C&CG algorithm terminates when the gap between upper bound and lower bound is smaller than 0.002 kW.

What's more,we also consider a scenario-based two-stage stochastic programming (STSP) model.The settings of uncertainty sets,deterministic model and parameters for STSP are listed in Table 7 and Table 8.For deterministic model,the uncertain parameter is replaced by average value of the varying range.For STSP,the selected uncertainty scenarios and probabilities are listed in Table 8.Similart with Gorissen et al.(2015),for the simulation we draw suction pressure and volumetric flow rate values uniformly from the box uncertainty.For each setting,the simulation is run for 500 times.The same input is used for deterministic,robust and stochastic models.

We record the average cost and worst-case cost during the simulation in Table 9.The infeasible ratio is the infeasible times divided by total simulation times.The results show that in simulation,the decision variable of deterministic model has an average 9.9% probability of being infeasible.While the robust model and stochastic model ensure feasibility in all the simulations.The robust model has lower worst-case cost compared with stochastic model and the average cost is also higher.In the first setting,with an increase less than 3% in average cost,a reduction of 5.8%-14.7% in worst-case cost is obtained.The results also imply that stochastic programming is highly dependent on the preset scenario probability,while robust model need not assuming any information about uncertainty information.

It is worth mentioning that if the customer demand suddenly changes and the decision variables are infeasible for the new demand,just as shown in the deterministic case in Table 9,we can only adjust the discharge pressure or the working status of the compressor.However,due to the long pipelines with large diameters,the changes will take time to react(Behrooz,2016).Thus during the adjusting time,failure to satisfy all the constraints may happen.These sudden changes of the operational conditions may result in bad economic outcomes.

Table 4 Computation time and optimality gaps with varying discretization intervals.

Table 5 Computation time and optimality gaps with varying uncertainty set sizes.

Table 6 Computation time and optimality gaps with varying compressor number.

Table 7 Values of the uncertain and deterministic sets.

Table 8 Settings of the STSP.

Table 9 Computation results of the robust models.

Fig.4.Selected Plots of Pn（Qn，ps，pd）which is convex in Qn

5.4.Simulation results under different uncertainty sizes

Next,we examine the effect of uncertainty set size on the results.We first fix the varying range of volumetric flow rate at[21.4，26.6 m3/s.The varying range of suction pressure is 1%-10% of 8.6 MPa in Fig.6.To examine the effect of volumetric uncertainty set size,we fix the varying range of suction pressure at [7.8，8.8 MPa.The varying range of volumetric flow rate is 1%-10% of 24m3/s in Fig.7.As the uncertainty size increases,the objective value increases.Discharge pressure also has a similar pattern,as the uncertainty size increases,we need more pressure buffer to ensure the feasibility in all possible realizations of the uncertainty sets.Accordingly,the robust optimization model provides results that in all simulated cases,loads imposed can be met while all the feasibility constraints are respected.

Next,we also compare the robust model results with deterministic model and stochastic model results under varying uncertainty sizes of suction pressure and volumetric flow rate in Fig.8 and Fig.9.We first present the ratio of times where the deterministic model result is infeasible in all simulation runs.And then,we compare the results of robust model with stochastic model from the perspective of average cost and worst-case cost.

The results show that robust model and stochastic model both are feasible in simulation.However,deterministic model results can be infeasible.What's more,suction pressure variability has a greater impact than volumetric flow rate variability.The infeasible ratio of deterministic model with varying suction pressure is between 15.8% and 45%while the same value for varying volumetric flow rate is between 6.6% and 9.2%.When the range of volumetric flow rate is small,robust model and stochastic model has same results.Generally,the average cost of robust model is slightly higher than that of stochastic model,but the worst-case cost is significantly reduced by up to 25%.

In summary,the robust model provides appropriate discharge pressure and working-status settings for a compressor station in uncertain cases.A reduction of 6.6%-45% in infeasible ratio is gained.What's more,robust model has lower worst-case cost with slightly increased average cost compared with stochastic model.

6.Conclusions

In this study,we develop an innovative two-stage robust modeling framework in treatment of the inherent natural gas demand uncertainty,considering technical details such as the stonewall line and surge line constraints.To the best of our knowledge,little work has been conducted on the topic of natural gas robust optimization with detailed model.

Fig.5.Selected Plots of Pn（Qn，ps，pd）of Type D compressor which is nonconvex in Qn

We prove that when input parameters satisfy some condition,the constraints are convex/concave,making it easier to calculate the feasible range.We examine the conditions under which the first-stage decision feasibility holds and form feasibility cuts.To solve the problem,we utilize the C&CG algorithm.The SP is solved by gradient-ascent algorithm combined with MILP.As for the master problem,accelerate algorithm for original linear approximation formulations gains considerable computation time improvement.

We carry out extensive experiments and verify the benefits and understand the limitations of this approach.It turns out that the robust model and stochastic model can significantly increase the feasibility of first-stage decision in our settings.In robust model,the worst-case cost can be reduced with the price of extra power cost compared with stochastic model,which is the cost we need to pay to gain robustness.The operational cost increases averagely when the uncertainty size increases.The uncertainty size can be used to“tune”the conservativeness when the compressor station is operated.A larger uncertainty size provides a more robust operation,but the price will increase as well.These management insights can help guide the trade-off between robustness and operational cost in practice.In order to be less conservativeness,improving the forecast accuracy will be effective.

Future research directions of this work include expanding the research object from a single compressor station to a pipeline network,considering distributed robust optimization methods to make better use of historical data information.What's more,the

current work considers only steady-state model,the transient-state model can also be studied in the future.

Fig.6.Robust model results with varying uncertainty sizes of suction pressure.

Fig.7.Robust model results with varying uncertainty sizes of volumetric flow rate.

Fig.8.Simulation results with varying uncertainty sizes of suction pressure.

Fig.9.Simulation results with varying uncertainty sizes of volumetric flow rate.

Acknowledgments

We thank Dr.Jianjun Liu from China National Petroleum Corporation.He gives advice on the selection of variables at each stage of the robust model,which helps improve this study.The authors also thank the editor and three anonymous referees for their constructive suggestions,which significantly improve this paper.Tianhu Deng acknowledges the support from the National Science Foundation of China (Grant 71822105).