Mixed-integer linear programming model for tree-like pipeline scheduling problem with intermediate due dates on demands

Multi-product pipelines are a significant and extensive mean of transporting petroleum based products from refineries to distribution centers. Previous contributions on tree-like pipeline scheduling problem have considered a simple structure with a single refinery connected to a mainline and some secondary lines only emerged from the mainline. In practice, however, a tree-like pipeline may also have several branches on a secondary line resulting in a complex structure, the so called multi-level tree-like pipeline. This paper addresses the short-term scheduling of multi-level tree-like pipelines with multiple refineries through a continuous time mixed-integer linear programming (MILP) model that considers multiple intermediate due dates for product demands. The objective is to satisfy product demands on time at the minimum operational costs, such as pumping, interface and backorder costs. The proposed model performance’s is shown by solving four examples.


Introduction
Supply chain can be described as an integrated process, in which suppliers, manufacturers, distributors and retailers are linked together in order to deliver the goods to the customer properly (Karakitsiou and Migdalas 2008;Zaghian and Mostafaei 2016). Product distribution is a challenging problem in the oil supply chain, in which several transportation modes (e.g., road, vessel and pipeline) can be used to convey oil products from production areas to markets. Pipelines are the most efficient, reliable and safe way of transporting the petroleum products from the refineries to the distribution centers. They are often multi-product systems conveying different types of oil products (e.g., gasoline, diesel, kerosene, liquefied petroleum gas (LPG) and jet fuel). The products are pumped back-to-back from refineries, often without any separation device in between, resulting in a product contamination at the interface material of adjacent products. The size of interface material depends on the product sequences in the pipeline.
The pipeline scheduling problem involves sequencing/sizing product batches and timing product injections/deliveries in order to meet product demands on time at minimum interface and pumping costs. Mixed-integer linear programming (MILP) and mixed-integer nonlinear programming (MINLP) models for pipeline scheduling problems have received growing attentions in the last decades. They are often divided into discrete or continuous-time approaches, with the former partitioning the time horizon into lots of fixed durations (Magatão et al. 2004;Rejowski and Pinto 2003Zyngier and Kelly 2009;Mirhassani and Alizadeh 2008;Herrán et al. 2010) while the latter relaxing such assumption Cerdá 2004, 2008;Relvas et al. 2006;Castro 2010;Mostafaei and Ghaffari-Hadigheh 2014;Ghaffari-Hadigheh and Mostafaei 2015). Mixed integer formulations usually achieve a better pipeline schedule compared to other optimization tools, such as knowledge-based heuristic techniques (Sasikumar et al. 1997), discrete event simulation (Mori et al. 2007;García-Sánchez et al. 2008) and decomposition approaches (Hane and Ratliff 1995;Lopes et al. 2010;Cafaro et al. 2011Cafaro et al. , 2015. Rejowski and Pinto (2003) developed the first discrete time MILP method for the scheduling of a straight pipeline connecting a unique refinery to several distribution caners. The same problem was tackled by Cafaro and Cerdá (2004), who offered a smaller continuous time model meeting the problem targets in a less computational time.
The scheduling of multiple refineries pipelines, also the subject of this work, involving additional operations such as inserting a new product or increasing the size of a batch in transit by mid refineries was also addressed by several authors Cerdá 2009, 2010;Mostafaei et al. 2015aMostafaei et al. , 2016Mostafaei and Castro 2017). In Cafaro and Cerdá (2009) and Mostafaei et al. (2015a), just a single refinery can inject material into the pipeline at any time while in Cafaro and Cerdá (2010), Mostafaei et al. (2016) and Mostafaei and Castro (2017) several refineries can simultaneously inject products into the pipeline during a pumping operation. Mostafaei and Castro (2017) are the first to consider interacting pumping runs for multi-refinery pipeline systems, in which a segment of the pipeline simultaneously receives products from its refinery and upstream segment. All these approaches considered the optimal scheduling of a straight multi-refinery pipeline with a unidirectional flow and with a single pipeline between two nodes (i.e., refinery or depot). Scheduling of complex pipeline structures featuring multiple lines between two nodes and bidirectional flows in pipeline segments can be also found in the literature (Herrán et al. 2010;Lopes et al. 2010;Boschetto et al. 2010;Stebel et al. 2012).
In the oil industry, tree-like pipeline structures are the most used way to connect refineries to depots. Some of them compose of a mainline and several secondary lines only emerging from the mainline (i.e., single level tree-like) while most of them also involve branches splitting from secondary lines, turning the system into a multi-level tree structure pipeline. Mirhassani and Jahromi (2011), Cafaro and Cerdá (2012) and Mostafaei et al. (2015b, c) developed continuous time MILP formulations to address the operational scheduling of a single level tree-like pipeline. In the approach of Mirhassani and Jahromi (2011), at any time just a single batch can be transferred to an active secondary line whereas the approach by Cafaro and Cerdá (2012) overcomes this limitation. Unlike both approaches (Mirhassani and Jahromi 2011;Cafaro and Cerdá 2012), the model used by Mostafaei et al. (2015c) is able to perform simultaneous deliveries to depots and to control flow rate limitations on pipeline segments. Table 1 gives MILP models for the short term scheduling of tree like and multi-refinery pipelines.
In this paper, we introduce a novel MILP formulation based on a continuous representation of both time and volume scales for the short term scheduling of multiple-refiner multi-level tree like pipelines. It overcomes the drawback of previous contributions that only consider the optimal scheduling of single level treelike pipelines featuring a unique refinery. Moreover, contrary to the previous researches on multi-refinery pipelines only assuming a single-period planning Challenges on the scheduling of multi-level tree like pipeline networks with multiple refineries comes from two issues: (a) tracking volume and position of product batches: product batches in the pipeline will no longer be sequenced in the same order that they are pumped. This is because mid-input nodes can inject new material between two products existed already in the pipeline or enlarge the volume of product batches, and (b) sequencing product in branching lines: some product batches in mainline may not be transferred to branching lines, and so some empty batches will move between two reals batches. This issue greatly complicate accounting for interface material and forbidden sequences. Considering both issues for multi-refinery multi-level tree like pipeline is a great challenge, requiring a completely major reformulation of previous approaches. Both issues will be considered in this paper.
The rest of this paper is structured as follows. Section 2 presents a brief description of the problem. Section 3 presents a continuous time mathematical formulation for the operational scheduling of multi-level tree structure pipelines. Section 4 provides four test problems to illustrate the applicability and efficacy of the proposed model, with Examples 1 and 2 featuring multi-level tree like pipeline with two refineries. Example 3 taken from Mirhassani and Jahromi (2011) addresses the scheduling of a single level tree-like pipeline and Example 4 taken from Cafaro and Cerdá (2010) deals with the scheduling of a straight pipeline including two refineries. In the last section, conclusions are provided.

Problem statement
This paper addresses the scheduling of multi-level tree structure pipelines with multiple refineries, depots and dual purpose stations (see Fig. 1). Such a pipeline composes of a mainline and several branching lines emerging at branch points in Depot s2n3 Depot s1n4 Fig. 1 A multi-level tree-like pipeline different levels. Input nodes (i.e., refineries) inject materials into the pipeline while output nodes (i.e., depots) receive oil products. Dual purpose nodes having both input and output facilities, both featuring the same location on the mainline, can inject/withdraw products to/from the pipeline simultaneously. We use the term ''secondary line'' for boughs emerged from the mainline, ''split line'' for branches appeared on the secondary lines, and ''branching line'' for all branches on the system. In Fig. 1, the pipeline system includes two input nodes (i.e., R1 and R2), eight output nodes, one dual purpose node and four branching lines with three secondary lines and one split line emerging from secondary line n2. Pipeline n0 is referred to the mainline. All lines are full with products at any time and so the injected volume into the pipeline system should be equal to the volume withdrawn from it.
Product batches injected into the mainline at input nodes can be diverted to mainline depots or branching lines. Each batch contains a single product (batch I1 conveys product P2 in Fig. 1). In this paper, we use the same method for defining batches as given by Mostafaei et al. (2015c). Batches of each pipeline are categorized in two types: old batches and new batches. Old batches of pipeline n are those exist in that pipeline at the start time of planning horizon (I old n ). The new products lots that will enter into pipeline n during the scheduling horizon are new batches ðI new n Þ. The set of old batches of mainline in Fig. 1 is I old n0 ¼ I1; I2; I3; I4 f g , with the old empty batch I3 has been reserved for a future injection of a new product at dual purpose node Cerdá 2009, 2010).
We assume that the cardinality of mainline's new batches ðjI new n0 jÞ is known beforehand. If two new batches I5 and I6 are injected into the mainline of Fig. 1 during the planning horizon, we will have I new n0 ¼ I5; I6 f g. The set of new batches that can reach secondary line n4 is I new n4 ¼ I2; I3; I4; I5; I6 f g . I ¼ I old n [ I old n is the set of batches to move in pipeline n during scheduling horizon. Sets I r contain batches that can receive material from refinery r, and I s n are those that can be diverted to depot s n . For example, batches I4, I5 and I6 can receive material from refinery R1 during the planning horizon while batch I1 cannot receive material, and so we have I R1 ¼ I4; I5; I6 f g . In Fig. 1, batches I2 to I6 can be discharged into depot s3 n0 during the planning horizon while batch I1 cannot be discharged, and so we have I s3 n0 ¼ I2; I3; I4; I5; I6 f g . The aim is to optimally determine the size and sequence of batches that will be injected into each pipeline during the planning horizon at the minimum interface and pumping costs considering the following conditions: (1) the pump rate at refineries should be kept in the feasible range, (2) product demands requested by local markets of each depot should be satisfied through planning horizon h max partitioned into some periods t 2 T, with a known demand for period t, and (3) the inventory level of depots should be kept in the acceptable range.
The remaining assumptions are as follows: a. All refineries are located on the mainline and can only inject material into the mainline. Note however, in practice a refinery placed on the mainline's branch point can directly inject material into a secondary line. This new problem feature will be considered in a next paper.
b. All pipeline segments stay full of product at any time and operate in a single flow direction. c. Sequence, content and volume of old batches in each pipeline are known. d. Interface volume between two pair of products is a known parameter. Due to a vast interface, some product sequences are forbidden. e. The length of the planning horizon is fixed. f. Product inventories at refinery tanks is known during planning horizon. g. The initial inventory of each product at depot tanks is also known at the start time of planning horizon. h. Some empty old batches can be considered in the mainline to account for new product injections at mid refineries. The location of these empty batches are determined by the operators, but their product content and volume are all model decisions.

Mathematical model
The basis of the following MILP model is the mathematical framework by Mostafaei et al. (2015c), and for this reason, the same nomenclatures will be used in similar sections. Tables 2, 3 and 4 define different indices/sets, parameters, and variables. The objective function is as follows:  Set of oil products indexed by p, p 0 CB p;s;n;t :Back p;s;n;t

Pump sequence constraints
In common with Mostafaei et al. (2016), we use the term ''composite run'' for sequencing pumping operations through the planning horizon. A composite run composes of a set of pumping operations accomplished at refineries within a time interval of the scheduling horizon (each composite run may involve one or more pumping operations, see the detailed description in the result section). Equation (1) states that the start time of composite pumping run k should be after finishing all runs belonging to the previous one k -1. Note that the start time of the first composite run is a known parameter (i.e., ST 1 = stf). The length of all composite runs should be lower than the length of planning horizon h max , as shown in Eq. (2).  Mixed-integer linear programming model for tree-like… 405

Tracing batches inside pipelines
To track batches inside the pipeline, we define continuous variable LPV i,k,n denoting the upper coordinate of batch i from the origin of pipeline n at the end of composite run k. This variable will be equal to the total volume of pipeline n from the origin to the farthest extreme of batch i, as imposed by Eq. (3). Note that LPV i,k,n -SPV i,k,n is the lower coordinate of batch i, and also the upper coordinate of batch i ? 1 (LPV i?1,k,n ), see Eq. (4). Since a batch i 2 I n always moves forward, its upper Back p,s,n,t Unsatisfied demand of product p in depot s n due at period t (m 3 ) IPV i,r,K Size of batch i 2 I r injected into the mainline from the refinery r during composite pumping run k (m 3 ) PPV i,p,r,k Size of batch i 2 I r conveying product p injected from the refinery r to mainline during composite pumping run k (m 3 ). coordinate will distance from the origin over the planning horizon, as shown in Eq. (5).
Each pumping operation at refinery r 2 R involves at most one batch injection, as shown in Eq. (6). Injection of batch i from refinery r during composite run k is possible (k i;r;k ¼ 1) if the upper coordinate of this batch at run k -1 reaches the input facility of the refinery (h r ) and its lower coordinate at run k -1 does not surpass h r . When k i;r;k ¼ 1, refinery r injects a positive volume of batch i into the mainline, as imposed by Eq. (9). By Eq. (10), the activity length of an active refinery r at run k is a function of the injected volume and pump rate belonging to [vr min r ; vr max r ]. The longest pump operation in composite run k defines the value of L k , as stated by Eq. (11). X i2I r In this research, we assume that the number of batches to be conveyed is known before hand, and so all new batches should be injected to the mainline during the planning horizon, as imposed by Eq. (12). Note that each batch conveys only one product and its content is determined by binary variable y i,p . If batch i is conveying product p, the size of continuous variable PPV i,p,r,k will be equal to IPV i,r,K ; otherwise, it is set to zero. Volume p entered from refinery r to the mainline during the horizon length should not surpass reft r,p , the inventory of product p at refinery tank r. X k X r k i;r;k ! 1; 8 i 2 I new n0 ð12Þ X p y i;p ¼ 1; 8 i 2 I ð13Þ X p PPV i;p;r;k ¼ IPV i;r;k ; 8 i 2 I r ; r 2 R; k 2 K ð14Þ X k PPV i;p;r;k K j j:IPV max n0 :y i;p ; 8 i 2 I r ; p 2 P; r 2 R ð15Þ X k2K X i2I r PPV i;p;r;k reft r;p ; 8 r 2 R; p 2 P ð16Þ

Discharge to depots
Equations (17)-(21) handle batch removal constraints at depots, and the similar equations can be found in Mostafaei et al. (2016) and Cafaro and Cerdá (2010). Discharge from batch i to depot s n during run k is possible (x i;s;k;n ¼ 1) if the upper coordinate of the batch at run k reaches the output facility of the depot (s s;n ) and the lower coordinate of the same batch at run k -1 has not surpassed s s;n , as imposed in Eqs. (17) and (18). If x i,s,k,n = 1, depot s n receives a positive volume of batch i, which is bounded by Eqs. (19) and (20). If batch i is conveying product p, the size of batch i containing product p discharged to depot s n computed through continuous variable PDPV i,p,s,k,n will be equal to DPV i,s,k,n ; otherwise, it is set to zero, as imposed by Eqs. (21) and (22). LPV i;k;n ! s s;n :x i;s;k;n ; 8 i 2 I n s ; k 2 K; s 2 S; n 2 N; ð17Þ LPV i;kÀ1;n À SPV i;kÀ1;n s s;n þ PV n À s s;n À Á 1 À x i;s;k;n À Á ; 8 i 2 I n s ; k 2 K; s 2 S n ; n 2 N ð18Þ x i;s;k;n DPV min s;n DPV i;s;k;n x i;s;k;n DPV max s;n ; 8 i 2 I n s ; k 2 K; s 2 S n ; n 2 N ð19Þ X s 0 s DPV i;s 0 ;k;n s s;n À LPV i;kÀ1;n À SPV i;kÀ1;n À Á þ IPV i;k;n þ PV n 1 À x i;s;k;n À Á ; 8 i 2 I n s ; 2 S n ; k 2 K; n 2 N ð20Þ X p PDPV i;p;s;k;n ¼ DPV i;s;k;n ; 8 i 2 I s n ; k 2 K; n 2 N; s 2 S n ð21Þ X k PDPV i;p;s;k;n K j jDPV max s;n :y i;p ; 8 i 2 I s n ; n 2 N; s 2 S n ð22Þ To avoid product contaminations, no material can be transferred from segment (s -1) n0 to s n0 when refinery r 2 R located at the start of segment s n0 injects materials into the mainline, as imposed by Eq. (23). When a refinery r 2 R injects material into the mainline, the volume transferred to downstream mainline's depots and secondary lines should be as large as the injected volume at the refinery r 2 R; as imposed by Eq. (24).

Delivering materials to branch lines
Through pumping run k, transfer from batch i 2 I n0 to secondary line n is possible (u i,k,n = 1) if the upper and lower coordinates of the batch satisfy LPV i;kÀ1;n0 À SPV i;kÀ1;n0 r n and LPV i;k;n0 ! r n , where r n is the volumetric coordinate of branch point n on the mainline, as imposed by Eqs. (25)-(26). When u i,k,n = 1, a portion of batch i in the mainline is delivered in the secondary line n during run k, as imposed by Eq. (27). LPV i;k;n0 ! r n u i;k;n ; 8 i 2 I n ; k 2 K; n 2 SP ð25Þ LPV i;kÀ1;n0 À SPV i;kÀ1;n0 r n þ PV n0 À r n ð Þ1 À u i;k;n À Á ; 8i 2 I n ; k 2 K; n 2 SP IRV min n u i;k;n IRV i;k;n IRV max n u i;k;n ; 8 i 2 I n ; k 2 K; n 2 SP Through pumping run k, transfer from batch i in secondary line n to split line n 0 is possible (ur i;k;n;n 0 ¼ 1) if the upper and lower coordinates of batch i in the secondary line satisfy LPV i;kÀ1;n À SPV i;kÀ1;n d n;n 0 and LPV i;k;n ! d n;n 0 , where d n;n 0 is the volumetric coordinate of branch point n 0 on the secondary line n, as shown in Eqs. (28) and (29). When ur i;k;n;n 0 ¼ 1, a portion of batch i in secondary line n is branched to split line n 0 during run k, as shown in Eq. (30). LPV i;k;n ! d n;n 0 ur ;k;n;n 0 ; 8 i 2 I n ; k 2 K; n 2 SP; n 0 2 SL n ð28Þ LPV i;kÀ1;n À SPV i;kÀ1;n d n;n 0 þ PV n À d n;n 0 À Á 1 À ur i;k;n;n 0 À Á ; 8 i 2 I n ; k 2 K; n 2 SP; n 0 2 SL n ð29Þ ur i;k;nn 0 IRV min n 0 IRV i;k;n 0 ur i;k;n 0 IRV max n 0 ; 8 i 2 I n ; k 2 K; n 2 SP; n 0 2 SL n ð30Þ

Interface volume and forbidden sequence
In multi-product pipelines, there is always a mixing volume (interface) between two successive batches conveying different products. In this paper, it is assumed that the mixing volume between products is independent of pipeline length, diameter and flow rate and is a known parameter MIX p;p 0 ;n : Eqs. (31) and (32) compute the interface volume between product batches inside the mainline and branching lines, respectively. In Eq. (32), binary variable z i,n denotes the existence of batch i in branching line n and its value satisfies Eqs. (33) and (34). Due to excessive interface, some products should not touch each other inside the pipeline.
Equations (35) and (36) avoid forbidden sequences inside the mainline and branching lines, respectively.
INFT i;p;p 0 ;n0 ! MIX p;p 0 ;n y iþ1;p þ y i;p 0 À 1 À Á ; 8 i 2 I n0 ; p; p 0 ð31Þ INFT i;p;p 0 ;n ! MIX p;p 0 ;n y i;p þ y i 0 ;p 0 þ z i;n þ z i 0 ;n þ X iÀ1 j ! i 0 þ1 z j;n À Touch p;p 0 À 2 ! ; 8 i; i 0 2 I n i [ i 0 ð Þ; p; p 0 ; n 6 ¼ n0 ð32Þ z i;n X k u i;k;n K j jz i;n ; 8 i 2 I new n ; n 2 N; n 6 ¼ n0 ð33Þ z i;n 0 X k ur i;k;n;n 0 K j jz i;n 0 ; 8 i 2 I new n ; n 2 SP; n 0 2 SL n ; ð34Þ y i;p þ y iþ1;p 0 Touch p;p 0 þ 1; 8 i 2 I; p; p 0 2 P ð35Þ z i;n þ z i 0 ;n X j2I n iÀ1 j ! i 0 þ1 z j;n À y i;p À y i 0 ;p 0 þ Touch p;p 0 þ 3; 8 i 2 I new n ; i 0 2 I n i [ i 0 ð Þ; p; p 0 2 P; n 6 ¼ n0 Note that the product sequencing problem for multi-product pipeline with minimum interface includes the Hamiltonian path which is a known NP complete problem (Milidiu and Liporace 2003;Mostafaei et al. 2016) and therefore Eqs. (32) and (36) will add additional complexity to the problem if a long term planning horizon is to be considered. To overcome this complexity, future works will involve developing heuristic methods for finding product sequences inside the pipeline.

Size of batch i in pipelines and mass balance
The continuous variable SPV i,k,n is the volume of batch i in pipeline n at the end of run k. This value will be equal to SPV i,k-1,n if there is no flow input or output to/ from batch i during run k. The size of batch i in the mainline can be increased by receiving material from one of refineries and decreased by sending materials to mainline's depots and secondary lines, as imposed by Eq. (37). Equations (38) and (39) compute the size of batch i in secondary and split lines. Note that the size of old batches is known parameters ISPV i,n and we have ISPV i;n ¼ SPV i;0;n : SPV i;k;n0 ¼ SPV i;kÀ1;n0 þ X r IPV i;r;k À X s2S n0 DPV i;s;k;n0 À X n2SP IRV i;k;n ; 8 i 2 I n0 ; k 2 K; k ! 1 ð37Þ SPV i;k;n ¼ SPV i;kÀ1;n þ IRV i;k;n À X s2S n DPV i;s;k;n À X n 0 2SL n IRV i;k;n 0 ; 8 i 2 I n ; k 2 K; n 2 SP ð38Þ SPV i;k;n ¼ SPV i;kÀ1;n þ IRV i;k;n À X s2S n DPV i;s;k;n ; When some materials are injected into pipeline n, the same volume should exit the pipeline. This is because the pipeline is always full with products, as imposed by Eq. (40). The volume balance in split lines, secondary lines and mainlines are controlled by Eqs. (41), (42) and (43), respectively. X i2I n SPV i;k;n ¼ PV n ; 8 n 2 N; k 2 K ð40Þ X i2I n IRV i;k;n ¼ X s2SL n X i2I n DPV i;s;k;n ; 8 k 2 K; n 2 N À n0 [ SP f g ð41Þ

Bound on materials from batch depots and branching lines
During run k, a batch i 2 I n0 may feed several secondary lines and depots on the mainline. In this case, the volume of material from the batch i 2 I n0 diverted to the secondary lines and depots is bounded by Eqs. (44) and (45). Equations (46) and (47) control the volume transferred from a batch i in secondary line n to its depots and its split lines. X s2S n0 d n ! s s;n0 DPV i;s;k;n0 þ X n n 0 2SP IRV i;k;n d n À LPV i;kÀ1;n0 À SPV i;kÀ1;n0 À Á þ X r h r d n IPV i;r;k þ PV n0 1 À u i;k;n À Á ; 8 i 2 I n0 ; k 2 K; n 2 SP ð44Þ X s s 0 2S n0 DPV i;s 0 ;k;n0 þ X n2SP d n s s;n0 IRV i;k;n s s;n0 À LPV i;kÀ1;n0 À SPV i;kÀ1;n0 À Á þ X r h r s s;n0 IPV i;r;K þ PV n0 1 À x i;s;k;n0 À Á ; 8 i 2 I n0 ; k 2 K; s 2 S n0 ð45Þ X s2S n d n;n 0 ! s s;n DPV i;s;k;n þ X n 0 n 00 2SP n IRV i;k;n 00 d n;n 0 À LPV i;kÀ1;n À SPV i;kÀ1;n À Á þ IRV i;k;n þ PV n 1 À ur i;k;n;n 0 À Á ; 8 i 2 I n ; k 2 K; n 2 SP; n 0 2 SL n ð46Þ X s 0 2S n s 0 s DPV i;s;k;n þ X n 0 2SP n d n;n 0 s s;n IRV i;k;n 0 s s;n À LPV i;kÀ1;n À SPV i;kÀ1;n À Á þ IRV i;k;n þ PV n 1 À x i;S;k;n À Á ; 8 i 2 I n ; s 2 S n ; k 2 K; n 2 SP; n 0 2 SL n ; ð47Þ

Transferring materials to local markets
Suppose that the planning horizon h max has been partitioned into several time periods t and E t is the end of period t. Binary variable v t,k is equal to one if run k is completed during period t. Every pumping operation k must be completed at one of the periods, as imposed by Eq. (48), and the start of this run must satisfy Eqs. (49) and (50). X Oil products are sent from depots to local markets after receiving from pipeline. The value sent from depot s n to local markets during run k computed through continuous variable LM s,p,k,n is restricted by Eq. (51), where lmr s,p,n is the maximum dispatch rate of product p from depot s n to local markets. Demand of local markets due at period t should be met on time, and unsatisfied demand will result backorder cost, Eq. (52). To reduce backorder cost due at period t, one pumping operation is at least required to be completed in this period, Eq. (53). LM s;p;k;n lmr s;p;n L k ð51Þ X 1 k 0 k LM s;p;k 0 ;n ! X 1 t 0 t Demand s;n;t 0 ;p v t;k À v t;kþ1 À Á À Back P;S;n;t À Back p;s;n;tÀ1 ; 8 p 2 P; s 2 S n ; n 2 N; k 2 K; t 2 T

Inventory at depots
During the execution of pumping run k, the inventory of product p at depot s n (W p,s,k,n ) will not change if there is no flow input or output to/from this depot.
Equations (54) and (55) compute inventory of product p at a pure output node and a dual purpose node during any pumping run, respectively. Due to operational rules, the value of variable W p,s,k,n should be placed in a feasible range, with W0 p,s,n is the initial volume, i.e., W p,s,0,n = W0 p,s,n . W p;s;k;n ¼ W p;s;kÀ1;n À LM p;s;k;n þ X i2I s n PDPV i;p;s;k;n ; 8 p 2 P; s 2 S n ; n 2 N; k 2 K may be injected into the mainline in consecutive runs, the number of pumping operations is a lower bound on the number of new batches ( I new   n0 K j j). On the other hand, since IPV min is the minimum volume to be injected from an active input node and product demands are deterministic data, we have K j j P p;s;n;t demand p;s;n;t IPV min . Mixed-integer linear programming model for tree-like… 413

Example 1
The pipeline structure for this example is depicted in Fig. 2. The pipeline network has three pipelines, two refineries (R1, R2), five depots (N1-N5) and a single dual purpose node containing refinery R2 and depot N2. The secondary line emerges from the mainline at point 7000 m 3 (r n1 = 7000) and split line leaves the secondary line at point 1000 m 3 (d n1,n2 = 1000). The pump rate can vary between 20 and 200 m 3 /h and the planning horizon h max = 240 h is partitioned into two periods T1 and T2, each containing five days (E T1 = 120 h and E T2 = 240 h). Demands of the periods should be satisfied on time; otherwise, backorder costs will arise. The unit backorder cost for each product is 200 $. Other data of this example can be found in Tables 5, 6 and 7. The maximum and minimum transferred volumes to each pipeline during each run are 4000 and 1000 m 3 , respectively. The same upper and lower bounds hold for the maximum and minimum transferred volumes from each pipeline to depots. The maximum and minimum inventory levels for the refined products in tanks are 10,000 and 2000 m 3 , respectively. Each depot sends products to its local markets at the maximum rate of 200 m 3 /h. Figure 2 shows the pipeline system of Example 1 at time stf = 0 h, which contains three old batches I1, I2 and I3. The aim is to sequence batch input and output operations in order to satisfy product demands at the minimum total cost, including pumping, interface and backorder costs.
The Optimal schedule for Example 1 depicted in Fig. 3 includes eight composite pumping runs with the following input and output operations: 1. From time 0.00 to time 14.00 h (see the second row of Fig. 3), composite pumping run 1 involving two pumping operations at input nodes R1 and R2 is performed. During this run, 2000 m 3 of product P4 from refinery R1 and 2800 m 3 of product from refinery R2 are injected into the mainline to direct 2000 m 3 of product P4 to depot N1, 1800 m 3 of P2 of batch I1 to depot N3 and 1000 m 3 of P1 of batch I2 to secondary line. When 1000 m 3 of batch I2 are injected in the secondary line, the same volume of materials is discharged into depot N4. 2. From time 14.00 to 30.00 h, composite run 2 with only a single pumping operation at refinery R1 is accomplished and inserts 3200 m 3 of product I4. Of the volume, 2000 m 3 of batch I2 are diverted to secondary line and the remaining volume are transferred to depot N3. In this operation, split line receives 1000 m 3 of batch I2 and discharges the same volume of batch I1 to depot N5. 3. From time 30.00 to 49.00 h, the injection of batch I4 is restarted, and 3800 m 3 of P4 are pumped to the mainline from refinery R1. During the execution of this pumping operation, the two following delivery operations to mainline depot are fulfilled: 2000 m 3 of batch I4 to depot N1 and 1800 m 3 from I1 to depot N3. 4. At time 49.00 h, composite run 4 with two pumping operations starts and two batches I5 and I3, both conveying product P4, are injected from refinery R1 and R2, respectively. During this operation, 2500 m 3 of batch I4 and 1000 m 3 of batch I3 are diverted to depots N1 and N5, simultaneously. No branching lines are active through composite run 4. 5. From time 61.50 h to 71.50 h, 2000 m 3 of batch I5 are injected from refinery R1 to divert 1000 m 3 of I4 to depot N2 and 1000 m 3 I3 to secondary line. No volume of products is transferred to split line and all volumes transferred to secondary line are discharged to depot N4. 6. From time 71.50 to 81.50 h, refinery R1 injects 2000 m 3 of batch I6 conveying product P3. During the executing of this pumping operation, split line receives 1000 m 3 of batch I3 originally diverted to secondary line from the mainline. Two depots N3 and N5 receive material from the pipeline during composite run 6. 7. At time 81.50 h, composite pumping run 7 involving two pumping operations at input nodes R1 and R2 starts. Refinery R1 injects 3000 m 3 of batch I6 into the mainline to send the same volume of material from the same batch to depot N1. Refinery R2 pumps 4000 m 3 of product P4 in batch I5. During the composite run, 2000 m 3 of batch I5 are diverted to secondary line and push 1000 m 3 of batch I5 to split line and 1000 m 3 of batch I2 to depot N4. 8. From time 235.50 to 240.00 h (last composite run 8), R8 injects 1000 m 3 of product P4 into batch I3 and the same volume leaves the mainline to the secondary line. During this operation, no volume is injected into a split line.
Like previous continuous time approaches (Cafaro and Cerdá 2009, 2010Mostafaei et al. 2015aMostafaei et al. , c, 2016Mostafaei and Castro 2017), we need an iterative procedure to determine the cardinality of set K at the optimum. To this end, we start with a few elements in set K and keep increasing its cardinality by one until no improvement is observed in the objective function. Table 8 gives the effect of the number of composite runs on solution quality and CPU time. When the number of composite runs is 5, an optimal solution of $359,215 is found in 33.743 s of CPU. However, some demands cannot be satisfied during the scheduling horizon. By increasing the number of runs to 6, demands are fully met and the optimal solution reduces to $55,800. By adding one more elements to the number of runs, the operational cost can be further reduced; however, the CPU time increases from 279.96 to 347.83 s. No improvement can be achieved by setting the number of composite pumping runs to 8, but the computational burden increases from 347.83 to 1179.3 s.

Example 2
This example is a variant of the real-world case study introduced by Mostafaei et al. (2015c) that incorporates a split line in one of the secondary lines and a dual purpose node on the mainline. Figure 4 shows the pipeline structure at time stf = 0 h. The pump rate at refineries should be kept between 300 and 800 m 3 /h. The length of planning horizon is h max = 192 h and is divided into two periods T1 and T2, each containing four days (E T1 = 96 h and E T2 = 192 h. Unit backorder cost for each product is 50 $. Other data can be found in Tables 9, 10 and 11. In Table 10, product sequences indicated with an 'x' are forbidden (Touch P1,P4 = 0). The maximum and minimum transferred volumes to each pipeline during each run are 12,000 and 500 m 3 , respectively. The maximum and minimum inventory level for product in tanks are 30,000 and 5000 m 3 , respectively. Distribution depots can dispatch their products to local markets at maximum rate of 800 m 3 /h. The aim satisfied and the objective function value significantly reduces. By increasing the number of runs from 9 to 10, a best solution worth $1,029,400 is found in 1061.61 s. A further increase in the number of composite runs produce no improvement in the objective function. Figure 5 shows the inventory level of products at different depots of Example 2. It can be observed that the proposed model rigorously controls the inventories and keeps product levels at depot tanks in acceptable limits at any time.

Comparison to previous formulations
Here we consider two examples to compare the performance of the proposed model with regards to previous approaches. The first one (i.e. Example 3) deals with the scheduling of a single level tree like pipeline connecting a single refinery to six depots. The second one (i.e. Example 4) addresses the scheduling of a real life straight pipeline connecting two refineries to three depots. Both examples assume a single period planning horizon.

Example 3
This example was first introduced by Mirhassani and Alizadeh (2008) and later considered by Mirhassani and Jahromi (2011) and Cafaro and Cerdá (2012). The pipeline system consists of a treelike pipeline with two secondary lines, connecting a single refinery to six depots (see Fig. 6). The product injection rate in refinery is restrained to the range of 800-1000 (m 3 /h) and the time horizon has a total length of 90 h. Oil products are dispatched to the local markets at a maximum rate of 400 m 3 / h. Other data can be found in Mirhassani and Jahromi (2011). The best schedule for Example 3 can be found in Fig. 6. Similar to the previous approaches, the proposed model needs three pumping operations to find the optimal solution, but it discovers the optimum in a much lower CPU time of 0.73 s (see Table 14). Such a reduction in CPU is mostly due to the smaller number of binary variables, which have been reduced by a factor of 1.2-13.  P  T1  T2  T1  T2  T1  T2  T1  T2  T1  T2  T1  T2  T1  T2   P1  65  85  70  -70  162  45  45  50  80  68  49  52  60   P2  68  49  58  39  72  48  12  34  30  34  20  -40  50   P3  68  89  31  79  53  162  40  65  24  56  53  11  30  40   P4  12  114  20  -89  178  38  18  40  30  34  -49  59  Table 12 Input to pipelines and depots of Example 2 during planning horizon in 10  This example is derived from a real world case study initially solved by Cafaro and Cerdá (2010). It concerns the distribution of five products from two refineries to three depots over a 10-day time horizon (see Fig. 7). The pipeline features a length of 1000 km and a total volume of 84,600 m 3 . The pump rate in refineries should be in range of 310-580 (m 3 /h). Oil products are dispatched to local markets at a maximum rate of 580 m 3 /h. Other data can be found in Cafaro and Cerdá (2010). The optimal schedule for Example 4 is given in Fig. 7 and was found in 245.21. It includes 7 composite runs, with runs 1, 3, 4 and 6 featuring simultaneous injections at refineries. From Table 15, one can observe that the propose model has 26% less constraints, 12% less continuous variables and 20% less binary variables with respect to the approach of Cafaro and Cerdá (2010). Therefore it is not surprising that the solution CPU time was reduced by a factor of 3.5 i.e., from 865.8 to 245.21 s. There is also a slight reduction in pipeline operation cost which is caused by saving in pumping costs.

Conclusions
In this paper, a continuous time MILP framework for short term scheduling of multi-level tree-like pipeline systems with multiple refineries has been presented. Such a pipeline network consists of a mainline and several branching lines emerging at branch points in different levels. The proposed model has been obtained by generalizing the approach of Mostafaei et al. (2015c) for the short term scheduling of single level tree-like pipeline with unique refinery. It rigorously tracks the position and volume of product batches in pipeline network and satisfies all operational constraints related to product sequencing, mass balances, loading/ unloading operations and capacity constraints at depots. Interface volumes and forbidden sequences between product batches in each pipeline of the network were accounted for by the proposed approaches.
Contrary to the previous approaches on multi-refinery pipelines, assuming a single-period planning horizon, the proposed model considered multiple intermediate due dates for product demands. The problem aimed to satisfy product demands Mixed-integer linear programming model for tree-like… 421 raised by local markets, by minimizing the total costs of pipeline including pumping, reprocessing of interface material and backordered demand. The proposed optimization framework was applied to four case studies. Example 1 and 2 considered a multi-level treelike pipeline with two refineries. Example 3 and 4 were to show the advantages of the proposed model based on both operational performance and computational cost. The results showed that the proposed MILP model is wider in scope and can find the optimal solution in a reasonable CPU time. Future work will involve extending the proposed approach with considering production scheduling at refineries and uncertainty on future demands of depots over a long time horizon.''