A MILP model for the detailed scheduling of multiproduct pipelines with the hydraulic constraints rigorously considered

Computers and Chemical Engineering 130 (2019) 106543

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

A MILP model for the detailed scheduling of multiproduct pipelines with the hydraulic constraints rigorously considered Xingyuan Zhou a, Yongtu Liang a, Xin Zhang a, Qi Liao a, Song Gao b, Wan Zhang a, Haoran Zhang c,∗ a

National Engineering Laboratory for Pipeline Safety/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 Viterbi School of Engineering, University of Southern California, 3650 McClintockAve, Los Angeles, CA 90089, United States c Center for Spatial Information Science, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8568, Japan b

a r t i c l e

i n f o

Article history: Received 13 December 2018 Revised 6 August 2019 Accepted 15 August 2019 Available online 7 September 2019 Keywords: Multiproduct pipeline Detailed scheduling Hydraulic constraints Mixed-integer nonlinear programming Continuous-time representation

a b s t r a c t Multiproduct pipelines are the most effective and important way to transport refined products from refineries to the downstream market. The detailed scheduling of a multiproduct pipeline with the hydraulic constraints considered plays a vital role in the pipeline’s safe and economic operation. Aiming at this issue, this paper proposes a continuous-time mixed-integer nonlinear programming (MINLP) model taking the sum of the pumps’ running and start/stop costs as the objective function. Based on a method proposed in a previous study for the rigorous description of the hydraulic loss changes in multiproduct pipelines, the pipeline hydraulic constraints and pump-related functions are considered directly and rigorously in the model. Other actual field processing constraints, such as the pipeline shutdown and contaminated oil flowrate, are also considered. To deal with the nonlinear relationship between the flowrate and pressure/volume, a piecewise linearization method is adopted, and the pipeline flowrate is graded according to the divided intervals. Finally, the application of the established model is successfully conducted on a real-world multiproduct pipeline through two case studies. A comparison with a previous discrete-time model is also performed to verify this model’s applicability, accuracy, and superiority. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. Background As the most commonly used method for refined products transportation, multiproduct pipelines play a significant role in the oil supply chain and are of vital importance for satisfying the downstream market demand (Laínez et al., 2009) and ensuring the energy supply security (Zhou et al., 2019a). The optimal scheduling of a multiproduct pipeline would have a great impact on transport efficiency and safety and is the critical issue for its operation and management. Due to the small and fluctuating demand for various refined products, a multiproduct pipeline generally delivers several kinds of refined products sequentially, which results in a complex pipeline operation. The adjacent batches, including different refined products, are pumped back-to-back into the pipeline and directly contact each other (Mostafaei et al., 2016), therefore, contaminated oil is unavoidably formed. To reduce the deprecia-

∗

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

https://doi.org/10.1016/j.compchemeng.2019.106543 0098-1354/© 2019 Elsevier Ltd. All rights reserved.

tion loss, certain refined products are forbidden to be adjacent during batch sequencing, minimum flowrate limits must be satisfied, and the shut-down operation of pipeline segments are not permitted if batch interfaces exist (Relvas et al., 2009). Meanwhile, as a pressure-tight system, the hydraulic characteristics of multiproduct pipelines change in a complicated manner during the pipeline operation (Qiu et al., 2018; Zhou et al., 2019b). More specifically, with the batch interfaces moving in the pipeline, the allocation of refined products changes gradually; with frequent injection/download operations, the pipeline flowrate also varies frequently. All of these result in the complex changes in a pipeline’s hydraulic loss over time, and thus influence the determination of the pump schedule with the minimum costs and carbon emissions (Rejowski and Pinto, 2008). Furthermore, with the continuous construction of multiproduct pipelines, longer pipelines with more delivery and pump stations and more complicated structures require the pipeline schedules to be more meticulous and reliable for the safe and economic operation (Kirschstein, 2018). With a given aggregate plan, the detailed scheduling of multiproduct pipelines determines the pipeline flowrate and injection/delivery flowrate as well as the accurate time for pump

2

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

λi

Nomenclature Sets and indices t ∈ T = {1, 2, . . . , tmax } i ∈ I = {1, 2, . . . , imax } j ∈ J = {1, 2, . . . , jmax } k ∈ Ki = {1, 2, . . . , kmax i } o ∈ O = {1, 2, . . . , omax } a ∈ A= {1, 2, . . . , amax }

Set of time windows. Set of stations along the pipeline. Set of batches. Set of pumps in station i. Set of oil types. Set of flowrate intervals for piecewise linearization.

Continuous parameters CE Industrial electricity price (CNY/kW•h). CU Unit costs for pump start/stop (CNY). VDEi, o Demand volume of delivery station i for oil o (m3 ). μ Maximum allowable deviation range between the actual download and the demand volume. τ The time horizon (h). g Gravitational acceleration (m/s2 ). Xi Volume coordinate of station i (m3 ). lZ The total length of the pipeline (m). qDmax i , qDmin i Maximum and minimum download flowrates of station i (m3 /h). qPmax i , qPmin i Maximum and minimum flowrates of pipeline segment (i, i + 1) (m3 /h). qXmini, j Minimum flowrate of pipeline segment (i, i + 1) when the upper coordinate of batch j exists (m3 /h). qJmax , qJmin Maximum and minimum injection flowrates (m3 /h). RFj Density of the refined product corresponding to batch j (kg/ m3 ). υF j Viscosity of the refined product corresponding to batch j (m2 /s). β, m Constant coefficients in the Liebenson (Лейбензон) equation, which are related to the flow regime (0 < m < 1). Hi, k Pumping head of pump k at station i (if station i is not a pump station, Ki = 0, therefore, there are no Hi, k ) (m). SAi Cross-sectional area of pipeline segment (i, i + 1) (m2 ). di Inner diameter of pipeline segment (i, i + 1) (m). vDmin i , vDmax i Maximum and minimum download volumes of station i in a time window (m3 ). vPmin i , vPmax i Maximum and minimum flowed volumes of pipeline segment (i, i + 1) in a time window (m3 ). vImin i , vImax i Maximum and minimum injection volumes in a time window (m3 ). wPmaxi.k Maximum output work of pump k at station i in a time window (kW). pFmaxi, j Maximum hydraulic loss of pipeline segment (i, i + 1) for batch j (MPa). pOmax i , pOmin i Maximum and minimum outlet pressures of station i (MPa). pImax i , pImin i Maximum and minimum inlet pressures of station i (MPa). qPUmax i, k , qPUmin i, k Maximum and minimum pump flowrates of pump k at pump station i (m3 /h).

qAmax a , qAmin a

The tilt coefficient of pipeline segment (i, i + 1). Upper and lower limits of flowrate interval a (m3 /h).

Binary parameters γ j, o If batch j corresponds to refined product o, γ j,o = 1; otherwise, γ j,o = 0. ψj If the shut-down operation is forbidden when the interface of batch j and j − 1 exists in the pipeline segments, ψ j = 1; otherwise, ψ j = 0. Continuous variables WPt, i.k Output work of pump k at station i in time window t (kW). TLt Start time of time window t (h). TAt Time length of time window t (h). COHt, j Upper coordinate of batch j at the start time of time window t (m3 ). VDOt, i, j Download volume of batch j at station i in time window t (m3 ). VNJt, j Injection volume of batch j in time window t (m3 /h). VPFt, i Flowed volume of pipeline segment (i, i + 1) in time window t (m3 ). QPFt, i Flowrate of pipeline segment (i, i + 1) in time window t (m3 ). LBt, i, j Length of batch j within pipeline segment (i, i + 1) at the start time of time window t (m). PFSJt, i, j Hydraulic loss of batch j within pipeline segment (i, i + 1) at the start time of time window t (MPa). PFEJt, i, j Hydraulic loss of batch j within pipeline segment (i, i + 1) at the end time of time window t (MPa). PFSt, i Hydraulic loss of pipeline segment (i, i + 1) at the start time of time window t (MPa). PFEt, i Hydraulic loss of pipeline segment (i, i + 1) at the end time of time window t (MPa). Pt, i, k Provided pressure of pump k at station i in time window t (MPa). PHt, i Provided pressure of station i in time window t (MPa). PINt, i , POUt, i Inlet and outlet pressure of station i at the start time of time window t (MPa). PIt, i , POt, i Inlet and outlet pressure of station i at the end time of time window t (MPa). Binary variables BDOt, i, j If station i is downloading batch j in time window t, BDOt,i, j = 1; otherwise, BDOt,i, j = 0. BFTt, i, j If batch j is flowing through station i in time window t, BFTt,i, j = 1; otherwise, BFTt,i, j = 0. BPAt, i, j If the upper coordinate of batch j has exceeded station i at the start time of time window t, BPAt,i, j = 1; otherwise, BPAt,i, j = 0. BFAt, i, j If the upper coordinate of batch j is within pipeline segment (i, i + 1) in time window t, BFAt,i, j = 1; otherwise, BFAt,i, j = 0. BIDt, i If pipeline segment (i, i + 1) is idle in time window t, BIDt,i = 1; otherwise, BIDt,i = 0. BSP t, i, k, j If batch j is flowing through station i and pump k is turned on in time window t, BSP t,i,k, j = 1; otherwise, BSP t,i,k, j = 0.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

BCP t, i, k

BAQ t, i, a

If the on/off status of pump k at pump station i in two adjacent time windows is different, BCP t,i,k = 1; otherwise, BCP t,i,k = 0. If the flowrate of pipeline segment (i, i + 1) is within interval a in time window t, BAQ t,i,a = 1; otherwise, BAQ t,i,a = 0.

start/stop and valve open/close on the basis of meeting the aggregate plan (Cafaro et al., 2012). If the hydraulic constraints are also rigorously and directly considered, the pressure distribution along the pipeline and the exact pump power consumption can also be determined. As for the relationship between the detailed scheduling and aggregated scheduling for a multiproduct pipeline, the detailed scheduling refines the aggregate plan and ensures that the formulated schedule can be executed directly by the pipeline operators. In other words, the detailed scheduling serves as the lower level of the aggregate scheduling. According to the market demands, refinery production capacity and the initial inventory, the aggregate plan including the injection batch sequence, size and the demands of delivery stations for each batch in a short term are solved (Cafaro et al., 2011; Liao et al., 2018a). All of these are input parameters of the detailed scheduling model, in which the constraints for contaminated oil reduction, hydraulic loss calculation, and other processing limits need to be considered. In this way, nonlinear items and complicated logical relationships between variables are usually ineluctable in the detailed scheduling model (Meira et al., 2018). 1.2. Related work The detailed scheduling of multiproduct pipelines is a common issue, and considerable efforts have been made based on discrete or continuous-time representation. The discrete-time representation divides a time period into several time windows with some known time nodes before model solving so that the establishment of the model can be easier. Rejowski and Pinto (2003) put forward two discrete-time models that divided the pipeline into packs. The oil depot’s inventory levels, product distributions, and optimal ordering can be generated through solving the models. Through the division of the planning horizon and products, Herrán et al. (2010) performed the short-term operational planning of multiproduct pipeline systems. Different from the discretetime representation, the continuous-time representation takes the actual times of events as the time nodes, and the model’s solving efficiency can be greatly enhanced. Therefore, most subsequent studies are based on continuous-time representation. Cafaro and Cerda (2008) established a continuous-time model for the dynamic scheduling of a multiproduct pipeline with the delivery due dates and the change in market environment considered. Taking the complex storage policy associated with the product flow into account, Castro (2010) proposed a continuous-time model for the short-term multiproduct scheduling based on the Resource-Task Network. Most previous studies were on the scheduling of a single pipeline with multiple depots (Cafaro et al., 2012; MirHassani et al., 2013; Zhang et al., 2018), which is the most common situation in the real-world. Some research studies on other pipeline structures such as the tree-like structure and mesh structure have also been conducted. Castro and Mostafaei (2019) introduced a continuous-time mixed-integer linear programming (MILP) model for the short-term scheduling of tree-like multiproduct pipelines, and the constraints of enforcing forbidden product sequences and ensuring no batch interface exits in the idle segments are considered. Aiming at the scheduling of a mesh-like multiproduct pipeline network, Cafaro and Cerda (2012) presented an MILP

3

model that solved the product sequence, timing of batch inputs, lot sizes, and flow to the terminals at the same time. As for the model solving strategy, the MILP method was the most widely used (MirHassani and Jahromi, 2011; Hossein Mostafaei et al., 2015; Relvas et al., 2013). Heuristic rules (Liao et al., 2018b; Relvas et al., 2009) and artificial intelligence algorithms (Chen et al., 2017a; Zhang et al., 2016) were also adopted for optimal or suboptimal solving results within an acceptable time. As an important sub issue of multiproduct pipeline scheduling, some studies have been done on pump scheduling with a given delivery schedule (Abbasi and Garousi, 2010; Ferreira and Vidal, 1984; Liang et al., 2004). Because of the complexity of the hydraulic pressure changes with time during the multi-batch sequential transportation process, most previous studies were based on the discrete-time representation method, and thus the relationship between time and other variables could be simplified. Liang et al. (2012a) and Liang and Liu (2011) established discrete-time hydraulic optimization models for a single multiproduct pipeline and for a pipeline network. These models included energy equilibrium constraints, pressure constraints, and so on. The variance in regional electricity prices was also considered, and the models were solved by the dynamic programming method. However, the influences of time on the pipeline conditions were not taken into account. Barreto et al. (2004) established a database for the pipeline pumping energy consumption, considering all possible pumping arrangements and viable flowrate ranges at discrete time nodes. Significant cost savings were obtained by applying the developed methodology to the ORBEL II pipeline in Brazil. Even though the discrete-time representation reduces the difficulty of model establishing and simplifies the coupling relationships of variables, with the shortening of the time steps, the model scale will increase exponentially, and the solving efficiency will be significantly reduced if the mathematical programming algorithms are still employed. To efficiently solve the established discrete-time pump scheduling model for large-scale multiproduct pipelines, Liang et al. (2017) proposed an integrated solving strategy that adopts dynamic programming to solve the pump schedule at the non-critical time points and uses the ant colony algorithm (ACO) to solve the pump schedule at the critical time points. The discrete-time models are limited in that all pump operations must be conducted at these discrete time nodes, different from that in the pipeline’s actual working conditions. To solve this problem fundamentally, Zhou et al. (2019c) proposed a method that rigorously describes the hydraulic loss changes in multiproduct pipelines under time continuity and established a hybrid-time MILP model to achieve the accurate control of the pumps and confirm the maximum energy-saving. Many studies have been done that aim at the multiproduct detailed scheduling issue considering the hydraulic constraints. Rejowski and Pinto (2005) and Rejowski and Pinto (2008) developed a hydraulic formulation for multiproduct pipelines and proposed an MINLP model for the pipeline scheduling problem with pumping yield rate constraints based on the continuous-time representation. Cafaro et al. (2013) introduced an innovative piecewise linear formulation for the calculation of energy loss varying with flowrate and set up an MILP model to minimize the total operating costs of the pipeline. Cafaro et al. (2015) rigorously calculated the energy consumption by introducing nonlinear equations into the established MINLP model for the detailed scheduling of single-source pipelines, which were solved by the GAMS– DICOPT algorithm. However, failing to consider the hydraulic characteristics of multiproduct pipelines, they assumed the constant physical properties of an average refined product in the pipeline for simplicity, which would lead to errors in calculating the hydraulic loss. Zhang et al. (2017a) established an MINLP model that was solved by a proposed ACO-SM algorithm with the consider-

4

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

ation of the nonlinear hydraulic related objective function; however, it did not take the hydraulic constraints directly into account. Liao et al. (2018c) proposed an MILP model for the multiproduct pipeline in pressure control mode, and the pressure and pump constraints are considered. However, it is based on the discrete-time representation and takes the average of the physical properties if there are two kinds of refined products in a pipeline segment simultaneously. Some other studies took the hydraulic constraints and pumping costs into consideration by employing heuristic rules. Chen et al. (2017b) demonstrated that the more stable the pumping rate of a pipeline segment was, the lower the frictional loss for pumping products was, and took the minimum summation of the pumping rate variations in each pipeline segment as the objective function. Based on the flowrate database, Liao et al. (2018a) related the pressure constraints under each pump scheme to the flowrate limits and saved them in the database. If a pump scheme is selected, its corresponding flowrate limits must be satisfied. It is a heuristic way to convert the pressure constraints to the flowrate constraints. Based on the literature review, most previous studies on the detailed scheduling of multiproduct pipelines did not establish continuous-time models that consider the hydraulic constraints directly, because it is difficult to rigorously describe the pipeline hydraulic loss changing with time and to efficiently and accurately solve the mathematical models with nonlinear items introduced by hydraulic constraints. Aiming to solve this issue, this paper establishes a continuoustime detailed scheduling model, in which the pipeline’s hydraulic characteristics are rigorously described under time continuity and the hydraulic loss is accurately calculated. The pipeline’s hydraulic constraints (i.e., the energy equilibrium equation, pressure limits, pipeline hydraulic loss, pumping head equation, and so on) are directly coupled with the common detailed scheduling constraints (i.e., batch tracking constraints, delivery constraints, flowrate limits, and so on). These are the most significant differences from the existing detailed scheduling models. 1.3. Contributions of this work The contributions of this paper are listed as follows: (1) A continuous-time MINLP model is proposed for the detailed scheduling of a multiproduct pipeline with the hydraulic constraints rigorously considered, and a piecewise linearization method that grades the pipeline flowrate is adopted to convert the MINLP model to a MILP model for efficient and accurate model solving. (2) In calculating the pipeline hydraulic loss, this paper rigorously considers the differences in the physical properties of different refined products and the length of each kind of refined product in the pipeline at different times. Due to conditional linear change in the pipeline hydraulic loss with time, the inlet/outlet pressure of each station at both the start and end times of a continuous time window are constrained to ensure the safe operation of the pipeline at all times. (3) Other actual field technical and operational constraints of the multiproduct pipeline, such as specified flowrate lower limits and no shut-down operation when certain contaminated products exist in the pipeline, frequent pump start/stop avoidance, and so on, are also considered in the model. (4) The established model is successfully applied to a real-world multiproduct pipeline in China, and its superiority and applicability are examined in detail by comparing it with a previous model.

1.4. Paper organization Section 2 introduces the composition and characteristics of a multiproduct pipeline, illustrates the problem to be solved, and provides the model requirements. Section 3 establishes the MINLP model for the detailed scheduling of multiproduct pipelines with the hydraulic constraints rigorously considered and gives the objective function and constraints of the model. Due to the difficulty of solving the MINLP model accurately and directly, a piecewise linearization method is adopted in this paper and introduced in Section 4. In Section 5, the application of the established model is performed on a real-world multiproduct pipeline in China. The solution results are compared with that of a previous study to prove the model’s practicability and superiority. Conclusions are provided in Section 6. 2. Methodology 2.1. Problem description 2.1.1. Multiproduct pipeline A typical multiproduct pipeline is shown in Fig. 1. Various refined products, shown in different colors, are transported sequentially in batches from the injection station (i.e., IS) to the delivery stations (i.e., DS1, DS2, DS3, and TS) through the pipeline system. The injection station connects with the upstream refineries or large-scale oil depots. The delivery stations connect with the downstream local markets, and the planned download volumes are determined according to the actual regional demand. Multiple pump stations (i.e., IS, P1, and P2) are set up along the pipeline to provide energy for the fluid moving in the pipeline. Pump stations can either coincide with delivery stations or be set between two delivery stations. The colored pumps indicate that this pump is running, and the corresponding refined products are flowing through. The green line below is the hydraulic grade line, which indicates the residual pressure of the fluid flowing inside the pipeline (i.e., hydrodynamic pressure). For the same refined product, if the pipeline flowrate does not change, the hydrodynamic pressure will decrease linearly as shown in Fig. 1. However, the download operations of the delivery stations will result in pipeline flowrate changes, and the simultaneous transportation of various refined products will lead to changes in the physical properties. Both of these result in the turning points in the hydraulic grade line. Meanwhile, with the batch interfaces moving in the pipeline, the start and end of the download/injection operations, and the changes in download/injection flowrate, the hydraulic condition of the pipeline will change complicatedly with time. 2.1.2. Detailed scheduling of multiproduct pipeline considering the hydraulic constraints According to the demand data of the regional market, the demand volume of each delivery station for various refined products is determined first. Based on the demand volume, an aggregate pipeline operating plan including the injection batch sequence and volume is calculated. With the above data, the detailed scheduling of the multiproduct pipeline decides the exact start and end times and the flowrate of injection/download operations to satisfy the downstream demands as much as possible under various technical and operational constraints. Because of the difficulty of describing the pipeline hydraulic condition changes accurately and the nonlinear relationship between flowrate and pressure, most previous studies on the detailed scheduling of multiproduct pipelines do not consider the hydraulic constraints rigorously and solve the pump operation schedule at the same time. On the contrary, the pipeline delivery and pump schedule are determined separately.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

5

Fig. 1. Multiproduct pipeline with multiple pump stations.

That is, without regard to the pressure and pump-related constraints, the delivery schedule is optimized first. Some heuristic rules are adopted to avoid the hydraulic constraints being directly coupled into the model. For example, because the pressure calculation is closely related to flowrate, the pressure limits are thus converted into the pipeline flowrate constraints according to the onsite operation experience. This means that if the flowrate constraints are satisfied, the pressure limits must be met, and the fluid can flow to the terminal safely. However, since the converted flowrate constraints usually cannot match the original pressure limits perfectly, the flowrate constraints are generally much more conservative, which may lead to wasted pump capacity. Meanwhile, the physical properties of different kinds of refined products vary greatly, and their hydraulic losses at the same flowrate are also quite different. The determined pipeline flowrate upper limit needs to ensure that all kinds of refined products with different physical properties can reach the terminal safely, and the underpressure accidents must not occur. In this way, the converted upper limit is usually much lower than that when taking the pressure constraints into consideration. Then, based on the solution results, the pump schedule is optimized in the second level using various algorithms. Although the relationships between the flowrate and pressure variables can be separated in this way and the solving efficiencies of the detailed and pump scheduling models can be greatly enhanced, the solved delivery and pump schedules do not couple well and cannot guide the pipeline to run efficiently and safely. This is because the pipeline delivery schedule in the first level can influence the pump schedule in the second level, but the reverse is not possible. Although the delivery schedule and pump schedule are the optimal solutions of the models, the combination of the two schedules may not be optimal for the entire detailed scheduling issue. The above problems can be effectively solved if the delivery schedule and pump schedule are optimized simultaneously through establishing a unified model. This model includes both the injection/download and batch moving constraints and the pressure and pump-related constraints. In this way, the solution results are globally optimal. The schematic diagram is shown as Fig. 2. As we can see from Fig. 2, the pipeline detailed scheduling issue with the hydraulic constraints rigorously considered includes three parts, namely, the normal pipeline detailed scheduling, the conversion of certain variables, and the pump scheduling. With the normal constraints, such as flowrate, delivery, and batch tracking considered, the pipeline delivery schedule can be solved, as the previous studies did. The calculation of the pipeline hydraulic loss is the basis of pump scheduling. To calculate the hydraulic loss accurately, some variable conversions (i.e., volume convert-

ing to flowrate, volume coordinate converting to length) are also needed. Then, the optimal pump schedule can be obtained based on the pipeline hydraulic conditions. In the previous studies, only the pipeline delivery schedule affects the pump schedule, and the pump schedule is based on the delivery schedule as mentioned above. With the establishment of the pipeline detailed scheduling model considering the hydraulic constraints rigorously, both the delivery schedule (i.e., flowrate, volume, etc.) and the pump schedule (i.e., pressure, output work, etc.) are the variables in the model. They can interact with each other. The black dashed line indicates that the pump schedule can also influence the determination of the pipeline delivery schedule. 2.1.3. Time representation method and nonlinear items of the model As mentioned above, time representation methods mainly include the discrete-time and continuous-time representations. The detailed scheduling model of multiproduct pipelines with the hydraulic constraints considered contains many constraints, such as delivery, flowrate, pump-related, and pressure constraints. If the model is based on the discrete-time representation, as in some previous studies, the model scale will grow sharply with the increase in pipeline scale, transport volume, and scheduling horizon, and a feasible solution cannot be found within a reasonable time. At the same time, differences between the model and the pipeline’s actual working condition will be caused due to the discretization of time. Furthermore, changes in the variables between two adjacent time nodes, such as batch interface locations, and the pipeline hydraulic loss cannot be considered. Even though this problem can be alleviated through shortening the time step, the model scale will increase exponentially. If the continuous-time representation is adopted, the above problems can be completely solved. However, this will cause other problems: (1) how to select time nodes (2) how to describe the relationship between the time and other variables under time continuity (3) how to deal with the nonlinear items caused by the unknown lengths of time. All these issues are considered in the model established in this paper and will be discussed below. The mathematical model for the detailed scheduling of a multiproduct pipeline coupled with hydraulic constraints can be expressed in the compact form as follows.

min

e,ta,p,q,v,l

f (e, ta, p, q, v, l )

s.t. ψL (ta, p, q, v, l ) ≤ 0 φL (ta, p, q) = 0 φNW (e, q ) = 0 φNP ( p, q ) = 0 φNV ( p, q ) = 0

(1)

6

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. 2. Schematic diagram of multiproduct pipeline detailed scheduling with the hydraulic constraints rigorously considered. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

where f represents the objective function of the model; ta is the vector of time variables (continuous-time representation), p is the vector of pressure variables, q is the vector of flowrate variables, v is the vector of volume variables, l is the vector of length variables, e is the vector of work variables; ψ L are the linear inequality constraints; φ L are the linear equality constraints; φ NW indicate the nonlinear relationship between the flowrate and the output work of pumps; φ NM indicate the nonlinear relationship between the flowrate and pressure; φ NV indicate the nonlinear relationship between the flowrate and volume; Eqs. (2) and (3) correspond to φ NM . Eq. (2) is the pump hydraulic equation, while Eq. (3) is the pipeline segment hydraulic equation. The pipeline flowrate Qt, i is a variable and its 2-m power results in Eq. (2) as a nonlinear equation. It is the same for Eq. (3), and the multiplication of Lt, i, j and Qt, i also leads to the nonlinearity. That is because Lt, i, j is also a variable and will change gradually with the batch interface moving in the pipeline. Eq. (4) corresponds to φ NW . Eq. (4) calculates the output work of the pumps and Eq. (5) indicates the relationship between the flowrate and volume. In addition to the 3-m power of Qt, i , the multiplication of Tt and Qt, i will also result in the nonlinearity, since Tt is also a variable due to the continuous-time representation.



 2−m

PHt,i,k = R j g ai,k − bi,k Qt,i

PFt,i =

 j

 R jg

2−m m β Qt,i υ j Lt,i, j

di5−m



Vt,i = Qt,i Tt

t ∈ T , i < imax

t ∈ T , i < imax , j ∈ J, k ∈ Ki

(4) (5)

where PHt, i, k is the provided pressure of pump k at station i in time window t; R j , υ j is the density and viscosity of the refined product corresponding to batch j; g is the gravity acceleration; Qt, i is the flowrate of pipeline segment (i, i + 1) in time window t; ai, k , bi, k is the pump head related coefficients; PFt, i is the hydraulic loss of pipeline segment (i, i + 1) at the start time of time window t; Lt, i, j is the length of batch j in pipeline segment (i, i + 1) at the start time of time window t; di is the inner diameter of pipeline segment (i, i + 1); λi is the tilt coefficient of pipeline segment (i, i + 1); Et, i.k is the output work of pump k at station i in time window t; Tt is the accumulated running time of pump k at station i; Vt, i is the flowed volume of pipeline segment (i, i + 1) in time window t. 2.2. Model requirements

(2)

To solve the detailed scheduling problem of multiproduct pipelines with the hydraulic constraints rigorously considered, an MILP model is set up. Given:

t ∈ T , i < imax , j ∈ J (3)

• Pipeline and station data: length and diameter of each pipeline segment, flowrate and pressure limits, horizontal coordinates and elevations of stations, pipeline initial conditions.

t ∈ T , i < imax , j ∈ J, k ∈ Ki

 + λi Lt,i, j



3−m Et,i.k = R j g ai,k Qt,i − bi,k Qt,i Tt

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

• Pump data: pump configuration of pump stations and pump performance parameters (i.e., pumping head and pump flowrate limits) • Injection and demand plan: injection batch sequence and volume, demand volume of each delivery station for various refined products. • Refined product data: physical properties of refined products (i.e., density and viscosity). • Cost data: electricity price, unit costs for pump start/stop. Determine: • Pipeline delivery schedule: injection oil type and flowrate in each time window, download flowrate of each delivery station, pipeline segment flowrates, and locations of different batch interfaces. • Pump schedule: exact time of each pump start/stop, frequency of pump start/stop, pressure distribution along the pipeline, and energy consumption of each pump. Objective: On the basis of meeting various operational and technical constraints, this model aims to solve a multiproduct pipeline detailed schedule with the minimal pumping costs, including pump running costs and start/stop costs. To establish and solve the model effectively, this paper makes the following reasonable simplifications and assumptions: • The mixed oil section between two adjacent batches of refined products is considered as an interface. • As the tankage profile will be considered in advance when the pipeline operators determine the overall goal of the short-term plan, inventory constraints are not considered in the model. • To highlight the research focus of this paper, the injection batch sequence and volume, which belong to the aggregate scheduling, are considered to be given beforehand (Cafaro et al., 2015). • The injection/delivery flowrate does not change in a time window, and an injection/delivery station can only inject/deliver one batch in a time window. • The fluid is incompressible, and the volume changes in the fluid when it is flowing through pumps and valves and moving in the pipeline can be ignored. The fluid’s physical properties are assumed to be constant. • The calculation of the pressure distribution along the pipeline is considered as a steady-state process, and the influence of fluid transient in the pipeline is not taken into account. • In the high efficiency working range of the pumps, the change in the pumping head with the flowrate is relatively small. Therefore, the pumping head is assumed to be a constant value as the operators on site did, and the nonlinear items in Eqs. (2) and (4) can be ignored. • The operation of the pump start/stop can be completed instantly and does not take up time. • Since the electricity price varies frequently in a day, which will increase the total number of time windows, the change in electricity prices over time is not considered to highlight the research focus of this paper. 3. Mathematical model In the previously proposed method for the rigorous description of the hydraulic loss changes in multiproduct pipelines (Zhou et al., 2019c), it was concluded that if the pipeline segment flowrate is unchanged and there is no batch interface flow in/out of the pipeline segment in a time interval, the hydraulic loss will change linearly with time. Based on that conclusion, this paper establishes a continuous-time model for the detailed scheduling of multiproduct pipelines with the hydraulic constraints rigorously considered.

7

In this model, all moments of pipeline flowrate changes and batch interfaces reaching stations are taken as time nodes. In this way, the above two requirements can be satisfied, and the pipeline hydraulic loss changes linearly in a time window. Furthermore, the moments in which the pump schedule changes in the above time windows are also taken as the time nodes, so that the pipeline flowrate, pump schedule, and refined products allocation are all the same in a time window, and the hydraulic constraints can be exactly satisfied in this time window. Another difference between the established model and the previous related models is that this paper tracks the batch interface locations accurately and calculates the pipeline hydraulic loss considering the allocations and the differences in the physical properties of different refined products strictly (i.e. Section 3.8), and this is the basis of solving a precise pipeline detailed schedule. 3.1. Objective function The pumping costs F include the pump running cost f1 and pump start/stop cost f2 .

min F = f1 + f2

(6)

The pump running cost is determined by the electricity price and the output work of the pump. The pump running costs can be calculated as follows.

f1 = CE



WPt,i.k

(7)

t∈T i∈I k∈Ki

The pump start/stop cost is equal to the product of the unit cost and the total times of pump start/stop during the scheduling horizon.

f2 = CU



BCPt,i,k

(8)

t∈T i∈I k∈Ki

3.2. Time node constraints The start time of time window t + 1 equals the start time of time window t plus its time length.

TLt+1 = TAt + TLt

t < tmax

(9)

The start time of the first time window equals zero.

TL1 = 0

(10)

3.3. Delivery constraints The actual delivered volume and the demand volume of oil o at delivery station i should be within a certain deviation range.



γ j,oVDOt,i, j ≥ (1 − μ )VDEi,o i ∈ I, o ∈ O

(11)

γ j,oVDOt,i, j ≤ (1 + μ )VDEi,o i ∈ I, o ∈ O

(12)

t∈T j∈J

 t∈T j∈J

The following three constraints ensure that if delivery station i downloads batch j in time window t (i.e., BDOt,i, j = 1), the downloaded volume should be within a certain range; otherwise, it is equal to 0.

BDOt,i, j vD min i ≤ VDOt,i, j ≤ BDOt,i, j vD max i 

j∈J VDOt,i, j

qD max i

 ≤ TAt ≤

j∈J VDOt,i, j

qD min i

t ∈ T , i ∈ I, j ∈ J

 +τ 1−



(13)

BDOt,i, j

t ∈ T, i ∈ I

j∈J

(14) Delivery station i can download batch j only when batch j is flowing through it.

BDOt,i, j ≤ BFTt,i, j

t ∈ T , i ∈ I, j ∈ J

(15)

8

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

3.4. Batch moving constraints BPAt, i, j is valued as 1 if the upper coordinate of batch j has exceeded the volume coordinate of station i at the start time of time window t, otherwise BPAt,i, j = 0. The following two constraints realize this definition. XI refers to the volume coordinate of station I (the terminal station), and COH1, j is the initial upper coordinate of batch j, which is known in advance.

COHt, j ≤ Xi + BPAt,i, j (XI − Xi ) t ∈ T , i ∈ I, j ∈ J





COHt, j ≥ Xi + BPAt,i, j − 1 Xi − COH1, j



(16)

If batch j has not been injected at the start time of time window t (i.e., BPAt,1, j = 0), its increased upper coordinate is equal to the injected volume in time window t.

COHt+1, j ≤ COHt, j + VPFt,1 + BPAt,1, j XI

t < tmax , j ∈ J

(28)

COHt+1, j ≥ COHt, j + VPFt,1 − BPAt,1, j XI

t < tmax , j ∈ J

(29)

If batch j has reached the terminal station, its upper coordinate no longer increases.









COHt+1, j ≤ COHt, j + 1 − BPAt,i max, j XI t ∈ T , i ∈ I, j ∈ J

(17)

For the same batch j, the time at which the upper coordinate reaches station i − 1 must be earlier than that of station i. Similarly, for the same station i, the arrival time of the upper coordinate of batch j − 1 must be earlier than that of batch j. If batch j has reached station i in time window t − 1, it must also reach the same station in time window t.

BPAt,i−1, j ≥ BPAt,i, j

t ∈ T , i > 1, j ∈ J

(18)

BPAt,i, j−1 ≥ BPAt,i, j

t ∈ T , i ∈ I, j > 1

(19)

BPAt,i, j ≥ BPAt−1,i, j

t > 1, i ∈ I, j ∈ J

(20)

COHt+1, j ≥ COHt, j + BPAt,i max, j − 1 XI

VPFt,i VPFt,i ≤ TAt ≤ + τ BIDt,i qPmaxi qPmini

BFTt,i, j = BPAt,i, j − BPAt,i, j+1 BFTt,i, j max = BPAt,i, j max

t ∈ T , i ∈ I, j < jmax

t ∈ T, i ∈ I

(22) (23)

There must be a batch flowing through station i in time window t.



BFTt,i, j = 1 t ∈ T , i ∈ I

(24)

j∈J

If the upper coordinate of batch j has exceeded station i at the start time of time window t, and it has not exceeded station i + 1, it is within pipeline segment (i, i + 1).

BFAt,i, j = BPAt,i, j − BPAt,i+1, j

t ∈ T , i < imax , j ∈ J

(25)

If the upper coordinate of batch j is within pipeline segment (i, i + 1) in time window t, its location at the start time of time window t +1 is equal to that at the start time of time window t plus the flowed volume of pipeline segment (i, i + 1) in time window t.





COHt+1, j ≤ COHt, j + VPFt,i + 1 − BFAt,i, j XI

t < tmax , i < imax , j ∈ J (26)





COHt+1, j ≥ COHt, j + VPFt,i + BFAt,i, j − 1 XI

t ∈ T , i < imax

(32)

If pipeline segment (i, i + 1) is idle, pipeline segment (i + 1, i + 2) must also be idle.

1 − ψ j BFAt,i, j ≥ BIDt,i

If the upper coordinate of batch j has exceeded station i at the start time of time window t, and the upper coordinate of batch j + 1 has not exceeded station i, batch j is flowing through station i. As for the last batch, if the upper coordinate has exceeded station i, it is flowing through the station.

(31)

(1 − BIDt,i )vP min i ≤ VPFt,i ≤ (1 − BIDt,i )vP max i t ∈ T , i < imax (33)

COHt+1, j ≤ Xi + 1 −BPAt+1,i, j + BPAt,i, j (XI − Xi ) t < tmax , i ∈ I, j ∈ J (21)

t < tmax , j ∈ J

If pipeline segment (i, i + 1) is idle, the flowed volume in time window t equals 0; otherwise, it must be within the required range.

BIDt,i+1 ≥ BIDt,i



(30)

3.5. Pipeline flowrate constraints

The following constraint ensures that the times the batch interfaces arrive at the stations must be the time nodes. That is, when BPAt+1,i, j − BPAt,i, j = 1, COHt+1, j must be equal to Xi in consideration of Eqs. (17) and (20).



t < tmax , j ∈ J

t ∈ T , i < imax − 1

(34)

The shut-down operation is forbidden when the upper coordinate of batch j exists in the pipeline segment (i, i + 1) (i.e., ψ j BFAt,i, j = 1).

t ∈ T , i < imax , j ∈ J

(35)

The flowed volume of pipeline segment (i, i + 1) equals the flowed volume of pipeline segment (i-1, i) minus the download volume of delivery station i.

VPFt,i = VPFt,i−1 −



VDOt,i, j

t ∈ T , 1 < i < imax , j ∈ J

(36)

j∈J

If pipeline segment (i, i + 1) is not shut down (i.e., BIDt,i = 0), and the upper coordinate of batch j exists in the pipeline segments (i.e., BFAt,i, j = 1), the flowed volume in time window t must satisfy the specific lower limit so that the contaminated products between two adjacent batches will not expand greatly.

TAt ≤

 VPFt,i + τ 1 − BFAt,i, j + BIDt,i t ∈ T , i < imax , j ∈ J qXmini, j

(37)

The pipeline flowed volume equals the pipeline flowrate multiplied by the time length of the time window.

VPFt,i = QPFt,i TAt

t ∈ T , i < imax

(38)

3.6. Injection constraints If BFTt,1, j = 1, that is, batch j is flowing through the initial station in time window t, batch j is being injected, and the injection volume must be within the required range; otherwise, if BFTt,1, j = 0, the injected volume of batch j equals zero.

 VNJt, j VNJt, j ≤ TAt ≤ + τ 1 − BFTt,1, j t ∈ T, j ∈ J qJ max qJ min

(39)

BFTt,1, j vI min ≤ VNJt, j ≤ BFTt,1, j vI max

(40)

t ∈ T, j ∈ J

The injection flowrate in time window t equals the flowrate of pipeline segment (1, 2).

t < tmax , i < imax , j ∈ J (27)

 j∈J

VNJt, j = VPFt,1

t∈T

(41)

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

3.7. Pump station constraints

If the upper coordinate of batch j is within pipeline segment

The following two equations constrain the pump-provided pressure. If batch j is flowing through station i and pump k is turned on in time window t(i.e., BSP t,i,k, j = 1), its provided pressure can be calculated as follows.

Pt,i,k =



t ∈ T , i ∈ I, k ∈ Ki

BSP t,i,k, j RF j gHi,k

9

(i, i + 1) and the lower coordinate of batch j has not exceeded station i, the length of batch j in pipeline segment (i, i + 1) can be calculated as follows.





COHt, j − Xi SAi

1 − BFA t,i, j lZ + LBt,i, j ≥ −BPA t,i, j+1 lZ + t ∈ T , i < imax , j < jmax

(42)

(52)

j∈J

If pump k in station i is not turned on in time window t (i.e., BSP t,i,k, j = 0), its output work must be 0. If it is turned on, the output work is related to the pumping head, pump running time, and flowrate and density of the refined product flowing through the pump. 



j∈J

+(1 − BSP t,i,k, j )vP max i WPt,i.k ≤

COHt, j − Xi SAi

t ∈ T , i < imax , j < jmax

(53)

If the lower coordinate of batch j is within pipeline segment

(i, i + 1) and the upper coordinate of batch j has exceeded station i + 1, the length of batch j in pipeline segment (i, i + 1) can be cal-

WPt,i.k WPt,i.k ≤ VPFt,i ≤ RF j gHi,k RF j gHi,k





BFA t,i, j − 1 lZ + LBt,i, j ≤ BPA t,i, j+1 lZ +

t ∈ T , i ∈ I, j ∈ J, k ∈ Ki

t ∈ T , i ∈ I, k ∈ Ki

BSP t,i,k, j wPmaxi.k

(43)

culated as follows.



(44)







1 − BPA t,i+1, j lZ + LBt,i, j ≥ BFA t,i, j+1 − 1 lZ +

Xi+1 − COHt, j+1 SAi

t ∈ T , i < imax , j < jmax

(54)

j∈J

Only when batch j is flowing through station i, can BSP t, i, k, j equal 1.

BSP t,i,k, j ≤ BFTt,i, j

t ∈ T , i < imax , j ∈ J, k ∈ Ki

(45)

If the on/off statuses of pump k at station i in two adjacent time windows are different, BCPt,i,k = 1; otherwise, it is not bounded.

BCPt,i,k ≥



BSP t,i,k, j −

j∈J

BCPt,i,k ≥





t ∈ T , i ∈ I, k ∈ Ki

BSP t−1,i,k, j

(46)

j∈J

BSP t−1,i,k, j −

j∈J



t ∈ T , i ∈ I, k ∈ Ki

BSP t,i,k, j



VPFt,i qPU min i,k

+τ 1−

 TAt ≥

VPFt,i qPU max i,k

+τ





j∈J



t ∈ T , i < imax , k ∈ Ki (48)

BSP t,i,k, j



BSP t,i,k, j − 1



Xi+1 − COHt, j+1 SAi (55)

If the lower coordinate of batch j has not exceeded station i and the upper coordinate of batch j has exceeded station i + 1, pipeline segment (i, i + 1) is filled with batch j.





1 − BPA t,i+1, j lZ + LBt,i, j ≥ −BPA t,i, j+1 lZ +

Xi+1 − Xi SAi

t ∈ T , i < imax , j < jmax

(56)



j∈J





COHt, j − COHt, j+1 SAi

t ∈ T , i < imax , j < jmax



BFAt,i, j − 1 lZ + LBt,i, j ≤ 1 − BFAt,i, j+1 t ∈ T , i < imax , j < jmax



BPA t,i+1, j − 1 lZ + LBt,i, j ≤ BPA t,i, j+1 lZ +

Xi+1 − Xi SAi

t ∈ T , i < imax , j < jmax

(57)

If neither the upper or the lower coordinates of batch j has exceeded station i or if both have exceeded station i + 1, the length of batch j in pipeline segment (i, i + 1) equals zero.

BPA t,i, j lZ + LBt,i, j ≥ −BPA t,i, j+1 lZ

t ∈ T , i < imax , j < jmax

(58)

−BPA t,i, j lZ + LBt,i, j ≤ BPA t,i, j+1 lZ

t ∈ T , i < imax , j < jmax

(59)



1 − BFA t,i, j lZ + LBt,i, j ≥ BFAt,i, j+1 − 1 lZ +





t ∈ T , i < imax , k ∈ Ki (49)

To calculate the pipeline hydraulic loss accurately, the allocation of the refined products within pipeline segment (i, i + 1) in time window t(i.e., LBt, i, j ) must be determined in advance. Since the time nodes t include all the time nodes of pipeline flowrate changes and batch interfaces reaching the stations, the location relationship between batch j and pipeline segment (i, i + 1) includes the following six cases. If both the upper and the lower coordinates of batch j are within pipeline segment (i, i + 1), the entire batch j is within pipeline segment (i, i + 1).





t ∈ T , i < imax , j < jmax

(47)

3.8. Hydraulic loss calculation constraints





BPA t,i+1, j − 1 lZ + LBt,i, j ≤ 1 − BFA t,i, j+1 lZ +

j∈J

To ensure the high-efficiency running of the pump and avoid possible accidents, the pump flowrate must be within the upper and lower bounds. Otherwise, the pump cannot be turned on in this time window.

TAt ≤



(50)









(60)





BPA t,i+1, j − 1 lZ + LBt,i, j ≤ 1 − BPA t,i+1, j+1 lZ t ∈ T , i < imax , j < jmax

(61)

The Liebenson (Лейбензон) equation is adopted to calculate the pipeline hydraulic loss in this paper. The Liebenson equation, which is derived from the Darcy-Weisbach equation, is widely used in the multiproduct pipelines in China (Chen et al., 2017b; Zhang et al., 2017a). The hydraulic losses of batch j in pipeline segment (i, i + 1) at the start and end times of time window t are expressed as follows.

PFSJt,i, j = RF j g (51)



t ∈ T , i < imax , j < jmax



COHt, j − COHt, j+1 lZ + SAi



1 − BPA t,i+1, j lZ + LBt,i, j ≥ BPA t,i+1, j+1 − 1 lZ

2−m m β QPF υ L t,i F j Bt,i, j

di5−m



+ λi LBt,i, j

t ∈ T , i < imax , j ∈ J (62)

10

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543



β

PFEJt,i, j = RF j g

2−m QPF t,i

υ

m L F j Bt+1,i, j 5−m di

way, the nonlinear function can be replaced by a series of fixed values, which is shown as follows.

 + λi LBt+1,i, j

t < tmax , i < imax , j ∈ J

(63)

The total hydraulic losses of pipeline segment (i, i + 1) at the start and end times of time window t can be calculated as follows.

PFSt,i =



PFSJt,i, j

t ∈ T , i < imax

(64)



t ∈ T , i < imax

PFEJt,i, j

(65)

j∈J

3.9. Pressure constraints The pressure provided by pump station i can be expressed as follows.

PHt,i =



Xmin 1 +Xmax 1 

Xmin 2 +2Xmax 2 

f f

2

.. .

 f Xmin a +2 Xmax a ⎪ ⎪ ⎪ .. ⎪ ⎪ ⎪ . ⎪   ⎪ ⎪ ⎩ Xmin amax +Xmax amax f

j∈J

PFEt,i =

f (x ) ≈

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨

Pt,i,k

t ∈ T, i ∈ I

(66)

k∈Ki

These two constraints are the energy balance constraints along the pipeline.

POUt,i = PINt,i + PHt,i

t ∈ T , i < imax

PINt,i+1 = POUt,i − PFSt,i

t ∈ T , i < imax

(67) (68)

The above two constraints are employed for the calculation of the inlet/outlet pressure of station i at the start time of the time window. However, with the linear and continuous change of the hydraulic loss of a pipeline segment in a time window, it is also necessary to calculate the inlet/outlet pressures at the end time of the time window (i.e., Eqs. (69) and (70)).

POt,i = PIt,i + PHt,i

t ∈ T , i < imax

PIt,i+1 = POt,i − PFEt,i

t ∈ T , i < imax

(69) (70)

The inlet/outlet pressure of each station along the pipeline at both the start and the end times of a time window should satisfy the upper and lower pressure limits.

pO min i ≤ POUt,i ≤ pO max i pI min i ≤ PINt,i ≤ pI max i pO min i ≤ POt,i ≤ pO max i

t ∈ T, i ∈ I t ∈ T, i ∈ I t ∈ T, i ∈ I

(71) (72)

2



t ∈ T, i ∈ I

(74)

4. Model solving Due to the nonlinear relationship of flowrate and pressure and the multiplication between variables (i.e., Eqs. (38), (62), and (63)), the mathematical model is an MINLP model. The piecewise linearization method (Zhang et al., 2017b) is adopted to deal with the nonlinear items. The MINLP model is converted into the MILP model, which can be solved by a commercial MILP solver. In this method, a continuous variable is divided into several intervals, and the nonlinear curve in each interval is approximately equal to the average value of the upper and lower bounds. In this

(75)

Xmin amax ≤ x ≤ Xmax amax

PFSJt,i, j ≤ RF j g

βυFmj LBt,i, j  qA min a + qA max a 2−m di5−m

+(1 − BAQ t,i,a ) pFmaxi, j

PFSJt,i, j ≥ RF j g

di5−m

 PFEJt,i, j ≤ RF j g

PFEJt,i, j ≥ RF j g



+ λi LBt,i, j

t ∈ T , i < imax , j ∈ J, a ∈ A

2

di5−m

+ λi LBt,i, j

t ∈ T , i < imax , j ∈ J, a ∈ A

2

+(BAQ t,i,a − 1 ) pFmaxi, j

2

(77)

 + λi LBt+1,i, j

t ∈ T , i < imax , j ∈ J, a ∈ A

βυFmj LBt+1,i, j  qA min a + qA max a 2−m di5−m

(76)



βυFmj LBt+1,i, j  qA min a + qA max a 2−m

+(1 − BAQ t,i,a ) pFmaxi, j



2

βυFmj LBt,i, j  qA min a + qA max a 2−m

+(BAQ t,i,a − 1 ) pFmaxi, j pI min i ≤ PIt,i ≤ pI max i

Xmin 2 ≤ x < Xmax 2 .. . Xmin a ≤ x < Xmax a .. .

Taking a nonlinear function as an example, Fig. 3 is presented for the illustration of this method. The number of divided intervals can be adjusted based on the model scale and the non-linearity degree of the equations. The denser the intervals are, the higher the accuracy of fitting the nonlinear function is, but the model scale will also increase sharply at the same time. On the contrary, the fewer the number of intervals is, the lower the accuracy is, but the model solving efficiency will be enhanced. In addition, there is no need to divide the intervals evenly, and the encryption intervals can be used where the extent of the nonlinearity is strong. Based on the above linearization method, the nonlinear equations are replaced by the following equations to convert the original MINLP model to the MILP model, and the set of flowrate intervals a is added. In more detail, Eqs. (62) and (63) are replaced by Eqs. (76)–(79), and Eq. (38) is replaced by Eqs. (80)–(82). In fact, the adoption of the piecewise linearization in this paper grades the pipeline flowrate and does not cause any deviation between the original and the converted equations, because if a flowrate interval a is selected, the pipeline flowrate QPFt, i is replaced by ( qA min a +2qA max a ) to plug into Eqs. (38), (62), and (63). This means that a limitation that QPFt, i can only equal a set of average values is made to deal with the nonlinear items in the model. This limitation is consistent with the onsite operation process of the pipeline, because the pipeline flowrate and download flowrate usually tend to be round numbers for the convenience of operation and measurement. If the flowrate of pipeline segment (i, i + 1) is within interval a in time window t, the hydraulic loss of pipeline segment (i, i + 1) at the start time and end time of time window t is calculated as follows.

 (73)

Xmin 1 ≤ x < Xmax 1

(78)

 + λi LBt+1,i, j

t ∈ T , i < imax , j ∈ J, a ∈ A

(79)

If the flowrate of pipeline segment (i, i + 1) is within interval a in time window t, the flowed volume in pipeline segment (i, i + 1)

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

11

Fig. 3. Schematic diagram of piecewise linear approximation.

Fig. 4. Schematic diagram for the studied multiproduct pipeline.

is calculated as follows.

VPFt,i ≤

q

A min a + qA max a



2

TAt + (1 − BAQ t,i,a )vP max i

t ∈ T , i < imax (80)

VPFt,i ≥

q

A min a + qA max a



2

TAt + (BAQ t,i,a − 1 )vP max i

t ∈ T , i < imax (81)



are shown in Table A-2. The pipeline segment flowrate and injection/download flowrate limits are shown in Tables A-3 and A-4, respectively. The pressure limits of each station are given in Table A5. The pipeline transports 92# gasoline, 95# gasoline and 0# diesel in sequence. The physical properties of each kind of refined product are shown in Table A-6. The electricity price at each station is set as 0.8 CNY/(kW•h) according to the industrial electricity price standards in China, and the flowrate interval is set as 50 m3 /h according to the pipeline’s actual working condition. Based on the operational experience on site, the unit costs for a pump start/stop is set as 500 CNY.

There must be a flowrate interval selected.

BAQ t,i,a = 1 t ∈ T , i < imax

(82)

a∈A

5. Results and discussion 5.1. Pipeline basic data The application of the established model was performed on a real-world multiproduct pipeline in Zhejiang, China. The total length of the pipeline is 375.5 km. There are 5 stations along the pipeline, including an injection station (IS) which is the initial station, and four delivery stations (DSs). Fig. 4 illustrates the pipeline system. Stations IS and DS1 are also pump stations (P1, P2). The basic parameters of the pipeline are shown in Table A-1. The pump configuration of each pump station and the related parameters

5.2. Case studies Two oil delivery tasks with different injection and demand plans in June and August 2018 were taken as two cases to verify the model’s practicability and accuracy for the detailed scheduling of multiproduct pipelines. A comparison was also conducted between the established model and Liao’s model (Liao et al., 2018c) to illustrate the superiority of this model. Liao’s model (hereafter “the compared model”) was based on the discrete-time representation. Even though the hydraulic constraints were directly considered, the model uses the average physical properties when calculating the hydraulic loss if there is more than one kind of refined product in a pipeline segment simultaneously. Furthermore, the model also requires that there should be no more than one batch interface in a pipeline segment during a time interval, which is inconsistent with the actual operation process of the pipeline. Based

12

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Table 1 Calculation results of the established model and the compared model. Case

Model

Cont. var.

Disc. var.

# of eq. con.

# of ineq. con.

The cardinality of set t

Pump running cost (CNY)

Pump start-stop cost (CNY)

1

Established model Compared model Established model Compared model

4956 2684 4956 2684

6540 32047 6540 32047

2081 6320 2081 6320

24433 77072 24433 77072

30 57 30 57

214920 258610 242400 273990

15000 10000 21000 7000

2

on the Intel Xeon E5-2630v3 (2.4 GHz) dual thread parallel processor, the MILP model was programmed with MATLAB R2015a and solved by GUROBI 8.1.0. The calculation results for the two cases are listed in Table 1. As shown in Table 1, for the established model, the model scales of Case 1 and Case 2 are the same. This is because the only difference between the two cases is the aggregate injection and demand plan. The total demands of all kinds of refined products in Case 2 are a little more than those in Case 1. In addition, the demands of the delivery stations farther from the initial station are also greater in Case 2. Therefore, more energy is needed to transport the oil, and the pumping costs in Case 2 are a bit higher. The solving time of the established model is concentrated in a thousand seconds or so and varies according to different oil supplies and demands. The scale of the compared model is much bigger than that of the established model due to the discrete time representation, and thus it takes much more time to be solved. In the compared model, the cardinality of set t is 57, while that of the established model is first set as 30. After the model is solved, it is shown that only 20 time windows are used in both cases. The other ten are invalid. Therefore, we find that more than half of the time slots can be reduced when moving from the discretetime representation to the continuous-time representation. At the same time, since both the cardinality of set t and the time interval of a time window are predetermined before model solving in the discrete-time representation, the entire scheduling horizon is also fixed, which is unfavorable for finding the optimal solution. The pump running cost of the compared model is a bit higher than that of the established model. However, the discretization of the scheduling horizon and the deviation in the hydraulic loss calculation in the compared model will lead to an inaccuracy in the pumping cost. Even though the pump start-stop cost in the compared model is lower than that in the established model, it is also unreliable due to the discrete-time representation. That is, if a pump is turned on at the start time of a discrete time window, it must continue until the start time of the next time window, which is unreasonable. More explanations are given in Section 5.3.3. 5.2.1. Case 1 The sequence and injection volume of the batches in Case 1 are shown in Table B-1, and the demand volume of each delivery station are shown in Table B-2. Based on the aggregate injection and demand plan and pipeline basic data, the optimal pipeline detailed schedule, including the delivery schedule (Hossein Mostafaei et al., 2015) (shown in Fig. 5) and pump schedule (shown in Fig. 6) can be solved. The actual downloads of each delivery station are shown in Table B-3, and the pipeline flowrate and download flowrate are shown in Figs. B-1 and B-2, respectively. The inlet/outlet pressures of each station are displayed in Fig. B-3. As we can see from Table B-3, the actual downloads of each delivery station can generally satisfy the demands. All the deviations are within ±5%, and the maximum is 998t (0# diesel for DS4). In this way, the delivery schedule solved by the established model can complete the oil transport task well and ensure the normal and continuous supply of the refined products.

Model solving time (s) 825.1 4013.4 971.4 3879.3

From Figs. B-1 and B-2, the pipeline flowrates and the download flowrates of the delivery stations are all within the upper and lower limits. Furthermore, all flowrates are multiples of 50, which is consistent with the flowrate interval division in Section 5.1. The variation in the flowrate lower limits in Fig. B-1 is owing to the existence of the batch interface. That is, if the batch interface exists in a pipeline segment, the flowrate lower limit must be higher than the normal one to avoid the formation of a large amount of contaminated oil. The partial disappearance of the lower limits in Fig. B-2 is because the delivery station does not download refined products in the corresponding time periods. As mentioned in Section 3, some of the time nodes (i.e., 28.15 h, 95.76 h, etc.) solved by the established model in Fig. 6 belong to the key nodes in the delivery schedule, namely, the time of pipeline flowrate changing or batch interfaces arriving at stations. This kind of time nodes can also be displayed in the batch migration chart. Other time nodes (i.e., 123.26 h, 205.22 h, etc.) are for the accurate control of pumps, which means that the pump scheme must change at that moment, otherwise the pressure limits will not be satisfied, and underpressure/overpressure accidents may be caused. As shown in Fig. 6, the minimum cumulative start/stop time in the pump schedule is 5.8 h (205.22h–211.02 h), and the continuity can meet the onsite requirements for the continuous and safe operation of the pumps. The pump schedule continuity can also be enhanced through adding pump start/stop time limit constraints, which are not considered in this paper to highlight the research focus of this paper and solve the established model efficiently. Fig. B-3 shows that the inlet/outlet pressure of all stations along the pipeline can meet the upper and lower limits. The inlet pressure lower limits of DS2 and DS3 decrease in certain time periods. This is because some pipeline segments shut down, and thus, the inlet pressure does not need to be constrained. Additionally, 0.5 MPa is the minimum allowable suction pressure of the pumps in P2 (DS1) to avoid cavitation. If no pump in P2 is turned on, this limit is not necessary. The disappearance of the blue lines (inlet pressure) indicate that the inlet pressure equals the outlet pressure. Except for the initial station, the inlet/outlet pressure of other four stations change linearly in some time periods, instead of being constant in a time window as in the previous studies. This is the result of considering the changes in pipeline hydraulic loss under time continuity. For instance, the inlet pressure of DS4 decreases from 0.47 MPa at 201.72 h to 0.1 MPa at 205.22 h, and if the pump provided pressure remains unchanged, the pressure lower limits will not be satisfied after this moment. To prevent this from happening, P2-3 turned on at 205.22 h. This is the reason why we define this kind of time nodes in the established model for the accurate control of the pumps. 5.2.2. Case 2 The injection and demand volumes of another oil delivery task are shown in Tables B-4 and B-5, respectively. According to these basic data, the pipeline delivery schedule and actual downloads solved by the established model are displayed in Fig. 7 and Table B-6.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. 5. Pipeline schedule for Case 1 solved by the established model.

Fig. 6. Pump schedule for Case 1 solved by the established model.

13

14

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. 7. Pipeline schedule for Case 2 solved by the established model.

As shown in Fig. 7, the pipeline is filled with 0# diesel at the end of the scheduling horizon, while the batch interface of 92# gasoline and 0# diesel locates at DS3 in Case 1. This is determined by the upstream supply and the downstream demand. If the supply of gasoline is more than the demand, some gasoline will inevitably be left in the pipeline. In this case, the batch interfaces should locate at stations to avoid the increase and facilitate the processing of the contaminated oil. Some related constraints have been added to during the solving of the established model. As we can see from Table B-6, the actual downloads are nearly equal to the demands. The maximum deviation is 600 t, 5% of the demand of DS2 for 0# diesel, and all other deviations are within ±5%. The pipeline flowrate and download flowrate solved by the established model are shown in Figs. B-4 and B-5, respectively. As we can see, all the flowrates fall in the required range. The pump schedule and inlet/outlet pressure are displayed in Figs. 8 and B-6. The case of avoiding underpressure/overpressure accidents has been illustrated in Case 1, and an example of the energy saving is

explained in this case. As shown in Fig. B-6, the inlet (outlet) pressure of DS2 increases to 1.91 MPa at 17.39 h, and if P1-3 is turned off at this time, the pressure will drop to 0.6 MPa, which equals the lower limit. Therefore, to save energy as much as possible, P13 is turned off from 17.39 h to 35.34 h. 5.2.3. Comparison According to the input parameters in the above two cases, the compared model was adopted to calculate the solutions for comparison. The discrete time interval is set as 5 h. The solving results for Case 1 are shown in Table C-1 (actual downloads), Fig. C1 (delivery schedule), Fig. C-3 (pipeline flowrate), Fig. C-4 (download flowrate), Fig. C-5 (pump schedule), and Fig. C-6 (inlet/outlet pressure), and those for Case 2 are shown in Table C-2, Fig. C-2, and Figs. C-7–C-10. Since there are so many (i.e., 57) discrete time intervals in the pipeline schedule solved by the compared model, they are not suitable to be expressed in the form of Case 1 and 2 (the figure will be too long), but in the form of a batch migra-

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

15

Fig. 8. Pump schedule for Case 2 solved by the established model.

tion chart, which is a type of Gantt chart adopted in many previous related studies and suitable for the expression of the pipeline schedule with many time intervals (Liang et al., 2012b; Liao et al., 2018b) . The vertical and horizontal axes in the figure represent the pipeline running time and volume coordinate, respectively. The colored block on the vertical axis represents the allocation of the corresponding refined products in the pipeline at the initial time. The horizontal colored blocks indicate the injection/download operations of the corresponding refined products (yellow for 95# gasoline, cyan for 92# gasoline, and blue for 0# diesel) at the injection/delivery station. The length of the horizontal blocks represents the operation duration, while the width represents the operation flowrate (the larger the flowrate, the wider the width). The black lines in the figure show the locations of the batch interfaces over time. Comparisons were conducted on the solving efficiency and the economy, applicability, and accuracy of the solved pipeline detailed schedule. The compared model is based on the discrete-time representation and the model scale is quite large. Therefore, for both oil delivery tasks, we choose to relax the deviation range parameter until the solution can be obtained within a reasonable time. The maximum deviation solved by the compared model is 30%, and the others are all within this range. Some of deviations are relatively small. For example, in Case 2, the difference between the demand and the actual download volume of 95# gasoline at DS3 is 63 t, 1.3% of the total demand. The number of discrete time windows in the compared model is set as 57, and the entire scheduling horizon of both pipeline delivery schedules is 285 h. This number has to be estimated according to the pipeline’s actual working conditions, and repetition tests are also needed before the model can be solved, which is another disadvantage of the compared model. As shown in the batch migration charts, the download operation of the delivery stations solved by the compared model usually starts at the time of the batch interface arriving, ends at the time of the batch interface leaving, and lasts for a long time. Meanwhile, the pipeline is always active, and no segments shut down during the scheduling horizon. In the delivery schedule solved by the established model (Figs. 5 and 7), the oil delivery tasks of the stations are completed centrally within one or several time periods, which could make the entire schedule more flexible. This is also beneficial for the management of stations and oil depots and can avoid the processing of the contaminated oil for the intermediate stations. Furthermore, as shown in Fig. C-1, the batch interface of 92# gasoline and 0# diesel stops between DS2 and DS3 at the end time, which will lead to the diffusion of contaminated oil.

As shown in Figs. C-3 and C-7, the pipeline flowrate lower limits when the batch interfaces exist cannot be considered in the compared model. Due to the close coupling relationship between the flowrate and pressure, the rigorous consideration of pressure constraints and pipeline flowrate lower limits for batch interfaces will decrease the feasible flowrate range. Therefore, based on the discretization of the scheduling horizon, a solution may not be available. Meanwhile, as shown in Figs. C-4 and C-8, the download flowrate solved by the compared model varies frequently during the scheduling horizon, which increase the difficulty for the download operation and is harmful to the stable and safe running of the pipeline. As shown in Figs. C-5 and C-9, the continuity of the pump schedule solved by the compared model is better than that of the established model. However, the compared model does not consider the accurate allocation of each kind of refined product in a pipeline segment; it simply uses the average physical properties of two kinds of refined products to calculate the hydraulic loss if a batch interface exists in a pipeline segment. In this way, there are only two cases (i.e., the average of 0# diesel and 92# gasoline, and the average of 92# gasoline and 0# diesel) for the hydraulic loss calculation. Due to the small range of variation in the pipeline flowrate, the hydraulic loss will not change widely, which eventually leads to a stable pump schedule. However, the method of hydraulic loss calculation in the compared model is quite rough. For example, as we can see from Table C-3, the difference between the actual hydraulic loss and that calculated by the above method is 0.5 MPa. Therefore, it will result in great differences with the pipeline’s actual working conditions, and thus, the detailed schedule solved by the compared model may not be able to guide the safe and economy operation of the pipeline. The pressure lower limits in Figs. C-6 and C-10 remains unchanged because the pipeline is always active and there is always a pump turned on in P2. The inlet/outlet pressure solved by the compared model remains unchanged in a time window. This is because the hydraulic loss changes between two adjacent discrete time nodes are not considered, and the pressure at the start time of a time window is regarded as that over the time window, which is obviously inaccurate. In addition, according to the rigorous calculation of the pipeline hydraulic loss in the established model, the solved inlet/outlet pressure can show the pressure distribution along the pipeline precisely. Furthermore, neither of the cases for energy saving and avoidance of underpressure/overpressure accidents as mentioned above can be considered in the compared model.

16

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. 9. Pipeline schedule for the testing case solved by the original model.

Concluding from the above comparisons, the compared model mainly has two disadvantages for the detailed scheduling issue of multiproduct pipelines considering the hydraulic constraints: (1) the model is based on the discrete-time representation, which will enlarge the model scale, greatly lower the solving efficiency, and reduce the applicability of the solved detailed schedule as mentioned above; (2) the hydraulic loss calculation in the model has a relatively large deviation, which will result in the solved pump schedule being unreliable to achieve the safe operation of the pipeline and pumps. Meanwhile, the simplification that no more than one batch interface is allowed in a pipeline segment during a time interval also narrows the scope of finding the optimal pipeline detailed schedule. The established model can overcome the above disadvantages and provide the dispatchers with a practical pipeline detailed schedule that can achieve the reliable completion of oil transport tasks, the safe and efficient operation of stations and pumps, and the maximum energy saving. 5.3. Piecewise linearization influence test For the efficient solving of the established MINLP model, the piecewise linearization method is adopted and limits the pipeline flowrate to a series of values. Even though the adoption of the method will not lead to errors and fits with the pipeline’s actual working conditions (the flowrates are usually multiples of 50 or 100 in the actual field for easy metering and operation), it will lead to a solution that may not be the globally optimal one of the original model. Therefore, the BARON solver of GAMS 25.1.3 was employed to solve a small-scale case directly and test how far the solution provided by the piecewise linearization method is with respect to the global optimum. The solving parameters are set according to (Tawarmalani and Sahinidis, 2002). There are only three batches to be injected, and the volume is relatively small compared with the normal cases. The batch sequence and injection volume as well as the demand volume of each delivery station are shown in Tables D-1 and D-2. The solving results of the original

model and the linearized model are shown in Table 2. The solved pipeline schedules are shown in Figs. 9 and 10. The actual downloads of each delivery station are shown in Tables D-3 and D-4. The pipeline flowrates are shown in Tables D-5 and D-6. As shown in the above results, the directly solved solution can satisfy the demands of the delivery station exactly, and there is almost no deviation. However, the solving time is much longer than that of the MILP model, even though the model scale is a little smaller. If the piecewise linearization method is adopted, download deviations will inevitably be caused, but all are within ±5%, which can be accepted by the actual field operators. Since more refined products are downloaded at the later stations (0# diesel at DS3, 0# diesel at DS4, and 92# gasoline at DS4), more energy is needed for fluid transportation, and the pump related costs solved by the MILP model are a bit higher than those solved directly. In summary, the adoption of the piecewise linearization method can reduce the solving time significantly, which will greatly improve the efficiency of formulating the onsite pipeline schedule. The testing case’s scale is relatively small compared with the normal ones. For the delivery tasks onsite, the model scale and solving time will increase exponentially, and it is possible that the solution will not be obtained within a tolerable time or that a solution will be obtained that is worse than that from the piecewise linearization (flowrate discretization). Even though small download deviations and the resulting slight increase of pump costs are caused, this will not influence the pipeline’s safe and efficient operation. The improvement in pipeline operation safety and the reduction of operation costs due to the strict consideration of pipeline hydraulic constraints are more significant than the download deviation caused by the piecewise linearization. In addition, as we can see from Table D-5, in the directly solved solution, most of the pipeline flowrates are equal to multiples of 50, which proves that the piecewise linearization (flowrate discretization) does not have a great influence on the optimal pipeline schedule determination. From another perspective, the adoption of the piecewise linearization method in this paper is equivalent to the discretization of

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

17

Table 2 Calculation results of the original model and the linearized model. Model

Cont. var.

Disc. var.

# of eq. con.

# of ineq. con.

The cardinality of set t

Pump running costs (CNY)

Pump start-stop costs (CNY)

Original model Linearized model

2544 2544

3429 3765

907 1027

11475 13815

12 12

103851 108505

7500 6000

Model solving time (s) 25687.1 384.2

Fig. 10. Pipeline schedule for the testing case solved by the linearized model.

the pipeline flowrate, which is the same as discretizing time using the discrete-time representation. However, the feasible range of the pipeline flowrate is much smaller than the scope of scheduling time, and the number of discrete time intervals is usually 5–6 times as much as the number of flowrate intervals. Since the established model is based on the continuous-time representation, the flexibility of the flowing volume determination in the model can still be ensured. 6. Conclusion In this paper, a continuous-time MILP model for the detailed scheduling of multiproduct pipelines is proposed to consider the hydraulic constraints directly and rigorously, and the optimal pipeline detailed schedule including the delivery schedule and pump schedule can be solved simultaneously. The allocation of the refined products in each time window is tracked rigorously and adopted in the accurate calculation of the hydraulic loss. The calculation of the inlet/outlet pressure of stations at both the start and the end times of a time window ensures that the hydraulic constraints can be satisfied exactly throughout the time window, so that the safe and economic pipeline operation can be guaranteed. Other operational and technical constraints, including delivery constraints, injection constraints, and so on, are also considered in the model. The adoption of the piecewise linearization method converts the original MINLP model into an MILP model from which the global optimal solution is easier to be obtained within a reasonable solving time range. It only limits the pipeline flowrate to a series of values (multiple of 50 m3 /h in case studies) instead of causing any deviation. The application of the established model to a real-world multiproduct pipeline was successfully conducted in two cases. The ap-

plicability and superiority of the established model was verified in detail through comparisons with a previous model that is based on the discrete time representation and makes some simplifications in the calculation of the hydraulic loss. The pipeline detailed schedule solved by established model is more practical to guide the safe and efficient operation of the pipeline, to ensure the reliable supply of refined products to the downstream stations, and to minimize the energy consumption and carbon emissions of the pipeline. Even though the piecewise linearization method does not lead to any deviation, it limits the value of the pipeline flowrate to a series of specific numbers. Therefore, in the future, the authors intend to seek effective methods to deal with the nonlinear relationships between flowrate/pressure and flowrate/volume. Other sub-problems in the scheduling of multiproduct pipelines, such as inventory management and batch sequencing, can also be incorporated. 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. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.compchemeng.2019. 106543.

18

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Appendix A. Pipeline basic parameters Tables A-1–A-6. Table A-1 Basic data of the pipeline. Pipeline segment

Pipeline length(km)

Inner diameter(mm)

Terminal elevation(m)

IS(P1)-DS1 DS1(P2)-DS2 DS2-DS3 DS3-TS

130.4 139.5 66.9 38.7

492.2 492.2 441.2 441.2

7 111 53 51

Table A-2 Basic data of pump stations. Pump station

Pump

Pumping head(m)

Flowrate range(m3 /h)

P1

P1-1 P1-2 P1-3 P2-1 P2-2 P2-3

290 280 180 300 290 180

500–1100 500–1100 500–1100 500–1100 500–1100 500–1100

P2

Table A-3 Download flowrate limits of each station. Station

Injection/download lower limits (m3 /h)

Injection/download upper limits (m3 /h)

IS(P1) DS1(P2) DS2 DS3 DS4(P3)

600 0 0 0 0

1100 950 600 600 600

Table A-4 Flowrate limits of each pipeline segment. Pipeline segment

Normal lower limits (m3 /h)

Lower limits (if the batch interface of 0# D and 92#G exists) (m3 /h)

Lower limits (if the batch interface of 92#G and 95#G exists) (m3 /h)

IS-DS1 DS1-DS2 DS2-DS3 DS3-DS4

600 100 100 100

650 650 550 550

550 550 400 400

1100 1100 1100 1100

Table A-5 Pressure limits of each station. Station

Inlet pressure lower limits (MPa)

Inlet pressure upper limits (MPa)

Outlet pressure upper limits (MPa)

IS(P1) DS1(P2) DS2 DS3 DS4(P3)

/ 0.5 (pump station) 0.6 (high-elevation point) 0.1 0.1

/ 8.5 8.5 8.5 2.2 (the terminal station)

8.5 8.5 8.5 8.5 8.5

Table A-6 Physical properties of each kind of refined product. Product

Density(kg/m3 )

Kinematic viscosity (× 10−6 m2 /s)

0# diesel 92# gasoline 95# gasoline

845 740 750

6.11 2.05 1.07

Appendix B. Calculation results of the established model Tables B-1–B-6, Figs. B-1–B-6.

Upper limits (m3 /h)

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Table B-1 Injection batch sequence and volume for Case 1. Batch

Injection volume(t)

Refined product type

BA1 BA2 BA3 BA4 BA5

0 38500 25000 40000 61000

0# diesel 92# gasoline 95# gasoline 92# gasoline 0# diesel

Table B-2 Refined products demands of each delivery station for Case 1. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

19000 7500 9000 26500

26500 15000 14000 17500

11000 6000 3500 4500

Table B-3 Refined products actual downloads of each delivery station for Case 1. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

19310 7500 8857 25502

27269 14991 14430 17431

11222 5853 3487 4438

Table B-4 Injection batch sequence and volume for Case 2. Batch

Injection volume(t)

Refined product type

BA1 BA2 BA3 BA4 BA5

0 30000 24000 40000 61000

0# diesel 92# gasoline 95# gasoline 92# gasoline 0# diesel

Table B-5 Refined products demands of each delivery station for Case 2. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

18000 12000 6500 25000

21000 22000 13000 13500

7500 6500 4500 5300

Table B-6 Refined products actual downloads of each delivery station for Case 2. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

18247 11400 6300 24603

21083 22419 13596 12901

7777 6792 4275 5156

19

20

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. B-1. Pipeline segment flowrate for Case 1 solved by the established model.

Fig. B-2. Download flowrate of each delivery station for Case 1 solved by the established model.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. B-3. Inlet/outlet pressure at each station for Case 1 solved by the established model.

21

22

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. B-4. Pipeline segment flowrate for Case 2 solved by the established model.

Fig. B-5. Download flowrate of each delivery station for Case 2 solved by the established mode.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. B-6. Inlet/outlet pressure at each station for Case 2 solved by the established model.

Appendix C. Calculation results of the compared model Tables C-1–C-3, Figs. C-1–C-10.

Table C-1 Refined products actual downloads solved by the compared model for Case 1. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

16318 6500 10400 34592

24849 15139 10368 18200

8125 7888 4550 4437

23

24

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Table C-2 Refined products actual downloads solved by the compared model for Case 2. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

18853 15600 7800 27458

19680 17884 12009 16588

7088 5849 4563 6500

Table C-3 Difference between the actual hydraulic loss and that calculated in the compared model. Pipeline segment

Type of batch interface

Location of batch interface (m3 )

Flowrate (m3 /h)

Actual hydraulic loss (MPa)

IS-DS1

92# gasoline (behind) −0# diesel(ahead)

2000

1000

2.9

Fig. C-1. Batch migration chart for Case 1 solved by the compared model.

Fig. C-2. Batch migration chart for Case 2 solved by the compared model.

Hydraulic loss calculated in the compared model (MPa) 3.4

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. C-3. Pipeline segment flowrate for Case 1 solved by the compared model.

Fig. C-4. Download flowrate of each delivery station for Case 1 solved by the compared model.

25

26

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. C-5. Pump schedule for Case 1 solved by the established model.

Fig. C-6. Inlet/outlet pressure at each station for Case 1 solved by the compared model.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. C-7. Pipeline segment flowrate for Case 2 solved by the compared model.

Fig. C-8. Download flowrate of each delivery station for Case 2 solved by the compared model.

27

28

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Fig. C-9. Pump schedule for Case 2 solved by the established model.

Fig. C-10. Inlet/outlet pressure at each station for Case 2 solved by the compared model.

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543

Appendix D. Calculation results of the testing case Tables D-1–D-6.

Table D-1 Injection batch sequence and volume for the testing case. Batch

Injection volume(t)

Refined product type

BA1 BA2 BA3

0 20000 57500

0# diesel 92# gasoline 0# diesel

Table D-2 Refined products demands of each delivery station for the testing case. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

19000 14500 7300 16700

6500 4000 4500 5000

0 0 0 0

Table D-3 Refined products actual downloads of each delivery station solved by the original model. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

19005 14497 7296 16702

6497 4004 4496 5003

0 0 0 0

Table D-4 Refined products actual downloads of each delivery station solved by the linearized model. Delivery station

0# diesel(t)

92# gasoline(t)

95# gasoline(t)

DS1 DS2 DS3 DS4

18197 14282 7605 17409

6543 3885 4393 5179

0 0 0 0

Table D-5 Pipeline flowrate solved by the original model. Pipeline flowrate (m3 /h) Time window 1 2 3 4 5 6 7 8 9 10 11

IS-DS1

DS1-DS2

DS2-DS3

DS3-DS4

820 650 650 1006 650 815 650 723 700 600 600

100 0 650 650 650 815 650 723 700 600 0

100 0 100 315 100 550 550 723 700 600 0

0 0 29 315 100 253 550 550 100 600 0

29

30

X. Zhou, Y. Liang and X. Zhang et al. / Computers and Chemical Engineering 130 (2019) 106543 Table D-6 Pipeline flowrate solved by the linearized model. Pipeline flowrate (m3 /h) Time window 1 2 3 4 5 6 7 8 9

IS-DS1

DS1-DS2

DS2-DS3

DS3-DS4

650 700 1000 750 700 1050 800 600 600

100 650 650 750 700 1050 800 550 600

0 100 350 400 550 550 800 550 600

0 0 0 400 550 550 600 100 600

References Abbasi, E., Garousi, V., 2010. An MILP-based formulation for minimizing pumping energy costs of oil pipelines: beneficial to both the environment and pipeline companies. Energy Syst. 1, 393–416. Barreto, C.V., Pires, L.F.G., Azevedo, L.F.A., 2004. Optimization of pump energy consumption in oil pipelines. In: International Pipeline Conference, pp. 23–27. Cafaro, D.C., Cerda, J., 2008. Dynamic scheduling of multiproduct pipelines with multiple delivery due dates. Comput. Chem. Eng. 32, 728–753. Cafaro, D.C., Cerda, J., 2012. Rigorous scheduling of mesh-structure refined petroleum pipeline networks. Comput. Chem. Eng. 38, 185–203. Cafaro, V.G., Cafaro, D.C., Cerdá, J., 2013. Improving the mathematical formulation for the detailed scheduling of refined products pipelines by accounting for flow rate dependent pumping costs. Iberoam. J. Ind. Eng. 5, 115–128. Cafaro, V.G., Cafaro, D.C., Mendez, C.A., Cerda, J., 2012. Detailed scheduling of single-source pipelines with simultaneous deliveries to multiple offtake stations. Ind. Eng. Chem. Res. 51, 6145–6165. Cafaro, V.G., Cafaro, D.C., Méndez, C.A., Cerdá, J., 2011. Detailed scheduling of operations in single-source refined products pipelines. Ind. Eng. Chem. Res. 50, 6240–6259. Cafaro, V.G., Cafaro, D.C., Méndez, C.A., Cerdá, J., 2015. MINLP model for the detailed scheduling of refined products pipelines with flow rate dependent pumping costs. Comput. Chem. Eng. 72, 210–221. Castro, P.M., 2010. Optimal scheduling of pipeline systems with a resource-task network continuous-time formulation. Ind. Eng. Chem. Res. 49, 11491–11505. Castro, P.M., Mostafaei, H., 2019. Batch-centric scheduling formulation for treelike pipeline systems with forbidden product sequences. Comput. Chem. Eng. 122, 2–18. Chen, H., Wu, C., Zuo, L., Diao, F., Wang, L., Wang, D., Song, B., 2017a. Optimization of detailed schedule for a multiproduct pipeline using a simulated annealing algorithm and heuristic rules. Ind. Eng. Chem. Res. 56, 5092–5106. Chen, H., Zuo, L., Wu, C., Wang, L., Diao, F., Chen, J., Huang, Y., 2017b. Optimizing detailed schedules of a multiproduct pipeline by a monolithic MILP formulation. J. Petrol. Sci. Eng. 159, 149–163. Ferreira, J.A.N.S., Vidal, R.V.V., 1984. Optimization of a pump-pipe system by dynamic programming. Eng. Optimiz. 7, 241–251. Herrán, A., Cruz, J.M.D.L., Andrés, B.D., 2010. A mathematical model for planning transportation of multiple petroleum products in a multi-pipeline system. Comput. Chem. Eng. 34, 401–413. Kirschstein, T., 2018. Planning of multi-product pipelines by economic lot scheduling models. Eur. J. Oper. Res. 264, 327–339. Laínez, J.M., Puigjaner, L., Reklaitis, G.V., 2009. Financial and financial engineering considerations in supply chain and product development pipeline management. Comput. Chem. Eng. 33, 1999–2011. Liang, Y., Li, M., Li, J.F., 2012a. Hydraulic model optimization of a multi-product pipeline. Petrol. Sci. 9, 521–526. Liang, Y., Liu, Z.Z., 2011. Optimal operation of multi-product pipeline network. J. China U. Petrol. Nat. Sci. 35, 115–118. Liang, Y., Ming, L., Ni, Z., 2012b. A study on optimizing delivering scheduling for a multiproduct pipeline. Comput. Chem. Eng. 44, 127–140. Liang, Y., Zhou, X.Y., Yan, X.H., Zhang, H.R., Liu, J., 2017. Optimization of pump start-up schemes for large-scale multiproduct pipelines. J. China U. Petrol. Nat. Sci. 41, 130–138. Liang, Y.T., Gong, J., Kang, Z.L., Yang, F.F., Wang, Y.H., 2004. Optimal operation of multi-product pipeline. J. China U. Petrol. Nat. Sci. 28, 97–101. Liao, Q., Liang, Y., Xu, N., Zhang, H., Wang, J., Zhou, X., 2018c. An MILP approach for detailed scheduling of multi-product pipeline in pressure control mode. Chem. Eng. Res. Des. 136, 620–637. Liao, Q., Zhang, H., Wang, Y., Zhang, W., Liang, Y., 2018b. Heuristic method for detailed scheduling of branched multiproduct pipeline networks. Chem. Eng. Res. Des. 140, 82–101.

Liao, Q., Zhang, H., Xu, N., Liang, Y., Wang, J., 2018a. A MILP model based on flowrate database for detailed scheduling of a multi-product pipeline with multiple pump stations. Comput. Chem. Eng. 117, 63–81. Meira, W.H.T., Magatão, L., Relvas, S., Barbosa-Póvoa, A.P., Neves, F., Arruda, L.V.R., 2018. A matheuristic decomposition approach for the scheduling of a single-source and multiple destinations pipeline system. Eur. J. Oper. Res. 268, 665–687. MirHassani, S.A., Abbasi, M., Moradi, S., 2013. Operational scheduling of refined product pipeline with dual purpose depots. Appl. Math. Model. 37, 5723–5742. MirHassani, S.A., Jahromi, H.F., 2011. Scheduling multi-product tree-structure pipelines. Comput. Chem. Eng. 35, 165–176. Mostafaei, H., Castro, P.M., Ghaffari-Hadigheh, A., 2015. A novel monolithic MILP framework for lot-sizing and scheduling of multiproduct treelike pipeline networks. Ind. Eng. Chem. Res. 54, 9202–9221. Mostafaei, H., Castro, P.M., Ghaffari-Hadigheh, A., 2016. Short-term scheduling of multiple source pipelines with simultaneous injections and deliveries. Comput. Oper. Res. 73, 27–42. Qiu, R., Zhang, H., Gao, X., Zhou, X., Guo, Z., Liao, Q., Liang, Y., 2018. A multi-scenario and multi-objective scheduling optimization model for liquefied light hydrocarbon pipeline system. Chem. Eng. Res. Des. 141, 566–579. Rejowski, R., Pinto, J.M., 2003. Scheduling of a multiproduct pipeline system. Comput. Chem. Eng. 27, 1229–1246. Rejowski, R., Pinto, J.M., 2005. A rigorous MINLP for the simultaneous scheduling and operation of multiproduct pipeline systems. Comput. Aided Chem. Eng. 20, 1063–1068. Rejowski, R., Pinto, J.M., 2008. A novel continuous time representation for the scheduling of pipeline systems with pumping yield rate constraints. Comput. Chem. Eng. 32, 1042–1066. Relvas, S., Barbosa-Povoa, A.P.F.D., Matos, H.A., 2009. Heuristic batch sequencing on a multiproduct oil distribution system. Comput. Chem. Eng. 33, 712–730. Relvas, S., Magatao, S.N.B., Barbosa-Povoa, A.P.F.D., Neves, F., 2013. Integrated scheduling and inventory management of an oil products distribution system. Omega-Int. J. Manage. S. 41, 955–968. Tawarmalani, M., Sahinidis, N.V., 2002. GAMS/BARON: a tutorial and empirical performance analysis. In: Convexification and Global Optimization in Continuous and Mixed-Integer Nonlinear Programming: Theory, Algorithms, Software, and Applications, Boston, MA: Springer US, pp. 313–401. Zhang, H., Liang, Y., Liao, Q., Shen, Y., Yan, X., 2018. A self-learning approach for optimal detailed scheduling of multi-product pipeline. J. Comput. Appl. Math. 327, 41–63. Zhang, H., Liang, Y., Liao, Q., Wu, M., Yan, X., 2017a. A hybrid computational approach for detailed scheduling of products in a pipeline with multiple pump stations. Energy 119, 612–628. Zhang, H., Liang, Y., Xiao, Q., Wu, M., Shao, Q., 2016. Supply-based optimal scheduling of oil product pipelines. Petrol. Sci. 13, 355–367. Zhang, H., Liang, Y., Zhou, X., Yan, X., Qian, C., Liao, Q., 2017b. Sensitivity analysis and optimal operation control for large-scale waterflooding pipeline network of oilfield. J. Petrol. Sci. Eng. 154, 38–48. Zhou, X., Liang, Y., Xin, S., Di, P., Yan, Y., Zhang, H., 2019b. A MINLP model for the optimal waterflooding strategy and operation control of surface waterflooding pipeline network considering reservoir characteristics. Comput. Chem. Eng. 129, 106512. Zhou, X., Zhang, H., Qiu, R., Liang, Y., Wu, G., Xiang, C., Yan, X., 2019c. A hybrid time MILP model for the pump scheduling of multi-product pipelines based on the rigorous description of the pipeline hydraulic loss changes. Comput. Chem. Eng. 121, 174–199. Zhou, X., Zhang, H., Qiu, R., Lv, M., Xiang, C., Long, Y., Liang, Y., 2019a. A two-stage stochastic programming model for the optimal planning of a coal-to-liquids supply chain under demand uncertainty. J. Clean. Prod. 228, 10–28.