Optimal design and operation for supply chain system of multi-state natural gas under uncertainties of demand and purchase price

Computers & Industrial Engineering 131 (2019) 115–130

Contents lists available at ScienceDirect

Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie

Optimal design and operation for supply chain system of multi-state natural gas under uncertainties of demand and purchase price

T

Haoran Zhanga, , Yongtu Liangb, Qi Liaob,c, Jinyu Chena, Wan Zhangb,d, Yin Longe, Chen Qianf ⁎

a

Center for Spatial Information Science, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8568, Japan Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum-Beijing, Fuxue Road No. 18, Changping District, Beijing 102249, PR China c Centro de Matemática Aplicações Fundamentais e Investigação Operacional, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal d Department of Human and Engineered Environmental Studies, Graduate School of Frontier Sciences, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8563, Japan e Centro de Estudos de Gestão do Instituto Superior Técnico (CEG-IST), Universidade de Lisboa, Avenida Rovisco Pais, 1049-001 Lisboa, Portugal f University of Calgary, 2500 University Drive, NW, Calgary, Alberta T2N 1N4, Canada b

ARTICLE INFO

ABSTRACT

Keywords: Supply chain system Multi-state natural gas Optimal design and operation Uncertainty Monte Carlo sampling

The transportation and storage mode of natural gas is various, and natural gas can be transformed into multiple states. Furthermore, it is realized that there are uncertainties in the natural gas purchase price and cities’ demands. Thus, the optimal design and operation for the supply chain of multi-state natural gas are revealed as a complex issue. Although many scholars have researched this issue, they usually focused on the design and operation optimization for one-state natural gas. To consider the uncertainties in natural gas purchase price and demand, this paper develops a multi-simulation MILP model based on Monte Carlo sampling, which is used for simulation runs to generate a deal of uncertain parameters. The model captures four seasons, three states of natural gas and four transport options. Finally, four cases are presented as an example to illustrate the solution that can satisfy the real requirements, testifying the application of the model. Simultaneously, it is confirmed that the price uncertainty and the demand uncertainty respectively play an important role in the entire operation cost and the construction scheme of the system.

1. Introduction 1.1. Background Recently, the application of natural gas (NG) has covered a broad field, ranging from electricity generation to households as an effective and clean energy. With the completion of national NG reserve strategy, NG industry and infrastructures get constantly developed. It contributes to a significant ongoing improvement in the NG supply situation (Touretzky, Mcguffin, Ziesmer, & Baldea, 2016). Currently, international NG supply chain is developing rapidly and encountering an opportunity as well as challenge. Specifically, on one hand, the NG transportation cost accounts for one-third of the total NG industry cost. A well-designed supply chain system can help energy companies to boost profits (Wang, Liang, Zheng, Yuan, & Zhang, 2018; Wang, Yuan, Zhang, Zhao, & Liang, 2018; Zhang et al., 2017). On the other hand, with the increase of the NG demand, the supply chain becomes more and more complex to manage. Any problem in the links, including exploration, storage, transportation, ⁎

and sales, will affect the operation and development of natural gas supply systems in a country or even the entire natural gas industry. To ensure the reliability of every stage in a NG supply chain, each link consisting of NG purchase, transportation, and usage in the chain must be planned (Alves, Souza, & Costa, 2016; Baccanelli, Langé, Rocco, Pellegrini, & Colombo, 2016). Hence, management of the natural gas supply chain will remain one of the main concerns in the near future. However, the planning, construction and operation management of an NG supply chain is a huge systematic project. The optimization of the natural gas supply chain must meet various transportation modes, flexible traveling routes, fluctuated demands, structure layouts and transportation schemes. Thus, it is very difficult to determine the optimal construction layout and work out the minimal cost supply scheme to ensure economic profit and stability of natural gas supply within the satisfaction of market supply situation (Kashani & Molaei, 2014; Lu, Su, Fath, Zhang, & Yan, 2016; Sanaye & Mahmoudimehr, 2013). Although many existing researches have focused on this challenge, there are still some limitations need to be improved.

Corresponding author. E-mail address: [email protected] (H. Zhang).

https://doi.org/10.1016/j.cie.2019.03.041 Received 4 April 2018; Received in revised form 6 February 2019; Accepted 20 March 2019 Available online 22 March 2019 0360-8352/ © 2019 Elsevier Ltd. All rights reserved.

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Nomenclature

VSLmaxk, i

Sets and indices

VSCmaxk, i

s S = {1, …NSCE} set of simulation numbers, NSCE is the maximum number of simulations j J = {1, 2, 3} set of transport options, where 1 denotes PNG transportation, 2 denotes LNG transportation, 3 denotes CNG transportation i, i I set of nodes k, k K set of seasons e Ej set of carriers under transport option j , where E1 contains pipeline, E2 contains LNG ships and road tankers while E3 only contains CNG road tankers

VQLmaxi VQGmaxi

VQCmaxi VQAmaxi

VSTLmaxi VSTLini

Continuous parameters

RB RP RCj, e

RS RT RV FSP FPC FPE FPM FCj, e FCMj, e FSGs, i, k FSLs, i, k FSCs, i, k

FST FSTTi FSP

FTL FLL FTG FLG FST FCGC FTG FCA FLSL FCC

GT

LT CT

G

A L C

G

A

VSGmaxk, i

QLPmini, i

depreciation rate of infrastructures depreciation rate of the pipeline depreciation rate of the carrier e under the transport option j depreciation rate of peaking underground storage depreciation rate of NG processing facilities depreciation rate of LNG, CNG storage facilities NG pipeline and station construction cost (CNY) unit price of NG pipeline construction cost (CNY) unit price of NG PNG cost (CNY) unite price of NG pipeline maintenance cost (CNY) unit price of the carrier e under the transport option j (CNY) unit transportation price of the carrier e under the transport option j (CNY) unit NG purchased price at the node i in the season k of the simulation s (CNY/m3) unit LNG purchased price at the node i in the season k of the simulation s (CNY/m3) unit CNG purchased price at the node i in the season k of the simulation s (CNY/m3) unit cost of peaking underground storage (CNY/m3) construction cost of peaking underground storage at node i (CNY) construction cost of liquefaction, gasification, compression and regulation stations (CNY) unit construction price of liquefaction facilities (CNY) unit price of the liquefaction process (CNY/m3) unit construction price of gasification facilities (CNY/m3) unit price of the gasification process (CNY/m3) unit construction price of compression facilities (CNY/m3) unit price of compression process (CNY/m3) unit construction price of regulation facilities (CNY) unit price of the regulation process (CNY/m3) LNG storage facility cost (CNY/m3) CNG storage facility cost (CNY/m3) NG pipeline design margin. periodical rate of LNG, CNG storage PNG transmission loss rate LNG transmission loss rate CNG transmission loss rate gasification process loss rate regulation process loss rate liquefaction process loss rate compression process loss rate volume ratio of LNG to NG volume ratio of CNG to NG maximum volume of purchasable NG at the node i during

QLPmaxi, i

QREj, i, i , e

TAj, e

RAj, e VTj, e NTj, e Li , i , j , e

the season k (m3) maximum volume of purchasable LNG at the node i during the season k (m3) maximum volume of purchasable CNG at the node i during the season k (m3) upper boundary of the liquefaction ability at the node i (m3) upper boundary of the gasification ability at the node i (m3) upper boundary of the compression ability at the node i (m3) upper boundary of the pressure regulation ability at the node i upper boundary of the storage ability at the node i (m3) initial storage volume of the peak shaving underground storage at the node i (m3) lower boundary of NG PNG capacity from the node i to the node i (m3) upper boundary of NG PNG capacity from the node i to the node i (m3) upper seasonal transportation difference boundary under the mode j from the node i to the node i by using the carrier e (m3) transportation, load and unload time under the mode j by using the carrier e (m3) average speed of the carrier e under the mode j (km/h) filling volume of the carrier e under the mode j (m3) average effective running time of the carrier e under the mode j required running distance of the carrier e under the mode j from the node i to the node i (km)

Discrete parameters

n maxj, i, e Dk

maximum number of the carrier e under the mode j when the node i is the origin the number of days in the season k (d)

Binary parameters

HBi

HSTi

HLi

HGi HCi HPi, i

HAi

ZNj, i, i , e

ZBi, i

ZSTi

116

binary parameters of existing pipeline stations, HBi = 1 if the node i is an NG pipeline station, otherwise, HBi = 0 binary parameters of existing peak shaving underground storage s, HSTi = 1 if the node i is a peak shaving underground storage, otherwise, HSTi = 0 binary parameters of existing liquefaction stations, HLi = 1 if the node i is a liquefaction station, otherwise, HLi = 0 binary parameters of existing gasification stations, HGi = 1 if the node i is a gasification station, otherwise, HGi = 0 binary parameters of existing compressor stations, HCi = 1 if the node i is a compressor station, otherwise, HCi = 0 binary parameters of the existing pipeline from the node i to the node i , HPi, i = 1; otherwise, HPi, i = 0 binary parameters of existing pressure regulation stations, HAi = 1 if the node i is a pressure regulation station, otherwise, HBi = 0 binary parameters of transportation limits, ZNj, i, i , e = 1 if it is possible to use the carrier e for transporting under the mode j from the node i to the node i , otherwise, ZNj, i, i , e = 0 binary parameters of adjacent nodes, ZBi, i = 1 if the node i is adjacent to the node i , otherwise, ZBi, i = 0 binary parameters of peak shaving underground storage limits, ZSTi = 1 if the node i is suitable for building a peak shaving underground storage here, otherwise, ZSTi = 0

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Continuous variables

VSTs, k, i

VCs, k, i

daily peak shaving volume at the node i of the simulation s . Positive when inputting gas and negative when outputting gas (m3)

VLmaxi VCmaxi

Positive continuous variables

Qs, k, j, i, i , e

Qmaxj, i, i , e

VSGs, i, k VSLs, i, k VSCs, i, k VSTs, k, i

VQLs, k, i VQLmaxi

VQGs, k, i VQGmaxi

VQCs, k, i VQCmaxi

VQAs, k, i VQAmaxi VGs, k, i VLs, k, i

VDGs, k, i

average daily transportation volume from the node i to the node i by using the carrier e under the mode j in the season k of the simulation s (m3) maximum average daily transportation volume from the node i to the node i by using the carrier e under the mode j during four seasons (m3) purchased NG volume at the node i in the season k of the simulation s (m3) purchased LNG volume at the node i in the season k of the simulation s (m3) purchased CNG volume at the node i in the season k of the simulation s (m3) peak shaving NG volume storage at the node i in the season k of the simulation s (m3) liquefaction NG volume at the node i in the season k of the simulation s (m3) maximum liquefaction NG volume at the node i during four seasons (m3) gasification NG volume at the node i in the season k of the simulation s (m3) maximum gasification NG volume at the node i during four seasons (m3) compressed NG volume at the node i in the season k of the simulation s (m3) maximum compressed NG volume at the node i during four seasons (m3) regulated NG volume at the node i in the season k of the simulation s (m3) maximum regulated NG volume at the node i during four seasons (m3) available NG volume at the node i in the season k of the simulation s (m3) available LNG volume at the node i in the season k of the simulation s (m3)

VDLs, k, i VDCs, k, i

available CNG volume at the node i in the season k of the simulation s (m3) maximum available LNG volume at the node i during four seasons (m3) maximum available CNG volume at the node i during four seasons (m3) required NG volume at the node i in the season k of the simulation s (m3) required LNG volume at the node i in the season k of the simulation s (m3) required CNG volume at the node i in the season k of the simulation s (m3)

Positive discrete variables

n s , k , j, i , e

n maxj, i, e

number of the carrier e under the option j in the season k of the simulation s when the node i is the origin maximum number of the carrier e under the option j during four seasons if the node i is the origin

Binary variables

CBi

SPi, i

CSTi

CLi

CGi

CCi

CAi

1.2. Related work

pipeline station construction binary variables. CBi = 1 if the node i is required to build an NG pipeline station, otherwise, CBi = 0 pipeline construction binary variables. SPi, i = 1 if the node i is required to build NG pipeline, otherwise, SPi, i = 0 peak shaving underground storage construction binary variables. CSTi = 1 if the node i is required to build a peak shaving underground storage, otherwise, CSTi = 0 liquefaction station construction binary variables. CLi = 1 if the node i is required to build a liquefaction station, otherwise, CLi = 0 gasification station construction binary variables. CGi = 1 if the node i is required to build a gasification station, otherwise, CGi = 0 compressor station construction binary variables. CCi = 1 if the node i is required to build a compressor station, otherwise, CCi = 0 pressure regulation station construction binary variables. CAi = 1 if the node i is required to build a pressure regulation station, otherwise, CAi = 0

operation, maintenance, the environmental in operational stage (Shi, Zhou, Xu, Zhou, & Shi, 2018). Hence, life-cycle cost analysis is presented for a more comprehensive and accurate calculation of cost (Ally & Pryor, 2016; Li, Ma, Liu, & Zhang, 2018; Zhang, Li, Yang, Li, & Qian, 2014). Kashani and Molaei (2014) developed techno-economic and life cycle analysis models to conduct a comparative assessment of the delivered costs and lifecycle greenhouse gas (GHG) emissions of the NG supply chain from production sites in Canada to the north and southwest Europe. With the increasing demand for LNG in recent years, the economic evaluation of the LNG supply chain gradually becomes the focus in research. For example, Kim, Seo, Chang, and Yan (2016) put forward an economic evaluation method to analyze the life-cycle cost (LCC) and the life cycle profit (LCP) of LNG supply chain and concluded that the maximum profit would be obtained when the pressure of the LNG product was approximately 7 bar. Meanwhile, they conducted sensitivity analysis for the various LNG price, LN2 cost, and flow rate of LNG product. Recently, some scholars have researched on overall optimization of NG or another energy supply chain (Chen, He, Li, & Zhang, 2018; Fernandes, Relvas, & Barbosa-Póvoa, 2013). Hamedi, Farahani, Husseini, and Esmaeilian (2009) took the minimum operation cost of

Generally, the optimization object of NG supply chain includes construction layout and detailed supply scheme. A large number of factors that involve various transport options, flexible transportation routes, fluctuated demands, structure layouts need to be taken into account to determine the optimal construction layout and work out the most cost-saving supply scheme (Kashani & Molaei, 2014; Lu et al., 2016; Sanaye & Mahmoudimehr, 2013). As the preliminary but indispensable work of NG supply chain optimization, the economic analysis of different NG transport options have been widely studied. Thomas and Dawe (2003) presented many of essential technical points and broad economic pointers on transport options. They include pipelines, LNG, CNG, gas to solids (GtS), gas to the wire (GtW) and gas to liquids (GtL). The preferable transport options must not only be determined by economic risks, but also the negative effects of possible terrorist activity, political changes as well as trade embargos over long periods of time. Later on, the object of economic assessment is converted from a single transport option to the supply chain during the whole life cycle. Any transport option requires not only a huge infrastructure investment in the early stage but also the 117

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

the overall NG supply chain system as the objective function and established MILP method that considers constraints such as volume, supply-demand balance, flowrate balance, etc. Based on the work of Hamedi et al. (2009), Azadeh, Raoofi, and Zarrin (2015) also paid attention to the impact of greenhouse gas on the environment during the production and transportation of NG. They also set up the fuzzy index for double-objective optimization on the economy and environmental protection. The two models referred above are established based on the PNG modes. Samsatli and Samsatli (2015) developed a common MILP method for energy system network consisting of transportation, storage, and infrastructures. The model could not only determine the optimal network structures but also consider short-term and long-term scheduling. Elia, Baliban, and Floudas (2013) put forward MILP method for GTL transport option to minimize daily operation cost by solving out the processing capacity, oil types and producing a proportion of all NG synthetic oil plants involved in the supply chain. However, their discussions were based on one single transport option and did not consider the direct transformation between multi-states of NG and the uncertainties of demands and prices. Resat and Turkay (2015) investigated the network operation with multimodal transportation of different transport options in Marmara, Turkey and proposed a MILP method to obtain the optimal operation scheme. They also took the least operation cost as the objective function. However, the NG supply chain is usually undertaken with a lot of uncertainties, especially in demand, price and profit. As mentioned in the background, the demand for NG during peak period is of great difference from that during valley period. If predictive statistic of the demand is not detailed enough, it will result in a great effect on operation and production (Ficco, Dell’Isola, Vigo, & Celenza, 2015; Kou, Liang, & Gao, 2016). Whereas, the demand for NG is usually volatile. If the uncertainty is ignored, there will be an imbalance between supply and demand, leading to an increase in investment and operation cost. As an important facility for adjusting the disequilibrium between supply and demand, the underground gas storage reservoir will be constructed near the downstream market to solve the gas storage occurred in winter. Moreover, the price of NG and LNG strongly influences the overall construction of systems (Choi, Lee, & Lee, 2016; Zhang, Li, Wang, & Li, 2016). Hence, recent works attempt to explicitly model the demand uncertainties in NG or other energy supply chains (Choi, Realff, & Lee, 2005; Nosratabadi, Hooshmand, & Gholipour, 2016; Verdejo, Awerkin, Saavedra, Kliemann, & Vargas, 2016; Zamarripa, Aguirre, Méndez, & Espuña, 2013). The typical approach for processing uncertain parameters of a model is Monte Carlo sampling, which can generate a series of parameters during each simulation according to the probable distribution of uncertain parameters (Lima, Relvas, & Barbosa-Póvoa, 2018; Marmidis, Lazarou, & Pyrgloti, 2008; Marseguerra, Zio, & Podofillini, 2002; Oliveira, Gupta, Hamacher, & Grossmann, 2013). Calvillo, Sanchez-Miralles, Villar, and Martin (2016) developed an optimization model with considering stochastic prices for the optimal planning and operation of aggregated DER systems. Halvorsen-Weare, Fagerholt, and Nnqvist (2013) addressed a real-life ship routing and scheduling problem from the LNG business considering uncertainties in ship travel time and production level. A simulation-optimization framework was proposed for evaluating solutions. But this paper only optimized the operation plan of NG supply chain but ignored the supply chain design. Li, et al. presented a series of stochastic pooling problem optimization methods to address product quality and uncertainty issues of NG production network design and operation (Li & Barton, 2015; Li, Armagan, Tomasgard, & Barton, 2011; Li, Tomasgard, & Barton, 2017). They took into account two gas states, viz., LNG and NG, but missed the impact of different transport options, the construction of gas storages and the seasonal changes on the supply chain. In summary, at present, few studies considered the factors that involve the direct transformation of NG multi-states, a variety of transport

options, construction of gas storages, seasonal changes, uncertainties of demands and prices to work out the optimal design as well as operation of complicated NG supply chain simultaneously. Aiming at these points, a multi-simulation MILP method coupling with Monte Carlo sampling is proposed in this paper. Inter-states transformation of NG, LNG, and CNG are also considered in the model. By solving the model, we can get the construction plan of NG supply chain infrastructure development and the detailed seasonal storage and transportation plan of NG supply chain. 1.3. Contributions of this work

• An uncertainty optimization approach to the design and operation of multi-state NG supply chain is presented. • The model considers three transport options (PNG, LNG, CNG) and four transport carriers. • Through this method, we can obtain detailed seasonal storage and transportation plan of NG supply chain. • This paper analyzes the impact of uncertainty parameters on the model. • A real example of China testifies the practicality of the method. 2. Methodology 2.1. Problem description The upstream of an NG supply chain includes NG production or import; the middle part includes transportation, storage and processing; and the downstream (consuming cities) includes the sales of NG products. Currently, NG products include a variety of forms, among which NG (Bilgin, 2009; Rodger, 2014; Wang, Yuan, et al., 2018), LNG (Lin, Zhang, & Gu, 2010; Tsatsaronis & Morosuk, 2010) and CNG (Boostani, Ghodsi, & Miab, 2010) are world-widely common forms for transportation or usage. CNG is made by compressing NG and is stored in hard containers at a pressure of 20–25 MPa, whereas LNG is made by cooling down NG to liquid form for ease and safety of non-pressurized storage or transport. The schematic diagram for the supply chain is depicted in Fig. 1. Upstream is characterized by various sources and big price difference. For instance, NG can be obtained from inland gas fields, import, LNG gasification or CNG regulation; LNG can be obtained by NG liquefaction or import, and CNG can be obtained by NG purification and compression. The downstream fluctuation of an NG supply chain reflects the uneven usage of time and area. There is a huge seasonal gas consumption difference, especially in winter and summer (Spoladore, Borelli, Devia, Mora, & Schenone, 2016). Thus, it is necessary to build underground gas storages near the market to ensure the NG supply in winter. Additionally, the local resources and economy differences in different areas causes a huge difference in city gas consumption (Li, Dong, Shangguan, & Hook, 2011; Soldo, 2012). Transportation is the main part of the midstream of an NG supply chain, which has a strong impact on the economy. According to the state of NG products, we firstly divide transport options into three types, namely NG, LNG, and CNG. Then based on investigation and economic analysis, the common carriers for NG, LNG, and CNG are pipelines, LNG ships/road tankers, and CNG road tankers, respectively. Considering downstream demand and transport economy, the mutual transformation among these three forms would take place before being transported or after reaching the destination. Hence, the processing technologies and devices for gasification, liquefaction, compression, and decompression need to be considered in midstream. Since NG can only be transported by pipeline in this paper, this transport option is defined as Pipeline Natural Gas (PNG). Different from LNG and CNG, PNG concerns about the design capacity of the carrier (pipeline) rather than the number of carriers and it has been 118

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Fig. 1. Supply chain system of multi-state NG.

2.2. Model requirements

mature enough after decades of development. More specifically, it requires less space, has a large transport capacity, suffers less friction, is free for any extreme weather conditions and can operate continuously and stable for a long time (Guandalini, Colbertaldo, & Campanari, 2017; Ríos-Mercado & Borraz-Sánchez, 2015). Meanwhile, PNG is limited to gas sources and investments, requiring a huge amount of initial outlay and infrastructures. The cost depends on the transportation capacity and distance. Large capacity leads to lower cost, and thereby PNG is suitable for large capacity, stable NG sources and inland transportation with long-term supply (Liu, Li, & Yang, 2014). LNG supply chain represents the process of NG supply in the form of LNG, including upstream, midstream and downstream stages. The upstream includes NG purification, separation, liquefaction and import; the midstream includes LNG road tankers or ships transportation; and the downstream includes terminal markets, LNG satellite stations, gas stations, and LNG filling stations (Xia, Cardin, Ranjbar-Bourani, & Caunhye, 2015). Since LNG storage is characterized by high density, low pressure, and strong safety, LNG transportation is widely applied to sea shipping and inland transportation (Burel, Taccani, & Zuliani, 2013). Recently, many countries have built LNG receiving stations along sea shores to alleviate the increasing demand for NG. LNG sea shipping is the only long-distance transoceanic mode, which transports LNG from liquefaction plants to LNG receiving stations. LNG inland transportation sends LNG from receiving stations to satellite stations using road tankers and LNG ships. The road tanker transportation is suitable for small capacity and sparse distribution situation, while the inland river shipping has a promising market future, which is suitable for large capacity, stable and flexible situation (Kuwahara, Bajay, & Castro, 2000). CNG transportation is a complementary mode as well (Khan et al., 2016). The lower pressure NG got compressed and transported to CNG gas stations or regulation stations of consuming places by road tankers with high compressed cylinders (Frick, Axhausen, Carle, & Wokaun, 2007). In the light of complex process and much cost and plant investment of LNG liquefaction, CNG transportation is easy to manipulate and less cost, which is applicable for small gas capacity and long distance situation (Demirbas, 2002).

Aiming at this issue, we proposed a MILP approach to find the optimal global results and work out the economic NG supply scheme. The objective function is to minimize the total construction and operation cost of LNG supply chain, considering purchase prices of NG and LNG and demand uncertainties of all-consuming cities. Since the system is complex, to establish and solve the model effectively, assumptions are made as follows:

• The gas field and the city governed area are set as one node (Kim et al., 2016). • Only the common transport options (PNG, LNG, and CNG) are involved, where LNG transportation includes LNG ships or road tankers and CNG transportation include CNG road tankers only.

The mean values and variances of the seasonal demands of each city and the NG purchase prices are the key inputs of the model which indicate the uncertainties of the NG market; Additionally, parameters of transportation, for instance, filling ability, load and unload speed and running speed are also needed to be given before the computation; The geographical information, such as locations of local gas fields, cities, and alongshore LNG receiving stations, is the key parameter of the model. Finally, the capacity limitations of transportation, processing, and storage should be given. After the computation, we could obtain the follows from the proposed model: the construction schemes of pipeline, stations (pipeline stations, liquefaction stations, gasification stations, compression stations and regulation stations) and underground gas storages; the seasonal local transportation scheme (including transport option, capacity, and carriers) of each simulation; the seasonal purchase amount of NG of each simulation; and the seasonal storage and peak shaving volume of each simulation. 3. Mathematical model In this section, a multi-simulation MILP model is developed based on Monte Carlo sampling. Specifically, Monte Carlo sampling is adapted 119

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

to generate a set of simulations according to the distribution features of uncertain demand and price. And then, these simulations are being modeled as stochastic parameters in the MILP model to find the optimal solution. Note that the best-found solution is not optimal in all the simulations, but is satisfactory in most of the simulations. That is to say, the obtained solution enables the NG supply chain running at relatively low cost even when demand or price fluctuations occur. Overall, the proposed mathematical model involves five sets: the set of simulation numbers (S ); the set of nodes (I ); the set of seasons (K ); the set of transport options ( J = {1, 2, 3} ), where 1 denotes PNG, 2 denotes LNG, 3 denotes CNG; the set of carriers under transport option j (Ej ), where E1 only contains pipeline, E2 contains LNG ships and road tankers while E3 only contains CNG road tankers. To better understand the meaning of each symbol, a Nomenclature section is added towards the back of the paper, where readers can refer from it.

2 2s

f1s = 1 1

=

+

1 3s

f3s = 3 1

+

RB (1

+

s

HBi ) CBi FSP

=

RP (1 i I i

1 3s

3 2s

1 4

HPi, i ) Qmax1, i, i , e Li, i ,1, e FPC

Dk Qs, k,1, i, i , e Li, i FPE I k K e E1

=

SPi, i Li, i FPM i I i

I

2 1

2 1

+

2 2s

=

s

S

RC2, e n max2, i, e FC2, e e E2 i I

(3a)

S

RC3, e n max3, i, e FC3, e

=

(3b)

Dk ns, k,3, i, e FCM3, e

Dk

(3c)

k K

Dk (VSGi, k, s FSGi, k, s + VSLi, k, s FSLi, k, s k K

i I

+ VSCi, k, s FSCi, k, s )

Dk

s

S

(4)

k K

NG storage cost for peak shaving includes the depreciation cost of underground storage and the operational cost for storing these redundant gas. The depreciation cost ( 15 ) is greater than zeros only when new underground storage is built. While the operational cost( 25s ) is linear with average daily stored volume.

f5s = 5 1

5 1

=

+

5 2s

s

RS (1

(5a)

S

HSTi ) CSTi FSSTi

(5b)

i I

5 2s

=

Dk VSTs, k, i FST

Dk

i I k K

(5c)

k K

Liquefaction and gasification cost includes the depreciation cost of liquefaction and gasification station construction ( 16 ), the depreciation cost of equipment ( 26 ), the depreciation cost of storage ( 36 ) and processing cost ( 46s ). Station construction cost does not correlate with flowrate or distance. The depreciation cost of equipment and storage depends on the maximum daily processing volume, while the processing cost depends on the average daily processing volume.

f6s = 6 1

6 1

=

+

6 2

+

[RB (1

6 3

+

6 4s

s

(6a)

S

HLi ) CLi FSP + RB (1

HGi ) CGi FSP]

i I 6 2

6 3

=

[RT VQLmaxi FTL + RT VQGmaxi FTG]

Dk k K

=

(6b) (6c)

RV VLmaxi FLSL

(6d)

i I

(1c)

=

Dk (VQLs, k, i FLL + VQGs, k, i FLG ) i I k K

(1d)

Dk k K

(6e)

Compression and regulation cost are similar to liquefaction and gasification cost.

(1e)

f6s =

LNG transportation cost includes depreciation cost of LNG carriers ( 12 ), transportation cost and maintenance cost ( 22s ): The former one depends on the maximum daily number of carriers, whereas the latter one depends on the average daily number of carriers.

f2s =

s

NG purchase cost is correlated with average daily purchase volume and price:

6 4s i I i

3 2s

e E3 i I k K

(1b)

I e E1

=

+

i I

i I 1 2

(2c)

e E2 i I

(1a)

S

3 1

=

f4s =

1 4

Dk k K

CNG transportation cost is similar to LNG transportation cost.

During each simulation, the uncertain parameters are determined by a group of certain ones. In this paper, the infrastructure and equipment investment are dynamically depreciated in the way of equivalent payment by installments. The objective function (minF = s (f1s + f2s + f3s + f4s + f5s + f6s + f7s )/ NSCE ) represents the mean value of all the costs per day of different simulations. The mean total cost per day mainly includes transportation cost for PNG as f1s , LNG as, CNG as f3s , purchase cost for NG as, peak adjustment for storage of NG as f5s , liquefaction and gasification cost as f6s and pressure adjustment cost as f7s . The transportation cost for PNG is given by Eq. (1), which is composed of four parts. The first term ( 11) stands for the total depreciation cost of initial station construction. Construction cost for the initial station is the basic civil construction cost, which is not related to flowrate or transportation distance. Note that this cost will be saved if a station has already been built up at the node i . Therefore, it is greater than zero only when a new initial station needs to be constructed in any node i (i.e., HBi = 0, CBi = 1). The second one ( 21) stands for the depreciation cost of pipeline installation, which is assumed to be linear with the transportation distance and the designed flow rate. Similar to the first term, this cost is greater than zero only when a new pipeline needs to be constructed from node i to node i (i.e., HPi, i = 0, Qs, k,1, i, i , e > 0 ). The third one ( 31s ) stands for transportation cost, which is linear with average flow rate and transportation distance. While the last one ( 41) represents the maintenance and inspection cost of pipelines, which is only considered to be related with the length of pipelines. (Zhang, Liang, Liao, Wu, & Yan, 2017). 1 2

Dk ns, k,2, i, e FCM2, e e E2 i I k K

3.1. Objective function

1 1

=

6 1

6 1

=

+

6 2

+

[RB (1

6 3

+

6 4s

s

(7a)

S

HCi ) CCi FSP + RB (1

HAi ) CAi FSP]

i I 6 2

(2a)

=

[RT VQCmaxi FTC + RT VQAmaxi FTA] i I

6 3

(2b)

=

RV VCmaxi FCC i I

120

(7b) (7c) (7d)

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

6 4s

=

Dk (VQCs, k, i FCGC + VQAs, k, i FCA ) i I k K

(

Dk

(7e)

k K

i I

e E1

Qs, k,1, i, i , e ).

VGs, k, i = VDGs, k, i + VSTs, k, i + VQLs, k, i + VQCs, k, i +

Qs, k,1, i, i , e

s

i I e E1

S, i

3.2. Transport option constraints As for PNG, pipeline flowrate must equal to zero if there is no pipeline between two nodes.

SPi, i M

Qs, k,1, i, i , e

i, i

I

(8)

s S k K e E1

ZNj, i, i , e ZBi, i M

s

S, j

J , i, i

I, k

K, e

VCs, k, i = VDCs, k, i + VQAs, k, i +

Qs, k,1, i , i, e + i

S, i

G

G VQGs, k , i

+

A A VQAs, k , i

Qs, k,2, i , i, e +

LT i

I e E2

L VQLs, k , i

s

S, i

I, k

s

VCs, k, i = VSCs, k, i +

I e E3

I, k

C VQCs, k , i

K

VLs, k, i

s

S, i

I, k

K

(16a)

VCmaxi

VCs, k, i

s

S, i

I, k

K

(16b)

VQLmaxi

VQLs, k, i

s

S, i

I, k

K

(17a)

VQGmaxi

VQGs, k, i

s

S, i

I, k

K

(17b)

VQCmaxi

VQCs, k, i

s

S, i

I, k

K

(17c)

VQAmaxi

VQAs, k, i

s

S, i

I, k

K

(17d)

VSTLmaxi

VSTLini +

NS VSTs, k, i

s

S, i

I, k

K

(18)

k=1

3.5. Flowrate constraints The determination of maximum flow rate is equal to Eq. (19).

Qmaxj, i, i , e

s

A

Qs, k, j, i, i , e

s

S, j

J , i, i

I, k

K, e

(19)

Ej

Flowrates must have boundaries for economic concerns.

(12)

K

I, k

k

K

G

) at this node.

Qs, k,3, i , i, e +

CT i

S, i

A

S, i

Peak shaving NG should be less than the maximum allowable storage volume of the underground gas storage.

Daily available CNG volume (VCs, k, i ) equals to the sum of the purchased CNG (VSCs, k, i ), transported CNG( CT i I e E3 Qs, k,3, i , i, e ) and C VQCs, k, i

s

VLmaxi

(11)

compressed NG volume (

Qs, k, j, i, i , e

The determination of the average maximum value of liquefaction, gasification, compression and regulation volume is equal to Eq. (19).

Similar with PNG, daily available LNG volume (VLs, k, i ) equals to the sum of the purchased LNG (VSLs, k, i ), transported LNG V ( LT i I e E2 Qs, k,2, i , i, e ) and liquefied NG volume ( L QLs,k ,i ) at this G node.

VLs, k, i = VSLs, k, i +

K

(15b)

(10)

K

I, k

The model always seeks for the minimum transportation cost. Additionally, the maximum available volume of each node has a positive linear correlation with the total cost. Thus the maximum available volume must be the minimum value under the constraints. In other words, this model only needs to determine the lower boundary of the maximum available volume to obtain the value of that in the objective function.

I e E1

I, k

S, i

i I e E3

In the season k under the simulation s , daily available NG volume at each node (VGs, k, i ) equals to the sum of purchased NG (VSGs, k, i ), transported NG from other nodes ( GT i I e E1 Qs, k,1, i , i, e ), gasified LNG ( G G VQGs, k, i ) and regulated CNG volume ( A A VQAs, k, i ) at this node. Here, GT, G and A represent the loss rate of PNG transportation, gasification, and regulation. G represents the volume ratio of LNG to NG, and A represents that of CNG to NG. GT

s

(15a)

3.3. NG state constraints

VGs, k, i = VSGs, k, i +

Qs, k,2, i, i , e i I e E2

(9)

Ej

(14)

K

VLs, k, i = VDLs, k, i + VQGs, k, i +

Based on path and equipment limitations, some area only allows special transport option and carrier to transport. Simultaneously, transportation is only allowed between two adjacent points to reduce useless calculation of the model. Let ZNj, i, i , e = 0 denotes the forbiddance of transporting under the option j from the node i to the node i by the carrier e , while ZBi, i = 1 denotes that the node i is adjacent to the node i . Therefore, the transportation volume Qs, k, j, i, i , e can be greater than zero only when ZNj, i, i , e = 1 and ZBi, i = 1, as stated in constraint (9).

Qs, k, j, i, i , e

I, k

Similarly, the volume of LNG and CNG is expressed in Eqs. (17) and (18). Taking LNG as an example, the available daily LNG volume at each node must equal to the total volume of consumption (VDGs, k, i ), gasification (VQGs, k, i VQLs, k, i ) and transfer under option LNG ( i I e E1 Qs, k,1, i, i , e ).

QLPmini, i

Qs, k, j, i, i , e

QLPmaxi, i

s

S, j

J , i, i

I, k

K, e

Ej (20)

3.4. Volume constraints

Flowrate variation among different seasons must result in a certain range for equipment and operation cost concerns.

Due to the production and import constraints, each node has available purchased NG, LNG and CNG volume range.

VSGs, k, i

VSGmaxk, i

s

S, i

I, k

K

(13a)

VSLs, k, i

VSLmaxk, i

s

S, i

I, k

K

(13b)

VSCs, k, i

VSCmaxk, i

s

S, i

I, k

K

(13c)

QREj, i, i , e

Qs, k, j, i, i , e K, e

Qs, k , j, i, i , e Ej

QREj, i, i , e

s

S, j

J , i, i

I, k (21)

3.6. Carrier constraints The number of tanker cars and carriers starting from each node shall be an integer. And the needed number of tanker cars and carriers is related to productive mooring time of transportation equipment, transportation route and speed, effective transport time, filling capacity

Considering volume conservation, in the season k under the simulation s , the available daily NG volume at each node must equal to the total volume of consumption (VDGs, k, i ), peak shaving (VSTs, k, i ), liquefaction (VQLs, k, i ), compression (VQCs, k, i ) and transfer under the option PNG 121

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

and average daily transportation flowrate.

CGi M

Qs, k, j, i, i , e (TAj, e + 2Li, i , j, e / RAj, e )

n s , k , j, i , e i

s

VTj, e NTj, e

I

K, e

S, j

{2, 3}, i

VTj, e NTj, e

i

S, j

I, k

CCi M

{2, 3}, i

I, k

K, e

CAi M +1 s

n maxj, i, e

s

S, j

{2, 3}, i

(23)

Ej

I, k

K, e

(Qs, k,1, i, i , e + Qs, k,1, i , i, e ) i I e E1

Ej

(24)

VSTs, k, i

i

I

s S k K

CSTi

ZSTi

i

I

(25)

(26a) (26b)

I

It will be required to build the corresponding stations if liquefaction, gasification, compression or regulation processes are needed at the node.

CLi M

VQLs, k, i s S k K

i

I

I

VQAs, k, i

i

I

(27c) (27d)

This paper optimizes an along-shore multi-state NG supply chain system. All the nodes and numbers are shown in Fig. 2. This block includes one LNG receiving station R1, ten NG consuming cities from C1 to C10, and two NG fields F1 and F2. Pipelines from F1 to C6 and from C6 to C7 have already existed, and there is a large NG liquefaction plant in C7. In this system, two NG fields and the LNG receiving stations are the sources and gasification stations have been built near the LNG receiving stations. Cities along the river can employ LNG ships for transportation, and all the LNG ships are equipped with 8 Forty-foot Equivalent Unit (FEU) containers. Inland cities can transport LNG or CNG by LNG road tankers. C2, C3, C7, and C10 can be chosen to build underground NG storages due to the construction and transportation conditions. A year is divided into four seasons with equal duration. Distance by land and water between each node are shown in Table 1 where the unit is km. Water distance is shown within the parentheses, and the land distance is shown at the outside. Depreciation rate of each equipment and carriers are shown in Table 2. Basic parameter of each carrier is shown in Table 3. Loss rates during transportation or process are shown in Table 4. The loss rate of liquefaction equipment is 2.1%, and when compared with the ratio of LNG production rate to input NG rate, the loss rate could be ignored. The same theory can apply to the gasification, regulation and compression loss rates.

It will be required to build a peak shaving underground storage if the node needs peak shaving. Only certain nodes can be built as peak shaving underground storages due to geologic limitations.

CSTi M

i

4.1. Basic data

It will be required to build a station if the node is the initial station of the whole transportation pipeline, and NG gets downloaded from the pipeline here. s S k K i

VQCs, k, i

(27b)

4. Example

3.7. Station construction constraints

CBi M

I

s S k K

The number of transportation equipment at each node should be lower than the upper boundary of this node.

n s , k , j, i , e

i

s S k K

(22)

Ej

Qs, k, j, i, i , e (TAj, e + 2Li, i , j, e / RAj, e )

n s , k , j, i , e

VQGs, k, i s S k K

(27a)

Fig. 2. Positions and numberings of nodes. 122

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Table 1 Distance between the nodes. Unit: km Node

F1

F2

C1

C2

C3

C4

C5

C6

C7

C8

C9

C10

R1 F1

1 276 (–) –

F2

–

890 (–) 896 (–) –

C1

–

–

173 (162) 1126 (–) 855 (–) –

C2

–

–

–

530 (542) 754 (–) 708 (–) 372 (383) –

C3

–

–

–

–

685 (664) 651 (–) 824 (–) 517 (502) 180 (119) –

C4

–

–

–

–

–

612 (–) 1 049 (–) 1 159 (–) 375 (–) 459 (–) 345 (–) –

C5

–

–

–

–

–

–

865 (–) 816 (–) 1 173 (–) 696 (–) 495 (–) 251 (–) 320 (–) –

C6

–

–

–

–

–

–

–

1 037 (1056) 252 (–) 805 (–) 882 (894) 510 (511) 399 (382) 804 (–) 598 (–) –

C7

–

–

–

–

–

–

–

–

920 (–) 361 (–) 640 (–) 778 (–) 321 (–) 363 (-) 790 (–) 649 (–) 171 (–) –

C8

–

–

–

–

–

–

–

–

–

606 (–) 858 (–) 307 (–) 552 (–) 428 (–) 574 (–) 864 (–) 917 (–) 689 (–) 521 (–) –

C9

–

–

–

–

–

–

–

–

–

–

236 (–) 1130 (–) 654 (–) 257 (–) 442 (–) 621(–) (–) 702 (–) 884 (–) 910 (–) 769 (–) 281 (–) –

669 (–) 675 (–) 413 (–) 565 (–) 306 (–) 411 (–) 765 (–) 760 (–) 491 (–) 260 (–) 200 (–) 477 (–)

Table 2 Depreciation rate.

Table 4 Loss rate.

Name

Depreciation Rate

Mode

Loss rate

Station Pipeline Ship/road tanker LNG/CNG storage infrastructure

14% 14% 17% 14%

LNG transportation LNG transportation CNG transportation Liquefaction Gasification Compression Regulation

0.1% 0.3% 0.3% 2.1% 0% 3% 0%

The purchase price of each source is shown in Table 5. Since the purchase price belongs to the interval uncertainty (Yang, Zhang, & Xiao, 2017), the price limits are listed here. NG fields and LNG receiving stations have the largest NG LNG supply volume respectively. Supposing the demands of consuming cities conform to a normal distribution, the mean demand of each city is shown in Fig. 3. The variance is 15% of mean value. The maximum of underground gas storage capacity is 5 × 108 m3 of which the base load is 30%. The volume ratio of LNG and CNG to NG is 620 and 300. The capacity of LNG liquefaction station per day is 1 × 106 m3, and the CNG compress station is 3 × 105 m3. And the rotation period of storage is 15 days.

Table 5 Basic parameters of each source.

4.2. Setting of cases

Source

Resource

Maximum supply volume (m3/d)

Price upper limit

Price lower limit

F1 F2 R1

NG NG LNG

2.70 × 107 3.50 × 107 –

0.70 CNY/m3 0.68 CNY/m3 5 410 CNY/t

0.63 CNY/m3 0.62 CNY/m3 5 250 CNY/t

double effect of the two factors on the model. Based on the discussion above, there will be four situations given for solution and discussion to pursue the effect of uncertainties on the model. Gurobi MILP solver is employed to solve all cases by Matlab R2014a. The CPU type we used is Intel (R) Core(TM) i5-7200U CPU 2.50 GHz.

In virtue of many kinds of uncertainties existing in the system, these uncertainties need to be further studied. Uncertainties in the system can be divided into three types: a) the effect of uncertain purchase prices on the model; b) the effect of uncertain demands on the model; and c) Table 3 Basic parameters of each carrier. Transportation facility

Filling capacity (m3)

Unit price (104 CNY)

Speed (km/h)

Filling time (h)

Effective running time (h)

LNG road tanker LNG ship CNG road tanker

47 864 15

150 8000 70

70 30 70

3.5 7 3.5

21 21 21

123

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Fig. 3. Demands of consuming cities.

Case 1: Determine the purchase prices and cities’ demands as their mean values, and the MILP method (i.e. s = 1) is applied for a solution. Case 2: Take uncertainty of purchase prices of inland NG and imported LNG into consideration and determine cities’ demands as the mean values. The uncertainty optimization model is applied for the solution. Case 3: Consider the uncertainty of cities’ demands and determine purchase prices as the mean values. The uncertainty optimization model is applied for the solution. Case 4: Consider the uncertainties of purchase prices and cities’ demands, and the uncertainty optimization model is applied for the solution. 5. Results and discussion

Fig. 4. Objective function of Case 4 for different numbers of simulations.

5.1. Results analysis for cases Simulation number has a tremendous impact on the model solving efficiency and effect. If the number is too small, the uncertainty effect on the model will not be sufficiently considered. If the number is too large, the model scale and solving time will get increased. To explore the influence of simulation number on the model solving effect, take Case 4 as an example and solve the model with the number gradually increasing. The implementation results are shown in Fig. 4. It is known the model solution tends towards stability and gets converged basically when the number is larger than 27. Hence, it is suggested that the

simulation number can be fixed as 30 which is easy to satisfy convergence. Table 6 displays the model scale and solving time of each case when the simulation number is 30. The supply chain construction is shown in Table 7 for different cases. The total cost is the lowest when ignoring the uncertainties. And the total cost is the highest when involving the uncertainties. Meanwhile, the price uncertainty has a greater influence on the entire cost while less on the construction scheme of the supply chain. On the contrary, the demand uncertainty has a greater influence on the

124

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Table 6 Calculation result.

Table 8 Purchase scheme of case 4.

Case

1

2

3

4

Unit: m3/d

Number of continuous variables Number of discrete variables Number of constraints CPU time (s)

4 010 3 003 33 735 211.3

84 630 3 003 719 355 18 326.5

84 630 3 003 719 355 23 658.2

84 630 3 003 719 355 16 589.6

Season

R1 LNG

F1 NG

F2 NG

Spring Summer Autumn Winter

13 10 11 14

2.59 × 107 2.65 × 107 2.67 × 107 2.69 × 107

3.48 × 107 3.47 × 107 3.49 × 107 3.49 × 107

Table 7 Supply chain construction of each case. Case

1

2

3

4

Length of pipeline (km) Number of liquefaction plant Number of the gasification plant Number of CNG plant Number of underground gas storage F(107 CNY)

2 917 1 1 4 2 12.773

2 917 1 1 4 2 13.224

3 097 1 1 4 2 12.981

3 097 1 2 4 2 13.355

228 929 446 305

to the nearby consuming cities C3, C5, C6, and C10. For the un-gasified LNG in R1, the part will be transported to C1, C2, and C3 along the river, while the other will be transported to inland cities with small demand through road tankers. When prices and demands fluctuate, there will be a small-scale gasification station constructed at C2, which is the location of geographical and pipeline network center, to satisfy the NG demand of its own as well as nearby cities. Taking the compression station cost, capacity and the demand volume and location of the city into consideration, it is suggested to build CNG compression station in C1, C3, C7, and C8. The NG is taken directly from pipelines at each compression station, then goes through desulfurization and dehydration process and finally is stored in the form of CNG. Each compression station must meet the local CNG demand as well as the transportation free CNG to nearby cities through road tankers. C2, C3, C7, and C10 are alternative locations of underground NG storages due to the geological limits. The optimal result suggests new underground NG storages in C3 and C10 for peak shaving in winter. From comparison with other cases, when considering the uncertainty of cities’ demands (comparing Fig. 5 with Fig. A1), a pipeline from C3 to C2 and a gasification plant at C2 will be added in the construction scheme. The pipeline connects the two pipeline networks (F1-C6-C7-C3-C5-C4 and F2-C8-C10-C2-C9-C1-R1) hence promotes the

construction scheme of the supply chain while less on the mean operation cost. The construction scheme of case 4 is shown in Fig. 5 and that of other cases are shown in Appendix. Gas fields F1 and F2 and gasification station R1 are the NG sources in this system. In the most simulations, F1 is the NG source for C3, C4, C5, C6 and C7, and F2 is the NG source for C1, C2, C8, C9, and C10. R1 gasifies most of the imported LNG and transports it to C1 through pipelines. In some special cases, it is possible for F1 to transport gas to C2 and C3. Therefore, there exists an NG pipeline constructed between the two places. In this way, all the nodes can be interconnected, thereby forming an NG pipeline network. Along-shore LNG receiving stations and liquefaction station C7 are the LNG sources in this system. C7 transports the liquefied NG from F1

Fig. 5. Construction scheme of case 4.

125

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Fig. 6. Transport volume of PNG in case 4.

Fig. 8. Transport volume of CNG in case 4.

flexibility of the integral transportation system. When considering the uncertainty of purchase prices (comparing Fig. 5 with Fig. A2), a gasification plant at C2 will be added in the construction scheme, which makes the state-transformation of NG in the supply chain more flexible. When prices fluctuate wildly, the supply chain manager can reduce its impact on the cost of supply through changing the purchasing ratios of NG, CNG, and LNG, then using state-transformation facilities to meet the NG demand of each state. 5.2. Results analysis for simulation To verify the model’s accuracy, we set a special simulation where the price of LNG and NG and demand of each city are substituted by their mean values. The optimal detailed transportation scheme in each case is worked out. Here, Case 4 is taken as an example. The purchase scheme is shown in Table 8 (The purchase schemes of other cases are shown in Appendix). The average transport volumes of the simulation by PNG, LNG, and CNG are shown in Figs. 6–8, and the average number of each kind of carrier is given in Fig. 9. The storage schemes of underground gas storages C3 and C10 are shown in Fig. 10. From the result, we find that there will be shifts of transport and storage plans in different seasons. Limited by the economical throughput, the throughput shifts of the pipeline in different seasons are relatively stable than that of LNG and CNG transportation modes. During the summer, C3 and C10 will storage amounts of natural gas and export during the winter. All the results can meet with constraints of volume conservation, downstream demand, and also the engineering process, which proves the practicability of the method.

Fig. 9. The number of carriers in case 4.

5.3. EOR analysis According to the determination method of economic operation range (EOR) for NG pipeline, and the LNG and CNG transport options referred in (Wan & Xie, 2015) (the detailed methodology is shown in Appendix), the economy of the transportation scheme is analyzed in this paper, shown as Fig. 11. Each point in the figure represents a transportation task in the result. The blue represents CNG transportation task and the green represents LNG transportation task. For

Fig. 7. Transport volume of LNG in case 4. 126

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

6. Conclusion This paper analyzes three common transport options, namely PNG, LNG, and CNG transportation, and summarizes the factors that affect the economy. This model takes the transportation cost as well as corresponding infrastructures and processing cost into consideration. PNG transportation cost includes initial station construction cost, pipeline installation cost, depreciation cost of compressors and other equipment, electricity cost and daily maintenance and inspection cost. LNG transportation cost includes depreciation cost of LNG road tankers and ships, transportation cost, daily operation and maintenance cost, and liquefaction and gasification cost. CNG transportation cost includes depreciation cost of CNG road tankers, transportation cost, daily operation and maintenance cost, and compression and regulation cost. This MILP method takes the minimum daily operation cost as the objective function based on the transportation, purchase, storage, construction and processing cost and the seasonal demand variation. The model is established considering constraints of transport options, volume, equipment and station limitations, and uncertainties of the purchase price and demand. And a MILP method is proposed and works out the optimal construction, purchase, transportation and storage scheme. A real case is presented to illustrate the application of the proposed model. Based on the calculation, this model carries out an optimal construction and the purchase, transportation, and storage scheme of each simulation. It is concluded that there will be the lowest total cost without involving uncertainties, while there will be the highest total cost with the consideration of uncertainties. And the price uncertainty has a greater influence on the entire cost while less on the construction scheme of the supply chain. On the contrary, the demand uncertainty has a greater influence on the construction scheme of the supply chain while less on the mean operation cost. The implementation results testify the accuracy, applicability, and efficiency of the model which can be further applied to the system construction of the NG supply chain. The approach proposed in this paper enables to consider uncertainties in the prices and demands of multi-state NG; however, the cities’ demand is related to NG price which is a crucial index to determine the mean value and variation of cities’ demands. Many factors decide the relationship between them, thereby contributing to complex analysis. Accordingly, we will proceed to further research on this issue in the future.

Fig. 10. Storage scheme of underground gas storages in case 4.

Fig. 11. EOR of different transport options.

example, there is a 100 m3/d CNG transportation task (equivalent to 2.5 × 104 m3/d natural gas) from C1 to C4 (375 km), hence a point with (375, 2.5) is drawn in the figure. The black line represents the margins of EORs. If the transportation task point falls into or close to the area of corresponding transportation mode EOR, the transportation plan could be regarded to be cost-effective. Since the PNG capacity is far greater than CNG and LNG and satisfies the pipeline EOR, there won’t be display. Most of the CNG and LNG can satisfy the EOR. In virtue of all cities demanding for CNG and LNG, the cost for separately constructing a single compression station and gasification plant will be higher, leading to some transportation tasks not belonging to but near the EOR. Hence, it can demonstrate the solution with good economy.

Acknowledgments This work was part of the Program of “Study on Optimization and Supply-side Reliability of Oil Product Supply Chain Logistics System” funded under the National Natural Science Foundation of China, grant number 51874325. The authors are grateful to all study participants.

Appendix A. Methodology of economic operation range determine The method we used is based on an existing research (Wan & Xie, 2015). In this paper, based on the empirical equation, unit costs for the three transportation modes (PNG, CNG, and LNG) are reviewed to obtain critical curves for economical transportation of these modes. The terms of transportation cost for each mode the paper proposed is the same with our paper (see Figs. A1 and A2). The unit volume transportation cost by PNG mode is concluded as: (A1)

f1s = 0.0734L/ qP , 4

3

where L is the transportation distance, km; qP is throughput by PNG transportation, 10 m /d. The unit volume transportation cost by LNG mode is concluded as: (A2)

f2s = 4.519 × 10 4L + 0.888/ qL + 0.926, 4

3

where qL is throughput by LNG transportation, 10 m /d (equivalent volume of natural gas).

127

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Fig. A1. Construction scheme of cases 1 and 2.

Fig. A2. Construction scheme of cases 3.

The unit volume transportation cost by CNG mode is concluded as: (A3)

f3s = 2.932 × 10 3L + 0.604/ qC + 0.361, 4

3

where qC is throughput by CNG transportation, 10 m /d (equivalent volume of natural gas) (see Tables A1–A3)

128

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al.

Table A1 Purchase scheme in case 1. Unit: m3/d Season

R1 LNG

F1 NG

F2 NG

1 2 3 4

13 10 11 13

2.55 × 107 2.59 × 107 2.64 × 107 2.64 × 107

3.50 × 107 3.50 × 107 3.50 × 107 3.50 × 107

275 810 059 748

Table A2 Purchase scheme in case 2. Unit: m3/d Season

R1 LNG

F1 NG

F2 NG

1 2 3 4

13 11 11 13

2.58 × 107 2.60 × 107 2.67 × 107 2.67 × 107

3.50 × 107 3.50 × 107 3.50 × 107 3.50 × 107

174 048 040 805

Table A3 Purchase scheme in case 3. Unit: m3/d Season

R1 LNG

F1 NG

F2 NG

1 2 3 4

13,175 10,817 11,037 14,628

2.57 × 107 2.62 × 107 2.68 × 107 2.68 × 107

3.50 × 107 3.50 × 107 3.50 × 107 3.50 × 107

Elia, J. A., Baliban, R. C., & Floudas, C. A. (2013). Nationwide, regional, and statewide energy supply chain optimization for natural gas to liquid transportation fuel (GTL) Systems. Industrial & Engineering Chemistry Research, 53(13), 5366–5397. Fernandes, L. J., Relvas, S., & Barbosa-Póvoa, A. P. (2013). Strategic network design of downstream petroleum supply chains: Single versus multi-entity participation. Chemical Engineering Research & Design, 91(8), 1557–1587. Ficco, G., Dell’Isola, M., Vigo, P., & Celenza, L. (2015). Uncertainty analysis of energy measurements in natural gas transmission networks. Flow Measurement & Instrumentation, 42(12), 58–68. Frick, M., Axhausen, K. W., Carle, G., & Wokaun, A. (2007). Optimization of the distribution of compressed natural gas (CNG) refueling stations: Swiss case studies. Transportation Research Part D Transport & Environment, 12(1), 10–22. Guandalini, G., Colbertaldo, P., & Campanari, S. (2017). Dynamic modeling of natural gas quality within transport pipelines in presence of hydrogen injections. Applied Energy, 185, 1712–1723. Halvorsen-Weare, E. E., Fagerholt, K., & Nnqvist, M. (2013). Vessel routing and scheduling under uncertainty in the liquefied natural gas business. Computers & Industrial Engineering, 64(1), 290–301. Hamedi, M., Farahani, R. Z., Husseini, M. M., & Esmaeilian, G. R. (2009). A distribution planning model for natural gas supply chain: A case study. Energy Policy, 37(3), 799–812. Kashani, A. H. A., & Molaei, R. (2014). Techno-economical and environmental optimization of natural gas network operation. Chemical Engineering Research & Design, 92(11), 2106–2122. Khan, M. I., Yasmeen, T., Shakoor, A., Khan, N. B., Wakeel, M., & Chen, B. (2016). Exploring the potential of compressed natural gas as a viable fuel option to sustainable transport: A bibliography (2001–2015). Journal of Natural Gas Science & Engineering, 31, 351–381. Kim, J., Seo, Y., Chang, D., & Yan, J. (2016). Economic evaluation of a new small-scale LNG supply chain using liquid nitrogen for natural-gas liquefaction. Applied Energy, 182, 154–163. Kou, P., Liang, D., & Gao, L. (2016). Distributed EMPC of multiple microgrids for coordinated stochastic energy management. Applied Energy, 185, 939–952. Kuwahara, N., Bajay, S. V., & Castro, L. N. (2000). Liquefied natural gas supply optimisation. Energy Conversion & Management, 41(2), 153–161. Li, X., Armagan, E., Tomasgard, A., & Barton, P. I. (2011). Stochastic pooling problem for natural gas production network design and operation under uncertainty. AIChE

References Ally, J., & Pryor, T. (2016). Life cycle costing of diesel, natural gas, hybrid and hydrogen fuel cell bus systems: An Australian case study. Energy Policy, 94, 285–294. Alves, F. D. S., Souza, J. N. M. D., & Costa, A. L. H. (2016). Multi-objective design optimization of natural gas transmission networks. Computers & Chemical Engineering, 93, 212–220. Azadeh, A., Raoofi, Z., & Zarrin, M. (2015). A multi-objective fuzzy linear programming model for optimization of natural gas supply chain through a greenhouse gas reduction approach. Journal of Natural Gas Science & Engineering, 26(3), 702–710. Baccanelli, M., Langé, S., Rocco, M. V., Pellegrini, L. A., & Colombo, E. (2016). Low temperature techniques for natural gas purification and LNG production: An energy and exergy analysis. Applied Energy, 180, 546–559. Bilgin, M. (2009). Geopolitics of European natural gas demand: Supplies from Russia, Caspian and the Middle East. Energy Policy, 37(11), 4482–4492. Boostani, A., Ghodsi, R., & Miab, A. K. (2010). Optimal location of compressed natural gas (CNG) Refueling station using the arc demand coverage model. Fourth Asia international conference on mathematical/analytical modelling and computer simulation (pp. 193–198). . Burel, F., Taccani, R., & Zuliani, N. (2013). Improving sustainability of maritime transport through utilization of Liquefied Natural Gas (LNG) for propulsion. Energy, 57(57), 412–420. Calvillo, C. F., Sanchez-Miralles, A., Villar, J., & Martin, F. (2016). Optimal planning and operation of aggregated distributed energy resources with market participation. Applied Energy, 182, 340–357. Chen, Y., He, L., Li, J., & Zhang, S. (2018). Multi-criteria design of shale-gas-water supply chains and production systems towards optimal life cycle economics and greenhouse gas emissions under uncertainty. Computers & Chemical Engineering, 109, 216–235. Choi, G. B., Lee, S. G., & Lee, J. M. (2016). Multi-period energy planning model under uncertainty in market prices and demands of energy resources: A case study of Korea power system. Chemical Engineering Research & Design, 114, 341–358. Choi, J., Realff, M. J., & Lee, J. H. (2005). Stochastic dynamic programming with localized cost-to-go approximators: Application to large scale supply chain management under demand uncertainty. Chemical Engineering Research & Design, 83(6), 752–758. Demirbas, A. (2002). Fuel properties of hydrogen, liquefied petroleum gas (LPG), and compressed natural gas (CNG) for transportation. Energy Sources, 24(7), 601–610.

129

Computers & Industrial Engineering 131 (2019) 115–130

H. Zhang, et al. Journal, 57(8), 2120–2135. Li, X., & Barton, P. I. (2015). Optimal design and operation of energy systems under uncertainty. Journal of Process Control, 30(1), 1–9. Li, J., Dong, X., Shangguan, J., & Hook, M. (2011). Forecasting the growth of China’s natural gas consumption. Energy, 36(3), 1380–1385. Li, J., Ma, X., Liu, H., & Zhang, X. (2018). Life cycle assessment and economic analysis of methanol production from coke oven gas compared with coal and natural gas routes. Journal of Cleaner Production, 185, 299–308. Li, X., Tomasgard, A., & Barton, P. I. (2017). Natural gas production network infrastructure development under uncertainty. Optimization and Engineering, 18(1), 35–62. Lima, C., Relvas, S., & Barbosa-Póvoa, A. (2018). Stochastic programming approach for the optimal tactical planning of the downstream oil supply chain. Computers & Chemical Engineering, 108, 314–336. Lin, W., Zhang, N., & Gu, A. (2010). LNG (liquefied natural gas): A necessary part in China's future energy infrastructure. Energy, 35(11), 4383–4391. Liu, E., Li, C., & Yang, Y. (2014). Optimal energy consumption analysis of natural gas pipeline. Scientific World Journal, 2014(4), 57–78. Lu, W., Su, M., Fath, B. D., Zhang, M., & Yan, H. (2016). A systematic method of evaluation of the Chinese natural gas supply security. Applied Energy, 165, 858–867. Marmidis, G., Lazarou, S., & Pyrgloti, E. (2008). Optimal placement of wind turbines in a wind park using Monte Carlo simulation. Renewable Energy, 33(7), 1455–1460. Marseguerra, M., Zio, E., & Podofillini, L. (2002). Condition-based maintenance optimization by means of genetic algorithms and Monte Carlo simulation. Reliability Engineering & System Safety, 77(2), 151–165. Nosratabadi, S. M., Hooshmand, R. A., & Gholipour, E. (2016). Stochastic profit-based scheduling of industrial virtual power plant using the best demand response strategy. Applied Energy, 164, 590–606. Oliveira, F., Gupta, V., Hamacher, S., & Grossmann, I. E. (2013). A Lagrangean decomposition approach for oil supply chain investment planning under uncertainty with risk considerations. Computers & Chemical Engineering, 50(50), 184–195. Resat, H. G., & Turkay, M. (2015). Design and operation of intermodal transportation network in the Marmara region of Turkey. Transportation Research Part E Logistics & Transportation Review, 83, 16–33. Ríos-Mercado, R. Z., & Borraz-Sánchez, C. (2015). Optimization problems in natural gas transportation systems: A state-of-the-art review. Applied Energy, 147, 536–555. Rodger, J. A. (2014). A fuzzy nearest neighbor neural network statistical model for predicting demand for natural gas and energy cost savings in public buildings. Expert Systems with Applications, 41(4), 1813–1829. Samsatli, S., & Samsatli, N. J. (2015). A general spatio-temporal model of energy systems with a detailed account of transport and storage. Computers & Chemical Engineering, 80, 155–176. Sanaye, S., & Mahmoudimehr, J. (2013). Optimal design of a natural gas transmission network layout. Chemical Engineering Research & Design, 91(12), 2465–2476. Shi, Y., Zhou, L., Xu, Y., Zhou, H., & Shi, L. (2018). Life cycle cost and environmental assessment for resource-oriented toilet systems. Journal of Cleaner Production, 196,

1188–1197. Soldo, B. (2012). Forecasting natural gas consumption. Applied Energy, 92(4), 26–37. Spoladore, A., Borelli, D., Devia, F., Mora, F., & Schenone, C. (2016). Model for forecasting residential heat demand based on natural gas consumption and energy performance indicators. Applied Energy, 182, 488–499. Thomas, S., & Dawe, R. A. (2003). Review of ways to transport natural gas energy from countries which do not need the gas for domestic use. Energy, 28(14), 1461–1477. Touretzky, C. R., Mcguffin, D. L., Ziesmer, J. C., & Baldea, M. (2016). The effect of distributed electricity generation using natural gas on the electric and natural gas grids. Applied Energy, 177, 500–514. Tsatsaronis, G., & Morosuk, T. (2010). Advanced exergetic analysis of a novel system for generating electricity and vaporizing liquefied natural gas. Energy, 35(2), 820–829. Verdejo, H., Awerkin, A., Saavedra, E., Kliemann, W., & Vargas, L. (2016). Stochastic modeling to represent wind power generation and demand in electric power system based on real data. Applied Energy, 173, 283–295. Wan, F., & Xie, D. (2015). Application scope and economical comparison of natural gas transportation. Oil & Gas Storage and Transportation, 34(7), 709–713. Wang, B., Liang, Y., Zheng, T., Yuan, M., & Zhang, H. (2018). Multi-objective site selection optimization of the gas-gathering station using NSGA-II. Process Safety and Environmental Protection, 119, 350–359. Wang, B., Yuan, M., Zhang, H., Zhao, W., & Liang, Y. (2018). An MILP model for optimal design of multi-period natural gas transmission network. Chemical Engineering Research and Design, 129, 122–131. Xia, T., Cardin, M. A., Ranjbar-Bourani, M., & Caunhye, M. A. (2015). A decision-making tool to design a flexible liquefied natural gas system under uncertainty. IEEE international conference on systems, man, and cybernetics (pp. 789–794). . Yang, Y., Zhang, S., & Xiao, Y. (2017). Optimal design of distributed energy resource systems based on two-stage stochastic programming. Applied Thermal Engineering, 110, 1358–1370. Zamarripa, M. A., Aguirre, A. M., Méndez, C. A., & Espuña, A. (2013). Mathematical programming and game theory optimization-based tool for supply chain planning in cooperative/competitive environments. Chemical Engineering Research & Design, 91(8), 1588–1600. Zhang, Q., Li, Z., Wang, G., & Li, H. (2016). Study on the impacts of natural gas supply cost on gas flow and infrastructure deployment in China. Applied Energy, 162, 1385–1398. Zhang, J., Li, H., Yang, S., Li, X., & Qian, Y. (2014). Investigation on life-cycle cost of coalbased synthetic natural gas (SNG). In J. J. Klemeš, P. S. Varbanov, & P. Y. Liew (Vol. Eds.), Computer aided chemical engineering: vol. 33, (pp. 991–996). Elsevier. Zhang, H. R., Liang, Y. T., Liao, Q., Wu, M. Y., & Yan, X. H. (2017). A hybrid computational approach for detailed scheduling of products in a pipeline with multiple pump stations. Energy, 119, 612–628. Zhang, H., Liang, Y., Liao, Q., Yan, X., Shen, Y., & Zhao, Y. (2017). A three-stage stochastic programming method for LNG supply system infrastructure development and inventory routing in demanding countries. Energy, 133, 424–442.

130