Benders decomposition with integer sub-problem applied to pipeline scheduling problem under flow rate uncertainty

Accepted Manuscript

Benders decomposition with integer sub-problem applied to pipeline scheduling problem under flow rate uncertainty Asl N. Beheshti , S.A. MirHassani PII: DOI: Reference:

S0098-1354(18)31120-7 https://doi.org/10.1016/j.compchemeng.2019.01.003 CACE 6321

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

25 October 2018 5 December 2018 2 January 2019

Please cite this article as: Asl N. Beheshti , S.A. MirHassani , Benders decomposition with integer subproblem applied to pipeline scheduling problem under flow rate uncertainty, Computers and Chemical Engineering (2019), doi: https://doi.org/10.1016/j.compchemeng.2019.01.003

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.

ACCEPTED MANUSCRIPT slhgilhgiH   

AC

CE

PT

ED

M

AN US

CR IP T



Introduced a new scenario-based two-stage MILP model for scheduling pipeline under flow rate uncertainty Proposed an accelerated Sample Average Approximation (SAA) method to overcome the computational challenges. In the developed SAA method, the enhanced BD approach is proposed to solve the MILP problems. Experiments conducted on generated large instances as well as real life cases.

1

ACCEPTED MANUSCRIPT Benders decomposition with integer sub-problem applied to pipeline scheduling problem under flow rate uncertainty Beheshti Asl N., and 1MirHassani S. A. Faculty of Mathematics and Computer Science, Amirkabir University of Technology, Tehran, Iran

AN US

CR IP T

Abstract Important issues in a pipeline system are energy efficiency, reliability and throughput flexibility. Practically conventional pumps are not capable of operating at the highest attainable efficiency for long running time. This deficiency has prevented a pipeline system to operate close to its predefined program. A possible remedy is to take into account uncertainty due to pumps operations. In this paper, a stochastic two-stage mixed integer programming (MILP) model is developed for the multiproduct pipeline-scheduling problem under flow rate uncertainty. The problem arises in a number of settings, and the real-world applicability discussed and demonstrated. The stochastic MILP model involves many discrete variables that make it intractable for real-life cases. As a solution method, the sample average approximation is combined with a three-step algorithm based on Benders decomposition. The modeling and solving approach is evaluated in some case studies including a real-life problem from NIOPTC.

AC

CE

PT

ED

M

Keywords: Pipeline scheduling; Uncertain flow rate; Two-stage stochastic model; Disruption; Benders decomposition method; Sample average approximation.

1

Corresponding author: [email protected]

2

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

1 Introduction In reality, petroleum products such as different grades of gasoline, kerosene, jet fuel and diesel can be transported through road, railroad and vessels as well as pipeline. However, large amounts of petroleum products are carried through pipelines from refineries to distribution terminals which is safer and less expensive than other modes of transportation. Therefore, pipelines are the most effective and reliable means for transporting refined petroleum products from major sources to distant distribution centers (DCs). Zhang et al. (2016) reported approximately 17.95 million barrels of oil products are imported and exported every day around the world, most of which are transported to different cities through pipelines. As a result, the scheduling and planning of transportation operation of oil products through pipelines are becoming more important. Furthermore, it becomes more complex when uncertainty assumption is considered, because multiproduct pipelines face a broad array of disruptions and uncertainties. The uncertainty causes disruptions to an operational plan, makes difficulties to the smooth execution of a perfect plan, and even causes it to fall apart, Yu and Qi (2004). The uncertainty of certain parameters such as flow rate due to unexpected events such as pumps failure has the most effect on pipeline operation. Muhlbauer (2004) reported that the flow rate may suddenly decrease in 17% of cases. This demonstrates the impotence of the uncertainties in the models. Indeed, without such a modeling tool, decision-makers might need to carry out significant modifications on their original plans during the recourse; which means loss of time and effort in saving human life and property. In spite of this fact, because solving the large-scale stochastic oil pipeline problem is an optimization challenge, disruption assumption is rarely addressed in the pipeline-scheduling problem. In this paper, we investigate the use of two-stage stochastic MILP for scheduling of multiproduct pipeline problems. The computational complexity of such problems has yet remained the major difficulty. This requirement led to the development of a large-scale model that makes the solution process more complicated. Polli et al. (2017) showed that the MILP problem is NP-hard, even in a deterministic context; Consequently, developing an efficient solution approach becomes essential. There are several exact methods for solving stochastic problems. However, these methods become impractical when the number of scenarios is large. Therefore, we herein focus on two types of methods: a sampling based method and a three-step solution approach based on Benders Decomposition (BD) method. We develop a three-step solution approach based on the BD method in order to solve the model appropriately. Many researchers have investigated the problem in different detailed levels and proposed models and solution methods taking into account technical and operational restrictions. As early started by Hane and Ratliff (1995), they proposed the discrete model for the problem of sequencing commodities in a single pipeline. Later on, Rejowski and Pinto (2003) developed a discrete time MILP model for a real case pipeline. In this system, a single pipeline transfers oil products from a refinery to several DCs. In their scheduling model, they focused on the convenience of irregular operations so that working during peak electricity periods could be avoided. Relvas et al. (2006) developed a mathematical formulation to account for pipeline scheduling and inventory management at the DCs. This model become more realistic by changing the fixed pumping rate to a variable rate. Later, they considered a special settling period for each product, as reported in Relvas et al. (2007) and Relvas et al. (2009). MirHassani and Ghorbanalizadeh (2008) proposed a MILP formulation for tree-structured pipeline problems and supposed that pipeline segments could be divided into two types, common and private. They defined a path between refinery and each DC and investigated the transfer of refined product in the relevant pathway.

3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Cafaro and Cerdá (2009, 2009) and Castro (2010) proposed the first mixed integer model based on continuous time and volume, which does not require time to be discrete or the pipeline to be seqmented. They formulated the system with multiple sources and destinations. Afterwards, Cafaro and Cerdá (2010) extended their approach to allow simultaneous batch injections at more source nodes at any time. Other works have emphasized the field of scheduling on pipeline networks with more complex structures. MirHassani and Jahromi (2010) presented a new continuous-time MILP formulation for the tree-structure singlesource pipeline operational planning problem. Cafaro and Cerdá (2012) considered a meshstructure pipeline network system that can be described as a set of interconnected pipelines with several inputs and distribution terminals. These approaches provide the set of batch stripping operations to be diverted to the distribution terminals during every pumping run; however, they specify neither the detailed sequence of individual cuts to be performed by the pipeline operator nor the flow rate profile at every pipeline segment. In response to this need, some approaches have been recently proposed for generating a detailed pipeline scheduling. Accordingly, Cafaro et al. (2015) provided an approach with parallel deliveries. Their optimization models can also handle a specific flow rate range for each pipeline segment. Recently some other researchers, Taherkhani et al. (2017), Castro and Mostafaei (2017) and Hitoshi et al. (2018) proposed optimization model with simultaneous deliveries assumption that make the model more realistic. Taking different limitations and assumptions into account can make the model more realistic and at the same time, make it more difficult to solve. They proposed hybrid approaches that integrates heuristic procedures and MILP models, to solve complex model. Furthermore, Cafaro & Cerdá (2014) and Moradi and MirHassani (2016) consider further assumptions reversible flow and demand uncertainty respectively that are two features that make the model more realistic. The main contributions of this paper are: (1) proposing a new scenario-based two-stage MILP model for scheduling pipeline under flow rate uncertainty and enhancing its capability for dealing with more realistic conditions; (2) introducing an accelerated Sample Average Approximation (SAA) method to overcome the computational challenges inherent in efficiently solving these stochastic problems. In the developed SAA method, the accelerated BD approach is proposed to solve the MILP problems. We use the SAA method and BD method together to overcome the computational complexity associated with the model structure and a large number of scenarios. The three-step solution scheme decomposed the original problem to the two MILP problems, called master problem (MP) and sub-problem (SP). Because of discrete variables in the SP, it is not possible to utilize the conventional BD method. Therefore, a three-step approach based on the BD method is proposed. The advantage of the three-step algorithm is fast convergence. The introduced algorithm leads to reduction in the number of binary variables and logic constrains and generates a favorable structure that can be exploited in the conventional BD method. This paper is organized as follows: Section 2 includes the problem definition and assumptions regarding scheduling, blending and inventory management. Section 3 presents the stochastic MILP model and detailed mathematical formulations of the two-stage model. A new accelerated SAA approach is combined with a three-step based on the BD algorithm introduced in Section 4. In Section 5, case studies and computational results for a real-life pipeline are presented, and the last section 6 gives the conclusions. 2 Problem statement Multiproduct pipelines permit to transport batch of different products through pipeline without any physical separation. Batches of oil-refined products are pumped back-to-back into pipeline. The products that are pumped into pipeline must be discharged and stored in dedicated tanks at the DCs, and the total inventory level update is done daily. The sequence 4

ACCEPTED MANUSCRIPT and volume of such batches of products should be carefully selected in order to satisfy DC’s demands at the promised dates. The problem under study consists of a straightforward multiproduct pipeline used to transfer various types of oil products from a single refinery to a single DC. The challenge is to specify when and how much of each product should be injected into the pipeline to satisfy daily demands with minimum delay while observing many operational restrictions. The pipeline scheduler develops a schedule of pipeline activities over a time horizon. This approach avoids the complex pipeline operational representation, by focusing on the end rather than the beginning, as reported in Relvas et al. (2013). Two sets of given information about pipelines are as follows:

M

AN US

CR IP T

Given: 1. The operational parameters such as tank capacities. Minimum/maximum storage capacity of each type of products at the distribution center. The number and different types of products. The pipeline volumetric capacity. The planning time horizon and whole number of days to be considered. The lower and upper limit of the size of each batch, the daily product demands at the DC during the scheduling horizon. The normal flow rate, disruption flow rate for each scenario, start/end time for each scenario. The matrix of forbidden sequences between pairs of different products in consequent batches (some sequences of products in the pipeline are forbidden due to excessive contamination) and interface volume between each pair. 2. The initial pipeline state such as sequence of batches inside the pipeline at the beginning of scheduling time horizon and their volumes (these products are separated in predefined batches, and must be preserved) and the initial inventory level of every refined product in terminal tanks.

AC

CE

PT

ED

2-1 Disruption and uncertainty assumption Uncertainty is hardly avoidable or preventable. Optimization models derived from real world problems (such as pipeline scheduling) regularly involve uncertainty. It often disrupts pipeline operations, causing significant adverse effects in the system; thus, it is essential to be considered as a fact in planning and scheduling of multiproduct pipelines. A typical method of generating an operational plan within an uncertain environment is to use models based on scenarios. Stochastic programming deals with optimization problems whose parameters take values from given probability distributions. This situation is closer to the actual industrial practice and reality. It is important to develop stochastic based optimization models to reduce or mitigate the impact of these uncertainties. One of the common troubles in oil pipelines is disruption at pumps or burst pipeline. As a result, flow rate is severely reduced or cut off suddenly. To minimize adverse effects, we provide a novel scenario-based MILP two-stage stochastic model for scheduling and managing pipeline under flow rate uncertainty. Future uncertainties can be described by the concept of scenarios, where each scenario is the outcome in which the system parameters are changed within certain ranges over a period of time. We considered the set as representing all the scenarios. In many cases, there is one scenario, , whose realization probability is by far the largest. We call this the normal scenario and all the other scenarios in , the disruption scenarios. When a scenario other than is realized, we say that a disruption has occurred, see Yu and Qi (2004). For different operating systems, types of disruptions vary. In general, these types include but are not limited to uncontrollable events, changes in availability of resources, and new external restrictions. In this paper, we consider the type of disruptions that cause the flow rate to decrease, such as pump failures, power outage, and pipeline burst. In such situations, the flow rate has to be reduced or vanish. 5

ACCEPTED MANUSCRIPT

CR IP T

In the present study, we consider the flow rate to be reduced deeply due to disruptions at different time-periods. When disruptions occur, for whatever reason, the normal flow rate is reduced for a while. The normal flow rate can vary during the time horizon, but to avoid the complex pipeline operation, schedulers usually practice an average flow rate that we show with ̅̅̅̅̅ . Disruption scenarios are defined based on three important sources of data: start and end time of disruption and the level of flow rate reduction. Let the set of all the scenarios be , a scenario represents three parameters: start time, end time and flow rate, shown by ̅̅̅̅̅̅̅ ̅̅̅̅̅ and ̅̅̅̅̅ , respectively. By the scenario (̅̅̅̅̅̅̅ ̅̅̅̅̅ ̅̅̅̅̅ ) we mean flow rate changes from normal ̅̅̅̅̅ to reduced one,̅̅̅̅̅ , during interval time ̅̅̅̅̅̅̅ ̅̅̅̅̅ . In Section 5, we explain how to generate different scenarios.

CE

PT

ED

M

AN US

2-2 Operational assumptions In order to provide our method, the following assumptions are considered. These assumptions usually have been considered in all previous models of scheduling multiproduct pipeline, see Cafaro and Cerdá (2008). 1. A single multiproduct pipeline with unidirectional flow rate is considered, conveying refined petroleum products from a refinery to one DC. 2. The pipeline remains full of products at any time. By assuming liquid incompressibility, as measure of the relative volume change of the fluid under pressure, the only way to discharge a volume of product into DC is by injecting an equal volume at the refinery. 3. The DC contains a tank farm, where all tanks for each product are considered as aggregated tanks. 4. The volume of each batch injected in the pipeline must lie between known lower and upper bounds corresponding to the product type, and no two adjacent batches could contain the same product. 5. Product injection takes place only in the refinery, and DC is not able to inject. Moreover, one tank at the refinery and one tank at DC are simultaneously connected to the pipeline during inject/discharge operations. 6. The refinery is capable of providing the required products. Therefore, there is no need to monitor product inventory levels in the refinery. 7. Daily demand for each product is determined in the DC throughout the time horizon and must be provided before the day ends. There will be penalty costs for back-orders. 8. It is assumed that product cannot be sent after another batch of the same product, i. e., every new batch has a new interface. 3

Mathematical formulation

AC

The problem goal is to establish the optimal sequence of batches, which is relatively respectable for all the scenarios. In other words, it is necessary to determine an efficient and appropriate schedule for pipelines, consisting of a sequence of products to be injected, batch size, discharge duration and inventory management for satisfying demands with minimum possible delay and minimum operational cost. In brief, the model is composed of four main tasks: a) Determining the type of batches that are injected into the pipeline during the time horizon (batch sequence). b) Determining the volume of each batch (batch size). c) Inventory management in DC. d) Calculating the overall operational cost. 6

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Here, the pipeline scheduling problem is formulated as a scenario-based, Two-Stage Stochastic Problem (TSSP) model to represent the randomness arising from disruption of pumping. In general, TSSP consist of two distinct components: structural components (firststage) that is fixed and free of any uncertainty; and status components (second-stage) that is affected by the uncertainty. Nomenclatures Sets Set of batches Set of scheduled days Set of refined petroleum products Set of scenarios Indicated by over lined symbols Parameters ̅̅̅ Minimum batch volume of the product (vu) ̅̅̅ Maximum batch volume of the product (vu) Cost or penalty of one unit shortage of the product at the day ̅ Cost of interface between consecutive batches ̅ ̅ Cost of transportation of product ̅̅̅̅̅̅ Total demand of the product to be delivered at the day (vu) ̅ End of time horizon (days) ̅̅̅̅̅ Minimum storage capacity of the product (vu) ̅̅̅̅̅ Maximum storage capacity of the product (vu) ̅̅̅̅̅ Initial volume of the product at DC (vu) ̅̅̅̅̅ Initial volume of the product inside the pipeline (vu) ̅̅̅̅̅ Sufficient large constant Sufficient small constant ̅̅̅̅ Minimum number of batches to be sent during scheduling horizon Probability of the disruption scenario ̅ The volume injected inside the pipeline since the beginning of time ̅̅̅̅ horizon until immediately before the start of the disruption scenario , ̅̅̅̅̅ ̅̅̅̅̅̅̅ (vu) which is obtained by ̅̅̅̅ The volume injected inside the pipeline since the beginning of time ̅̅̅̅ horizon until immediately after the end of the disruption , which is ̅̅̅̅ ̅̅̅̅̅ ̅̅̅̅̅̅̅ ̅̅̅̅̅ (vu) obtained by ̅̅̅̅ The permitted sequences inside the pipeline. If ̅̅̅̅̅ , the product ̅̅̅̅̅ can be injected after ; otherwise, it is forbidden. ̅̅̅̅̅̅̅ Start of the disruption scenario (days) ̅̅̅̅̅ End of the disruption scenario (days) ̅̅̅̅̅ The flow rate under the disruption scenario (vu/day) ̅̅̅̅̅ The flow rate under the normal scenario (vu/day) Continuous variables Shortage of the product at the day under scenario (vu) The portion of volume of the batch that is discharged during interval time ̅̅̅̅̅̅̅ ̅̅̅̅̅ under scenario (vu) Inventory of the product at the day under scenario (vu) Duration of discharging time of the batch under scenario (days) Complete discharging time of the batch under normal flow rate (days) 7

ACCEPTED MANUSCRIPT

AN US

CR IP T

Duration of discharging time of the batch under normal flow rate (days) Volume of the batch containing the product (vu) A continuous variable which is positive if batch completely discharged before the start of disruption under scenario (vu) A continuous variable which is positive if batch completely discharged after the start of disruption under scenario (vu) A continuous variable which is positive if batch completely discharged before the end of disruption under scenario (vu) A continuous variable which is positive if batch completely discharged after the end of disruption under scenario (vu) Complete discharging time of the batch at DC under scenario (days) Volume of the product in the batch completely discharged from the pipeline between day and day under scenario (vu) Cumulative delay of discharging time until batch under scenario (days) Delay of discharging time of the batch under scenario (days) Binary Variables 1 if batch is discharged between day and day under scenario 1 if batch contains the product 1 if batch is discharged after the start of disruption under scenario 1 if batch is discharged after the end of disruption under scenario In the following, the constraints of the two-stage model are classified.

M

3-1 Batch-defining constraints Each new batch might contain at most one of the refined products. Thus ∑

(1)

∑

PT

∑

ED

If a batch does not contain any product, it is a fictitious batch and never gets pumped into the pipeline. The fictitious batches are confined to the end of the pumping run sequence through the following constraint: (2)

AC

CE

Because of product contamination, some product sequences in the pipeline are forbidden. For a pair of forbidden sequence the parameter is zero (̅̅̅̅̅ ) and it would be excluded from batch sequences. ̅̅̅̅̅ =0 (3) Because of some operational restrictions, the volume of each batch, based on product type, must be placed in a specific range. Therefore, we have ̅̅̅ ̅̅̅ (4) 3-2 Inventory management The daily demand of distribution centers dictates that the time horizon be discretized in several time-period (days). Under each scenario the total inventory of the product at the end of day is equal to the inventory level at the previous day plus the discharged volume of product and the shortage at day , minus the volume of product used to satisfy demand in the day and the shortage of product at day . ̅̅̅̅̅̅ ∑ (5) 8

ACCEPTED MANUSCRIPT

Inventory profiles is limited by the predefined values. ̅̅̅̅̅ ̅̅̅̅̅

(6) When the batch is completely discharged under the scenario in the day ( ),) then the discharged volume of the batch under the scenario in the day must be equal to the volume of batch and conversely. In other words, if and only if . These ̅̅̅ ̅̅̅̅̅ logical constraints are defined as follows where . ̅̅̅̅̅ ̅̅̅̅̅

(8) (9)

is equal to 1 if and only if, the receiving ∑

AN US

3-3 Scheduling constraints The constraints indicate that the binary variable of batch is completed at day . ̅ ∑ ∑

CR IP T

(7)

(10)

The discharging of batch will be completed in one day of time horizon, for each scenario . ∑

(11)

ED

M

The normal discharging time of batch is equal to the volume of batch divided by normal flow rate. ∑ (12) ̅̅̅̅̅ Complete discharging time of the batch without any disruption is equal to the sum of discharging time of all previous batches (from 1 to ). (13)

PT

∑

CE

Complete discharging time of batch under the scenario is equal to the complete discharging time of the batch under the normal flow rate plus cumulative delay of complete discharging time until the batch under the scenarios . (14a) Note that, in summary can be stated as follows:

AC

̅̅̅̅̅ ∑ ∑

̅̅̅̅̅ ∑

(14b)

Cumulative delay of complete discharging time until the batch under the scenarios is equal to the cumulative delay of complete discharging time until the batch under the scenario plus the delay discharging time for the batch under the scenario . (15) (16) ̅̅̅̅̅ In order to calculate the discharging time delay under scenario , we need to know the DCs that cope with the batch during the interval time ̅̅̅̅̅̅̅ ̅̅̅̅̅ . Therefore, to calculate , 9

ACCEPTED MANUSCRIPT we need to know the status of the batch . There are three cases for discharging time as is depicted in Figure 1: I. The discharge of batch is completed before the disruption ( see batch 1 in Figure 1. II. The discharge of batch is completed after the disruption ( see batch 5 in Figure 1. III. The discharge of batch is effected by disruption. See batch 2 ( batch 3 ) and batch 4 ) in Figure 1. In brief, we face with the following statuses under scenario . Status of batch vs disruption

1

The discharge of batch 1 is started before and is completed before the disruption. The discharge of batch 2 is started before and is completed during the disruption. The discharge of batch 3 is started during and is completed during the disruption. The discharge of batch 4 is started during and is completed after the disruption. The discharge of batch 5 is completed disruption

2

,

disruption disruption

disruption

4 5

Batch 5 Refinery

M

AN US

3

),

CR IP T

Batch #

),

Batch

ED ̅̅̅̅̅̅ 𝑠𝑡𝑎𝑟𝑡

Time

𝑠

after the

Flow direction Batch

Batch

4

disruption

Batch 1 2

3

̅̅̅̅̅ 𝑒𝑛𝑑

DC

𝑠

PT

Figure 1:Status of batches in ouccurance of disruption.

AC

CE

In addition, we need to determine how much of the batch is discharged in DC with normal flow and disrupted flow rate. Based on definition, ̅̅̅̅ is the total volume injected into the pipeline prior to the initiation of the scenario . If the sum of volume from the first batch to the batch is less than ̅̅̅̅ , then the indicator variable is 0; Otherwise is 1 ( ). The same explanation is applied to variables . These equations show the status of the batch at the start/end of the scenario . ̅̅̅̅

∑

∑

(17)

̅̅̅̅

∑

∑

(18)

̅̅̅̅̅ (19) ̅̅̅̅̅ (20) ̅̅̅̅̅ (21) ̅̅̅̅̅ (22) These equations determine how much of the batches is discharged in DC under flow rate of the scenario . 10

ACCEPTED MANUSCRIPT ̅̅̅̅̅

̅̅̅̅̅

(23) If one part of the batch is discharged prior to the start of the scenario and the other parts are discharged after the initiation of the scenario (see batch 2, Figure 1), then the can be calculated based on the following equation: ̅̅̅̅̅ ̅̅̅̅̅ (24) If the total volume of the batch is discharged during the scenario (see batch 3, Figure 1), then the will be calculated according to the following equation: ̅̅̅̅̅ ̅̅̅̅̅ ∑ (25) and another part is can be calculated

CR IP T

If one part of the batch is discharged before the end of the scenario discharged after the end of the scenario (see batch 4, Figure 1), then based on the following equation: ̅̅̅̅̅ ̅̅̅̅̅ ∑

AN US

If the total volume of the batch is discharged after the end of the scenario Figure 1), then the can be obtained from the following constraint: ̅̅̅̅̅ ̅̅̅̅̅

(26)

(see batch 5, (27)

AC

CE

PT

ED

M

3-4 Tightening up the feasible region The major objective here is to add linear constraints that cut off parts of the feasible region in the Linear Programming (LP) relaxation of the problem, without cutting off any integer feasible solutions. Therefore, the feasible region for the LP relaxation will be cut down and close to the convex hull of the integer feasible points. In other words, the following are some valid inequalities based on the nature of the problem to tighten the search space. If batch is discharged after the start of the scenario ( ), then batch will be discharged after the start of the scenario too ( ). (28) If batch is discharged after the end of the scenario ( ), then batch will be discharged after the end of the scenario too ( ). (29) If batch is discharged after the end of the scenario ( ), then it has to be discharged after the start of the scenario ( ). (30) If batch is discharged after the start of the scenario ( ) then the discharging day of ̅̅̅̅̅̅̅ and it must not be before the ̅̅̅̅̅̅̅ . In other words, for , we have . ̅̅̅̅̅̅̅ (31) If batch is discharged after the end of the scenario ( ) then the discharging day of ̅̅̅̅̅ and it must not be before the ̅̅̅̅̅ . In other words, for , we have . ̅̅̅̅̅ (32) If a batch is not a fictitious batch, in other words, it contains product (∑ ), then it has to be discharged at most one day for each scenario (∑ ).

∑

∑ ∑

11

(33) (34)

ACCEPTED MANUSCRIPT If a batch is a fictitious batch, in other words, it does not contain any product (∑ ), then it does not discharge after the start/end of any scenario ( and ). ∑

(35)

∑

(36)

CR IP T

3-5 Objective function The objective function consists of the first-stage decision costs and the expected value of the second-stage recourse costs including the penalty of unsatisfying demands. It minimizes the following cost functions: a) The operational cost of the interface between consecutive batches: ∑ ̅

(37)

AN US

b) The product transportation cost: ∑̅

(38)

c) Expected value of penalty of shortage imposed by the unsatisfied demands under each scenario: (39)

M

∑ ̅ ∑ ̅

AC

CE

PT

ED

4 Solution strategy The proposed TSSP model consists of the configuration decisions (batches sequence and volume) as a first stage decisions and the second stage involves a large number of recourse decisions. The proposed MILP formulation clearly yields a large and complex model that cannot be solved directly. As shown in the constraints section, the proposed MILP model involves logical implications modeled by Big-M (BM) coefficient and is among the most complicated problems to solve. The SAA method is a well-known approach to solve stochastic optimization problems and to reduce a scenario set to a manageable size. Regarding the difficulties posed to solve the MILP problem using the SAA method, an improved stochastic solution method is introduced. In the proposed approach, a three-step algorithm is developed based on the BD to find the lower bound on the original two-stage stochastic problem. Then, the deterministic equivalent of the SAA problem is solved to find the upper bound. The procedure repeats until the optimality gap gets small enough.

4-1 Sample average approximation In section 3, we introduced the TSSP model as: (TSSP)

∑ ̅

̅

∑

̅ ∑ ̅

(40)

For the real-world optimization problems, TSSP model (40) is to solve all scenarios of real conditions. In these cases, the total number of the scenarios is extremely large or even infinite. Therefore, the model has an infinite number of variables and constraints, which are 12

ACCEPTED MANUSCRIPT

CR IP T

impossible to be solved. The SAA method is a common approach to breakdown a scenario set to a manageable size. This method was proposed to solve stochastic linear programs with finite set of scenarios ( ). Many researchers such as Shapiro (1996), Shapiro and Homem de Mello (1998) and Mak et al. (1999) employed this method. The SAA proceeds by repeatedly solving a stochastic problem over M independent samples of the size N. The procedure is repeated several times until the stopping criterion, which is the optimality gap, gets small enough. To illustrate the point, the outlines of the SAA scheme are reported according to Kleywegt et al. (2001). The SAA approach consists of three main steps: The first two steps calculate statistical lower and upper bounds of the problem (40) and the last step is to estimate the optimality gap. Step1: Lower bound estimation Generate Independent and Identical Distributed (IID) samples of size , i.e., . For , solve the corresponding SAA problem as below: ̅̅̅̅̅

(

)

AN US

∑

where

(41)

̅ ∑ ( ) ∑ ̅ ̅ and . In addition, let ̅̅̅̅̅ be the corresponding optimal objective value. Then one statistical lower bound and its sample variance can be estimated as follows: ∑ ̅̅̅̅̅

̅̅̅̅̅

(̂

̂

(42) (43)

)

(44)

PT

∑

̅̅̅̅̅

ED

∑ ̅̅̅̅̅

̅̅̅̅̅

M

̅̅̅̅̅

AC

CE

Step 2: Upper bound estimation Let ̂ be the minimum ̅̅̅̅̅ corresponding to the best-computed ̂ of the first stage variables. Generate an IID sample solution ̂ with sufficient large size . The statistical upper bound can be computed as: and the corresponding variance can be estimated as follows: ̅̅̅̅̅

∑

(̂

̂

)

̅̅̅̅̅

(45)

Step 3: Optimality gap estimation Compute the optimality gap between the lower and upper bound. (46) and (47) then give the estimated variance of the above gap, respectively. The procedure repeats until the optimality gap gets small enough ( ̅̅̅̅̅ ). ̅̅̅̅̅ ̅̅̅̅̅ (46) ̅̅̅̅̅ ̅̅̅̅̅ ̅̅̅̅̅ (47) Although the SAA method is effectively used in many applications, it is not practical due to a significant computational effort required for solving NP-hard problem (MILP problems). 13

ACCEPTED MANUSCRIPT Even after approximating the expectation using the SAA technique, the computational burden of solving the resulting MILP problems may still be constraining. In other words, the first step is computationally the most difficult step and solving the M stochastic programs with N scenarios in each replication is extremely time consuming, for more detail see Mohammadi Bidhandi and Patrick (2017). To overcome the difficulty of solving MILP problems, we developed a solution method. The proposed method is a three-step algorithm based on the BD method to overcome this difficulty.

ED

M

AN US

CR IP T

4-2 Benders decomposition Benders decomposition is a solution method fixing complicating variables temporarily and consequently making a problem significantly easier to handle, Benders (1962). It has become a popular technique to solve certain types of complex problems, including MILP problems and stochastic programming problems. In the BD method, the original problem is decomposed into a master problem (MP) and a sub problem (SP) problem. (48) (49) (50) (51) The MP/SP is solved iteratively by utilizing the solution of one for the other. The MP contains a set of the complicating (integer or continuous) variables and their related constraints. The SP is obtained by temporarily fixing complicating variables in the original problem (OP). The SP solution generates the Benders cut for each iteration repeatedly added back to the MP. The algorithm iteratively solves the MP and SP until an optimal solution is found. To confirm the convergence of algorithm, the optimality gap is calculated at each iteration. The objective function of the MP and the SP proposes a valid lower bound and an upper bound based on the optimal cost. In each iteration, , MP is a relaxation of the original problem and its solution provide a lower bound. On the other hand, adding contribution of MP’s variables to the objective value of the SP yields a valid upper bound. (52) ̅ ̅ ̅ (53) ̅ ̅ ̅

PT

In general, the MP is considered as MILP of the form:

AC

CE

(54) (55) ̅ (56) ̅ ( (57) ̅ ̅ ) ̅ (58) ̅ ̅ (59) The parameter ̅ is a lower bound for continuous free variable . Let the solution vector of MP at iteration be ̅ . If ̅ is feasible for SP then the constraints (57) must be added to MP and if ̅ is infeasible for SP then the constraints (58) must be add to MP. Constraints (57) and (58) are referred to as optimality and feasibility Benders cuts, respectively, as demonstrated Van Slyke and Wets, (1969), where ̅ is representing the value of dual variables associated with the constraint (62). The sets and were defined as sets of indices for the optimality cut and feasibility cut respectively. The SP is considered as: (60) (61) (62) ̅ (63) 14

ACCEPTED MANUSCRIPT

CR IP T

Many successful applications of the BD are found in many diverse fields including planning and scheduling; Hooker (2007), transportation; Galereh et al. (2015), resource management; Zhang and Ponnambalam (2006). Benders, (1962) initially proposed for the MILP programming for which integer variables were considered to be complicating variables and appear in the MP. When some of the SP variables are required to be integer, a different theoretical framework to develop the Benders cut is needed. A common strategy is based on solving the linear relaxation of the SP to generate cut; Papadakos (2008). Many extensions of the method have developed to address an integer SP. Hooker (2011) and Codato and Fischetti (2006) introduced a new extension known as logic-based BD in which the MP is a 0-1 integer problem and the SP is a feasibility problem. These methods do not have a standard template to develop valid optimality cut and are dependent upon the problem structure. The major disadvantage of the BD approach is slow convergence. Many extensions (e.g. Hooker and Ottosson ( 2003), Geoffrion (1972) and Magnenti and Wong (1981)) have been developed to apply the algorithm and improve its convergence. In our case, the proposed TSSP decomposed to the two MILP problems as follows: ̅

(64)

̅ ̅

∑̅

̅

(

constraints (1) – (4) And

)

M

∑ ̅ ∑ ̅

AN US

∑ ̅

ED

̅

(65) (66)

(67) (68)

AC

CE

PT

4-3 Discrete sub-problem and sub-problem reduction procedure In this work, both MP and SP are MILP problems; hence, the standard BD cannot be applied to generate the Benders cut. On the other hand, other methods that are dealing with the discrete SP are very slow, especially when the number of binary variables is large. Moreover, there are many logical constraints in SP, which are not directly related to MP solution and keep big M coefficients in SP while MP variables are fixed. Therefore, a modified approach was presented to overcome this difficulty and solve the problem efficiency. The macro structure of original problem is: (69) (70) (71) (72) (73) Where, are vectors of both continuous and discrete variables, respectively and represent a set of nonnegative variables and . By fixing the master variables to appropriate value; Having in mind that and considering the model structure, allowed us to fix by using the second set of constraints. This process was applied as a SP reduction (SPR) procedure to remove all binary variables and logical constraints from the SP. The real sub problem that follows the Benders requirement is obtained by fixing in the third set 15

ACCEPTED MANUSCRIPT of constraints. Now, this is an LP problem to generate the Benders cut and guarantee convergence of BD algorithm. Considering the notation of the TSSP model, and represent sets of integer and continuous variables { }, and { } respectively and finally represents set of continuous variables { }. Also and represent ̅ ∑ ̅ and ∑ ̅ ∑ ̅ receptively and .

M

AN US

CR IP T

4-4 The three-step approach As a result, we proposed a three-step algorithm based on the BD method. In the first step, MP is solved to obtain a trial value for the variables and .Then, the SPR procedure is run in the second step, by considering the MP solution, to obtain the value of the following variables: , Finally, a linear SP is solved for continuous variables in third step, (See Figure 2).

Figure 2: Three-step approach based on BD method.

ED

First step: Solve the following MP problem. ̅ ∑ ̅

PT

̅ ̅

∑̅

(

(64) (65)

̅

)

(66)

AC

CE

constraints (1) – (4) Where ̅ represents the value of dual variables associated with the SP constraint (89). It is worth mentioning that this is a complete recourse model and any solution proposed based on the MP is feasible for the SP. Therefore, no need to add feasibility cut. Second step: Run the SPR procedure. By considering ( ̅ ) as a part of the MP solution in the following logical equations, we are able to evaluate the value of several continuous/binary SP variables. The solution obtained at each step is fixed on the next step and is indicated by over lined symbols. ̅̅̅̅ ̅̅̅̅ ∑ ∑ ̅ ∑ ∑ ̅ 1) { (74) ̅̅̅̅ ∑ ∑ ̅ 2)

16

{

̅̅̅̅

∑

∑ ̅

̅̅̅̅

∑

̅̅̅̅

∑

∑ ̅ ∑ ̅

(75)

ACCEPTED MANUSCRIPT

3)

{

4)

{

̅

(76) ̅

̅̅̅

(77)

̅̅̅ ̅ ̅ ̅

∑ ̅ ∑ ̅

5)

̅̅̅

̅ ̅

̅

̅ ̅

{ ̅ ̅̅̅̅̅̅

7)

∑ ̅ ̅̅̅̅̅̅

8)

∑

9)

̅̅̅̅ ̅̅̅̅

11)

{

(80)

̅̅̅̅

(81)

̅̅̅̅ ̅̅̅̅ ̅

̅

(82) (83) (84) (85)

M

̅

AN US

10)

̅

(78)

(79)

6)

12)

̅ ̅

CR IP T

̅

CE

PT

ED

After fixing some variables in second step (SPR procedure), the SP will be a convex optimization problem (linear programming), therefore the convergence of BD algorithm will be guaranteed; Benders (1962).

Third step: Solve the following linear SP.

AC

∑ ̅ ∑ ̅

̅̅̅̅̅

(86) ∑ ̅

̅̅̅̅̅̅

̅̅̅̅̅

̅ Note that the SP is bounded for any arbitrary MP solution. The flowchart of three-step algorithm with as the threshold is shown in Figure 3.

17

(87) (88) (89)

AN US

CR IP T

ACCEPTED MANUSCRIPT

M

Figure 3: The flowchart of three-step decomposition algorithm.

AC

CE

PT

ED

The SAA method combined with three-step technique has a modular structure as shown in Figure 4.

18

ACCEPTED MANUSCRIPT Figure 4: Schematic diagA

Start Generate 𝑀 independent sample sets of size 𝑁

CR IP T

Solve the SAAP model by three-step method for all samples

Obtain ̅̅̅̅̅ 𝑜𝑏𝑗𝑙𝑜𝑤𝑒𝑟 𝑣𝑎𝑟 ̅̅̅̅̅𝑙𝑜𝑤𝑒𝑟

Evaluate the best solution by (44)

Enlarge 𝑁

No

AN US

Obtain ̅̅̅̅̅ 𝑜𝑏𝑗𝑢𝑝𝑝𝑒𝑟 𝑣𝑎𝑟 ̅̅̅̅̅𝑢𝑝𝑝𝑒𝑟

Converge? Yes

ED

M

Set the current solution as optimal

End

Figure 4: Schematic diagram of the SAA method.

AC

CE

PT

4-5 Improving master problem The iterative solving of the MP and SP is a major computational bottleneck. In particular, the MP as a MILP formulation often lacks a special structure and its size continuously grows that makes it more difficult to be solved. It has often been reported that more than 90% of the total execution time of the BD method is spent on solving the MP (Magnenti & Wong, 1981). Some strategies are proposed to solve it more efficiently. Three approaches are proposed to improve the quality of the solution or to solve it more quickly: (1) using alternative formulations, (2) improving the MP formulation, and (3) using heuristic strategies to independently generate solutions or to improve those already found. In order to improve the MP in initial iterations and strengthen the MP, adding valid inequalities (VIs) can be of benefit as demonstrated Naoum Sawaya and Elhedhli (2010). Saharidis, et al. (2011) applied the BD method to a refinery system, adding two groups of VIs to the initial MP and also Tang et al. (2010) used VIs to improve the initial lower bound. We added some valid cuts to the MP to provide information about a part of the original objective function being removed from the MP. In this section, two speed-up constraints are presented to add the MP to reduce the search space. We introduced a new parameter ̅ as "the ratio of the total pumped volume and total volume of demand during time horizon". The initial added cut is as follows: 19

ACCEPTED MANUSCRIPT ̅ ∑ ̅̅̅̅̅̅

∑

(90)

The total volume of the product pumped through the pipeline during time horizon is often smaller than the total product demand. The added cuts forced the model to obtain value of variable in a way that the difference between total pumped volume and total demand reduced. In other words, this cut relates the MP solution to the SP objective function so the MP obtains better solution in each iteration and the algorithm converges more quickly. The value of ̅ can be changed at each iteration as ∑

̅

∑

̅

(91)

̅̅̅̅̅̅̅

AN US

CR IP T

Where ̅ is the solution of MP at the iteration. The initial value of the parameter ̅ is obtained by using some experimental information and it should usually be a small value at the beginning. When the original problem is decomposed into two problems (MP and SP), some dependent information may be destroyed and this would result in extremely weak solutions in early iterations. Thus, the number of the iteration to convergence would be extremely high. For example, the solution will be zero for large numbers of initial iterations, because the MP objective function minimizes the penalty for interface between each of the two successive batches while there is no constraint on the MP. In other words, the algorithm with many initial iterations manages to find a sequence of batches. Hence, a useful initial cut should be added to MP: ̅̅̅̅ ∑ (92)

∑

̅̅̅̅̅̅ {̅̅̅ }

(93)

PT

̅̅̅̅

ED

M

This cut prevents the MP from setting the value of all to zero at early iterations and ensures that the sequence has at least ̅̅̅̅ batches. The ̅̅̅̅ value is lower-bounded by dividing the total demand into the maximum batch size. Hence,

AC

CE

5 Computational result In this section, the new proposed model is solved by proposed three-step approach based on the BD method. The algorithm performance was compared to the results obtained from exact direct method on the same set of instances. The results were studied in terms of both CPU time and solution quality. The AIMMS 3.12 modeling software, using CPLEX 12.6.1 as a solver on a PC platform with a Pentium IV 2.4 GHz processor and 1 GB RAM, is used to implement all experiments. In a majority of previous studies, the scheduling time horizon was considered at most one month and this duration is ideal.

5-1 Performance of the three-step approach In order to verify the performance of the three-step solution method, a set of computational experiments were conducted with regard to the stochastic scheduling pipeline problems. We apply the proposed formulation for three case studies: CS1, CS2, and CS3 with time horizons of 10, 20, and 30 days, respectively. All the three test problems solved in this section are case studies based on the real world NIOPTC scenarios. It involves the transportation of five different oil products from refinery to the depot (Rey–Esfahan). Refinery is in Esfahan and is 400-km away from Ray. The pipeline capacity is 59,000 volumetric units; more details is given in MirHassani & BeheshtiAsl (2013). At the beginning, the pipeline is filled with the 20

ACCEPTED MANUSCRIPT product P1 (detailed data are shown in Table 1 and the forbidden sequences are represented by (*)). We consider a constant transportation cost which is independent of the product type. The flow rate can vary from 11,800 to 17,280 vu/day. We consider an average flow rate in the model. In presented model, the uncertainty is associated with scenarios and a solution is sought which is immunized against all possible scenarios. A scenario is a description of a future situation and the hypotheses used to create a proper scenario must simultaneously be relevant, consistent, and plausible. Product

̅̅̅̅̅̅̅̅

̅̅̅̅̅̅̅̅

̅̅̅̅̅

̅̅̅̅̅

P1 P2 P3 P4 P5

(uv) 17600 6300 5200 3400 6700

(uv) 26000 12000 16,000 16000 20000

(uv) 5100 3400 1000 1100 1300

(uv) 54136 38828 38760 39276 56771

Forbidden sequences P1 * *

P2 * * * *

P3 * * -

P4 * * -

P5 * * *

CR IP T

Table 1: Common data used in all case studies

̅̅̅̅̅

̅̅̅̅̅̅̅

(95) (96)

CE

̅̅̅̅̅

̅̅̅̅̅

PT

ED

M

AN US

The first scenario in each case study, called normal scenario, has a higher probability and is used for deterministic models. We then generate a large number of instances (100 different scenarios for each case). The number of occurrence of disruption in the pipeline is an event with Poisson process; more details had been shown in Dawotola et al. (2011); therefore in order to generate scenarios, we use the Gamma and Exponential distribution to have ̅̅̅̅̅̅̅ and by (94)-(96) formulae. These stochastic parameters are determined by using the Inverse Transform Method described in Graham and Talay (2013). The parameter is defined as the duration time of the system reforming. The parameters of these statistical distributions (Gamma and Exponential) are based on knowledge obtained from well-known expertise in oil industry. In all cases, the parameters of Gamma and Exponential distributions are considered as and respectively; see Dawotola et al. (2011) and Kent Muhlbauer (2004). The probability of each scenario is determined by using the techniques described in Graham and Talay (2013). (94) , ̅̅̅̅̅̅̅ ,

AC

Table 2 shows how to calculate the size of models and SPR procedure. Moreover, Table 3 summarizes the characteristics of three models: MP, SP and original problem (OP) for all case studies (namely CS1, Cs2 and CS3) with 100 scenarios. Table 2: The model size for different problem MP SPR Var. | || | | || | | || | | |

Int. Con.

| || | || | | | |

| || | | | | |

0

||

| || | || | || | ||

SP || | | 0 || | |

| || || | | |

OP || | | | | | || | | | | || |

| || | | | | || | | | | |( | | | || | | | | | || | |

)

The number of variables and constraints of the OP model in three case studies revealed the higher degree of complexity of the stochastic model. The sizes of MP is not dependent on and SP is nearly | || || | while the size of original problem strongly depends on the time 21

ACCEPTED MANUSCRIPT horizon and is the major reason of high computational solution time . The number of integer variable of the OP for CS3 is 661 times greater than the number of integer variable of MP. Moreover, when the time horizon ranges from 10 days to 31 days, the number of integer variable of the MP grows less than 3 times; however, it grows more than 7 times in OP.

Var.

OP Int. var.

Con.

109568 374024 780056

16870 59535 125590

293892 1053652 1896531

CR IP T

Table 3: The size of MP, SPR, SP and OP for each case MP SPR SP Case Var. Int. Con. Var. Int. Var. Con. var. var. CS1(10days) 142 70 474 99355 16800 10071 10070 CS2(20days) 272 135 895 353617 59400 20136 20135 CS3(31days) 382 190 1256 748483 125400 31191 31190

AN US

Table 4 compares the results obtained from three-step algorithm with the findings of the direct method for the original two-stage model. As observed, only the first case, CS1 reached the optimal solution in several minutes by using direct method. In this case, the required CPU time to find the optimal solution is 2471.42s. For other cases, it was not possible to obtain optimal solution in a limited time (5 hours). On the other hand, all instances reached the optimal solution in several minutes by using the three-step algorithm and it converges quickly after a few iterations. This could be due to improving MP by adding VIs (see Section 4-5) that push the MP solution toward optimality. Table 4: Comparing the result of the three-step algorithm with the findings of the direct method

CS1 CS2 CS3

684.35 782.43 831.12

23 28 37

Relative gap (%)

M

Time (sec)

112352 154231 333157

ED

Prob.

Three Steps Method Number of Objective Iteration value

0.04 0.01 0.03

Time (sec)

2471.42 -

Direct Method Objective Relative value gap (%) 112304 197137 2788672

0.00 66.75 94.13

AC

CE

PT

It can be seen from computational results that the algorithm performed faster than the direct method. The experimental results show that, for a larger number of scenarios, the proposed algorithm can perform faster on average than the direct method on CS1. In addition, increasing time horizon resulted in growing the problem size as well as the solution time. For the direct method, the size of the model and solution time depended strongly on the number of days within the time horizon. Note, the size of MP and SP as well as the solution time is not greatly affected by the length of the time horizon and it is possible to be extended to longer time horizon. In all the cases, initial value of ̅ was considered 0.5. The value of ̅ cannot be set equal to 1 at first, because the optimal solution may be obtained in a lower level due to the initial inventory and the volume of product in the pipeline at the beginning of the time horizon. The initial value of ̅ can be calculated by

∑

̅̅̅̅̅ ∑

̅̅̅̅̅ ̅̅̅̅̅̅̅

.

The value of ̅ was updated in next iterations by (91). The variation of ̅ for the case study CS3 and problem bounds (lower/upper) for CS3 at 20 iterations are shown in Figure 5. It can be seen that, by increasing the number of iterations, the value of ̅ increased. Moreover, with the growth of the value of ̅ , the rate of change of the lower bound on the upper bound increased to 1 during the algorithm convergence. The MP becomes stronger with this cut, a better lower bound is obtained, and the algorithm converges faster.

22

ACCEPTED MANUSCRIPT

1.2

Volume

1 0.8 0.6 0.4 0.2 0 2

3

4

5

6

7

8

9

10

11

12

Iteration μ

LB/UB

13

14

15

16

17

18

19

20

CR IP T

1

Figure 5: Variation of ̅ and objective function bounds in the case CS3.

ED

M

AN US

Hence, TSSP has the reputation of being computationally difficult. Solving a simplified version like the single scenario problem (SSP) is a natural temptation when faced with reallife problems. In order to show the importance of using TSSP model, we solve CS1 by using the model TSSP considering all scenarios as well as single common scenario (normal scenario). Let ̅ and ̅ be the first stage decisions corresponding to the normal scenario. Then, we solve all 100 scenarios problems individually by fixing the first stage variables to ̅ and ̅ . Finally, we calculate the average shortage for each day. Total volume of shortage in each day for both the TSSP and SSP problems is addressed in Figure 6. It is obvious that, in almost all the days, the volume of shortage of TSSP is less than the average shortage of the SSP problems. The average shortage is 7130 for TSSP, while this value is 10197 for SSP problems. shortage(TSSP)

1.8

PT

1.6 1.4

CE

1.2 1

0.8 0.6

AC

VOLUME (*10000)

shortage(SSP)

0.4 0.2

0

1

2

3

4

5

6

7

8

9

DAY

Figure 6: Comparison of the TSSP and SSP performance in terms of cumulative shortage.

23

10

ACCEPTED MANUSCRIPT 5-2 Computational result of the SAA method

Table 5: The result of the SAA method for the three cases

Time (sec) 2341 2892 3781

AN US

Prob. CS1 CS2 CS3

SAA Method Objective value 111459 155732 332145

CR IP T

The computational study discussed in the previous section shows that the synergy of the improved BD approach expedites the computational process significantly; thus, allowing larger, more realistically sized problems to be solved. In this section, the results of the cases CS1, CS2 and CS3 represent the execution of the SAA method. We choose to solve all the cases using , SAA problems (SSAP model) with sample sizes of . The CPU-time is increasing in the number of scenarios in SAA problem. Moreover, the growth of the computational process increases substantially with the size of the problem. Table 5 gives objective values and optimality gaps, as previously observed. Our computational results indicate that the proposed method is successful in solving problems with up to 1000 scenarios within a maximum 4% of optimality gap. The value of upper and lower bound for CS1 is shown in Figure 7 for each iteration.

Relative gap (%) 3.81 3.78 3.96

M

As a result, we are able to find provably near-optimal solutions to these difficult stochastic programs using only a moderate amount of computation time. Integrating this scheme within the SAA method allows us to solve problems with a large number of scenarios.

ED

120

110

CE

100

PT

Volume(*1000)

130

90

200

500

800

1000

1100

Number of scenarios

AC

50

Figure 7: Upper and lower bound of the objective function in the SAA method for CS1

6

Conclusion

In this paper, we have addressed the straightforward multiproduct scheduling problem under uncertain flow rate. Our main goal was to yield a new formulation and new solution approach to obtain an optimal solution in a reasonable time. The two-stage stochastic MILP model developed takes into account the inventory management, integrates it into the final solution and represents the real world synergies between the tank farm product needs and pipeline operation. 24

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

We developed a practical solution method to solve two-stage stochastic programming pipeline problems. Considering the computational challenges inherent in solving the twostage stochastic program using the standard sample average approximation, an accelerated stochastic solution method (a three-step approach based on BD method) is developed to solve the two-stage stochastic program. The three-step approach takes advantage of problem structure. This reduces the number of binary variables and logical constraints and generates a favorable structure that can be exploited in the BD algorithm. Finally, a real-world case study for three time horizons is tackled using the proposed algorithm and the optimal detailed schedule is found at reasonable CPU time. We expect that the proposed methodology can be extended to even more complex configurations, like pipeline networks with a tree structure. However, developing an efficient solution algorithm for solving this complicated problem is still a big challenge for future research. Moreover, the three-step approach can be considered for solving other optimization problems and applications that follow the same structure. In practice, we face with cases in which more than one disruption occurs during the time horizon. Therefore, as a feature work, it is possible to take into account scenarios with more than one disruption.

25

ACCEPTED MANUSCRIPT 7 Bibliography Benders, J. (1962). Partitioning procedures for solving mixed-variables programming problems. Numerical Mathematic, 4(1), 238-252. Cafaro, D. C., & Cerdá, J. (2008). Efficient Tool for the Scheduling of Multiproduct Pipeline and Terminal Operations. Industrial and Engineering Chemical Research, 47(24), 9941-9956. Cafaro, D., & Cerdá, J. (2009). Optimal Scheduling of refined products Pipelines with multiple sources. Industrial and Engineering chemistry Research, 48, 6675-6689.

CR IP T

Cafaro, D. C., & Cerdá, J. (2010). Operational Scheduling of Refined Products Pipeline Networks with Simultaneous Batch Injections. Compters and Chemical Engineering, 34, 1687-1704. Cafaro, D., & Cerdá, J. (2012). Rigorous scheduling of mesh-structure refined petroleum pipeline networks. Computers and Chemical Engineering, 38, 185-203.

AN US

Cafaro, D., & Cerdá, J. (2014). Rigorous formulation for the scheduling of reversible-flow multiproduct pipelines. Computers & Chemical Engineering, 61, 59-76. Cafaro, V., Cafaro, D., Méndez, C., & Cerdá, J. (2015). MINLP model for the detailed scheduling of refined products pipelines with flow rate dependent pumping costs. Computers and Chemical Engineering, 72, 210-221.

M

Castro, P. M. (2010). Optimal Scheduling of Pipeline Systems with a Resource-Task Network. Industrial & Engineering Chemistry Research , 49(22), 11491-11505.

ED

Castro, P. M., & Mostafaei, H. (2017). Product-centric continuous-time formulation for pipeline scheduling. Computers & Chemical Engineering, 104, 283-295.

PT

Codato, G., & Fischetti, M. (2006). Combinatorial benders cuts for mixed-integer linear programming. Operations Research, 54(4), 756-766.

CE

Dawotola, A., Gelder, P., Charima, J., & Vrijling, J. (2011). Estimation of failure rates of crude product pipelines. Applications of Statistics and Probability in Civil Engineering. London.

AC

Galereh, S., Moneni, R., & Nichel, S. (2015). Multi-priod hub location problems in transpoetation. Transportation research Part E: Logistics and Transportation Review, 75, 6794. Geoffrion , A. (1972). Generalized Benders Decommposition. Journal of Optimization Theory and Applications, 10(4), 237-260. Graham, C., & Talay, D. (2013). Stochastic Simulation and Monte carlo Methods. New York: Springer. Hane, C., & Ratliff, H. (1995). Schedulinginputs to multi-commodity pipeline. Annals of Operations Research, 57(1), 73-101.

26

ACCEPTED MANUSCRIPT Hitoshi Tsunoda Meira, W., Magatão, L., Relvas, S., Barbosa-Póvoa, A., Neves Jr, F., & Lúcia, V. (2018). A matheuristic decomposition approach for the scheduling of a singlesource and multiple destinations pipeline system. European Journal of Operational Research, 268(2), 665-687. Hooker, J. N. (2007). Planning and scheduling by logic-based Benders decomposition. Operations research, 55(3), 588-602. Hooker, J. (2011). logic-based method for optimization: combining optimization and constraint satisfaction. Canada: John Wiley and Sons.

CR IP T

Hooker, J. N., & Ottosson, G. (2003). Logic-based Benders decomposition. Mathematical programming, 96(1), 33-60. Kent Muhlbauer, W. (2004). Pipeline Risk Management Manual: Ideas, Techniques and Resources. USA: Elsevier.

AN US

Kleywegt, A. J., Shaprio, A., & Homem De Kello, T. (2001). The sample average approximation method for stochastic discrete optimization. Society for Industrial and Applied Mathematic, 12, 479-502. Magnenti, T., & Wong, R. (1981). Accelerating Benders decomposition: Algorithm enhancment and model selection criteria. Operations Research , 29(3), 464-484.

M

Mak, W., Morton, D., & Wood, R. (1999). Monte Carlo bounding techniques for determining solution quality in stochastic programs. Operations Research Letters , 24, 47–56.

ED

MirHassani, S., & BeheshtiAsl, N. (2013). A heuristic batch sequencing for multiproduct pipelines. Computers and Chemical Engineering, 56, 58-67. MirHassani, S., & Ghorbanalizadeh, M. (2008). The Multi-Product Pipeline Scheduling System. Computers and Mathematics with Applications, 56(4), 891-897.

PT

MirHassani, S., & Jahromi, H. (2010). Scheduling Multi-product Tree-Structure Pipelines. Computers and Chemical Engineering, 35(1), 165-176.

CE

Mohammadi Bidhandi, H., & Patrick, J. (2017). Accelerated sample average approximation method fortwo-stage stochastic programming with binary first-stage variables. Applied Mathematical Modelling, 41, 582–595.

AC

Moradi, S., & MirHassani, S. A. (2016). Robust scheduling for multi-product pipelines under demand uncertainty. The International Journal of Advanced Manufacturing Technology, 87(9), 2541-2549. Naoum Sawaya, J., & Elhedhli, S. (2010). A nested Benders decomposition approach for telecommunication network planning. Naval Research Logistics (NRT), 57(6), 519-539. Papadakos, N. (2008). Practical enhancment to the magnenti-wong method. Operation Research letters, 36(4), 444-449. Polli, H. L., Magatão, S. N., & Neves, J. (2017). Collaborative Approach Based on Heuristic Algorithm and MILP Model To Assignment and Sequencing of Oil Derivative 27

ACCEPTED MANUSCRIPT Batches in Pipeline Networks. Industrial Engineering & Chemistry Research, 56(9), 24922514. Rejowski, R., & Pinto, J. (2003). Scheduling of A Multiproduct Pipeline System. Computers and Chemical Engineering, 27(8-9 ), 1229-1246. Relvas, S., Barbosa-Póvoa, A., & Matos, H. (2009). Heuristic Batch Sequencing on A Multiproduct Oil Distribution System. Compters and Chemical Engineering, 33, 712-730.

CR IP T

Relvas, S., Boschetto Magatão, S., Barbosa-Póvoa, A., & Neves Jr, F. v. (2013). Integrated scheduling and inventory management of an oil products distribution system. Omega, 41, 955-968. Relvas, S., Matos, H., Barbosa-Póvoa, A., Fialho, J., & Pinheiro, A. (2006). Pipeline Scheduling and Inventory Management of A Multiproduct Distribution Oil System. Industrial & Engineering Chemistry Research, 45, 7841-7855.

AN US

Relvas, S., Matos, H., Barbosa-Póvoa, A., Fialho, J., & Pinheiro, A. (2007). Reactive Scheduling Framwork for A Multi-Product Pipeline with Inventory Management. Industrial & Engineering Chemistry Research, 46, 5659-5672. Saharidis, G., Boile, M., & Theofanis, S. (2011). Initialization of the benders master problem using valid inequalities applied to fixed-charge network problems. Expert Systems with Applications, 38(6), 6627-6636.

M

Shapiro, A. (1996). Simulation-based optimization: Convergence analysis and statistical inference. Communications in Statistics. Stochastic Models, 12, 425–454.

ED

Shapiro, A., & Homem de Mello, T. (1998). A simulation-based approach to two-stage stochastic programming with recourse. Mathematical Programming, 81, 301–325.

PT

Taherkhani, M., Tavakkoli-Moghaddam, R., Seifbarg, M., & Fattahi, P. (2017). Detailed Scheduling of Tree-like Pipeline Networks with Multiple Refineries. International Journal of Engineering, 30(12), 1870-1878.

CE

Tang , L., Jiang, W., & Saharidis, G. (2010). An improved Benders decomposition algorithm for the logistics facility location problem with capacity expansions. Annals of Operation Research , 210(1), 165-190.

AC

Van Slyke, R. M., & Wets, R. (1969). L-Shaped linear programs with applications to optimal control and stochastic programming. Society for Industrial and Applied Mathematics: SIAM , 17(4), 638-663. Yu, G., & Qi, X. (2004). Disruption managment Framework, Models and Applications. USA: World Scientific Publishing. Zhang, H., Liang, Y., Xiao, Q., Wu, M. Y., & Shao, Q. (2016). Supply-based Optimal Plan Making for Scheduling of Oil Product Pipeline. Petroleum Science, 13(2), 355-367. Zhang, J., & Ponnambalam, K. (2006). Hydro energy optimization in a deregulated electricity market. Optimization and Engineering, 7(1), 47-61.

28