New batch-centric model for detailed scheduling and inventory management of mesh pipeline networks

Computers and Chemical Engineering 130 (2019) 106568

Contents lists available at ScienceDirect

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

New batch-centric model for detailed scheduling and inventory management of mesh pipeline networks Qi Liao a,b, Pedro M. Castro b,∗, Yongtu Liang a, Haoran Zhang c a

China University of Petroleum-Beijing, Beijing Key Laboratory of Urban oil and Gas Distribution Technology, Fuxue Road No.18, Changping District, Beijing 102249, PR China b Centro de Matemática Aplicações Fundamentais e Investigação Operacional, Faculdade de Ciências, Universidade de Lisboa, Lisboa 1749-016, Portugal c Center for Spatial Information Science, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwashi, Chiba 277-8568, Japan

a r t i c l e

i n f o

Article history: Received 17 April 2019 Revised 31 July 2019 Accepted 3 September 2019 Available online 4 September 2019

a b s t r a c t The scheduling of a multi-source and multi-destination pipeline with unidirectional flow is a challenging task since the flexibility of the network allows for products to take alternative routes. To tackle this issue efficiently, this paper develops a generic continuous-time batch-centric model that permits the injection, by an input node, and delivery, to an output node, of multiple batches over a time slot. The objective function is to minimize the operation cost related to pumps, contamination, inventory, backorder and capacity utilization, while satisfying forbidden product sequences and inventory limits. The model is validated through the solution of four benchmark problems, featuring a mesh network with growing complexity in topological structure. Compared to the literature, we are now able to generate detailed schedules with lower operating costs, due to a better tradeoff between interface (different product sequences) and inventory cost. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Pipelines are the common long-distance carrier of crude oil to refineries and of petroleum products to consumers. In China, the total mileage of crude oil and refined multiproduct pipelines reached 27 and 21 thousand kilometers in 2015, and will rise to 32 and 33 thousand kilometers during the “Thirteenth Five-Year Plan” period (2016–2020) (Commission National Development and Reform, 2017). Because crude oil is not a final product and will be further processed at refineries, crude oil scheduling mainly focuses on the optimal operation of the tanks feeding the crude distillation units rather than the pipeline (Hamisu et al., 2013; Panda and Ramteke, 2019). In contrast, pipelines for refined products often convey various grades of gasoline and fuel oil from one or more refineries to several depots for sales. Their scheduling should be determined after a multilateral discussion among suppliers (refineries), pipeline operators and dealers (depots). Given the fluctuant market demand and the product contamination generated during the transportation, there is a need for pipeline operations to ensure on-time delivery of sufficient amounts of qualified products to their respective destinations.

∗

Corresponding author. E-mail address: [email protected] (P.M. Castro).

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

In early practice, the decision-making of pipeline operation relied on heavy manual calculations, statistics, verification and proofreading work, involving a great deal of manpower and time. Operation profit is influenced by multiple factors that include energy consumption, throughput, inventory management and contamination treatment, factors often conflicting in practice, such as energy consumption versus throughput (Cafaro and Cerdá, 2012; MirHassani and Fani Jahromi, 2011; Mostafaei et al., 2016). As a result, it is difficult to find an optimal compromise satisfying all facility and operational constraints just by manual work. With the aid of computer-based methods, operations are now gradually creating high-quality schedules faster. Due to the diversity of pipeline structures and operational requirements, a growing number of novel mathematical models with different emphasis have been proposed during the past decade. Based on the representation of time and volume, these models can be divided into continuous-time (Cafaro et al., 2015a, 2012; MirHassani and Fani Jahromi, 2011; Mostafaei et al., 2015), discrete-time (Chen et al., 2017; Liao et al., 2018a), discrete-volume (De Souza Filho et al., 2013; Herrán et al., 2010), and continuousvolume, with the latter subdivided into product-centric (Castro, 2017; Castro and Mostafaei, 2017) or batch-centric approaches (Castro and Mostafaei, 2019; Liao et al., 2019a, 2019b). Each representation has advantages and limitations, and the preference depends upon the objectives and operation requirements involved.

2

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Nomenclature

max N,min vN, /vn,p n,p

Sets/Indices I/i, i batches Ilold old batches in line l Il batches that can appear in line l Il, n batches in line l that can appear in node n L/l, l lines in the pipeline network N/n, n nodes Nl Nodes along line l ND depot nodes NlD depot nodes along line l

vLl 0 vL, l,i

JI Nl JO Nl R N

input nodes of junctions along line l

NlR NlT P/p, p Pl Pn Pp S/s, s Sl SnD SUn T/t

output nodes of junctions along junction j refinery nodes refinery nodes along line l terminal node of line l products products that can appear in line l products that can appear in node n products that cannot enter a segment after p pipeline segments segments in line l downstream segments to node n upstream segments to node n event points of the time grid

Parameters cbn, p c f p,p cidn, p cpl, n, p cmk D,max D,min fl,n / fl,n J,max

fl,n

J,min

/ fl,n

flL,max / flL,min N,end fn,p R,max R,min fl,n / fl,n

f pP,max fsS,max / fsS,min h L,0 lcl,i L,0 rcl,i 0 vN, n,p

vN,end n,p

backorder cost for 1 of product p at node n ($/m3 ) cost of reprocessing the interface p − p ($) cost of storing 1 m3 of product p at node n ($/m3 ) cost of transporting 1 m3 of product p though line l to node n ($/m3 ) hourly cost for operating the pipeline ($/h) maximum/minimum volume transferred from line l to depot node n during a slot (m3 ) maximum/minimum volume transferred from/to line l through junction node n during a slot (m3 ) maximum/minimum volume of one batch transported in line l during a slot (m3 ) demand at the end of scheduling horizon for node n and product p (m3 ) maximum/minimum volume transferred from refinery node n to line l during a slot (m3 ) maximum volume of product p that can be transferred during the time horizon (m3 ) maximum/minimum volume into segment s during a slot (m3 ) time horizon (h) left coordinate of batch i in line l at the start of the time horizon (m3 ) right coordinate of batch i in line l at the start of the time horizon (m3 ) initial volume of product p in storage system of refinery or depot node n (m3 ) minimum inventory of product p in node n at the end of scheduling horizon (m3 )

initial volume of batch i in line l (m3 ) equal to one when old batch i in line l contains product p volumetric coordinate of node n in line l (m3 )

yl, i, p

σ l, n D,max D,min ρl,n /ρl,n maximum/minimum flowrate that can be LM,max ρn,p R,max R,min ρl,n /ρl,n

ρsS,max /ρsS,min

transferred from line l to depot node n (m3 /h) maximum flowrate that can be transferred from depot node n to local market (m3 /h) maximum/minimum flowrate that can be transferred from refinery node n to line l (m3 /h) maximum/minimum flowrate in segment s (m3 /h)

Variables L Xl,i,t binary variable indicating batch i is transported in line l during slot t L,d Xl,i,p,t binary variable indicating batch i containing prodD Xl,n,i,t D,idle Xl,n,t

m3

maximum/minimum inventory of product p in node n (m3 ) volume of line l (m3 )

D,1 Xl,n,i,t

D,2 Xl,n,i,t

J

Xl,n,i,t J,idle

Xl,n,t

R Xl,n,i,t R,idle Xl,n,t S Xs,t I Xl,i,i  ,p,p

Yl, i, p Bn, p LCl, i, t Dt D Fl,n,i,t D,d Fl,n,i,p,t J Fl,n,i,t J,d

Fl,n,i,p,t

uct p is transported in line l during slot t binary variable indicating node n is receiving batch i from line l during slot t binary variable indicating node n is not receiving batches from line l during slot t binary variable indicating batch i is the first to be transferred from line l to depot node n during slot t binary variable indicating batch i is not the first to be transferred from line l to depot node n during slot t binary variable indicating node n is transferring batch i from/into line l during slot t binary variable indicating node n is not transferring batches from/into line l during slot t binary variable indicating node n is pumping batch i into line l during slot t binary variable indicating node n is not pumping batches into line l during slot t binary variable indicating segment s is active during slot t binary variable indicating batch i of product p is being conveyed after batch i of product p in line l during slot t binary variable indicating batch i in line l contains product p backorder of product p at node n left coordinate of batch i in line l at event point t (m3 ) duration of time slot t (h) volume of batch i transferred from line l to depot node n during slot t (m3 ) volume of batch i of product p transferred from line l to depot node n during slot t (m3 ) volume of batch i transferred from/to line l through junction node n during slot t (m3 ) volume of batch i of product p transferred from/to line l through junction node n during slot t (m3 )

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

LM Fn,p,t

volume of product p sent to local market during

R Fl,n,i,t

from node n during slot t (m3 ) volume of batch i transferred from refinery node n

R,d Fl,n,i,p,t

to line l during slot t (m3 ) volume of batch i of product p transferred from re-

Tt N Vn,p,t

finery node n to line l during slot t (m3 ) volume transported in segment s during slot t (m3 ) the interface between batches i and i in line l ($) right coordinate of batch i in line l at event point t (m3 ) time of event point t (h) inventory of product p in node n at event point t

L Vl,i,t

(m3 ) volume of batch i existing in line l at event point t

S Fs,t

MIXl,i,i RCl, i, t

(m3 )

Normally, the continuous model is superior to its discrete counterpart for the smaller model scale, except in the cases of timevarying cost, demand and flowrate-related constraints. In terms of continuous-volume models, the batch-centric approach allows multiple batches (products) to be processed by a node over a slot, resulting in a need for fewer slots compared to the productcentric approach that allows at most one (Castro and Mostafaei, 2019; Liao et al., 2019b). Considering the impact of the number of time slots on problem size and computational performance, the batch-centric approach is an overall better option in most cases, but it may perform worse than its product-centric counterpart in a multi-source or a reversible system with short-distance segments. This is because a product-centric approach needs not to define an initial batch sequence neither the total number of batches, which are known to influence solution quality (Castro, 2017; Liao et al., 2019a). Besides improving computational performance, recent work is gradually focusing on the integration of pipeline scheduling with pump scheduling (Cafaro et al., 2015b), inventory management (Dimas et al., 2018; Moradi and MirHassani, 2014; Relvas et al., 2013) and contamination reprocessing (Cafaro and Cerdá, 2012; Liao et al., 2018b; Mostafaei et al., 2016, 2015). It is also considering more practical factors like hydraulic calculation, pump switching, time-varying electricity cost for pumps, time-dependent and volume-dependent inventory cost, forbidden product sequences, treatment fee for different types of contamination, forbidden delivery operations of contamination at intermediate depots and forbidden idle periods for segments with contamination inside. As the problem complexity grows, commercial optimizers may fail at finding a satisfactory or even a feasible solution within acceptable time. Therefore, a series of improved algorithm and decomposition strategies were put forward in the literature (Herrán et al., 2012; Relvas et al., 2009; Zhang et al., 2018, 2017) to save computation time without severely compromising solution quality. Although several papers have been published in the field of multiproduct pipeline scheduling, few of them dealt with a multisource and multi-destination network. To avoid the high computational burden from solving a detailed model for such a complex arrangement, some authors adopted decomposition strategies. Lopes et al. (2010) proposed a two-stage hybrid approach for a large-scale network as an extension of their earlier work (Moura et al., 2008a, 2008b). It consists of a planning phase that generates delivery orders (product, volume, origin and destination tanks) through heuristics and a scheduling phase that determines the start and ending times for each order with a Constraint Programming (CP) model. Boschetto et al. (2010) presented a hierarchical approach consisting of four sequential blocks: assignment,

3

sequencing, simulation, and timing; which rely on a continuoustime mixed-integer linear programming (MILP) model and a set of heuristics. It was later improved by using two new MILP models for the assignment and sequencing blocks (Magatão et al., 2012), as well as another for the lower-level timing block (Magatão et al., 2015). The approach is capable of finding feasible solutions in reasonable computational time while satisfying operational constraints related to production/consumption requirements, pipeline stoppages, batch tracking, use of preferential routes to avoid contamination losses, avoiding pumping on peak electricity prices, reversible flows, surge tank operations, etc. However, it cannot guarantee global optimality. Cafaro and Cerdá (2012) were the firsts to develop a continuous-time batch-centric MILP model that can generate an approximate schedule for a mesh-structure pipeline network in a single step. Afterwards, Castro (2017) presented a product-centric formulation for detailed scheduling that allowed for reversible flow. However, it assumes that each segment contains at most one batch per product, which may degrade solution quality or even compromise feasibility. To further enhance computational performance and solution quality, Chen et al. (2019) recently proposed a batch-centric model that rigorously enforces forbidden sequences. This paper extends our efficient batch-centric model (Liao et al., 2019a) from straight pipelines to generic pipeline networks with unidirectional flow. Assignment, sequencing and timing decisions are obtained for every node all at once. Compared to the work of Castro (2017) and Chen et al. (2019), which have minimized the makespan, the objective of this new model is to minimize an overall operation cost related to pumps, contamination, inventory, backorder and capacity utilization. More importantly, by allowing the operation of multiple batches at a node over a slot, the proposed model can find a better detailed schedule using fewer time slots and in shorter computational time. The remainder of the paper is organized as follows. The problem statement is given in Section 2 where we provide instructions on how to divide a mesh network into a system of lines, refineries, depots and junction nodes. We also reproduce the problem data for the benchmark instances, taken from Cafaro and Cerdá (2012). Basic concepts of batch-centric formulations, like the one being proposed, are the subject of Section 3, where we discuss the main reason for choosing batch numbering per line. The continuoustime representation is introduced in Section 4, while the remaining constraints of the model are given in Section 5. Section 6 presents the computational results for the four benchmark instances, analyzing the optimal schedules and inventory profiles. Finally, the conclusions are left for Section 7. 2. Problem statement This paper focuses on pipeline networks that comprise multiple refineries and depots. Throughout the paper, we consider the mesh-structure pipeline network in Cafaro and Cerdá (2012) and substructures of it. This is given in Fig. 1. The system conveys products from three refineries (N1, N2 and N4) to six depots (N3-N8) through six lines (L1-L6). Each line l ∈ L has a single source, located at its origin, and possibly multiple destinations, with one necessarily located at the end of the line. The source node can be a refinJI ery n ∈ NlR and/or an input node of a junction n ∈ Nl . The destiD nation nodes can be a depot n ∈ Nl and/or an output node of a JO

junction n ∈ Nl . In the specific case of Fig. 1 with no branches, all sources are refinery nodes and all destinations are depot nodes. The major difference between refineries and input nodes of junctions (or depots and output nodes of junctions) is that junctions have no storage tank to receive/send material from/to the pipeline, only transferring products from their connected upstream segments s ∈ SUn to downstream segments s ∈ SnD .

4

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Fig. 1. Mesh-structure pipeline network considered in this work.

Fig. 2. Possible operating modes for a junction node with a refinery and depot. While there may be multiple input and output lines, the illustration considers just two lines and two products for simplification.

In Fig. 1, central node N4 is the most interesting, since it acts as a triple-purpose node. As seen in Fig. 2, it can: (A) behave as a depot to receive products from lines L1 and L2; (B) behave as a refinery to pump products into lines L3, L4 and L6; (C) behave as a junction to transfer products from lines L1 and L2 to lines L3, L4 and L6. These are the simplest operating modes, but as discussed in the literature (Mostafaei et al., 2016; Mostafaei and Castro, 2017), it is desirable to simultaneously execute two of the above operations for a shorter makespan. If the simultaneous operation involves the transfer between two lines, the crossing product must be the same as the one being received by (D) or injected from the storage tanks (E). Otherwise, the injected product can be different from the received one (F). With two upstream lines and three downstream lines, several operating modes can be explored for node N4. In Fig. 3, node N4 is receiving two different products, P3 from line L1 and P1 from line L2. P3 is simultaneously being transferred to line L3 while P1 goes entirely into the storage system. Two other prod-

Fig. 3. The operation of a junction node with a refinery and depot connecting several lines can be very flexible.

ucts leave the tanks into the pipeline, P2 into line L4 and P4 into line L6. Another noteworthy aspect is that products can be delivered to their destinations from multiple refineries and through different routes. For each product, the preferential source and route to a

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

5

Table 1 Product inventory and demand data (102 m3 ). Product

Node

N1

N2

N3

N4

N5

N6

N7

N8

P1

Min Max Initial Demand Min Max Initial Demand Min Max Initial Demand Min Max Initial Demand

100 1000 1000 – 100 1000 1000 – 20 200 200 – – – – –

– – – – 100 1000 1000 – – – – – 40 400 200 –

50 500 50 150 50 500 50 150 10 100 30 10 – – – –

50 500 250 – 100 1000 700 – 10 200 100 – 20 300 100 –

50 500 50 800 100 1000 100 900 – – – – – – – –

– – – – – – – – – – – – 20 200 20 200

50 500 50 200 100 1000 100 400 10 100 10 250 – – – –

– – – – – – – – – – – – 10 100 10 160

P2

P3

P4

Table 2 Interface and inventory costs. Interface cost (102 $)

P1 P2 P3 P4

Inventory cost ($/m3 )

P1

P2

P3

P4

N1

N2

N3

N4

N5

N6

N7

– 184 340 235

184 – 250 413

340 250 – ×

235 413 × –

0.19 0.16 0.24 –

– 0.16 – 0.24

0.24 0.20 0.32 –

0.19 0.16 0.24 0.24

0.24 0.20 – –

– – – 0.32

0.24 0.20 0.32 –

certain depot is not fixed and may change over time. As an example, there are two alternative routes between N1 and N5. The direct route is through line L5 while the indirect route goes through line L1, intermediate node N4 and line L3. Nodes N5-N8 located in lines L3, L4 and L6 can receive products from the closer refinery N4 as well as from refineries N1-N2, located farther away. Shipping from N4 is a more recommendable source to reduce pumping costs in lines L1-L2. Overall, the preferential source should be determined based on a tradeoff amongst real-time inventory level, inventory cost and pumping cost. The data for product inventory and demand are listed in Table 1 and there is a total of four products involved (P1-P4). Product backorders are not allowed in this system, so all depots should dispatch the required demand to local markets at a maximum rate of 700 m3 /h. Table 2 presents the cost of reprocessing any interface and storing one m3 of any product. Note that the inventory of products P1-P4 in intermediate node N4 at the end of the schedule, is required to be higher than 250 0 0, 50 0 0 0, 50 0 0 and 10 0 0 0 m3 , respectively, so as to not compromise the next run. Note that P3-P4 and P4-P3 are forbidden product sequences. Table 3 shows the cost of pumping one m3 of any product, which depends on the origin, destination and route. The hourly cost for operating the pipeline is set to $10 0 0/h. 3. Batch-centric model basic concepts In a batch-centric model like the one proposed in this paper, products are transformed into batches when entering the pipeline. As discussed in Liao et al. (2019a), there are two possibilities for numbering batches, global numbering and line numbering. Global numbering was shown more appropriate for straight pipelines but this work deals with a mesh network for which it is quite difficult to define an appropriate initial batch sequence. This is because some lines (L3, L4 and L6) no longer have a single source, leading to far more product sequences when considering the full length of the time horizon. It is an issue since batch-centric models work under the assumption that batch i < i always precedes batch i in the sequence.

Table 3 Pumping costs ($/m3 ). Line

Origin

Product

L1

N1

L2

N2

L3

N4

L4

N4

L5

N1

L6

N4

P1 P2 P3 P2 P4 P1 P2 P1 P2 P3 P4 P1 P2 P4

Destination N3

N4

N5

N6

N7

N8

0.70 0.72 0.96 – – – – – – – – – – –

1.05 1.08 1.44 0.90 0.93 – – – – – – – – –

– – – – – 0.88 0.90 – – – – 1.22 1.25 –

– – – – – – – – – – 0.56 – – –

– – – – – – – 1.05 1.08 1.44 – – – –

– – – – – – – – – – – – – 0.61

To avoid this drawback, our mathematical model uses batch numbering per line, which also has the advantage of avoiding the definition of empty batches when characterizing the initial status of the pipeline (Liao et al., 2019a). That is to say, we need both batch (i) and line (l) indices to identify a batch of product p in the network. This is done through binary variables Yl, i, p , where each pair needs to be assigned to exactly one product, as stated by Eq. (1). The products initially on each line are defined as old batches (i ∈ Ilold ) using parameters yl, i, p , which are generated assuming that the numbering is from the end to the origin of the line. Eq. (2) then fixes the corresponding binary variables.



Yl,i,p = 1 ∀l, i ∈ Il

(1)

p∈Pl

Yl,i,p = yl,i,p ∀l, i ∈ Ilold , p ∈ Pl

(2)

With line numbering, the batch number can be reset after crossing a junction, thus making it possible for the batch (product) sequence to change. As seen in Fig. 4, the batches being transferred to line L4 can come from either line L1 or line L2, leading to

6

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Fig. 4. The model uses line batch numbering to facilitate considering all possible product sequences in networks featuring lines with multiple sources.

duration Dt .

T1 = 0

(3)

T|T | ≤ h

(4)

Tt+1 = Tt + Dt Fig. 5. The model relies on a continuous-time representation with a single grid.

different batch-product assignments in line L4. Notice that in scenario (A), batches (I1,L1) and (I2,L1) are renumbered to (I3,L4) and (I4,L4) after crossing N4, whereas in scenario (B), batches (I1,L2) and (I2,L2) become (I4,L4) and (I5,L4), respectively. 4. Time representation The model described in the next section uses a continuoustime representation with a single reference grid. As illustrated in Fig. 5, a series of event points t = 1, 2, …, |T| divide the planning horizon h into a total of |T | − 1 time slots with not necessarily equal durations. Let continuous variable Tt represent the absolute time of event point t. The absolute time of the first event point is 0 (Eq. (3)), while that of the last event point, indicating the makespan, cannot be larger than h (Eq. (4)). Eq. (5) states that the difference between event points t and t + 1 corresponds to the slot

∀t < |T |

(5)

One drawback of continuous-time models is that it is difficult to predict the required number of event points. Depending on the form of the objective function, a low value leads to an infeasible problem while a large value can lead to an intractable problem or to a solution with dummy slots (Dt = 0). The latter is undesirable since it affects the estimation of the inventory cost (check later Eq. (52)). The search for the optimal solution thus involves single increments in the value of |T|. 5. Main mathematical formulation The continuous-time mathematical model for mesh pipeline networks is derived from our recent work for multi-source straight pipelines (Liao et al., 2019a). The new batch numbering per line model is a generalization in the following ways: (i) a junction in the system now can connect to multiple upstream lines and multiple downstream lines, making the model more modular and applicable to pipeline networks with any structure; (ii) the objective is

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

data. L Vl,i,t =

 L vL,l,i0 t=1 + Vl,i,t−1 + 

−

J Fl,n,i,t−1



7



R Fl,n,i,t−1 +

n∈NlR :i∈Il,n

−

J Fl,n,i,t−1

n∈NlJI :i∈Il,n



D Fl,n,i,t−1

∀l, i ∈ Il , t

(9)

n∈NlD :i∈Il,n

n∈NlJO :i∈Il,n

5.2. Node constraints

Fig. 6. Forbidden product sequences may appear between non-consecutive batches.

to find the minimum operational cost that includes pumping, backorder, interface, inventory and the utilization cost, rather than just makespan; (iii) new sets of constraints are proposed to estimate the interface cost generated during transportation. The main assumptions of the model are: • Unidirectional flow along every line. • A node can process the incoming/outgoing batches from/to several lines per time slot. • During a slot, an intermediate depot is allowed to receive multiple batches only when its downstream segment is idle. • Refineries and terminal depots that are not junctions are allowed to inject/receive multiple batches per time slot. If such nodes are also junctions, this can happen only if no batch is crossing. This assumption and the previous improve computational performance by reducing the number of slots. Note that schedules are not removed from the solution space since one can always increase the number of slots (Liao et al., 2019a, 2019b). • Considering the required time to monitor product quality, products received at a dual-purposed node should wait until the next time slot to be injected into downstream lines (Cafaro and Cerdá, 2012). • The contamination between two batches is modelled as an interface cost that depends on the adjacent products but not on the operating flowrate.

i ∈ Il : i ≥ i L Vl,i,t = RCl,i,t − LCl,i,t



L Vl,i,t = vLl

∀l, t

Vl,iL  ,t

R to be an inequality since variable Xl,n,i,t has an index for batch i R,idle while Xl,n,t does not. Supposing that only batch i is being pumped R,idle R R by the refinery (i.e. Xl,n,t = 0, Xl,n,i,t = 1 and Xl,n,i  ,t = 0 for any

i = i), the left side of Eq. (10) for any i = i will be zero and not one. R,idle R Xl,n,i,t + Xl,n,t ≤ 1 ∀l, n ∈ NlR , i ∈ Il,n , t < |T |

∀l, i ∈ Il , t

∀l, i ∈ Il , t

(6)

(7) (8)

i∈Il

Eq. (9) represents the multiperiod volumetric balance of batch i in line l. Compared to the previous event point, the volume increases through injection from a refinery or incoming flow from an input junction and decreases from outgoing flows to output junc0 tions and deliveries to depots. The initial batch volume vL, is given l,i

(10)

If the node is working as a refinery, the pumped volume should R,min R,max be within range [ fl,n , fl,n ], as stated by Eq. (11). The pumping  R rate, obtained by dividing the total pumped volume i∈I Fl,n,i,t by l,n

the slot duration Dt , should also satisfy the given flowrate limit R,min R,max [ρl,n , ρl,n ], as stated by Eq. (12). Since the refinery is located at the origin of a line, only batches with zero left coordinate at event point t can be injected during slot t (Eq. (13)). R,min R R,max R R fl,n Xl,n,i,t ≤ Fl,n,i,t ≤ fl,n Xl,n,i,t



R i∈Il,n Fl,n,i,t R,max l,n

ρ

∀l, n ∈ NlR , i ∈ Il,n , t < |T | (11)



R i∈Il,n Fl,n,i,t R,min l,n

≤ Dt ≤

ρ

R,idle + h · Xn,t

∀l, n ∈ NlR , t < |T | (12)

R LCl,i,t ≤ vLl · 1 − Xl,n,i,t

L , RC Continuous variables Vl,i,t l, i, t and LCl, i, t are used to specify the volume, right and left coordinates of batch i in line l at event point t. Eq. (6) states that the right coordinate for batch i in line l is obtained by summing the volume of all batches i ≥ i. The batch volume inside the line is computed by the difference between its right and left coordinates, as stated by Eq. (7). Since the liquid products are assumed incompressible, the volume of all batches inside line l must be equal to the line capacity vLl .



5.2.1. Refinery R Let binary variable Xl,n,i,t = 1 indicate that refinery n is pumping batch i into line l during time slot t. Eq. (10) also accounts for the R,idle refinery being idle, a case with Xl,n,t = 1. Note that Eq. (10) needs



5.1. Batch tracking in lines

RCl,i,t =

A node can have up to three functions: refinery, depot and junction. The constraints in this section are shared with our previous work (Liao et al., 2019a).



∀l, n ∈ NlR , i ∈ Il,n , t < |T |

(13)

5.2.2. Depot Some of the constraints for batch delivery operations to a depot are similar, as seen in Eqs. (14)–(16), with the variables and parameters reflecting the change of superscript from refinery R to depot D. D,idle D Xl,n,i,t + Xl,n,t ≤ 1 ∀l, n ∈ NlD , i ∈ Il,n , t < |T | D,min D D,max D D fl,n Xl,n,i,t ≤ Fl,n,i,t ≤ fl,n Xl,n,i,t



D i∈Il,n Fl,n,i,t D,max l,n

ρ

(14)

∀l, n ∈ NlD , i ∈ Il,n , t < |T | (15)

 ≤ Dt ≤

D i∈Il,n Fl,n,i,t D,min l,n

ρ

D,idle + h · Xn,t

∀l, n ∈ NlD , t < |T | (16) NlT

It is now important to distinguish between terminal and intermediate NlD \NlT depots of line l. Eq. (17) states that batch i can be delivered to terminal node n during slot t only if its right coordinate at the end of the slot is equal to the line volume vLl . D RCl,i,t+1 ≥ vLl Xl,n,i,t

∀l, n ∈ NlD ∩ NlT , i ∈ Il,n , t < |T |

(17)

The constraints for an intermediate depot are more complex. During slot t, the node can receive batches with right coordinate

8

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Fig. 7. Optimal solution for Ex1 (|T| = 8).

at event point t+1 greater or equal than the node location σ l, n and left coordinate less or equal than σ l, n , as seen in Eqs. (18)– (19). The right coordinate at the start of the slot is also required to exceed the node location σ l, n if batch i is the first to be reD,1 ceived (Xl,n,i,t = 1), Eq. (20). Eq. (21) then ensures that there can be at most one batch first-received. If the batch is delivered to D,2 the intermediate depot, it is either the first or non-first (Xl,n,i,t = 1), as stated by Eq. (22). Finally, multiple deliveries are allowed only

S = 0). when its adjacent downstream segment s is idle (Xs,t D RCl,i,t+1 ≥ σl,n Xn,i,t

∀l, n ∈ NlD \NlT , i ∈ Il,n , t < |T | 

D D LCl,i,t+1 ≤ σl,n Xl,n,i,t + vLl 1 − Xl,n,i,t



(18)

∀l, n ∈ NlO \NlT , i ∈ Il,n , t < |T | (19)

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

9

of batch i to occur. More specifically, the value of the right coordinate at event point t must be greater than the node location σ l, n while the left coordinate at event point t+1 must not exceed σ l, n . J,min J J J,max J fl,n Xl,n,i,t ≤ Fl,n,i,t ≤ fl,n Xl,n,i,t

∀l, n ∈ NlJI ∪ NlJO , i ∈ Il,n , t < |T | (27)

J RCl,i,t ≥ σl,n Xl,n,i,t

∀l, n ∈ NlJI ∪ NlJO , i ∈ Il,n , t < |T | 

J J LCl,i,t+1 ≤ σl,n Xl,n,i,t + vLl 1 −Xl,n,i,t



(28)

∀l, n ∈ NlJI ∪ NlJO , i ∈ Il,n , t < |T | (29)

Fig. 8. Inventory profile for N4 in Ex1.

D,1 RCl,i,t ≥ σl,n Xl,n,i,t



∀l, n ∈ NlD \NlT , i ∈ Il,n , t < |T |

D,1 Xl,n,i,t ≤ 1 ∀l, n ∈ NlD \NlT , t < |T |

(20) (21)

Nodes belonging to both refinery and junction sets have the flexibility to simultaneously receive products from the refinery and upstream line(s). However, the injected batch should be in accordance with the crossing one (recall Fig. 2(E)). This is accomplished by Eq. (30). Similarly, junction nodes simultaneously dispatching products to a depot and downstream line(s) should satisfy Eq. (31), as illustrated in Fig. 2(D).





i ∈Il,n ,i =i

∀l, n ∈ NlD \NlT , i ∈ Il,n , t < |T |





(22)

i ∈Il,n ,i =i D,2 S Xl,n,i,t + Xs,t ≤ 1 ∀l, n ∈ NlD \NlT , i ∈ Il,n , s ∈ Sl ∩ SnD , t < |T |

J

J

FL1,N4,I2,t = FL4,N4,I4,t , but these constraints are not known beforehand since the batch-line matchings are only determined at run time. The best we can do is to restrict the transfer to at most one product per time slot, ensuring that the volume entering the junction is equal to the volume leaving. Notice that Eq. (25) features J,d disaggregated variables by product Fl,n,i,p,t that can be related to J

Fl,n,i,t by Eq. (26).



J Xl,n,i,t ≤ 1 ∀l, n ∈ NlJI ∪ NlJO , t < |T |

(24)





J,d Fl,n,i,p,t =

l:n∈NlJI ,p∈Pl i∈Il,n





i∈Il  ,n l  :n∈N JO  ,p∈Pl 

FlJ,d  ,n,i,p,t

p ∈ Pn , t < |T | J Fl,n,i,t

=



(31)

In the batch-centric models of Mostafaei and Castro (2017) and Castro and Mostafaei (2019), segments are a building block of the network and an important model entity. This is no longer the case in this work, where lines replace segments as a building block. Note that a line may comprise one or more segments, with depots dividing consecutive segments. Given that intermediate depots may be responsible for a large portion of the line’s product demand, it is not uncommon for the segments to have different diameters (Castro and Mostafaei, 2019), which is translated into different intervals [ρsS,min , ρsS,max ] for the operating flowrate. S give the volume crossing segment s Let continuous variable Fs,t during slot t. The volume entering node n belonging to line l from its upstream segment s ∈ Sl ∩ SUn , refinery or input junction, must be equal to the volume leaving node n to its downstream segment s ∈ Sl ∩ SnD , depots and output junctions. Eq. (32) thus relates variS with the other volumetric variables. Eq. (33) states that ables Fs,t the former variables are different than zero only if the segment S = 1). The flowrate constraints are then enforced by is active (Xs,t Eq. (34).

=

l

 |Il,n | − 1 ∀l, n ∈ NlD ∩ NlJO ,

5.3. Segment constraints

S Fs,t +

s∈Sl ∩SUn

∀n ∈ NlJI ∪ NlJO ,

(30)

i ∈ Il,n , t < |T |



i∈Il,n

 

J D Xl,n,i  ,t ≤ 1 − Xl,n,i,t ·

(23)

5.2.3. Junction The main difference between junction nodes and refinery and depot nodes, is that at most one batch can be allowed to cross the junction per time slot (Liao et al., 2019a), translating into Eq. (24). This is a direct consequence of considering batch numbering per line, which causes some modeling difficulties. Let continuous variJ able Fl,n,i,t specify the volume of batch i entering/leaving line l through junction node n during slot t and recall the example in Fig. 4(A). The first two batches in line L1 are crossing junction N4 into line L4, being transformed in batches (I3,L4) and (I4,L4). To J J verify the volume balance, we must have FL1,N4,I1,t = FL4,N4,I3,t and

 |Il,n | − 1 ∀l, n ∈ NlR ∩ NlJI ,

i ∈ Il,n , t < |T |

i∈Il,n D,1 D,2 D Xl,n,i,t = Xl,n,i,t + Xl,n,i,t

 

J R Xl,n,i  ,t ≤ 1 − Xl,n,i,t ·



 

l:n∈NlR i∈Il,n

FsS ,t +

s ∈Sl ∩SnD

(25)

R Fl,n,i,t +

 

 

J Fl,n,i,t

l:n∈NlJI i∈Il,n D Fl,n,i,t +

l:n∈NlD i∈Il,n

 

J Fl,n,i,t ∀l, n ∈ Nl , t < |T |

l:n∈NlJO i∈Il,n

(32) J,d Fl,n,i,p,t

∀l, n ∈

NlJI

∪

NlJO , i

∈ Il,n , t < |T |

(26)

p∈Pl ∩Pn

The rest of the constraints are similar to before. Eq. (27) places minimum and maximum limits on the volume crossing the junction, while Eqs. (28) and (29) define the conditions for the transfer

S S S fsS,min Xs,t ≤ Fs,t ≤ fsS,max Xs,t

ρ

S Fs,t S,max s

≤ Dt ≤

S Fs,t

ρsS,min

∀s, t < |T | 

S + h · 1 − Xs,t



∀s, t < |T |

(33)

(34)

10

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Fig. 9. Optimal solution for Ex2 (|T| = 5).

5.4. Inventory management Every refinery or depot is linked to a storage system with one or more dedicated tanks for every product. Since the movement inside the pipeline is modelled through batches, we need to disaggregate the variables giving the volumes leaving the refinery and being delivered to the depots, in order to control product inventory. This is done through Eqs. (35) and (36), with Eq. (37) forcing all disaggregated variables to be zero when product p is not assigned to batch i of line l (Yl, i, p = 0). Notice that the disaggregated junction variables defined in Eq. (26) are also included in Eq. (37).



R Fl,n,i,t =





n∈NlJI ∪NlJO :i∈Il,n ,p∈Pn

t=|T |

+

R,d Fl,n,i,p,t

∀l, n ∈ NlR , i ∈ Il,n , t < |T |

D,d Fl,n,i,p,t

∀l, n ∈

(35)



=

NlD , i

∈ Il,n , t < |T |

p∈Pl ∩Pn





n∈NlR :i∈Il,n ,p∈Pn t=|T |

R,d Fl,n,i,p,t +





n∈NlD :i∈Il,n ,p∈Pn t=|T |

D,d Fl,n,i,p,t

(36)

(37)

N ) should be equal to the Product inventory at event point t (Vn,p,t N inventory at the previous event point (Vn,p,t−1 ) minus the volume   R,d being injected to downstream lines ( i l Fl,n,i,p,t−1 ) and the local LM market during slot t − 1 (Fn,p,t−1 ), plus the volume being received   D,d from upstream lines during that same slot ( i l Fl,n,i,p,t−1 ). Note 0 that the inventory at the first event point vN, n,p in Eq. (38) is given beforehand. N Vn,p,t =

  0 N vN, n,p t=1 + Vn,p,t−1 −

p∈Pl ∩Pn D Fl,n,i,t

J,d Fl,n,i,p,t ≤ f pP,maxYl,i,p ∀l, i ∈ Il , p ∈ Pl

+







i∈Il,n l:n∈N R ,p∈Pl l

D,d Fl,n,i,p,t−1



R,d LM  Fl,n,i,p,t−1 − Fn,p,t−1

∀n ∈ NR ∪ ND , p ∈ Pn , t

n∈N D

(38)

i∈Il,n l:n∈N D ,p∈Pl l

Eq. (39) takes into account the given upper bound for the delivLM,max ery rate from the depot to the local market, ρn,p . The total sent N,end volume should also meet the given demand fn,p but this is not

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

11

enforced as a hard constraint. Let nonnegative continuous variable Bn, p represent the backorder of product p at node n. By Eq. (40), it takes a positive value whenever the total volume sent to the local  LM < f N,end . market is not enough, i.e. t<|T | fn,p,t n,p LM LM,max Fn,p,t ≤ ρn,p Dt



∀n ∈ ND , p ∈ Pn , t < |T |

LM N,end fn,p,t + Bn,p = fn,p

(39)

∀n ∈ ND , p ∈ Pn

(40)

t<|T |

At any time, product inventory should be within the feasible min N,max range [vN, ]. n,p , vn,p min N max vN, ≤ Vn,p,t ≤ vN, ∀n ∈ NR ∪ ND , p ∈ Pn , t n,p n,p

(41)

At the end of the scheduling horizon, there should be enough inventory vN,end to not compromise the optimization run for the n,p next planning period. N R D vN,end ≤ Vn,p, n,p |T | ∀n ∈ N ∪ N , p ∈ Pn

(42)

According to the fifth assumption in Section 5, the injected volume of product p by refinery node n during slot t should be less than the inventory at the start of that slot minus the minimum level.





R,d N min Fl,n,i,p,t ≤ Vn,p,t − vN, n,p

∀n ∈ NR ∩ ND , p ∈ Pn , t

(43)

i∈Il,n l:n∈N R ,p∈Pl l

5.5. Forbidden product sequences The transportation of multiple products in the pipeline leads to contamination at the interface that should be minimized. It can be done through sequence-dependent costs (see Section 5.6) and by enforcing forbidden sequences for incompatible products. Let subset Pp include products p that cannot follow product p. Assuming product p has been assigned to batch i in line l (Yl, i, p = 1), no prod uct p ∈ Pp can be assigned to batch i + 1 ( p ∈P ∩Pp Yl,i+1,p = 0).

Yl,i,p +



l

Yl,i+1,p ≤ 1 ∀l, i ∈ Il , p ∈ Pl

(44)

p ∈Pl ∩Pp

Eq. (44) is however not sufficient due to the possible appearance of empty batches following a complete batch removal at an intermediate depot, as seen in the last row of Fig. 6. Let binary L,d variable Xl,i,p,t = 1 if batch i containing product p is already inside L > 0), or if it has entered the line line l at the start of the slot t (Vl,i,t  R from a refinery ( n∈NR :i I Fl,n,i,t > 0) or an input junction node l,n l  J during this slot ( Fl,n,i,t > 0). Such binary variable is acJI n∈N :i Il,n l

tivated by Eq. (45), with Eq. (46) ensuring that the batch-product assignment is respected. Compared to Liao et al. (2019b) that features a model with global instead of line batch numbering, Eq. (46) is now written for every line l. Forbidden sequences are then enforced through Eq. (47).

flL,min



L,d Xl,i,p,t

≤

L Vl,i,t

+



R Fl,n,i,t

n∈NlR :i Il,n

p∈Pl

≤ flL,max



L,d Xl,i,p,t



+

J Fl,n,i,t

n∈NlJI :i Il,n

∀l, i ∈ Il , t < |T |

(45)

p∈Pl

 t<|T |

L,d Xl,i,p,t ≤ Yl,i,p · (|T | − 1 ) ∀l, i ∈ Il , p ∈ Pl





i
L,d L,d Xl,i  ,p ,t ≥ Xl,i,p,t +

 p ∈Pl ∩Pp

(46)

L,d Xl,i  ,p ,t

− 1 ∀l, i , i ∈ Il , i < i, p ∈ Pl , t < |T |

Fig. 10. Inventory profile for N4 in Ex2.

5.6. Interface cost The contamination volume between any two adjacent batches is assumed a constant value in this paper and so the interface cost only depends on the number of interfaces and on the cost of each interface. Most previous works (Mostafaei and Ghaffari Hadigheh, 2014; Zhang et al., 2017) calculate interface cost based on pump sequences, with the interface cost of consecutives batches i and i+1 in line l, MIXl,i,i+1 , being computed by Eq. (48), where c f p,p is the given cost of reprocessing interface p − p .



∀l,

i < |Il |, p ∈ Pl , p ∈ Pl \Pp , t < |T |

(48)

However, as explained in the previous section, this method is not accurate since a new interface may be generated when a batch is fully dispatched to an intermediate node. In the example of Fig. 6, it can be observed that interface P2-P1 was not created at the injection point but still appears in the line at event point t+2, being received later by the terminal node. To compute a rigorous interface cost, we need to track all interfaces generated over the planning horizon. I Let auxiliary variable Xl,i,i  ,p,p = 1 if the interface between

batches i and i , featuring respectively products p and p , appears inside line l. Interface p − p is generated during slot L,d t, if the two batches appear inside the line (Xl,i,p,t = 1 and L,d Xl,i  ,p ,t = 1) and there are no other non-empty batches in between   L,d ( i
Eq. (49). L,d L,d L,d Xl,i,p,t ∧i
∀l, i ∈ Il , i, i < i , p ∈ Pl , p ∈ Pl \Pp , t < |T |

(49)

Eq. (49) can then be reformulated into Eq. (50), and its righthand side used in Eq. (51). Note that the reprocessing cost of an existing interface between batches i and i in line l is never lower than the unit cost c f p,p , being actually equal to this value since MIXl,i,i appears as a positive term in the objective function that is being minimized. If the interface does not exist, then the term inside the brackets in Eq. (51) will be less than or equal to zero and variable MIXl,i,i will assume the value of its lower bound (zero). L,d L,d I Xl,i,i  ,p,p ≥ Xl,i ,p ,t + Xl,i,p,t −

(47)



MIXl,i,i+1 ≥ c f p,p Yl,i,p + Yl,i+1,p − 1

 

i
L,d Xl,i  ,p ,t

− 1 ∀l, i ∈ Il , i, i < i , p ∈ Pl , p ∈ Pl \Pp , t < |T |

(50)

12

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Fig. 11. Optimal solution for Ex3 (|T| = 6).

 L,d L,d MIXl,i,i ≥ c f p,p Xl,i  ,p ,t + Xl,i,p,t −

  i
 L,d Xl,i  ,p ,t − 1

i ∈ Il , i ∈ Il , i < i , p ∈ Pl , p ∈ Pl \Pp , t < |T |

∀l, (51)

5.7. Objective function The objective function in Eq. (52) is to minimize the total operating cost that comprises five terms. The first term is the pumping cost, where the unit cost for taking a product from an origin to a destination is multiplied by the volume and added over all time

slots. Parameter cpl, n, p gives the cost for pumping one m3 of product p from the origin of line l to a depot or output junction node. In case of shipping a product from N2 to N7, the cost will be divided into two parts: (i) reaching node N4 from line L2, through J,d variables FL2,N4,i,p,t ; (ii) transportation from N4 to N7 in line L4, through variables FLD,d . 4,N7,i,p,t The second term gives the backorder cost, where the unit cost cbn, p is multiplied by the backorder volume computed by Eq. (40). The third term is the interface cost, independent of the interface volume that is assumed constant, as previously explained. Compared to Cafaro and Cerdá (2012), this term is more accurate since it considers the interfaces created for all touching batches (new

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

13

Table 4 Results for increasing batch number. Ex

|T|

|Il |b

CPUs

Pumping cost (102 $)

Inventory cost (102 $)

Interface cost (102 $)

Total cost (102 $)

Make-span (h)

Back-order (m3 )

Ex1

5

5,3,4,6 6c ,4c ,4,7 5c ,4c ,5c ,8c 5,3,4c ,6 6c ,4c ,4c ,7 5,3,4c ,6 6c ,4c ,4,7 5,3,4c ,6 6c ,4c ,4c ,7c 5,3,4c ,6 6,2,2,6,2 7c ,3c ,3c ,7,3c 7c ,3c ,3c ,8c ,3c 6,2,2,6,2 7c ,3c ,3c ,7,3c 7c ,3c ,3c ,8c ,3c 6,4,2,6,2,1 7c ,5c ,3c ,7c ,3c ,1 6,4,2,6,2,1 7c ,5c ,3c ,7,3c ,1 6,4,2,6,2,1 7c ,5c ,3c ,7,3,1 6,4,2,6,2,1 7c ,5c ,3c ,7,3c ,1 6,4,2,6,2,1

25.1 781 3908 152 7200 189 7200 5315 7200 7200 10.2 283 840 56.8 3617 7200 15.1 801 238 7200 138 7200 2855 7200 7200

– 5378.7 5378.7 5417.7 5446.5 5417.7 5431.2 5417.7 5417.7 5417.7 5023.5 5019.8 5019.8 5023.5 5019.8 5019.8 5263.9 5263.9 5263.9 5263.9 5308.9 5380.6 5323.0 5324.2 5324.5

– 586.0 586.0 652.8 648.1 610.1 519.8 586.5 617.4 592.2 569.1 556.4 556.4 532.5 522.4 522.4 556.7 556.7 515.4 516.1 560.7 577.4 504.5 583.9 551.9

– 3936 3936 3686 3523 3686 3876 3686 3686 3686 3529 3351 3351 3529 3351 3351 4355 4355 4355 4177 4355.0 4967.0 4171.0 4177.0 4077.0

– 11640.7 11640.7 11606.5 11567.6 11529.5 11617.0 11480.2 11511.1 11485.9 10621.6 10427.2 10427.2 10585.0 10393.1 10393.1 11675.6 11675.6 11634.3 11457.0 26749.6 16400.9 25573.5 25585.1 25528.4

– 174 174 185 195 181.6 179 179 179 179 150 150 150 150 150 150 150 150 150 150 152.5 147.6 157.5 150.0 157.5

– – – – – – – – – – – – – – – – – – – – 7500 2000 7000 7000 7000

6 7 8

Ex2

9a 5

6a

Ex3

5

Ex4

6a 6 5 6 7

a b c

One of the slots is dummy at the optimum. Values of |Il | in each line are separated by commas. At least one of the batches is fictitious at the optimum.

interfaces may be generated between old batches due to the appearance of empty batches), instead of only considering the interface cost between new batches and their immediate predecessors. The fourth term gives an estimation of the inventory cost by dividing the amount in storage by the number of slots in the grid and then multiplying by the unit cost cidn, p . As mentioned in Section 4, we should check for the appearance of dummy slots in the objective function, which by increasing the denominator of the inventory cost term, are an artificial way of reducing the total cost. It should also be noted that inventory costs should be computed in the trapezoidal way (Lee et al., 1996) by multiplying the average tank inventory by the duration of the time slot. However, this would lead to the appearance of nonconvex bilinear terms, as seen in Castro (2016), turning the model into a mixed-integer nonlinear problem, much more difficult to solve to global optimality. This is no longer an issue for discrete-time models for which the duration of every time slot is known a priori. Finally, the last term gives a capacity utilization cost, aimed at reducing pipeline stoppages. It is computed by multiplying the hourly utilization cost cmk by the makespan, given by the absolute time of the last event point in the grid, T|T| .

min







c pl,n,p

l  L n∈N D ∪N JO p∈Pl ∩Pn l l

+



cbn,p Bn,p +





n∈N R ∪N D p∈Pn

t<|T | i∈Il,n

 



D,d  Fl,n,i,p,t

n∈NlD





J,d  + Fl,n,i,p,t

n∈NlJO

MIXl,i,i

l  L i∈Il i ∈Il ,i >i

n∈N D p∈Pn

+

 

cidn,p

N  Vn,p,t t>1

|T | − 1

+cmk · T|T |

(52)

6. Computational results The model was implemented in GAMS 25.1.3 and solved by CPLEX 12.8.0 running in parallel deterministic mode using up to eight threads. The termination criteria were a relative optimality tolerance of 10−6 or a maximum wall time limit of 7200 CPUs. The

hardware consisted of a Windows 10, 64-bit desktop with an Intel i7-6700 K (4.0 GHz) processor and 16 GB of RAM. The proposed model is validated by solving three benchmark instances from Cafaro and Cerdá (2012), Ex1-Ex3, with growing complexity in topological structure. To observe the impact of product backorder on the schedules, we also provide another instance, Ex4, derived from Ex3 but with higher demands. More specifically, the demands of P1 and P2 at depot N5 were increased by 50 0 0 m3 for a unit backorder cost of $200/m3 . As discussed in our previous work (Liao et al., 2019a), two tuning parameters greatly affect the performance of a batch-centric model with line numbering: |T| and |Il |. We first initialize parameters |Il | with the values from the best-found schedules in Cafaro and Cerdá (2012). Then, we increase the number of event points |T| as well as the number of batches in each line |Il | one by one until: (i) solution quality is not improved within 7,200 CPUs; (ii) a fictitious slot appears, leading to a poor estimation of average inventory cost (values in italic in Table 4); (iii) a fictitious batch appears in every line (two consecutive batches with the same product assignment). From Table 4, one can see that for the same number of events, a single increase in the number of batches for every line leads to an order of magnitude increase in computational time, at least. While there is the guarantee that solution quality cannot degrade, there are three cases where CPLEX is unable to find as good a solution inside the 7200 CPUs (Ex1 for |T| = 7 and |Il | = (5,3,4,6) vs. (6,4,4,7), and for |T| = 8 and |Il | = (5,3,4,6) vs (6,4,4,7), and Ex4 for |T| = 6 and |Il | = (6,4,2,6,2,1) vs. (7,5,3,7,3,1)). Another interesting result occurs for Ex3 and |T| = 6. For |Il | = (6,4,2,6,2,1) there is a dummy slot, but it disappears when increasing the number of batches to (7,5,3,7,3,1). Overall, an increase in the number of batches and events tends to decrease most terms in the objective function. Still, it is not uncommon to find an increase in inventory or interface cost despite the reduction in total cost. Therefore, it is crucial to aim at the integrated optimization of lot-sizing, batch sequence and storage infrastructure, based on predictable market demand.

14

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568

Fig. 12. Optimal solution for Ex4 (|T| = 5).

According to the results in Table 4, the optimal parameters |T| and |Il | should be set as 8 and (5,3,3,6) for Ex1, 5 and (6,2,2,7,2) for Ex2, 6 and (6,4,2,7,2,1) for Ex3, and 5 and (6,4,2,7,3,1) for Ex4. Next, we focus on the comparison of Ex1-Ex3 between the results obtained by the proposed model for such settings and the ones reported by Cafaro and Cerdá (2012). We also study the impact of backorder costs in Ex4. Since our model computes the cost for all interfaces that appear in the system rather than those generated during injection, and to have a fair comparison, we have added the interface cost from old and empty batches to the values in Cafaro and Cerdá (2012). 6.1. Example 1 It can be seen from Table 5 that the optimal solution is found for |T| = 8, two slots more than in Cafaro and Cerdá (2012). Our schedule corresponds to a total cost of $1,148,020, 1.19% lower than the previously reported value. The success is mainly due to the $19,0 0 0 reduction in interface cost, reflecting a lower number of

interfaces. Comparing the optimal schedule in Fig. 7 with the one reported in Cafaro and Cerdá (2012), the difference can be observed in the product sequences of lines L3 (P1-P2-P1 vs P1-P2-P1P2) and L4 (P2-P1-P3-P2-P4-P2 vs P2-P1-P3-P2-P4-P1). In particular, the elimination of one of the P1-P2 interfaces in L3, accounts for $18,400. As for the last interface in L4, our P4(I5)-P2(I6) choice does lead to a $(41,30 0–23,50 0) higher cost compared to the interface P4(I5)-P1(I6) but then, it has the advantage of not incurring on additional cost after fully delivering P4(I5) to N6 (L4 ends full with P2 instead of P2-P1). The change in pumping sequence also alters the final products in lines L3 and L4. As a consequence, the quantities of P1 and P2 transferred from N1 to N4 change from 79,0 0 0 m3 to 95,0 0 0 m3 and from 60,0 0 0 m3 to 44,0 0 0 m3 , respectively. From Table 3, the unit pumping cost for P2 is $(1.08–1.05) higher than that for P1, which explains the $(1.08–1.05) × 16,0 0 0 = 480 reduction in pumping costs in Table 5. In the inventory profiles for N4 in Fig. 8, note that all products reach their lower limits towards the end of the time

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568 Table 5 Model sizes and computational results for Ex1.

|T| |Il |a Equations Continuous variables Binary variables CPUs Total cost (102 $) Interface cost (102 $) Pumping cost (102 $) Inventory cost (102 $) Makespan (h)

This work

Cafaro and Cerdá (2012)

8 5,3,3,6 5369 3503 1300 651 11480.2 3686 5417.7 586.5 179

6 5,3,4,6 3464 1785 330 1607 11618.7b 3876b 5422.5 530.2 179

15

steps. Consider the latter case as an example, where 17,0 0 0 m3 of P1 enter line L4. One could dispatch 13,0 0 0 m3 of P4 to N6, then dispatch 4,0 0 0 m3 of P2 to N7, both at feasible rates, stopping L4 for a minimum of 7.75 h. The issue is that dividing one aggregate operation into three, changes the number of event points in the continuous-time grid as well as the average inventory level, affecting the inventory cost (recall Eq. (52)). In other words, the inventory cost reported in Cafaro and Cerdá (2012) is not as accurate. As shown in Fig. 7, a total of 179,0 0 0 m3 products are pumped by N1 at the maximum rate of 10 0 0 m3 /h, leading to the same 179 h makespan, as in the schedule of Cafaro and Cerdá (2012). 6.2. Example 2

a

Values of |Il | in each line are separated by commas. b Interface cost caused by old and empty batches, equal to $1,82,900, added to the values reported in Cafaro and Cerdá (2012).

Table 6 Model sizes and computational results for Ex2.

|T| |Il |a Eqs. Continuous variables Binary variables CPUs Total cost (102 $) Interface cost (102 $) Pumping cost (102 $) Inventory cost (102 $) Makespan (h)

This work

Cafaro and Cerdá (2012)

5 6,2,2,7,2 3645 2268 820 60.0 10427.2 3351 5019.8 556.4 150

6 6,2,2,6,2 3160 1591 390 59.3 10601.4b 3529b 5020.5 551.9 150

a

Values of |Il | in each line are separated by commas. Interface cost caused by old and empty batches, equal to $2,01,300, added to the values reported in Cafaro and Cerdá (2012). b

Table 7 Model sizes and computational results for Ex3.

|T| |Il |a Eqs. Continuous variables Binary variables CPUs Total cost (102 $) Interface cost (102 $) Pumping cost (102 $) Inventory cost (102 $) Makespan (h)

This work

Cafaro and Cerdá (2012)

6 6,4,2,7,2,1 5082 3143 560 579 11457.0 4177 5263.9 516.1 150

6 6,4,2,6,2,1 3606 1779 360 174 11659.0b 4355b 5265 539.0 150

a

Values of |Il | in each line are separated by commas. Interface cost caused by old and empty batches, equal to $2,01,300, added to the values reported in Cafaro and Cerdá (2012). b

horizon, except P3, which returns to its initial value (10,0 0 0 m3 ). This is because the initial volume of P3 inside the system equals the 25,0 0 0 m3 demand at N7 (check Table 1). Although we have the same final inventory as Cafaro and Cerdá (2012), our inventory cost is a bit worse. This is because our model generates a detailed instead of an approximate schedule. As explained in Castro (2017), the schedule represented in Fig. 7 in Cafaro and Cerdá (2012) sometimes violates the lower flowrate limit (800 m3 /h). This is the case of line L3 during interval [133, 150] h (700 m3 /h) and line L4 during intervals [0, 60] h (166.7 m3 /h), [105, 133] h (714.3 m3 /h) and [150, 179] h (586.2 m3 /h). It might be that the approximate schedule in Cafaro and Cerdá (2012) is actually implemented in multiple

From the optimal schedule for Ex1, we know that operation is limited by the pumping rate in line L1, connecting refinery N1 to node N4. To alleviate the load of line L1, Ex2 features an alternative supply line (L5) from N1 to N5. Not surprisingly, the makespan is reduced to 150 h (Table 6). Apart from the 350 0 0 m3 of P1 initially inside line L3, the remaining demand for P1 at N5 is supplied by line L5 at the maximum rate of 40 0 m3 /h, i.e. (80,0 0 0– 35,0 0 0)=45,0 0 0 m3 . Consequently, there is no need to carry P1 through a much longer route (lines L1 and L3), leading to substantial savings in pumping and interface costs compared to Ex1. The optimal schedule in Fig. 9 corresponds to a total operating cost of $1,042,720, 1.64% below the results from Cafaro and Cerdá (2012). It uses one less time slot, being obtained in about the same computational time. The only key performance indicator that is slightly worse in our schedule is the inventory cost. Similarly to Ex1, the inventory of products P1, P2 and P4 at node N4 reaches the minimum levels towards the end (see Fig. 10), indirectly indicating that the pumping costs at refineries N1-N2 have been reduced as much as possible. As for the reduction in interface cost, it is due to line L4. By observing Fig. 10 in Cafaro and Cerdá (2012), we find that three new batches I4(P2)-I5(P4)-I6(P1) are successively injected to the right of P3(I3) and a new interface P2(I4)-P1(I6) is also generated during the last interval. Interestingly, by allowing one more batch in line L4, our I4(P2)-I5(P1)-I6(P4)-I7(P1) pumping sequence in Fig. 9 creates no new interface after fully dispatching batch I6(P4) to N6. In spite of the (1840 0+2350 0–4130 0)=$60 0 caused by the additional batch I5(P1), our schedule avoids paying the extra $18400 for the P2-P1 interface appearing in Cafaro and Cerdá (2012). 6.3. Example 3 Ex3 adds an additional line, L6, connecting N4 with a new depot, N8. From Fig. 11, one can see that L6 is filled with 120 0 0 m3 of P4 at the start, carrying 160 0 0 m3 of P4 over the planning horizon. Despite the new transportation task, the makespan is not affected compared to Ex2. However, both the pumping and interface costs increase. The latter is not due to L6, since its interface cost is zero (single product carried). It can be explained by the additional injection of P4 from node N2 into line L2 (recall that in the optimal solution for Ex2, the stock of P4 at N4 is depleted to the minimum, so P4 needs to come from somewhere else). We also notice another batch of P2 (I4) being pumped after P4. Although this operation raises the interface cost (P4-P2), it is necessary since the initial stock of P4 at N2 (20,0 0 0 m3 ) is much less than the capacity of line L2 (50,0 0 0 m3 ). Thus, another product needs to be injected to push enough P4 out of line L2. It is also important to note that the inventory at N8 can always be kept at the minimum level. This is because the maximum rate to N8, 400 m3 /h, is lower than the one to the local market, 700 m3 /h. The decrease in total cost compared to Cafaro and Cerdá (2012) is now equal to 1.73%. From Table 7, our schedule is

16

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568 Table 8 Model sizes and computational results for Ex4. |T|

|Il |a

Equations

Cont. variables

Binary variables

CPUs

5 Total cost (102 $) 16367.7

6,4,2,7,3,1 Interface cost (102 $) 4879

3942 Pumping cost (102 $) 5346

2500 Inventory cost (102 $) 567.7

894 Makespan (h) 157.5

2147 Backorder (m3 ) 2000

a

Values of |Il | in each line are separated by commas.

obtained in 579 CPUs and performs better in all indicators, even inventory cost. The main contributor to total cost reduction is still the interface cost, being the reason for the sharp decline similar to the one previously discussed.

posed model by accounting for time- and volume-dependent inventory cost.

6.4. Example 4

Financial support from the China Scholarship Council (201806440102) during a visit of Qi Liao to the University of Lisbon, from the National Natural Science Foundation of China, Grant no. 51874325, and from Fundação para a Ciência e Tecnologia through projects CEECIND/00730/2017 and UID/MAT/04561/2019.

From Table 8, we can see that by using |T| = 5 and |Il | = (6,4,2,7,3,1), we are able to improve the best-reported solution in Table 4 for Ex4, to $1,636,770. Still, there is a significant increase in total operation cost ($491,070) compared to Ex3, mainly due to the 20 0 0 m3 of backorder of P1 at N5. Note that for Ex3, the final inventories at the refinery nodes indicate a margin of 30 0 0 m3 of P2 (N1) and 50 0 0 m3 of P3 (N4) and so, we knew that increasing the demand by 10 0 0 0 m3 (Ex4), would lead to at least 20 0 0 m3 of product backorder. Since the unit backorder cost is much higher than the pumping cost, the optimization decides to inject all available products into the system. Compared to Fig. 11, node N1 in Fig. 12 pumps one batch of P2 into line L5, to push the additional 30 0 0 m3 of P1 to N5 at a rate of 400 m3 /h, thus prolonging the makespan by 7.5 h. There is also an important growth in the interface cost, which is caused by the change in the pumping sequence of line L4. Due to the additional demand of P2 at N5, the transported volume of P2 in line L3 is increased by 50 0 0 m3 , a volume that in Ex3 went to line L4 as a filler batch. The alternative to fill line L4 is for node N4 to inject its only surplus product, P3, thus creating a new interface. 7. Conclusion This paper has presented a batch-centric MILP formulation for the detailed scheduling of a generic pipeline network, aiming to minimize the operating cost related to pumps, contamination, inventory, backorder and pipeline utilization. This formulation can be viewed as an extension of our recent work that allows the injection/delivery of multiple batches by/to a refinery/depot node over a slot, leading to a reduction in the computational effort required to solve the resulting mathematical problems. Other novel aspects are the rigorous constraints for enforcing forbidden product sequences when defining batch numbering for every line and a more accurate definition of the interface cost between consecutive products inside the pipeline. The superiority of our model has been demonstrated by solving four benchmark problems involving a complex mesh-structure pipeline network, featuring a central triple-purpose node connecting four or five lines, depending on the instance. Compared with the work of Cafaro and Cerdá (2012), our model can now generate detailed schedules, accurately compute interface costs, and is able to reduce operating cost by over 1.1%. Problem size is larger and the computational time is similar. Inventory cost has been computed as an average inventory level over all event points in the grid. This is a rough estimation because slot duration is not taken into account and because the inventory may change at a variable rate during a slot, whenever multiple injections/deliveries take place. Future work will improve the pro-

Acknowledgements

References Boschetto, S.N., Magatão, L., Brondani, W.M., Neves, F., Arruda, L.V.R., Barbosa-Poévoa, A.P.F.D., Relvas, S., 2010. An operational scheduling model to product distribution through a pipeline network. Ind. Eng. Chem. Res. 49, 5661–5682. Cafaro, D.C., Cerdá, J., 2012. Rigorous scheduling of mesh-structure refined petroleum pipeline networks. Comput. Chem. Eng. 38, 185–203. Cafaro, V.G., Cafaro, D.C., Méndez, C.A., Cerdá, J., 2015a. Optimization model for the detailed scheduling of multi-source pipelines. Comput. Ind. Eng. 88, 395–409. Cafaro, V.G., Cafaro, D.C., Méndez, C.A., Cerdá, J., 2015b. MINLP model for the detailed scheduling of refined products pipelines with flow rate dependent pumping costs. Comput. Chem. Eng. 72, 210–221. Cafaro, V.G., Cafaro, D.C., Méndez, C.A., Cerdá, J., 2012. Detailed scheduling of single-source pipelines with simultaneous deliveries to multiple offtake stations. In: Industrial and Engineering Chemistry Research, pp. 6145–6165. Castro, P.M., 2017. Optimal scheduling of multiproduct pipelines in networks with reversible flow. Ind. Eng. Chem. Res. 56, 9638–9656. Castro, P.M., 2016. Source-based discrete and continuous-time formulations for the crude oil pooling problem. Comput. Chem. Eng. 93, 382–401. Castro, P.M., Mostafaei, H., 2019. Batch-centric scheduling formulation for treelike pipeline systems with forbidden product sequences. Comput. Chem. Eng. 122, 2–18. Castro, P.M., Mostafaei, H., 2017. Product-centric continuous-time formulation for pipeline scheduling. Comput. Chem. Eng. 104, 283–295. Chen, H., Zuo, L., Wu, C., Li, Q., 2019. An MILP formulation for optimizing detailed schedules of a multiproduct pipeline network. Transp. Res. Part E Logist. Transp. Rev. 123, 142–164. Chen, H., Zuo, L., Wu, C., Wang, L., Diao, F., Chen, J., Huang, Y., 2017. Optimizing detailed schedules of a multiproduct pipeline by a monolithic MILP formulation. J. Pet. Sci. Eng. 159, 148–163. Commission National Development and Reform, 2017. Long-term oil and gas pipeline network planning. De Souza Filho, E.M., Bahiense, L., Ferreira Filho, V.J.M., 2013. Scheduling a multi-product pipeline network. Comput. Chem. Eng. 53, 55–69. Dimas, D., Murata, V.V., Neiro, S.M.S., Relvas, S., Barbosa-Póvoa, A.P., 2018. Multiproduct pipeline scheduling integrating for inbound and outbound inventory management. Comput. Chem. Eng. 115, 377–396. Hamisu, A.A., Kabantiok, S., Wang, M., 2013. Refinery scheduling of crude oil unloading with tank inventory management. Comput. Chem. Eng. 55, 134–147. Herrán, A., de la Cruz, J.M., de Andrés, B., 2012. Global search metaheuristics for planning transportation of multiple petroleum products in a multi-pipeline system. Comput. Chem. Eng. 37, 248–261. Herrán, A., de la Cruz, J.M., de Andrés, B., 2010. A mathematical model for planning transportation of multiple petroleum products in a multi-pipeline system. Comput. Chem. Eng. 34, 401–413. Lee, H., Pinto, J.M., Grossmann, I.E., Park, S., 1996. Mixed-integer linear programming model for refinery short-term scheduling of crude oil unloading with inventory management. Ind. Eng. Chem. Res. 35, 1630–1641. Liao, Q., Castro, P.M., Liang, Y., Zhang, H., 2019a. Batch-centric model for scheduling straight multi-source pipelines. AIChE J. 65, e16712. https://doi.org/10.1002/aic. 16712. Liao, Q., Castro, P.M., Liang, Y., Zhang, H., 2019b. Computationally efficient milp model for scheduling a branched multiproduct pipeline system. Ind. Eng. Chem. Res. 58, 5236–5251. Liao, Q., Liang, Y., Xu, N., Zhang, H., Wang, J., Zhou, X., 2018a. An MILP approach for detailed scheduling of multi-product pipeline in pressure control mode. Chem. Eng. Res. Des. 136, 620–637.

Q. Liao, P.M. Castro and Y. Liang et al. / Computers and Chemical Engineering 130 (2019) 106568 Liao, Q., Zhang, H., Xu, N., Liang, Y., Wang, J., 2018b. 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. Lopes, T.M.T., Ciré, A.A., De Souza, C.C., Moura, A.V., 2010. A hybrid model for a multiproduct pipeline planning and scheduling problem. In: Constraints, pp. 151–189. Magatão, S.N.B., Magatão, L., Luis Polli, H., Neves, F., De Arruda, L.V.R., Relvas, S., Barbosa-Póvoa, A.P.F.D., 2012. Planning and sequencing product distribution in a real-world pipeline network: an MILP decomposition approach. Ind. Eng. Chem. Res. 51, 4591–4609. Magatão, S.N.B., Magatão, L., Neves, F., Arruda, L.V.R., 2015. Novel MILP decomposition approach for scheduling product distribution through a pipeline network. Ind. Eng. Chem. Res. 54, 5077–5095. MirHassani, S.A., Fani Jahromi, H., 2011. Scheduling multi-product tree-structure pipelines. Comput. Chem. Eng. 35, 165–176. Moradi, S., MirHassani, S.A., 2014. Transportation planning for petroleum products and integrated inventory management. Appl. Math. Model. 39, 7630–7642. Mostafaei, H., Castro, P.M., 2017. Continuous-time scheduling formulation for straight pipelines. AIChE J. 63, 1923–1936. 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. 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.

17

Mostafaei, H., Ghaffari Hadigheh, A., 2014. A general modeling framework for the long-term scheduling of multiproduct pipelines with delivery constraints. Ind. Eng. Chem. Res. 53, 7029–7042. Moura, A.V., De Souza, C.C., Cire, A.A., Lopes, T.M.T., 2008a. Heuristics and constraint programming hybridizations for a real pipeline planning and scheduling problem. In: Proceedings - 2008 IEEE 11th International Conference on Computational Science and Engineering, CSE 2008, pp. 455–462. Moura, A.V., De Souza, C.C., Cire, A.A., Lopes, T.M.T., 2008b. Planning and scheduling the operation of a very large oil pipeline network. In: Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), pp. 36–51. Panda, D., Ramteke, M., 2019. Preventive crude oil scheduling under demand uncertainty using structure adapted genetic algorithm. Appl. Energy 235, 68–82. Relvas, S., Barbosa-Póvoa, 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., Boschetto Magatão, S.N., Barbosa-Póvoa, A.P.F.D., Neves, F., 2013. Integrated scheduling and inventory management of an oil products distribution system. Omega 41, 955–968. 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., 2017. A hybrid computational approach for detailed scheduling of products in a pipeline with multiple pump stations. Energy 119, 612–628.