Transportation planning for petroleum products and integrated inventory management

Transportation planning for petroleum products and integrated inventory management

Accepted Manuscript Transportation planning of petroleum products and integrated inventory management S. Moradi, S.A. MirHassani PII: DOI: Reference: ...

469KB Sizes 0 Downloads 83 Views

Accepted Manuscript Transportation planning of petroleum products and integrated inventory management S. Moradi, S.A. MirHassani PII: DOI: Reference:

S0307-904X(15)00284-X http://dx.doi.org/10.1016/j.apm.2015.04.023 APM 10564

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

4 June 2014 8 March 2015 17 April 2015

Please cite this article as: S. Moradi, S.A. MirHassani, Transportation planning of petroleum products and integrated inventory management, Appl. Math. Modelling (2015), doi: http://dx.doi.org/10.1016/j.apm.2015.04.023

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Transportation planning of petroleum products and integrated inventory management S. Moradi and S. A. MirHassani∗ Faculty of Mathematics and Computer Science, Amirkabir University of Technology, Tehran, Iran

Abstract This paper addresses the scheduling of an oil derivatives transportation system by a multiproduct pipeline that connects a single refinery to a distribution center (DC) and meets the customer demands on a daily basis considering integrated inventory management at DC and refinery. A batch sequencing is presented which leads to a mixed integer linear programming (MILP) formulation in order to determine the optimal pipeline batch sequence and timing of product deliveries to the DC. Some operational restrictions such as continuous or discrete flow rate and high peak electricity periods are considered. Experimental results for several examples including a real-world multiproduct pipeline, considering different scenarios, show that the proposed model leads to an optimal pipeline schedule with less computational cost.

Keywords: Pipeline scheduling, integrated inventory management, peak hours, continuous and discrete flow rate

∗ Corresponding author: [email protected].

1

1. Introduction Planning products transfer from industrial units to distribution centers (DCs) is an important task, which often becomes a complicated problem. In the oil industry, products are transferred by different means. However, the use of pipelines is affordable and has a smaller impact on the environment and traffic than other means. Therefore, high percentage of petroleum products is typically transported by pipelines from refineries or ports to DCs and then is delivered to customers by trucks. Although pipeline is the most economical mode of transferring petroleum products, its scheduling is a complex managerial problem in the oil industry. A simple pipeline consists of a unidirectional line. In this system, a source is connected to some DCs. Different batches of oil derivatives, such as gasoline, gas oil, and kerosene, are pumped into the pipeline and discharged at DCs successively. At the beginning of each period, a planner takes into account some restrictions such as maximum and minimum flow rate, capacity of pipeline, storage tanks in source and destination, and time between loading and unloading of a batch. Then, s/he should decide on the sequence of injected batches, their volume, start and completion time of each injection, and discharging. All the related decisions should be made in a way that the customer demand is satisfied with minimum delay, minimum cost of pumping, and least interference between the injected products. In these systems, however, the most important decision is how to determine the sequence and schedule of injection and discharging in batches. Background of the problem: In the last two decades, many researchers have studied this problem from different points of view and most of them have provided MILP models as well as exact or heuristic algorithms for solving this problem. In some models, time and volume are considered as discrete variables ( [1]), while in others, they are treated as continuous ones ( [2]). In addition, in the continuous models, there are two different methods of batch sequencing. The first is based on the injecting time of batches ( [2]), whereas in the second method, their discharging time is considered ( [3]). Both methods lead to models with different structures and sizes. Rejowski and Pinto [1] studied the scheduling of a real pipeline system consisting of a multiproduct pipeline connecting one source to multiple DCs. In their study, a discrete MILP model minimized the cost of inventory, pumping, interface, and charge of high peak electricity periods. The process of obtaining the model's optimum solution was time consuming; therefore, adding some efficient cuts could somewhat improve it [4]. Magatao et al. [5] presented another discrete model for scheduling a multiproduct pipeline in a limited period, which was divided into two parts and solved in two steps. Cafaro and Cerda [2] provided the first MILP continuous formulation for scheduling a pipeline with one source and multiple outputs, whose number of binary variables and running time was less than that of discrete models. MirHassani and Ghorbanalizadeh [6] as well as MirHassani and Jahromi [7] presented discrete and continuous models for scheduling tree structure pipelines, respectively. Inventory management in DC, daily due dates, and settling period for quality 2

control were considered by Relvas et al. [8] and [9]. Later, Cafaro and Cerda [10] improved the size, running time, and performance of these models. In most of the above cases, the running time of the presented models for long-term scheduling is high; therefore, some heuristic methods ( [11], [12]) and algorithms ( [13]) are presented to address this shortcoming. MirHassani et al. [13] generalized Cafaro's model to the pipeline system with one source and multiple DCs by considering the daily-based demand. Also, MirHassani et al. [14] presented an MILP model for scheduling a unidirectional pipeline with one source and multiple DCs, some of which were able to receive and inject the products. Cafaro et al. [15] provided a new model for scheduling a pipeline with single source and many DCs, which could receive batches simultaneously. Recently, the pipeline systems with more complexity have been studied. Boschetto et al. [16] studied a real-world pipeline network and proposed an optimization model in order to address the network scheduling activities. Lopes et al. [17] presented a compact network flow model for scheduling a network of pipelines whose nodes were source, destination, or junction. Relvas et al. [3] have recently proposed a model for scheduling a pipeline connecting the refinery to a distribution center. Based on the volume of the batches, they offered two separate models. In the first model, the volume of each batch, based on the product type, was determined among several options (discrete volume). However, in the second one, the volume of each batch, based on its product type, was determined continuously within an interval (continuous volume). The main variable for scheduling in the previous models was based on starting and finishing the time of batch injection; but, in this new model, the completion time of discharging batches was the main variable. The method led to a more efficient formulation and reduced the number of variables and constraints as well as model solution time. Another new aspect of this work was a preprocessing method for removing some binary variables from the search space, which could reduce the CPU time for long-term scheduling. However, there was a trade-off between tighter bounds and CPU time. Cafaro and Cerda [18] recently proposed a model for the shortterm planning of a reversible multiproduct pipeline, whereby one DC received products from opposite directions. Based on Relvas' model [3], this paper addresses the problem of scheduling a single transmission pipeline carrying refined products from an oil refinery to a unique DC to meet the customer demands on a daily basis. In this paper, discharging process of batches is obtained, similar to Relvas' model [3] and injecting time of the lots is calculated. In addition, some important features, like peak electricity periods, continuous and discrete flow rate, and integrated inventory management at DC and refinery are considered. The new features, which are considered in this paper, are four dimensions: Initially, providing a new attitude for batch sequencing by considering some batches for each day, instead of total time horizon. This point reduces the number of binary variables and constraints and could reach the optimal solution in a reasonable time without preprocessing (fixing or deleting some variables before running the model). In addition, by increasing the number of total days, the 3

size of model does not considerably increase. Second, the additional cost of pumping in high peak electricity periods is added. Third, discrete and continuous flow rates are considered separately. Fourth, inventory management at both refinery and DC is considered simultaneously. The rest of this paper is organized as follows: Problem definition and assumptions are described in Section 2. The problem formulation is described in Section 3. Finally, in Section 4, the computational experiments and analysis of the results are presented. 2. Problem definition and assumptions This problem deals with scheduling a straight pipeline that connects a refinery to a DC. There are various tanks for each product at refinery and DC, each product is stored in its own storage tank, and DC must meet clients' daily demand. 2.1 Consider 1) Number of derivatives; 2) Time horizon including some days and last hour of each day; 3) Maximum number of batch injection per day; 4) Daily demand of clients; 5) Maximum and minimum allowed volume for any discharged batch considering type of product; 6) Initial conditions of pipeline and inventory level at DC and refinery; 7) Volume of interfaces between each pair of products and list of forbidden sequences; 8) Number of peak hours per day; 9) Allowable interval of flow rate (continuous flow rate) or available options for flow rate (discrete flow rate); 10) Settling time for each product (daily based); 11) Capacity of each tank at DC and refinery; and 12) Daily input of each product to the refinery tanks. On the other hand, the following should be obtained: 1) Sequence of product batches to be injected at refinery and discharged at DC; 2) Batch sizes and the start/end discharging/injecting times; and 3) Product inventory management at refinery and DC so that capacity restrictions are not violated and daily demands are provided.

2.2 Model assumptions The following assumptions are considered: • • • • •

The pipeline is always full and volumes of injection and discharging are equal at any time. Product batches are injected in a row in the pipeline and discharged at DC Interface volume between each pair of products, regardless of batch size, is a known constant. Sum of storage capacity of all tanks, assigned to each product, at DC or refinery is considered as an aggregated tank. Each batch should take a settling time before entering the storage tank at DC.

4

• • • • • • •

If the settling time of a batch is not completed, it is not used to supply the customer demand. Daily product demands are given at the beginning of the time horizon. The flow rate is continuous which is at a determined interval or constant and is selected from among a few options. The volume and discharging time of old batches, which are already inside the pipeline, are known. The type and schedule of some batches that are injected on the last day and remain in the pipeline at the end of the time horizon are not considered. On some days, a few hours are considered as peak hours and cost of electricity energy during these hours is greater than that of others. For refinery inventory management, daily production of each product is assumed to be determined. In addition, when the injection of a batch starts, total volume of this batch is subtracted from the inventory at refinery.

Nomenclature Sets ,  ∈ : Discharging batches during a day (, ) ∈ ⊆  × : Old batches that are injected in the previous scheduling period

∈ : Discrete flow rate options ,  ∈ : Refined petroleum products ,  ∈ : Daily periods of time horizon Parameters  : Last hour of time period t , : Demand of product p on day t  : Initial inventory level of product p at DC  − ,, : Volume of initial batches ((, ) ∈ )  − , : Discharging time of initial batches ((, ) ∈ )  !"#$%" ,  !"#$%&' : Minimum and maximum allowed storage capacity of clients for product p at DC (& $%" , (& $%&' : Minimum and maximum allowed total storage capacity for product p at DC , )* : Minimum and maximum allowed volume for every discharged batch of product p +: A big number )  : Peak hours of time period t -), : Produced volume of product p on day t /: Maximum allowable interference 5

0 : Initial inventory level of product p at refinery 0%" , 0%&' : Minimum and maximum allowed storage capacity for product p at refinery 1 : Duration of settling period for a batch of product p 23 : The kth option for discrete flow rate 2%" , 2%&' : Minimum and maximum allowed continuous flow rate : Pipeline volume /)1,4 : Contamination volume between consecutive batches of products p and q Continuous variables 52,, : Cumulative volume of a sequence of successive batches of product p leading to the ith batch on day t 6, : Shortage volume of product p at DC to meet demand of day t ,, : Volumetric distance between the ith batch of day t and first batch of day tt 77: Difference between the total discharged volume and total demand of all the products 7, : Contaminated volume of the ith batch of day t ,,, : Size of the ith discharged batch of day t including product p injected on day tt  !"# , : Inventory level of product p at DC on day t ready to supply customer demand (& , : Total inventory level of product p at DC on day t : Lowest total inventory among all products on the last day ℎ : Used peak hours of day t 9,,3 : Size of the ith batch of day t discharged at flow rate 23 0, : Inventory level of product p at refinery on day t :# , : Discharging time of the ith batch of day t ;,, : Size of the ith batch on day t including product p

Binary variables 7-,,3 = 1; then, the ith batch of day t is discharged at flow rate 23 . *,, = 1 if injection of the ith batch of day t starts on day tt; otherwise, 0. >,, = 1 if the ith batch of day t contains product p; otherwise, 0. 3. Mathematical formulation To formulate the pipeline scheduling, batches and days are shown by the indices i and t, respectively, and sorted in an ascending order. In this model, the set I denotes the batches which could be discharged during a single day and | | is independent from the number of days. By increasing the total number of days (| |), the value of | | remains constant and size of the model does not change considerably. In contrast, in the previous similar work ([3]), the set I denotes the batches that could be discharged during time horizon and, by increasing the total number of days, 6

the value of | | increases proportionately. For example, if time horizon includes 15 or 31 days, in the present model, in both cases, | | is considered as 3 batches, while in the previous work, for 15 and 31 days, | | is considered as 35 and 71 batches, respectively (see more details in Section 4.1). To obtain the proper number for the cardinality of I, similar to [14], we can start with an initial guess (provided by historical data) and then increase it to a state, in which no better solution is obtainable. Few groups of constraints are used to formulate the problem: 3.1. Allocating products to the batches For every day, a new discharged batch is one of the products or is empty. @ >,, B 1, ∀ ∈ ,  ∈ (1)

∈A

The empty batches should be left at the end of the sequence of batches for every day, i.e. if one of the considered batches contains no product, all the subsequence batches of the day will be considered empty. @ >DE,, B @ >,, , ∀ ∈ ,  F | |,  ∈ (2)

∈A

∈A

Also, if on day t, no batch is to be discharged (∑ >IE,, = 0), on the subsequent days, no batches will be discharged. @ >,DE, B @ >,, , ∀ ∈ ,  = 1,  ∈ ,  F | | (3)

∈A

∈A

Daily batch sequencing is shown in Fig. 1. Based on Fig. 1, on the first day, two batches including P2 and P3 with 12 and 6 volumetric unit sizes are discharged, respectively. On the second day, two batches including P4 and P1 with 8 and 9 volumetric unit sizes are discharged, respectively, and on the third day, just one batch including P5 with 18 volumetric unit sizes is discharged.

7

Fig. 1: An example of batch sequencing

3.2. Lot-sizing

Size of a batch of product p should be in a specific range ([ , )* ]). In addition, the volume of fictitious batches is zero. Therefore: >,, ×  B ;,, B >,, × )* , ∀ ∈ ,  ∈ ,  ∈  (4)

If we consider only the above constraint, it is possible for some consecutive nonempty batches to include product p and constitute a single lot of p; then, maximum batch size may be violated. Suppose 52,, is the cumulative volume of a sequence of successive batches of product p that leads to the ith batch on day t. This variable is determined by Eqs. (5) and (6) and should be lower than the maximum batch size. 52,, ≥ 52$E,, + ;,, − )* O @ >,,4 R, ∀ ∈ ,  > 1,  ∈ ,  ∈  (5) 4∈A$PQ

52,, ≥ 52U,$E, + ;,, − )* O @ >,,4 R, ∀,  ∈ ,  = 1,  = ||,  ∈ ,  > 1,  ∈  (6) 4∈A$PQ

52,, B )* , ∀ ∈ ,  ∈ ,  ∈  (7)

The nonnegative variable diff is defined equal to the difference between total discharged volume and total customer demand in time horizon and is minimized in the objective function. 77 = @ @ @ ;,, − @ @ , (8) ∈Y ∈X ∈A

∈X ∈A

3.3. Product inventories at DC At the end of each day, the inventory level of each product at DC is calculated. Each batch should go through settling time, for quality control, before entering the storage tank at DC (one or more days) and these batches will be available for client demands after settling periods. Therefore, there are two types of inventory level: total inventory level and available inventory for client demands. (& The total inventory level of product p at the end of day t (, ) is equal to the sum of inventory of the prior day and any inventory discharged during day t (∑ ;,, ), which backorders the prior day (6$E, ) and the demand of day t (, ) will be excluded. In

8

addition, it is possible for DC to be faced with inventory shortage as much as 6, . Thus, for the first day, we have: (& , =  + @ ;,, − , + 6," , ∀ ∈ ,  = 1,  ∈  (9) ∈Y

For other days, we have: (& (& , = $E, + @ ;,, − , − 6$E," + 6," , ∀ ∈ ,  > 1,  ∈  (10) ∈Y

It is supposed that the settling period is daily-based; therefore, if the settling period of product p is equal to 1 days, each batch of this product will be available on 1 days after discharging at DC. The available inventory level of product p on day t for meeting the client  !"# demand (, ) is equal to the sum of the available inventory level of product p from the prior day and volume of product p that is available on day t (∑ ;,$#!\ , ). Thus, for the first day, we have:

 !"# , =  + @ ;,, − , + 6," , ∀ = 1,  =  − 1(),  ∈  (11)  ∈

And for other days, we have:

 !"#  !"# , = $E, + @ ;,, − , − 6$E," + 6," , ∀ > 1,  =  − 1(),  ∈  (12)  ∈

The total and available inventory levels of each product on each day should be in a specified range: (&

$%"

(& B , B (&

$%&'

, ∀ ∈ ,  ∈  (13)

 !"#  !"#$%" B , B  !"#$%&' , ∀ ∈ ,  ∈  (14)

3.4. Interface volume between consecutive batches In Relvas' model, [3] interface between consecutive batches is ignored and the optimal solution may include a sequence with high contamination that must be avoided. If the ith batch of day t contains product p and the next batch contains product q, the contaminated volume of the ith batch of day t is supposed to equal /)1,4 (see Fig. 2). Continuous variable 7, represents the amount of contaminated volume of the ith batch on day t and is obtained using the next three constraints:

9

Fig. 2: Interface material between two consecutive batches

On each day, a batch is mixed with the next batch. If batch i is not the last batch of day t, its interface can occur with the i+1th batch of day t. 7, ≥ /)1,4 ]>,, + >DE,,4 − 1^, ∀ ∈ ,  F | |,  ∈ , ,  ∈  (15)

If batch i is the last batch of day t, its interface can occur with the first batch of day t+1.

7, ≥ /)1,4 ]>,, + >U,DE,4 − 1^, ∀,  ∈ ,  = ||,  = 1,  ∈ ,  F | |, ,  ∈  (16)

On day t , if the next batch after batch i is a fictitious batch, interface of batch i occurs with the first batch of day t+1: 7, ≥ /)1,4 O>,, + >U,DE,4 − 1 − @ >DE,,́ R , ∀,  ∈ ,  F ||,  = 1,  ∈ , ,  ∈  (17) ́ ∈A

For example, in Fig. 1, the first batch of day t=3 is product P5, the second batch of day t=3 is empty, the first batch of day t=4 is product P3, and therefore batches (i=1, t=3) and (i=1, t=4) have interface. Then, the volume of the interface between those consecutive batches will never be less than /)1,4 . Such a parameter denotes the volume of the trans-mix between products p and q, just in case two consecutive batches contain products p and q, respectively. Otherwise, the constraint will become redundant. If the value of variable 7, is minimized in the objective function, its upper bound will not be extreme. If the sequence of p and q is a forbidden sequence, the value of /)1,4 should be considered as / + 1 = max Pwasteǵ h́ |sequence ṕ and q́ is permissibleQ + 1 and the following constraint would be added to the model. 7, B /, ∀ ∈ ,  ∈ (18)

3.5. Discharging time During each day, one or more batches are discharged at DC. Discharging duration of each batch depends on the flow rate. There are many scenarios for flow rate. Constant and continuous flow rates have been considered in previous works; however, in some real systems, the scheduler must choose the flow rate from among several options. For example, 10

in National Iranian Oil Products Distribution Company, a scheduler decides that one pump, with a constant flow rate, or two or more pumps should work. Modeling of discharging time constraints, based on the flow rate scenarios, is different. Eqs. (19) and (20) are related to the continuous flow rate, whereas Eqs. (21) – (25) are related to the discrete flow rate. First, suppose that flow rate is continuous and should belong to a specified interval ([2%" , 2%&' ]). Discharging the first batch of day t starts from the earliest hour of day ($E ) and finishes at time $E + (tu) × ∑ ;,, , in which 2 is discharging flow rate E

of batch and ∑ ;,, is volume of this batch; therefore, (tu) × ∑ ;,, is equal to the E

length of discharging. Since 2 should be satisfied in a specified range ([2%" , 2%&' ]), we have: $E + v

1 1 :# w @ ;,, B , B $E + v w @ ;,, , ∀ ∈ ,  = 1,  ∈ (19) 2%&' 2%"  ∈

∈

x! :# Discharging of the other batch starts after prior batch ( $E, ) and finishes at time $E, +

(tu) × ∑ ;,, . Thus: E

:# $E, +v

1 1 :# :# w @ ;,, B , B $E, +v w @ ;,, , ∀ ∈ ,  > 1,  ∈ (20) 2%&' 2%"  ∈

∈

Next, suppose that the flow rate is discrete and the scheduler decides that pumping operation should be done by which flow rate. In this case, the length of discharging batch of the ith batch on day t is equal to ytu { ∑ ;,, if its flow rate is 23 . For the ith batch on day t, suppose E

z

that if binary variable 7-,,3 is equal to 1, this batch is discharged by flow rate 23 and continuous variable 9,,3 is equal to the volume of this batch; otherwise, it is equal to zero. Thus, the discharging length of the ith batch on day t is equal to ∑3

|},~,z tuz

; so, one of the 9,,3 is

equal to the volume of batch and others are equal to zero. Thus, Eqs. (21) – (23) are added. @ ;,, − +]1 − 7-,,3 ^ B 9,,3 B +7-,,3 , ∀ ∈ ,  ∈ , ∈  (21) 

9,,3 B @ ;,, , ∀ ∈ ,  ∈ , ∈  (22) ∈A

@ 7-,,3 B 1 , ∀ ∈ ,  ∈ (23)

∈A

In addition, Eqs. (19) and (20) are replaced with the two following equations: 11

:# , = $E + @

3∈

:# :# , = $E, +@

3∈

9,,3 , ∀ ∈ ,  = 1,  ∈ (24) 23

9,,3 , ∀ ∈ ,  > 1,  ∈ (25) 23

3.6. Peak hours For some companies, the cost of electricity energy varies at peak hours and schedulers prefer to shut down the pipeline as long as possible during these hours. Since the inventory level is calculated at the last hour of the day, without compromising generality, it is assumed that the last hours of the day are peak hours. If peak hours are intermediate hours of the day, the obtained schedule can be shifted to the right time. For example, if the interval (7, 15) of (0, 24) is the peak :# hours, we move this period forward to the last 8 h of the day, i.e. (16, 24). If IE,E = 6,

:# :# I€,E = 12, I,E = 16, and ℎE = 0 are obtained from running the model, the real discharging time of these batches is equal to time 6, 12+8=20, and 16+8=24, respectively (discharging time of the batches prior to the peak period remains unchanged and others are moved forward).

Thus, discharging the last batch on day t should be finished before peak hours, unless the required pumping is performed during ℎ hours of the peak hours. :# , B  − )  + ℎ , ∀ ∈ ,  = | | ∈ (26)

However, the used peak hours are less than or equal to the total peak hours of day t and these variables are minimized in the objective function. ℎ B )  , ∀ ∈ (27)

It should be noted that, if a forbidden period for pumping exists on day t, by choosing appropriate weight for peak hours in the objective function, the value of ℎ will move to zero. 3.7. Inventory management at refinery Based on our knowledge, in the previous works, inventory management at refinery or DC has been disregarded and integrated inventory management at refinery and DC is not considered. However, in this paper, inventory management at refinery and DC is considered simultaneously. For calculating inventory level at refinery, it is necessary to determine the day on which each batch is injected. Thus, variable ,, is defined equal to volumetric distance between the ith discharged batch of day t and the first batch of day tt ( B ) and obtained from Eqs. (28) - (30). Volumetric distance between the i-1th discharged batch of day

12

t and first batch of day tt ($E,, ) plus the volume of this batch is equal to the volumetric distance between the next and first batches of day tt. ,, = 0 , ∀ ∈ ,  = 1,  ∈ (28)

,, = U,$E, + @ ;U,$E, , ∀,  ∈ ,  = 1,  = | |, ,  ∈ ,  F  (29) ∈A

,, = $E,, + @ ;$E, , , ∀ ∈ ,  > 1, ,  ∈ ,  B  (30) ∈A

If volumetric distance between the ith batch of day t and the first batch of day tt lies at interval [vp, vp+∑∈ ∑∈ ;,, ], then injecting this batch causes discharging some batches of day tt and therefore its discharging starts at day tt and binary variable *,, will equal 1. Eqs. (31) (33) provide these conditions. ,, ≥  × *,, , ∀ ∈ , ,  ∈ ,  B  (31)

,, B  + @ @ ;, , + +(1 − *,, ) , ∀ ∈ , ,  ∈ ,  B  (32) ∈Y ∈A

@

∈ ,B

*,, = 1 , ∀ ∈ ,  ∈ ((, ) ∉ ) (33)

For example, in Fig. 1, DE,E,E = 0, D€,E,E = 12, DE,€,E = 18, D€,€,E = 26, DE,,E = 35, D€,,E = 53, DE,„,E = 53; if pipeline volume is equal to 18 volumetric unit, then the injection of batches (1,2), (2,2), and (1,3) starts at day tt=1 (DE,€,E , D€,€,E , DE,,E ∈ …,  + ∑∈Y ∑∈A ;,E, † = [18, 36]) and these injections cause discharging the batches (1,1) and (2,1). If injecting of the ith batch on day t starts during day tt, the injected volume must be equal to ;,, ; otherwise, it must be equal to zero; thus: ,,, B +*,, , ∀ ∈ , ∀,  ∈ ,  B , ∀ ∈  (34)

;,, − +]1 − *,, ^ B ,,, B ;,, , ∀ ∈ , ∀,  ∈ ,  B , ∀ ∈  (35)

Inventory level of product p at refinery on day tt is equal to the sum of inventory level on the prior day plus the production volume of product p on day tt minus the injected volume. For calculating the inventory level at refinery, it is assumed that the daily production of each product is known (-), ). Also, when a batch injection starts on day tt, total volume of this batch is subtracted from the refinery inventory of day tt.

13

0, = 0 + @ -), − @ @ ,,, , ∀  ∈ ,  = 1 (36) 7∈ˆ

∈Y ∈X

0, = 0$E, + @ -), − @ @ , ,, , ∀  ∈ ,  > 1 (37) 7∈ˆ

∈Y ∈X

0%" B 0, B 0%&' , ∀  ∈ ,  ∈ (38)

3.8. Initial conditions At the first hour of the time horizon, the pipeline is full of the products that are injected into the previous scheduling period. It is assumed that the type and discharging time of these batches are given. Thus, the following constraints are added as the initial conditions: ;,, =  − ,, , ∀(, ) ∈ ,  ∈  (39) :# , =  − , , ∀(, ) ∈ ,  ∈  (40)

3.9. Objective function

The objective function may be operationally or economically oriented. The orientation in the most often published works is economic. However, in real-world scenarios, operational objectives are frequently used. In the present work, the objective function is composed of both types by normalizing the objective function terms, including: (a) Minimizing the differences between total discharged volume and total demands of customers: min ‰E =

77 (41) ∑∈X ∑∈A ,

(b) Maximizing the minimum inventory among all the products on the last day of time horizon () which should be obtained through a minmax strategy so as to maximize the inventory level of the product with the lowest final inventory level, as stated in Eq. (43): max ‰€ =  (42) (& |X|,  B , ∀ ∈  (43) %&' (c) Minimizing the underutilizing pipeline capacity within total hours of time horizon:

14

:#  − |Y|, 1 min ‰ = @( ) (44) | |  ∈X

(d) Minimizing the total used peak hours: min ‰„ =

1 ℎ @( ) (45) | | )  ∈X

(e) Minimizing the volume interfaces occurring between the injected batches: 1 min ‰Š = @ @ 7, (46) | || |max,4 (/)1,4 ) ∈Y ∈X

(f) Minimizing the delays in supplying the demand in due dates: 1 min ‰‹ = @ @ 6, (47) ||| |6%&' ∈A ∈X

Thus, the objective function is:

min ‰ = /E ‰E − /€ ‰€ + / ‰ + /„ ‰„ + /Š ‰Š + /‹ ‰‹ (48)

The first three components are operational and the last three ones are economical. Also, the terms ‰E and ‰ mean that the minimum volume of the products should be discharged during the maximum time. Therefore, the rate of flow is indirectly minimized. / is the weight of term l and indicates the importance and priority of this term in the objective function from planners point of view. For example, in some cases, at first, satisfying all demands by minimum delay is more important than other terms; second, discharging the minimum volume during the maximum time (minimizing the flow rate) is important. 4. Results and discussion To maintain results and comparison, a real-world case study introduced by Relvas et al. [3] was implemented here. This case study includes a multiproduct pipeline that transports six products from an oil refinery to a single DC. The batch size for each product can be selected from an interval [ , )* ] (see Table 1). The minimum settling period for each product and maximum storage capacity for every product are given in Table 1. This table also represents the volume of interface between each pair of products. Some sequences of products are not allowed inside the pipeline. Interface volume of such sequences is shown in Table 1; but, in the model, these parameters are set as: pw+1= max{wasteg,h |sequence p, q is permisable}+1=39. The minimum inventory capacity considered in all the case studies is zero. Table 1: Common data used in all the scenarios

invgŽ(vu)

minlot g (vu) maxlot g (vu) setlg (day)

15

Interface volume (vu)

P1 P2 P3 P4 P5 P6

81500 32000 24000 27800 10320 13120

8000 8000 8000 3800 800 3100

21800 16000 16000 11800 3200 6200

1 2 1 1 1 1

P1 0 30 37 35 − −

P2 30 0 − − − −

P3 37 − 0 38 30 −

P4 35 − 38 0 25 −

P5 − − 30 25 0 30

P6 − − − − 30 0

For the case of continuous flow rate, similar to Relvas et al. [3], the flow rate can span from 450 to 650 vu/h. Monthly delivery requirements and the initial inventory for every product at DC are shown in Table 2. It is assumed that, on each day, the maximum discharged batches are equal to 3 (| | = 3) and the daily demands are removed from the available inventories at the end of each day. The monthly demand reported in Table 2 is supposed to be uniformly distributed all over the month. Moreover, the pipeline is initially full of product P1 and the initial product is discharged with minimum flow rate at time

E‘’’’ „Š’

= 40. Therefore, it is

discharged as 2 batches at time 24 and 40 and ;E,E,E = 24 × 450 = 10800, ;€,E, = :# :# ;,E, = 0, ∀ ∈ , ;E,€,E ≥ 7200 and E,E = 24, E,€ ≥ 40 are fixed as initial conditions.

All the models were run with the stopping criterion of 120 min of computation or the determination of the optimal solution. In addition, the proposed mathematical model was solved all at once without any decomposition or preprocessing approach on an Intel 2 GHz processor and 2 GB RAM with AIMMS software version 3.12 and solver CPLEX 12.4 Table 2 : Monthly demand and initial inventory for each product P1

P2

P3

Monthly demand (vu)

198043

64800

14642

Initial inventory (vu)

52397

17565

18569

P4

68244 19888

P5

10934 10027

P6 16955 7309

4.1. Comparing the proposed model with the previous work In this section, our model is compared with the previous work (Relvas et. al. [3]). Since in this work, continuous batch volume was considered, in order to compare the performance of the pipeline scheduling approach with that of the previous work, the instances previously defined by Relvas et al. [3] were considered (the VBs model). In the first step, peak hours, discrete flow rate, and refinery inventory management were ignored (Eq. 21-39 of constraints and term ‰„ of objective function were ignored). Similar to Relvas et al. [3], three case studies were presented, in which the time horizon was the major difference between them. All the data of these scenarios were based on the real month of operations at CLC. So, the following three case studies were considered: Case study 1 (CS1) with the time horizon of 1 week (7 days); 16

Case study 2 (CS2) with the time horizon of 2 weeks (15 days); Case study 3 (CS3) with the time horizon of 1 month (31 days) The VBs model in Relvas et al.'s paper had just one set of integer variables (that were binary variables) and was shown by 1,,3 . This variable set was equivalent to >,, in the present model. In Relvas' model [3], for each scenario, the number of all products was equal to 6, the number of all batches was equal to 16, 35, and 71, and the number of binary variables, without preprocessing, was 16× 7 ×6=672, 35×15×6=3150, and 71×31×6=13206, respectively. Then, after preprocessing and deleting some integer variables from search space, the number of these variables was reduced to 156, 516, and 1890, respectively. However, in the present work, the cardinality of set I was equal to 3 batches and, for each scenario, without any preprocessing, the number of binary variables was equal to 3× 7 ×6=126, 3×15×6=270, and 3×31×6=558, respectively, representing that the number of binary variables was reduced compared with Relvas' model (see Table 3). In all the scenarios, the weight of the objective function terms was set equal to 5, 4, 1, 2, 0.5, and 6, respectively. Computational results shown in Table 3 include the final status of running models, number of constraints, continuous/binary variables in the problem formulation, required CPU time, number of nodes explored, and iterations. For the sake of comparison, Table 3 also presents the results reported by Relvas et al. [3] in two rows without and with preprocessing. In the present formulation, for all the scenarios and without preprocessing, the number of problem constraints dropped by the factor of 6.7 and the number of binary variables fell by 2.1 times. Moreover, the problem was solved to optimality and the solution time compared with Relvas et al. [3] diminished by 3.3 times. It is worth mentioning that, in the present model, the interface constraints were taken into account. Note that, in none of the case studies of Relvas et al., delay and backorders were considered; also, in the present model, , backorder was equal to 0 in all the case studies.

17

Table 3: Comparing performance of models for examples 1-3 Relvas et al.'s model (2013) Scenario CS1 Final status (no pre.pro) Optimal Final status (after pre.pro) Optimal Constraints (no pre.pro) 2228 Constraints (after pre.pro) 2230 Continuous variables 1447 Binary var. (no pre.pro) 672 Binary var. (after pre.pro) 156 CPU time (no pre.pro) 0.71 min CPU time (after pre.pro) <0.01min Nodes (no pre.pro) 34644 Nodes (after pre.pro) 82 Iteration (no pre.pro) 138156 Iteration (after pre.pro) 6545 Note: pre.pro =Preprocessing.

CS2 Optimal Optimal 8220 8222 6518 3150 515 30.42 min 0.57 min 91842 8464 21656045 1107044

CS3 No solution found Optimal 30316 30318 26858 13206 1890 120 min 2.17 min 630 1718 2100 1115185

Proposed model CS1 Optimal 2075 654 126 3s 524 4728 -

CS2 Optimal 4539 1398 270 7.4 s 861 38116 -

CS3 Optimal 9467 2886 558 47.3 s 3421 126467 -

The other important point in this work was the improvement of the pipeline operation and final inventory. To compare the utilization level of pipeline, medium flow rate, minimum final inventory, and minimum inventory along time horizon, the model result was compared with the one reported by Relvas et al. [3] for CS3. Table 4 shows that the pipeline utilization level was 100% in the proposed model and differences between the discharged volume and total demand was zero, which was similar to Relvas' model. However, the medium flow rate in the present model fell to 20.53 unit and the final inventory level increased by 21.02% for all the products and by 33.57% for the minimum final inventory (see Fig. 4). Discharging sequence and volumes for CS3 obtained from the new model is presented in Table 5.

18

5 4.5 4 3.5

P6

3

P5

2.5

P4

2

P3

1.5

P2

1

P1

0.5 0 1

3

5

7

9

11 13 15 17 19 21 23 25 27 29 31

Fig.3: Daily inventory level (%) for CS3 obtained from the proposed model

Table 5: Discharged volumes and injecting day for CS3 obtained from the proposed model (batch, day)

Product

volume

Discharging time (h)

Injecting day

(1,1) (1,2) (1,3) (2,3) (1,4) (1,5) (1,6) (1,7) (1,8) (1,9) (1,10) (2,10) (3,10) (1,11) (2,11) (3,11) (1,12) (1,13) (1,14) (1,15) (1,16) (1,17) (1,18) (1,19)

P1 P1 P4 P1 P2 P1 P2 P1 P1 P4 P1 P4 P5 P6 P5 P4 P1 P2 P1 P4 P1 P1 P2 P1

10800 10800 3800 8000 12807 10800 10800 10800 10800 10800 8000 3800 800 4703 3200 3800 10800 10800 10800 11800 10800 10800 10800 10800

24 48 54.22 72 96 120 144 168 192 216 230 238.2 240 250.45 257.56 264 288 312 336 360 384 408 432 456

0 0

1

1 2 3 4 5 6 7 8 9 9 9 9 10 10 11 12 13 14 15 16 17 19

(batch, day)

product

Volume

Discharging time

Injecting day

(1,20) (2,20) (1,21) (2,21) (3,21) (1,22) (2,22) (3,22) (1,23) (2,23) (3,23) (1,24) (1,25) (1,26) (2,26) (1,27) (2,27) (1,28) (1,29) (1,30) (1,31) (2,31) (3,31)

P4 P1 P1 P4 P5 P6 P5 P6 P5 P4 P1 P2 P1 P1 P4 P4 P1 P2 P1 P3 P4 P5 P6

7600 8000 8000 5831 800 6200 2367 5637 800 3800 8000 15600 10800 8000 7600 3800 8000 10800 14922 14353 6900 800 3100

467.7 480 493.2 502.2 504 513.5 517.2 528 529.2 535.1 552 576 600 612.3 624 630.23 648 672 696 720 735.33 737.11 744

18 18 19 20 20 20 21 21 21 21 21 22 23 24 24 25 26 26 27 28 29 30 30

4.2. Studying additional features In this section, the results of adding peak hours, discrete flow rate, and refinery inventory management are studied separately. So, by adding Eqs. (26) - (27) and term ‰„ to the terms of constraints and objective function, respectively, the effect of peak hours on the schedule was studied. In addition, by adding Eqs. (21) - (23) to the constraints and replacing Eqs. (19) and (20) with Eqs. (24) and (25), the effect of discrete flow rate on the schedule was investigated. Moreover, by adding Eqs. (28) - (38) to the constraints, the effect of inventory management at refinery was examined. After the above changes, the peak hours of each day was set to 8 h and discrete flow rate options were to 450 and 650 (vu/h); then, these data were added to CS3 data. These results are shown in Table 6; although the above limitations were added, the size and performance of the proposed model did not change significantly. In the case that the flow rate was continuous and peak hour constraints were added, the final minimum inventory was equal to 76% and just 79.07 h of total 248 h peak hours was used. When discrete flow rate was considered, just 81.94 h of total 248 h peak hours was used. It should be noted that, if the weight of the second and fourth terms of the objective function changed to 2 and 4, then the used peak hours and final minimum inventory would diminish to 8 h (because of the initial condition, not to 0) and 51%, respectively. In addition, if Eqs. (28) - (35) were just added to the base model, injecting day of each batch would be obtained (this result is shown in the fifth and tenth columns of Table 5).

Table 6: Model performance for C3 considering peak hours, discrete flow rate, and refinery inventory management Discrete Peak hours & Refinery inventory Peak hours flow rate discrete flow rate management Final status Optimal Optimal Optimal Optimal Constraints 9529 9777 9808 15567 Continuous var. 2917 3086 3117 4876 Binary var. 558 651 651 831 CPU time 56s 63s 141s 219s Nodes 375 1825 2131 2562 Iteration 60311 120842 161217 239954 Total used peak hours 79.07h 81.94h Diff. (vu) 0 0 0 0 Minimum final inv. (%) 76 76 76 76

5. Conclusion Scheduling of a multiproduct pipeline connecting an oil refinery to a single DC was studied and led to an MILP model. The proposed model was an adaptation of the previous one provided by [3] which aimed to optimally schedule multiproduct pipelines. In the previous work, a unique set of batches for the whole time horizon was handled and the model decided when a batch would be 20

delivered (i.e., to which day it would be assigned). In the new representation, set I denoted the batches which could be discharged during a single day and, by increasing the number of total days, the number of batches would remain constant. Also, reductions in the model size and the computational effort would be yielded. A real-world example was studied that showed better performance in terms of solution time. Controlling the pipeline operation during peak electricity hours, continuous and discrete flow rates, and integrated inventory management at DC and refinery, which made the model closer to the real cases, were other features that was included in the model. Moreover, by introducing the additional restrictions (peak hours, discrete flow rate, and refinery inventory management), the model's running time did not significantly change. In future works, the proposed model could be extended to a case consisting of multi-DCs or the current formulation could be improved by adding more details. As an example, uncertain demands may be considered.

21

References [1] R. Jr. Rejowski, J. M. Pinto., Scheduling of a multiproduct pipeline system, Computers and Chemical Engineering, 27 (2003) 1229-1246. [2] D. C. Cafaro, J. Cerda, Optimal scheduling of multiproduct pipeline systems using a non-discrete MILP formulation, Computers and Chemical Engineering, 28 (2004) 2053-2068. [3] S. Relvas, S. N. B. Magatão, A. P. F. D. Barbosa-Póvoa, F. Neves., Integrated scheduling and inventory management of an oil products distribution system, Omega, 41 (2013) 955-968. [4] R. Jr. Rejowski, J. M. Pinto., Efficient MILP formulations and valid cuts for multiproduct pipeline scheduling, Computers and Chemical Engineering, 28 (8) (2004) 1511-1528. [5] L. Magatao, L. V. R. Arruda, F. Jr. Neves. , A mixed integer programming approach for scheduling commodities in a pipeline, Computers and Chemical Engineering, 28 (2004) 171-185. [6] S. A. MirHassani, M. Ghorbanalizadeh., The Multi-product Pipeline Scheduling System, Computers and Mathematics with Applications, 54 (4) (2008) 891-897. [7] S. A. MirHassani, H. F. Jahromi., Scheduling Multi-Product Tree-Structure Pipelines, Computers and Chemical Engineering, 35 (2010) 165-176. [8] S. Relvas, H. A. Matos, A. P. F. D. Barbosa-Po´voa, and J. O. Fialho, A. N. S. Pinheiro, Pipeline Scheduling and Inventory Management of a Multiproduct Distribution Oil System, Industrial & Engineering Chemistry Research, 45 (2006) 7841-7855. [9] S. Relvas, H. A. Matos, A. P. F. D. Barbosa-Po´voa, and J. O. Fialho, Reactive Scheduling Framework for a Multiproduct Pipeline with Inventory Management, Industrial & Engineering Chemistry Research, 46 (2007) 5659-5672. [10] D. C. Cafaro, J. Cerda., Efficient Tool for the Scheduling of Multiproduct Pipelines and Terminal Operations, Industrial & Engineering Chemistry Research, 47 (2008) 9941-9956. [11] S. Relvas, A. P. F. D. Barbosa-Po´voa, H. A. Matos., Heuristic batch sequencing on a multiproduct oil distribution system, Computers & Chemical Engineering, 33 (2009) 712-730. [12] S. A. Mirhassani, N. Beheshti Asl., A heuristic batch sequencing for multiproduct pipelines, Computers and Chemical Engineering, 56 (2013) 58– 67. [13] S. A. MirHassani, S. Moradi, N. Taghinezhad., Algorithm for long-term scheduling of multi-product pipelines, Industrial & Engineering Chemistry Research, 50 (2011) 13899-13910.

22

[14] S. A. MirHassani, M. Abbasi, S. Moradi., Operational scheduling of refined product pipeline with dual purpose depots, Applied mathematical modeling, 37 (2013) 5723-5742. [15] V. G. Cafaro, D. C. Cafaro, C. A. Mendez, J. Cerda., Detailed Scheduling of Single-Source Pipelines with Simultaneous Deliveries to Multiple Offtake Stations, Industrial & Engineering Chemistry Research, 51 (2012) 6145-6165. [16] S. N. Boschetto, L. Magatao, W. M. Brondani, F. Jr. Neves, L. V. R. Arruda, A. F. D Barbosa-Povoa, S. Relvas., An Operational Scheduling Model to Product Distribution through a Pipeline Network, Industrial & Engineering Chemistry Research, 49 (2010) 5661-5682. [17] T. M. T. Lopes, A. V. Moura, C. C. De Souza, and A. A. Cire, Planning the operation of a large realworld oil pipeline, Computers and Chemical Engineering, 46 (2012) 17-28. [18] D. C. Cafaro, J. Cerda., Rigorous formulation for the scheduling of reversible-flow multiproduct pipelines, Computers and Chemical Engineering, 61 (2014) 59-76.

23

This table should be placed on page 18 line 41.

Table 4: Operational indicators Impletion Medium flow rate (vu/h) diff (vu) Pipeline usage (%) Final inventory level (%) Minimum inventory (%) Minimum final inventory (%)

Relvas model 522.7 0 100 55.14 11.56 (P4) 42.59 (P1)

Proposed model 502.17 0 100 76.16 11.79 (P4) 76.16 (P1-P6)

24