Collaborative Approach Based on Heuristic Algorithm and MILP Model

Subscriber access provided by University of Newcastle, Australia

Article

Collaborative Approach based on Heuristic Algorithm and MILP Model to Assignment and Sequencing of Oil Derivative Batches in a Pipeline Network Helton Luis Polli, Leandro Magatão, Suelen Neves Boschetto Magatão, Flavio Neves, and Lúcia Valéria Ramos de Arruda Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b03516 • Publication Date (Web): 03 Feb 2017 Downloaded from http://pubs.acs.org on February 7, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Collaborative Approach based on Heuristic Algorithm and MILP Model to Assignment and Sequencing of Oil Derivative Batches in Pipeline Networks Helton L. Polli,† Leandro Magatão,∗,‡ Suelen N. B. Magatão,† Flávio Neves-Jr,‡ and Lucia V. R. Arruda‡ Sollve Optimization and Technology, and Federal Univesity of Technology - Paran´ a (UTFPR), Graduate Program in Electrical and Computer Engineering (CPGEI) E-mail: [email protected] Phone: +55 41 33104701

Abstract This paper presents a collaborative approach to the assignment and sequencing of batches in pipeline networks. The approach is based on the integration of a heuristic algorithm with a Mixed Integer Linear Programming (MILP) model. The pipeline scheduling problem is solved using a hierarchical decomposition [Ind. Eng. Chem. Res. 2015, 54, 5077], 1 but a new collaborative approach is proposed for Assignment/Sequencing tasks. At a first moment, the proposed heuristic algorithm (Assignment module) determines priorities for sending batches in order to respect deadlines. ∗

To whom correspondence should be addressed Sollve Optimization and Technology, Curitiba, PR, Brazil, 80220-060 ‡ Federal Univesity of Technology - Paraná (UTFPR), Graduate Program in Electrical and Computer Engineering (CPGEI), Curitiba, PR, Brazil, 80230-901 †

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The algorithm encompasses an analysis of production and demand plans, inventories, input and output of products in terminals, trying to use resources, namely tanks and pipelines, in an optimized form. This algorithm is used in cooperation with a proposed MILP Sequencing model, which allows overcoming computation difficulties previously indicated by a traditional scheduling approach that tried to aggregate into the same monolithic MILP model assignment and sequencing decisions [Ind. Eng. Chem. Res. 2012, 51, 4591]. 2 The proposed assignment/sequencing collaborative approach can be used to define operational batches with their volumes and routes in pipeline networks. Thus, the lot-sizing problem of batches in pipeline networks is addressed within the proposed paper. Tests were made in pipeline networks of different topologies. First, a small, but representative pipeline network is proposed and a data set for this network is made available for reproducibility purposes. Secondly, tests are made in a real-world pipeline network and results have been attained in computational times from seconds to few minutes.

1

Introduction

1.1

Introduction to the Cluster of Articles

Pipeline networks are an efficient way to oil derivatives and gas transport 3 and making good use of this transportation is crucial to the petroleum industry, where environmentally safe and economical transport solutions are sought. Thus, the optimized operational management of oil derivatives transport in pipeline networks is a decisive factor: a small set of operational decisions about pipeline operations can imply a significant transportation cost variation. The transportation decisions are influenced by the refineries’ production and consumption, variability in depots’ storage levels, pipeline usage, among other factors that depend on product distribution. Nowadays, the operational decision making for pipelines is based on the experience of experts and systems that only check the operations’ consistency. 4 However, a growing interest 2 ACS Paragon Plus Environment

Page 2 of 63

Page 3 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the oil industry has motivated the development of tools that are able to aid the decision making process, particularly the ones that apply optimization techniques in order to use the resources in an efficient, safe, and profitable way. Some optimization techniques to short term scheduling of pipeline activities can be found in literature as in refs 5–12. In particular, in ref 1 a taxonomy for pipeline-scheduling papers is given; important contributions from 2004 to 2014 are indicated. More recent contributions are presented by refs 13–16. In ref 15, for instance, the authors proposed a continuous time monolithic model which is able to handling the input/output operations of nodes that can send/receive material into/from pipelines. References 17 and 18 applied a single MILP model to solve hypothetical pipeline networks. In fact, short term pipeline-scheduling problems are difficult to be solved and ref 19 proved that they are NP-hard combinatorial optimization problems. In this way, many authors have used heuristic procedures/decomposition techniques 2,4,20–22 intending to circumvent the problem combinatorial complexity, and enabling the application for pipeline networks. Reference 23 also proposed a heuristic methodology. Reference 22 used a MILP model combined with post-processing heuristic to solve the scheduling problem to transport heavy oil products in a real-world pipeline network. A real-world pipeline network to transport oil derivatives was also (partially) solved with different techniques, such as: Heuristics and CP; 21 single MILP model for the planning phase; 24 and, Micro-genetic algorithm and MILP model. 4

1.2

The Main Contributions of the Paper

In this work a hierarchical decomposition approach 1 is used to solve pipeline-scheduling problems. Three main blocks are used for the entire solving strategy: the Planning block, 2 the Assignment/Sequencing block, and the Timing block. 1 The herein proposed paper addresses the solution of the Assignment/Sequencing block, which plays a fundamental role to obtain the final pipeline network scheduling. The new collaborative approach here proposed 3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to the Assignment/Sequencing is composed of the integration of a heuristic algorithm with a MILP model. Reference 2 formerly proposed a monolithic MILP model to solve in a unified framework the assignment-sequencing features for a real-world pipeline-scheduling problem. The proposed monolithic model spent several hours to present results with an elevated gap for pipeline-network scenarios with only a week of horizon. Thus, the referenced work indicated that to obtain the assignment-sequencing solution for pipeline-networks with a unique MILP model can be a hard computational task. The addressed scheduling problem is combinatorial and the decomposition approach herein presented allows overcoming computation difficulties indicated by the former work, 2 dealing with the assignment problem by a heuristic procedure and, in a following step, addressing the sequencing problem with a MILP model. To address the pipeline network scheduling in a non-prohibitive computational time (e.g., minutes) is a fundamental issue, mainly for real-world applications, where scheduling decisions of different nodes/pipelines can be highly interconnected 1 and the company needs a consensus scheduling solution for the entire pipeline network management. Based on the highlighted topics, the paper novelty can be summarized as follows: (i ) A new collaborative approach is proposed to solve the assignment and the sequencing tasks for pipeline-scheduling networks, a hard combinatorial problem; (ii ) The novel proposed approach was able to solve hypothetical and real-world instances that were not solved by a former literature approach based on a monolithic MILP model. 2 The remainder of the manuscript is organized as follows. The next section describes the main features of pipeline network scheduling problems and the proposed hierarchical structure to solve this class of problem. The heuristic algorithm for the assignment module is explained. Following, the proposed base MILP sequencing model and its extensions are detailed. Results and discussions include two case studies. The first one involves the solution of a small pipeline network in order to give insights to the reader about the problem main features. A set of data for this small case is provided in the Supporting Information for


Page 4 of 63

Page 5 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reproducibility purposes. A second case study involves the solution of a real-world pipeline network. In this case, results about each proposed block are presented and discussed and a comparison with a previous literature pipeline-sequencing MILP monolithic model for this network 2 is made. The final conclusions are then presented.

2

Problem Description

The proposed approach is applied to the solution of pipelines networks. Figure 1 presents an example of a pipeline network that includes six pipelines (D1 to D6) linking five nodes (N1 to N5). The nodes include two depots (N1 and N2), two refineries (N3 and N4), and a harbor (N5). Eight oil derivatives and ethanol can be transported in this network. Different colors indicate different products. At this example, in-transit batches are illustrated into the pipelines D1 to D6, namely, batches 27 (originally sent from N4 to N2), 28 (from N4 to N2), 63 (from N3 to N1), 64 (from N2 to N1), 65 (from N3 to N2), 66 (from N3 to N2), 73 (from N2 to N1), 74 (from N5 to N1), and 75 (from N5 to N2). Batch 74, that was originally sent from N5 to N1, is split in pipelines D6/D1, and batch 63 (from N3 to N1) is split in pipelines D3/D2. Within this small network, pipelines D4 and D6 can have the flow direction reverted.

Figure 1: Example of a Pipeline Network



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Pipeline network scheduling problems generally present the following characteristics: 20 • Each pipeline has a particular volume and always operates filled; • A batch remains in a pipeline until another batch “pushes” it; • Batches that are in the pipeline at the beginning of the scheduling horizon must be routed to their previously determined destination (in-transit batches); • Each origin area has upper and lower flow rate limits to pump a batch; • Each pipeline has upper and lower flow rate limits to transport a batch; • Each area has a tank farm that stores different products; • Each product should be stored in a different tank (or set of tanks) during the scheduling horizon; • The upper and lower tankage limits to the overall product inventory in each node should be respected; • Volumes of batches can be smaller or bigger than the volume of pipelines, but should respect the tankage limits on origin and destination areas; • The flow direction of predefined pipelines can be reverted, according to operational necessities; • Auxiliary batches can be used to manage reversions in flow directions of pipelines, since batches have to be dislocated (pushed) to the destination node before the direction of flow changes; 1 • Plug products can be used between two products that should be transported in sequence but are incompatible (in order to avoid prohibited interfaces); • The harbor can receive or send products, according to the vessels scheduling; 6 ACS Paragon Plus Environment

Page 6 of 63

Page 7 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


• The production/consumption rates can be approximated by a constant value during the entire horizon, depending on the product/area. However, stepwise production/consumption functions during smaller periods inside the horizon can be more suitable to represent production/consumption rates that vary along the scheduling horizon. For instance, campaigns in refineries or a vessel arriving in a harbor that generate a specific demand during a predefined period, typically one day; • Batches can be sent from an origin node to a destination node “passing through” intermediate nodes. In each intermediate node, a batch can be “aligned” to another pipeline and pumped. At this case, no intermediary tankage is demanded. Alternatively, intermediate storage can be used and a “surge tank operation” defined. This operation occurs when a batch is received in a tank of an intermediate area and has to be sent to another area at a different flow rate; • A limited set of pumps and manifolds exist in each area. Thus, only a restricted number of batches are allowed to be sent/received at the same time from/into each area. This fact characterizes “local constraints”; • The electrical energy has different costs, according to season and period of day. Therefore, at on-peak demand hours, the pumping should be periodically interrupted in some operational areas during the scheduling horizon; • Each batch can be linked with a specific route inside the network; • Each route is composed by a sequence of areas intercalated by pipelines; For instance, N5→D6→N2→D1→N1 (alternatively, N5→6→N2→1→N1); thus, the origin is node N5 and the final destination is N1; • Preferential routes for each product can be used.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.1

Page 8 of 63

Problem Hierarchical Decomposition

Pipeline scheduling problems are difficult to solve and can be categorized as NP-Hard problems. 19 In order to address the computational issues of pipeline networks, ref 20 proposed decomposition strategies to address a real-world pipeline-scheduling problem. More recently, these decomposition strategies have been modified. 1,2 Within the proposed paper, a decomposition approach presented by ref 1 is improved, as indicated in Figure 2. Planning Block (Magatão et al., 2012)

 Choice of routes and global volumes  Minimization of pipelines in reverse flow  Evaluation of usage rate of pipelines

Assignment / Sequencing Block Assignment module

MILP sequencing model

Definition of batches and time-windows

Definition of batch ordering

Simulation block (Boschetto et al., 2010)

Timing Block (Magatão et al., 2015)

MILP timing model - 2 phases

MILP timing model - single phase Contain all the constraints

MLC

MST

Final Solution

Figure 2: Hierarchical Decomposition Approach (Adapted from ref 1) The blocks’ functionalities have been already detailed in ref 1. In essence, the planning block 2 uses input data based on forecasting of production and consumption areas and it analyses the different products comprised within a scenario of, typically, one month. This block suggests the products’ transportation characteristics from producing areas with the goal of supplying the depots’ demand, final clients, and harbors. After running, the planning 8 ACS Paragon Plus Environment

Page 9 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model provides, as a final result, the routes to be used and the global volume to be transferred, while considering the problem operational constraints. These results are fed into the next block, the assignment/sequencing block. The herein proposed assignment/sequencing block is a procedure that determines parameters for each pumped batch and it indicates a sequence that batches should be pumped from each area. In a first step (Assignment module), the heuristic assignment algorithm is used to compute batch volumes, taking into account typical values for the range of operational batches, production and consumption information, and storage levels for each area. The algorithm also determines time-windows for each batch in origin/destination areas. Then (Sequencing model), the MILP sequencing uses such temporal intervals to determine the optimized sequence (order) of operational batches. In section 3 the execution flow for the proposed assignment/sequencing block is detailed (e.g., Figure 3). Thus, a list of batches, including the route, batch volume, and the time-windows for each batch is computed by this collaborative approach. As a result, pumping orders suggested by the assignment/sequencing block are respected by the following blocks. It is important to emphasize that the presented paper improves this specific block within the proposed methodology. In ref 20 just heuristic procedures were adopted. In ref 2 a monolithic MILP model was tested for integrating assignment/sequencing decisions, but the computational load was a concern. The presented paper address the previously detected computational difficulties 2 by the cooperation of an assignment module and a MILP sequencing model. Thus, lot sizing and sequencing decisions for pipeline network scheduling problems are determined by the proposed collaborative approach. The simulation block uses information derived from the assignment/sequencing block and evaluates, step by step, the batch that is being pumped. This block also analyses the influence of this pumping operation on the receiving of other batches. As a result, the batch movement along the pipeline network is analyzed, and “volumetric parts” (portions of this batch or sub-batches) are identified. 20 In this way, the simulation block provides valuable



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

information, used as input parameters to the MILP timing model. All the parameters calculated by the previous blocks are used in a continuous time MILP model, the timing block. 1 This block determines the operational short-term scheduling for the entire pipeline network using the parameters (e.g., time-windows and sub-batch volumes) determined by previous blocks. The blocks proposed in Figure 2 are integrated in a collaborative way and are sequentially executed. Thus, the output of a block can be used as an input to the next block. Additionally, the proposed approach of Figure 2 is embedded in a customized computational tool that allows visualizing reports, inventory graphics, and Gantt charts related to proposed solutions. It is important to emphasize that the proposed article presents a collaborative approach to the assignment/sequencing of batches based on the integration of a heuristic algorithm with a MILP model. This cooperation allows overcoming computation difficulties previously indicated. 1,2 Thus, the hierarchical decomposition approach previously presented 1 is herein improved and modified, as indicated in Figure 2, and can be used to solve pipeline network scheduling problems.

3

The assignment methodology

The proposed assignment methodology determines the batches that are pumped based on production and consumption forecasts and network transport characteristics. The assignment module is developed considering the time-windows’ propagation concept, and it is implemented in Java programming language. The module is responsible for the creation of batches to be transported by the pipeline network, determining the following characteristics: the used route, the batch volume for each product, and the calculated time-windows. As a consequence, the module also suggests an initial order for pumping batches. Figure 3 illustrates the three main steps of this module: pre-assignment, batches’ generation, and post-assignment/sequencing. The MILP


Page 10 of 63

Page 11 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sequencing model is executed between batches’ generation and post-assignment/sequencing steps. As hereafter detailed, the post-assignment/sequencing should be made after the MILP sequencing model because the ordering suggested by the model influences the number and position of batches generated in this step. Figure 3 is an overview diagram for the entire assignment/sequencing proposed approach. The main inputs of blocks are mentioned in this figure. Each block is detailed in the next subsections.

Pre-assignment

In pipeline batches Batches generation

Batches with time-windows

Batches sequencing

Sequenced batches Post-assignment/sequencing

Figure 3: Execution flow for assignment/sequencing blocks

The batch is computed based on storage management issues in each depot. This procedure is done separately for each product, assuming that the pipeline network is dedicated to this product. The storage situation in each depot is periodically verified, considering the analyzed product. This is done by monitoring a storage curve, using production and demand data. This abstraction allows to infer the product storage level during the studied scenario and to represent stepwise production and demand behaviors. In harbors, for instance, stepwise demand/production “curves” are quite common, where a vessel can arrive to be loaded or unloaded. Thus, the production/consumption occurs through curve peaks (steps), and not in a constant rate along the month. The algorithm works during all considered time intervals with an accurate representation of storage conditions in each depot. This fact results in a batch portfolio adherent to the reality. Time-windows are calculated by the assignment methodology for each one of the three storage levels, as illustrated in Figure 4. The first one, capacity level (“CAP ”), is comprised between the operational zero and the tankage capacity storage. The second, Min-Max storage level (“MnMx ”), is comprised between the minimum and maximum storage levels. Finally,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 63

the third, the target storage level (“Target”), is the storage level comprised between the minimum target storage and the maximum target storage. In this way, the assignment block generates three time-windows for each batch and node. Each time-window refers to a particular storage level within the considered node. Figure 4 also indicates typical volumetric bounds used for time-windows’ generation. These values were suggested by ref 2. In fact, such bounds can be tuned according to each product and area in order to better suit operational conditions. Volume (v.u.)

Capacity (Cap.) Maximum (~80% to 95% of Cap.)

Maximum target (~75% of Cap.) Target level (“Target”)

Min-Max level (“MnMx”)

Capacity level (“CAP”)

Minimum target (~25% of Cap.) Minimum (~5% to 20% of Cap.) Operational zero

Time (t.u.)

Figure 4: Storage levels definition In a following step (Sequencing model), the proposed MILP sequencing uses time-windows determined by the assignment methodology. So that, the proposed model manages the time-windows’ violations in origin and destination areas. Figure 5a illustrates the sending time-window generation on an origin area (for Min-Max level), comprising the time interval between TED (available sending time) and TEC (critical sending time). TED is determined as the time instant in which an origin area reaches a sufficient storage level to send a complete batch; TEC, as the time instant that the tankage in origin area reaches the maximum allowed volume. Similarly, the receiving time-window in destination area, illustrated by Figure 5b, is comprised between TRD (available receiving time) and TRC (critical receiving time). TRD is defined as the time instant that the tankage in destination area reaches the allowed volume to store a complete batch; TRC, as the time instant that the tankage in destination area reaches the minimum allowed volume. Therefore, the algorithm creates a new time-window 12 ACS Paragon Plus Environment

Page 13 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in origin area between TED and TEC, indicating the time interval wherein a product sending can and must occur. Another time-window is created in destination area between TRD and TRC, indicating the time interval wherein a product receiving can and must occur. Thus, to determine the batch assignment in origin/destination areas, key time values are derived from time-windows.

(b) Destination area

(a) Origin area

Figure 5: Time-windows calculus 20

3.1

Assignment pre-processing (Pre-assignment)

As indicated in Figure 3, this block is just a pre-processing step for the batch calculus, where the key decisions about origin-destination choices and volume allocations are, in fact, made. Figure 6 illustrates the assignment pre-processing, or pre-assignment execution flow. Initially, the ordering of “starting batches” is determined, that is, batches already in pipeline at the scheduling horizon beginning. This ordering occurs according to the position and the route previously linked to each starting batch. Following, the initialization of each pair “area-product” is carried out. Each area-product pair is obtained through the initial storage in each area and based on information provided by the Planning Block of Figure 2. For instance, consider an “A” area with seven tanks: two tanks with initial storage of P1 product, three tanks with P2 product, and two tanks with P3 product. We have the following relations: A-P1, A-P2, and A-P3, with aggregated capacity for each pair. The aggregated capacity is the sum of each tank capacity and the aggregate initial storage is the sum of each



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tank opening stock. If, for instance, P4 and P5 products do not exist within a considered scenario, the pair A-P4 and A-P5 are not generated because there is neither a tank in A area with initial storage of P4 nor P5. In this data generation methodology, the undesirable indexes are not created and, consequently, the execution time is reduced. The pre-assignment also comprises the generation of initialization batches. The procedure involves all starting batches and “scheduled batches” in the format required by the MILP model, with the time-windows’ information. In practice, there is the possibility of existing already scheduled batches at the horizon beginning. This fact is addressed within the initialization batch generation procedure. To finalize, the assignment pairs generation procedure makes the entire set of date consistent with the format required by the batch generation algorithm, the following block (Figure 7). In pipeline batch ordering

Starting pairs area-products

Assignment pairs generation

Initialization batches generation

Figure 6: Pre-assignment execution flow

3.2

Batch generation algorithm

The batch calculus algorithm was created based on the specialists’ working process. Figure 7 illustrates the algorithm flow. Basically, it represents a cyclic process were each cycle generates a new batch to attend a necessity, for instance, a demand of an area or a production pumping, until satisfying all necessities of each considered area. In each cycle, four main processes are executed: computing critical times, choosing a necessity area, choosing a pair, and batch generation.


Page 14 of 63

Page 15 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Assignment pairs generation

Critical times calculus Critical times Storage necessity is attempted?

yes

no Necessity area select Necessity area Choose pair Assignment of pair Batch generation Batch Assignment pairs generation

Figure 7: Batch generation algorithm flow 3.2.1

Computing critical times

The critical time computing process, illustrated in Figure 8, identifies the time limit in which each area needs to receive a product (in demand necessity case), or to send a product (in case of production pump necessity). Thus, a necessity area can be either a sending or a receiving area. However, an area can receive from more than one origin, or send to more than one destination, and it is necessary to check the feasibility of critical times. Figure 8 illustrates that the critical time computing is a cyclic process in which all pairs, of a considered area, are evaluated. The choice of a practical volume inside an operational range is executed within this process, based on time-windows propagation. It is verified if the time-windows are viable, considering the time shift from origin area until destination area. This time lag is estimated by the route volume by the maximum flowrate. If the time-window is feasible, then it is verified if the computed critical time is lower than the previously stored critical



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

time. It is important emphasizing that the area is started with a high enough critical time value to the first analyzed pair to be chosen. At the end of all pairs’ rating, the in analysis area has the volume of pumping chosen and the critical times recorded. An operational volume is chosen based on the route’s maximum and minimum viable volumes. According to ref 2, volume of batches should be compatible with storage constraints. Thus, an operational range for lot-sizing purposes has to be respected. The cyclical process starts with the highest possible volume inside the range that can be assigned to a batch. Therefore, the first tested volume in the “Choose volume in a range” block is the maximum limit. In sequence, for this volume, critical times are calculated and it is verified whether they are practicable or not, according to time-windows propagation. If so, the volume is chosen and critical times are recorded; otherwise, the tested volume is decremented to other reference in the viable volume range (e.g., a reduction of one thousand cubic meters is supposed and tested). Those verifications occur following a cyclical process involving the choice of pairs/volumes.

3.2.2

Choosing the necessity area

The area with the lowest critical time previously computed is selected as the in analysis “necessity area”. Either a sending area or a receiving area can be a candidate to the considered necessity area. The next step is to choose a pair that can better attend the necessity of the chosen area.

3.2.3

Choosing the pair

Figure 9 illustrates the execution flow process for choosing a node to attend the “necessity area”. Firstly, it is identified the necessity kind – receiving or sending. In case of receiving, a search for a product supplier area is made; otherwise, a search for a product receiver area is made. An average dislocation time (or delay time) from the origin to the destination of each evaluated pair is considered. Again, we can notice that the node choice is a cyclic process,


Page 16 of 63

Page 17 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Critical times reset

All pairs evaluated?

yes

no Choose pair

yes

Assignment bigger than planned volume?

no Choose volume in a range

Time-windows calculus

TED calculus

TRD calculus

TEC calculus

TRC calculus

Viable time-windows?

Update volume limits

no

yes Critical times record

Figure 8: Critical times computing process because all pairs are evaluated to choose an option. The main goal is to choose a suitable node to attend the necessity area considering time-windows’ values.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Identify necessity

Choose pair

Evaluate choosing criteria yes no

More options? no

Best choice until now? yes

Choose pair

Figure 9: Choosing a node to attend a necessity area 3.2.4

Batch generation

This process step, in fact, computes the batch characteristics. All information about batches was already defined, such as: origin, destination, product, volume, route, and time-windows. From this moment, the solution process does not consider the aggregated tankage capacity, but it evaluates the time-windows associated with each batch. The storage managing is hereafter made in a temporal form analysis.

3.3

A brief comment on the assignment methodology

The proposed assignment methodology (heuristic algorithm) allows determining all batches to be sent/received and their parameters. Information of each batch, such as: product, route, volume, and time-windows for each storage level indicated in Figure 4 are determined. The algorithm calculates sending priorities, looking for attendance of deadlines (time-windows). Thus, a preliminary order of batches is naturally given by the assignment module as batches


Page 18 of 63

Page 19 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are generated by the methodology. In a complementary step, the proposed MILP model starts the search procedure with no initial sequencing solution loaded and the model optimizes the sequencing of batches in each involved pipeline within the network, trying to address storage issues and minimization of plugs. The model is guided by temporal bounds of time-windows and is hereafter detailed in section 4.

4

The MILP sequencing

The sequencing for pumping batches made by means of MILP aims at an optimized inventory management and an efficient usage of pipelines. The modeling and optimization tool IBM ILOG CPLEX Optimization Studio 12.5 was used for computational purposes. Sparse sets are adopted to define the MILP model as presented, for instance, in Table 15 in Nomenclature section. This approach allows handling only valid indexes for variables/constraints. All parameters, variables, and sparse sets used in the MILP set of equations/inequalities are presented within the Nomenclature section, namely Table 13, Table 14, Table 15, and Table 16. For readability purposes, the notation of parameters and sets starts with an uppercase letter. On the other hand, all variables start with a lowercase letter. In particular the binary variables start with a “b” (e.g., bM Ob,b0 ,d , bRevb,b0 ,d , bBBN db,b0 ,d ). We present two MILP sequencing formulations: (i) Base model and, (ii) MI model, where MI stands for Minimization of Interfaces model. The MI model is, in fact, an extension of the Base model, which takes into account the additional goal of minimizing contamination due to undesired interfaces formed by two incompatible products that are transported one after another into the same pipeline.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.1

The MILP sequencing: Base model

The objective function and constraints of the MILP Base model are presented in the following subsections.

4.1.1

Base model objective function

The Base model objective function (expression 1) minimizes the costs of storage levels’ violations. As indicated in Figure 4, three main levels for violations are considered (CAP , M nM x, and T arget), which are concatenated by the index f x. The obtained values for violation variables adOb,n,f x , atOb,n,f x , adDb,n0 ,f x , atDb,n0 ,f x , adDdb,n0 ,f x , and atDdb,n0 ,f x are directly influenced by to the pumping order defined by the model. Within origin areas (BOf x), the storage levels in f x for pumping in advance (adOb,n,f x ) a batch b are considered, weighted by the cost factor Ktof x . In a similar way, the pumping delay (atOb,n,f x ) is also weighted by the cost factor Ktof x . Within the receiving areas, a differentiation is made according to the presence (BDf xDin) or absence (BDf x) of dynamic time-windows. So, if an area does not have the dynamic windows (BDf x), the receiving in advance (adDb,n0 ,f x ) or delay (atDb,n0 ,f x ) are penalized by the cost factor Ktdf x . In an area with dynamic windows (BDf xDin), the receiving in advance (adDdb,n0 ,f x ) or delay (atDdb,n0 ,f x ) are penalized by the cost factor KtdDf x , with an order of magnitude higher than Ktdf x . The order of magnitude for cost factors should be as following: (KtoCAP , KtdCAP , KtdDCAP ) >> (KtoM nM x , KtdM nM x , KtdDM nM x ) >> (KtoT arget , KtdT arget , KtdDT arget ).


Page 20 of 63

Page 21 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Minimize X

adOb,n,f x + atOb,n,f x ∗ Ktof x

{b,n,f x}∈BOf x

+

adDb,n0 ,f x + atDb,n0 ,f x ∗ Ktdf x

X

(1)

{b,n0 ,f x}∈BDf x

+

adDdb,n0 ,f x + atDdb,n0 ,f x ∗ KtdDf x

X {b,n0 ,f x}∈BDf xDin

4.1.2

Base model constraints

Groups of constraints with functional similarities are arranged into subsections.

4.1.3

Constraints for time-windows’ violations

Within batch receiving, constraint 2 restricts the start time to receive a batch b in destination area n0 (irb,n,n0 ,d ) to be greater than its available receiving time (T RDb,f x ). Violations can occur and variable adDb,n0 ,f x is penalized in the objective function. irb,n,n0 ,d ≥ T RDb,f x − adDb,n0 ,f x (2) 0

0

∀ {b, n, n , d} ∈ BN N D, {b, n , f x} ∈ BDf x In a similar way, constraint 3 restricts the final receiving time of a batch b in destination area n0 (f rb,n,n0 ,d ) to its critical receiving time (T RCb,f x ). Violations can occur (atDb,n0 ,f x ) and they are penalized in objective function. f rb,n,n0 ,d ≤ T RCb,f x + atDb,n0 ,f x (3) 0

0

∀ {b, n, n , d} ∈ BN N D, {b, n , f x} ∈ BDf x Equation 4 indicates that the advance in time may not occur in the pumping from the source node, when analyzing the physical storage level (CAP ). In practical terms, this



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 63

reflects an attempt to send the product without having it available on the pumping origin area. For violations of other storage levels, penalties are imposed.

adOb,n,f x = 0 ∀ {b, n, f x} ∈ BOf x | f x = CAP

4.1.4

(4)

Constraints for dynamic time-windows’ violations

In the studied real-world pipeline network the supply of a product in a destination node can be done by different source nodes. For instance, node N2 of Figure 1 can be supplied by nodes N3 and N4. Within the considered assignment algorithm, each pumping batch is generated and linked with a destination time-window to be supplied. However, at the MILP sequencing, the model can decide to change the batch that supplies a specific destination time-window. For doing this, the concept of dynamic time-windows is herein proposed. In essence, the model decides which pumping batch is the most appropriate to supply a time-window associated with a specific destination node product. The batch precedence is identified considering the receiving start time in destination area. Constraints 5 and 6 define the binary variable value (bBBN db,b0 ,n0 ), for the receiving start time of batches b and b0 . In constraint 5, bBBN db,b0 ,n0 = 1 if irb,n,n0 ,d is greater than irb0 ,nx,n0 ,d0 , then batch b0 precedes batch b. Otherwise, bBBN db,b0 ,n0 = 0 (constraint 6), batch b precedes batch b0 . irb0 ,nx,n0 ,d0 − irb,n,n0 ,d ≤ U ∗ (1 − bBBN db,b0 ,n0 ) (5) 0

0

0

0

0

0

0

0

∀ {b, n, n , d} ∈ BN N D, {b , nx, n , d } ∈ BN N D, {b, b , n } ∈ BBN din

irb,n,n0 ,d − irb0 ,nx,n0 ,d0 ≤ U ∗ bBBN db,b0 ,n0 (6) 0

0

0

0

∀ {b, n, n , d} ∈ BN N D, {b , nx, n , d } ∈ BN N D, {b, b , n } ∈ BBN din Constraint 7 calculates the batch ordering to receive a specific batch in a destination area n0 . The receiving ordering value of batch b in destination area n0 of product pr is attributed


Page 23 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to the integer variable iRecBatn0 ,pr,b . The (+1) term in the left part of this constraint avoids the occurrence of two batches with the same receiving order. iRecBatn0 ,pr,b0 − iRecBatn0 ,pr,b + 1 ≤ U ∗ (1 − bBBN db,b0 ,n0 ) (7) 0

0

0

0

0

∀ {n , pr, b } ∈ BDP din, {n , pr, b} ∈ BDP din, {b, b , n } ∈ BBN din Constraints 8 to 12 are used to identify the time-window that corresponds to a batch b receiving order (iRecBatn0 ,pr,b ). In case bBN P Ib,n0 ,pr,i = 1, batch b is attributed to timewindow with index i in destination area n0 . i − iRecBatn0 ,pr,b ≤ U ∗ (1 − bBN P I1b,n0 ,pr,i ) (8) 0

∀ {b, n , pr, i} ∈ BN P ind

i − iRecBatn0 ,pr,b ≥ (L − e) ∗ (bBN P I1b,n0 ,pr,i ) + e (9) 0

∀ {b, n , pr, i} ∈ BN P ind i − iRecBatn0 ,pr,b ≥ L ∗ (1 − bBN P I2b,n0 ,pr,i ) (10) 0

∀ {b, n , pr, i} ∈ BN P ind i − iRecBatn0 ,pr,b ≤ (U + e) ∗ (bBN P I2b,n0 ,pr,i ) − e (11) ∀ {b, n0 , pr, i} ∈ BN P ind bBN P Ib,n0 ,pr,i = bBN P I1b,n0 ,pr,i + bBN P I2b,n0 ,pr,i − 1 (12) 0

∀ {b, n , pr, i} ∈ BN P ind Equation 13 limits each batch b to only one time-window i. Similarly, equation 14 limits each time-window i to only one batch b. X

bBN P Ib,n0 ,pr,i = 1 ∀ {n0 , pr, b} ∈ BDP din

{b,n0 ,pr,i}∈BN P ind


(13)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

X

bBN P Ib,n0 ,pr,i = 1 ∀ {bx, n0 , pr, i} ∈ BN P ind

Page 24 of 63

(14)

{b,n0 ,pr,i}∈BN P ind

Constraints 15 and 16 calculate dynamic time-windows’ violations. Constraint (15) limits the receiving start time of batch b in destination area n0 (irb,n,n0 ,d ) to be greater than trd value of time-window i, for storage level f x, if i is attributed for batch b (bBN P Ib,n0 ,pr,i = 1). Violations can occur and are penalized in the objective function (adDdb,n0 ,f x ). Constraint 16 limits the receiving final time of batch b in destination area n0 (f rb,n,n0 ,d ) to be less than trc value of time-window i, for storage level f x, if i is attributed for batch b (bBN P Ib,n0 ,pr,i = 1). Violations can occur and are penalized in objective function (atDdb,n0 ,f x ). irb,n,n0 ,d − (trd − adDdb,n0 ,f x ) ≥ L ∗ (1 − bBN P Ib,n0 ,pr,i ) (15) 0

0

0

∀ {b, n, n , d} ∈ BN N D, {b, n , pr, i} ∈ BN P ind, {n , pr, i, f x, trd, trc} ∈ N P IJan

f rb,n,n0 ,d − (trc + atDdb,n0 ,f x ) ≤ U ∗ (1 − bBN P Ib,n0 ,pr,i ) (16) 0

0

0

∀ {b, n, n , d} ∈ BN N D, {b, n , pr, i} ∈ BN P ind, {n , pr, i, f x, trd, trc} ∈ N P IJan 4.1.5

Constraints for batches’ flow into the network

Equation 17 defines the pumping final time of batch b in pipeline d from n to n0 (f bb,n,n0 ,d ) as the start pumping time (ibb,n,n0 ,d ) of b added to the pumping time of this batch in the pipeline (T BBatb,d ). This set of constraints is dynamically generated for pipelines were batches will be transported. For the initialization batches, the final of pumping is considered as being equal to the start of pumping (constraint 18), since batches were already pumped. In an analogous way of equations 17 and 18, constraints 19 and 20 calculate the final receiving time (f rb,n,n0 ,d ) of batches. f bb,n,n0 ,d = ibb,n,n0 ,d + T BBatb,d ∀ {b, n, n0 , d} ∈ BN N Dbomb

(17)

f bb,n,n0 ,d = ibb,n,n0 ,d ∀ {b, n, n0 , d} ∈ BN N DEstDuto

(18)


Page 25 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


f rb,n,n0 ,d = irb,n,n0 ,d + T BBatb,d ∀ {b, n, n0 , d} ∈ BN N Dbomb

(19)

f rb,n,n0 ,d = irb,n,n0 ,d ∀ {b, n, n0 , d} ∈ BN N DEstDuto

(20)

Equation 21 indicates that the receiving start time (irb,n,n0 ,d ) of a batch b in pipeline d is equal to the start pumping time (ibb,n,n0 ,d ) added to the time to dislocate the batch b through this pipeline. This constraint is a simplification for the detailed pumping process. In this case, it is considered that a batch is moved at an average time. In fact, for a batch to be moved through pipelines other batch (or batches) has (have) to push it. The detailed pumping is addressed in the Timing Block 1 of Figure 2. For sequencing purposes this feature is addressed in a simpler manner.

irb,n,n0 ,d = ibb,n,n0 ,d + T DBatb,d ∀ {b, n, n0 , d} ∈ BN N D

(21)

Equation 22 is responsible for the batch b alignment between pipelines d and d0 connected to node n0 . So, the start pumping time should be equal to the start receiving time. ibb,n0 ,nx,d0 = irb,n,n0 ,d ∀ {b, n0 , nx, d0 } ∈ BN N Dbomb, {b, n, n0 , d} ∈ BN N Dbomb

(22)

| n 6= n0 ∧ n0 6= nx ∧ d 6= d0 4.1.6

Constraints for order change

Constraint 23 limits the option of changing some batch orders. This equation defines the value 1 for binary variable bM Ob,b0 ,d within the domain of set BBDrestringe. This is done to keep an initial order at specific cases, thus, batch b necessarily precedes batch b0 in pipeline d. The set BBDrestringe is constructed substantially by initialization batches, which are inside the pipeline and cannot have the order changed.

bM Ob,b0 ,d = 1 ∀ {b, b0 , d} ∈ BBDrestringe 25 ACS Paragon Plus Environment

(23)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 63

Constraints 24 and 25 define the batch precedence of b and b0 in pipeline d. If the initial pumping time of batch b0 (ibb0 ,nx,nx0 ,d ) is greater than the final pumping time of batch b (f bb,n,n0 ,d ) in pipeline d so, the binary variable bM Ob,b0 ,d is set to 1. Otherwise, bM Ob,b0 ,d is 0 and batch b0 precedes batch b in pipeline d. f bb,n,n0 ,d − ibb0 ,nx,nx0 ,d ≤ U ∗ (1 − bM Ob,b0 ,d ) (24) 0

0

0

0

∀ {b, n, n , d} ∈ BN N D, {b , nx, nx , d} ∈ BN N D, {b, b , d} ∈ BBDtotal

f bb0 ,nx,nx0 ,d − ibb,n,n0 ,d ≤ U ∗ (bM Ob,b0 ,d ) (25) 0

0

0

0

∀ {b, n, n , d} ∈ BN N D, {b , nx, nx , d} ∈ BN N D, {b, b , d} ∈ BBDtotal Constraints 26 and 27 are cuts to the solution search space. Let’s consider that two batches b and b0 are transported into the same pipeline d. After a junction node n0 , however, the batches are transported by distinct pipelines, respectively, d0 and dx0 . Let’s also consider that bM Ob,b0 ,d = 1. Thus, we know that batch b precedes b0 in pipeline d and the start pumping time of b in d0 will be necessarily smaller than the start pumping of batch b0 in dx0 (constraint 26). If bM Ob,b0 ,d = 0, the start pumping time of b0 in dx0 will be smaller than the start pumping of batch b in d0 (constraint 27). ibb,n0 ,nx,d0 − ibb0 ,n0 ,nx0 ,dx0 ≤ U ∗ (1 − bM Ob,b0 ,d ) ∀ {b, n0 , nx, d0 } ∈ BN N D, {b0 , n0 , nx0 , dx0 } ∈ BN N D, {b, b0 , d} ∈ BBDtotal

(26)

| b 6= b0 ∧ nx 6= nx0 ∧ d0 6= dx0 ∧ d 6= d0

ibb0 ,n0 ,nx0 ,dx0 − ibb,n0 ,nx,d0 ≤ U ∗ bM Ob,b0 ,d ∀ {b, n0 , nx, d0 } ∈ BN N D, {b0 , n0 , nx0 , dx0 } ∈ BN N D, {b, b0 , d} ∈ BBDtotal | b 6= b0 ∧ nx 6= nx0 ∧ d0 6= dx0 ∧ d 6= d0


(27)

Page 27 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.1.7

Surge tank constraints

For sequencing purposes, two cases are considered for surge tank operations. In the Timing Block 1 of Figure 2 further details of surge tank operations are addressed. (i) Flow rate input smaller than the flow rate output of surge tank node: in this case (constraint 28), the synchronism between received and sent batches in surge tank node is made by the final times. The pumping start of b0 is delayed until the tank accumulates enough volume to make viable to equal the final receiving time of batch b (f rb,n,n0 ,d ) with the final pumping of batch b0 (f bb0 ,n0 ,nx0 ,d0 ). (ii) Flow rate input greater than flow rate output of surge tank node: in this case (constraint 29) the synchronism between received and sent batches in surge tank node is made by the start times. The pumping start of batch b0 (ibb0 ,n0 ,nx0 ,d0 ) must be equal to the receiving start of batch b (irb,n,n0 ,d ). f bb0 ,n0 ,nx0 ,d0 = f rb,n,n0 ,d ∀ {b, b0 , n, n0 , nx0 , d, d0 } ∈ P U LM AO (28) | V Bombb,d ≤ V Bombb0 ,d0 ibb0 ,n0 ,nx0 ,d0 = irb,n,n0 ,d ∀ {b, b0 , n, n0 , nx0 , d, d0 } ∈ P U LM AO (29) | V Bombb,d > V Bombb0 ,d0 4.1.8

Reversion flow constraints

The BBDrev set has all viable combinations of batches b and b0 (b 6= b0 ) to be considered for the reversion flow in the pipeline d. Additionally, binary variable bRevb,b0 ,d identifies the precedence between batches b and b0 in pipeline d. For sequencing purposes, a simplified approach to evaluate the reversion flow procedure is adopted. A more detailed formulation is presented in the Timing block. 1 Constraint 30 adds to the pumping start of batch b (ibb,n,n0 ,d ) an estimated time of reversion flow operation (T DBatb,d ) in relation to pumping end of batch b0 (f bb0 ,nx,nx0 ,d ). 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 63

The goal is to represent the spent time for reversion flow operation if variable bRevb,b0 ,d assumes a value equal to 1. Constraint 31 limits only one variable bRevb,b0 ,d to be equal to 1. The additional time to reversion flow is added to b0 pumping time if bRevb,b0 ,d = 1, otherwise, is added to b pumping time. ibb,n,n0 ,d ≥ f bb0 ,nx,nx0 ,d + T DBatb,d − (1 − bRevb,b0 ,d ) ∗ U ∀ {b, n, n0 , d} ∈ BN N Dbomb, {b0 , nx, nx0 , d} ∈ BN N Dbomb, {b, b0 , d} ∈ BBDrev

(30)

| n 6= n0 ∧ n0 6= nx0 bRevb0 ,b,d + bRevb,b0 ,d = 1 ∀ {b, b0 , d} ∈ BBDrev

4.2 4.2.1

(31)

The MILP sequencing: MI model MI model objective function

The expression (32) is the MI model objective function. Comparing with expression (1) a term is aggregated considering the nP lugd variable. The total number of plugs in each pipeline d is penalized by a cost factor Kplug. So, the addition of plug products is avoided, unless the need of maintaining the inventory levels between targets induces the adjacent transportation of incompatible products, which demands a plug inclusion. In this way, the obtained solution balances storage violations and incompatibility occurrences. The main goal is to obtain inventory profiles without violations, mainly in the capacity levels, but plug inclusions are avoided, whenever possible. Afterwards, the post-assignment/sequencing procedure chooses an adequate product, according to adjacent batches and the pipeline characteristics, to each plug necessity detected by the MILP sequencing model.


Page 29 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Minimize X

adOb,n,f x + atOb,n,f x ∗ Ktof x

{b,n,f x}∈BOf x

+

adDb,n0 ,f x + atDb,n0 ,f x ∗ Ktdf x

X {b,n0 ,f x}∈BDf x

+

(32)

adDdb,n0 ,f x + atDdb,n0 ,f x ∗ KtdDf x

X {b,n0 ,f x}∈BDf xDin

+

nP lugd ∗ Kplug

X d∈DutosIncomp

4.2.2

MI model additional constraints

Constraints from (2) to (31) previously defined to the Base model are also used for the MI model. In addition, constraints from (33) to (39) are created to determine plug necessities. The total number of plugs is a factor minimized in the MI model objective function defined in (32). In order to address plug insertions it is necessary to determine if two specific batches are consecutively pumped in a pipeline, that is, they are adjacent, pumped one after another. It is not enough to determine that a batch b precedes a batch b0 in a pipeline d (or b0 precedes b), a function already mapped by the binary variable bM Ob,b0 ,d . In this section, the exact ordering of pumped batches in each pipeline d (first batch, second batch, third...) is identified and, therefore, the plug necessity between consecutive batches is set. The ordering of batch b in pipeline d is assigned to variable ordBDd,b . This is done by considering the values obtained by binary variable bM Ob,b0 ,d to pipeline d. Constraint (33) defines the ordering of each batch b by the number of batches that are transported in pipeline d, minus the number of batches that succeeds batch b. It is important to highlight two conditions for a batch b0 to succeed a batch b in pipeline d: (i) if bM Ob,b0 ,d = 1 then b0 succeeds batch b; (ii) if bM Ob0 ,b,d = 0 then b0 also succeeds batch b. These two conditions



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 63

are mapped in equation (33). X

ordBDd,b = N Batd −

bM Ob,b0 ,d

b,b0 ,d∈BBDtotal

−

X

(1 − bM Ob0 ,b,d )

(33)

b0 ,b,d∈BBDtotal

∀ {d, b} ∈ DBO Constraint (34) imposes that two batches b and b0 cannot have the same ordering for the same pipeline d. ordBDd,b + 1 ≤ ordBDd,b0 (34) ∀ {b, b0 , d} ∈ BBDrestringe, d ∈ Dutos Since we know the batch ordering in each pipeline, constraints (35) and (36) are used to identify consecutive batches in pipeline d. These constraints define the value for the binary variable bEmSeqb,b0 ,d . If the binary variable is equal to 1, batches b and b0 are one after another in pipeline d, 0 otherwise.

ordBDd,b0 − ordBDd,b − 1 ≤ (N Batd + 1) ∗ (1 − bEmSeqb,b0 ,d ) (35) 0

∀ {b, b , d} ∈ BBDseq

ordBDd,b0 − ordBDd,b − 1 ≥ −(N Batd + 1) ∗ (1 − bEmSeqb,b0 ,d ) (36) ∀ {b, b0 , d} ∈ BBDseq Constraint (37) limits the search space, since only one binary variable considering bEmSeqb,b0 ,d and bEmSeqb0 ,b,d can assume a value equal to 1. bEmSeqb,b0 ,d + bEmSeqb0 ,b,d ≤ 1 ∀ {b, b0 , d} ∈ BBDseq


(37)

Page 31 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Constraint (38) limits the number of binary variables bEmSeqb,b0 ,d that can assume a value equal to 1 on pipeline d. X

bEmSeqb,b0 ,d = N Batd − 1 ∀ d ∈ Dutos

(38)

{b,b0 ,d}∈BBDseq

Equation (39) calculates the number of incompatibilities in pipeline d by the sum of consecutive batches (bEmSeqb,b0 ,d = 1) that are incompatible ({b, b0 , d} ∈ BBDincomp). In cases of consecutive incompatible batches, a plug product has to be inserted to feasible the transportation.

nP lugd =

X

bEmSeqb,b0 ,d ∀ d ∈ DutosIncomp

(39)

{b,b0 ,d}∈BBDincomp

4.3

Post-assignment/sequencing

After the batch generation and the ordering changes suggested by the MILP model, a module called post-assignment/sequencing, which is illustrated in Figure 10, gives the output of Assignment/Sequencing Block of Figure 2. The post-assignment/sequencing process has the goal of handling some operational characteristics, as described: (i) Reversal flow batches procedure: Evaluates each pipeline and verifies whether pipeline flow reversion procedures are necessary or not. This analysis verifies if two batches in sequence, into the same pipeline, have distinct flow directions. In this case, it is generated a new batch with reversal flow route to ensure that both batches will be delivered in their respective destinations. (ii) Plug batches procedure: Evaluates each pipeline verifying prohibited interfaces between two consecutive batches, namely bx and by . If an undesired interface is identified, a new batch is generated (bplug ), with a minimum operational plug cost, and composed of a product that makes the sequence (bx − bplug − by ) with compatible interfaces. It is important to highlight that this post-processing procedure is not able to reduce the 31 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

number of plug insertions, which was already defined by the MILP model (variable nP lugd ). The module takes into account the involved batches (bx , by ) and chooses an adequate product/volume to form the plug batch (bplug ) in order to minimize, in each evaluated case, the interface cost. The “local” cost of each necessary plug is reduced by the plug batches procedure. Reversal flow Plug batches batches Reversal flow batches insertion Plug batches insertion

Figure 10: Post-assignment/sequencing execution flow

4.4

Brief comments on the proposed MILP sequencing

Two MILP formulations were proposed to address the sequencing for pumping batches. The model seeks for an optimized inventory management and an efficient usage of pipelines. The Base model comprises the objective function (1) and constraints from (2) to (31). The Minimization of Interfaces model, or MI model, involves the objective function defined in (32) and constraints from (2) to (31) and from (33) to (39). In fact, the MI model is an extension of the Base model, which takes into account the additional goal of minimizing contamination due to undesired interfaces formed by two incompatible products that are transported, one after another, into the same pipeline. Sparse sets were used for the implementation of models, a procedure that reduced the indexes’ generation to just viable choices. As before mentioned, simplified conditions are adopted within the MILP formulations for: surge tank operations, reversion flow procedures, and the batches’ flow propagation into the network. However, it is important to note that this simplified approach is applied at a relative early sequencing phase highlighted in Figure 2. Thus, an already optimized batch ordering is given to the MILP timing model, which indeed contributes to the problem solving in a reasonable computational time, since few ordering decisions remained to the timing phase. 1 Furthermore, the dynamic time-windows consideration added an important 32 ACS Paragon Plus Environment

Page 32 of 63

Page 33 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


flexibility for optimizing orders previously suggested by the heuristic algorithm during the assignment procedure. The following results section indicates the computational gain of the collaborative approach in relation to a former monolithic MILP model 2 for the assignment/sequencing block. It is highlighted that the proposed assignment/sequencing collaborative approach can be used to define operational batches with their volumes and routes in pipeline networks, besides providing an optimized ordering of batches from source nodes to destination nodes, in computational times from seconds to few minutes. The proposed approach is applied in pipeline networks of different topologies. Within case study 1, a small example is studied; case study 2 focus on the scheduling solution of a real-world pipeline network.

5

Results and discussion: Case study 1

This section presents results for an instance considering the pipeline network of Figure 1. The network has six pipelines linking five areas and it transports a set of nine products. The tests were conducted in a computer with Intel core i7, 2.1GHz of clock (octa-core), 8GB of RAM, operational system Windows 8 – 64 bits, and Java version 1.6.

5.1

Case study 1: Proposed decomposition approach

The used data for the herein proposed Assignment module and MILP sequencing model is available in the Supporting Information. A scenario involving an one month period was used in order to perform the computational test. This set of data can be used for reproducibility purposes of the proposed paper core task: the assignment and sequencing of batches in pipeline networks. This case study 1 is relatively small, but is intended to maintain the main features presented in the operational scheduling of pipeline networks. The set of data includes, for instance, a strategical planning output that specifies the total volume of product that should be sent from an origin node to a destination node by a specific



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

route during the in-analysis 30-day horizon (Table S1). The given total volume has to be broken in operational volumes of batches, a task defined by the proposed Assignment module (lot-sizing) based on time-windows propagation. Production and demand data, including stepwise functions, are given in the supporting information data. Tankage limits for each area and product are also indicated (Table S2). Tables S1 and S2 are input parameters to the Assignment module. The output given by the Assignment module is detailed in Table S3: 75 batches were generated in order to attend the planning given in Table S1. For each batch the algorithm indicates the batch number (in fact, a label), the pumping route, the involved product, the operational volume (lot-sizing) of batches, the estimated average pumping flow rate, and the calculated time-windows for sending/receiving purposes. Relatively small deviations from the total planned volume (given in Table S1) can occur due to recommendations of operational (minimum/maximum) volumes of batches. For instance, if we take the first line of Table S1 it is indicated that a total volume of 15943 vu of product 1 should be sent by route N4→5→N2→4→N3. If we take Table S3, we can notice that the Assignment module suggested two batches of 8000 vu each (labeled as batches 33 and 34), that sum up 16000 vu of product 1 to be sent by the specified route. The number of operational batches generated by the Assignment module varied according to product and route. For instance, for product 2, route N3→4→N2→1→N1, the total volume of 112175 vu indicated in Table S1, line 2, was broken in 7 batches with volumes, in the great majority, of 18000 vu. The Assignment module took a fraction of a second to generate the set of 75 batches with the associated time-windows that are input parameters to the next step, the sequencing model. The proposed MILP sequencing model also takes into account a series of parameters given in the Supporting Information, for instance, Tables S3 to S5, the list of incompatible products given in Table S6, in-transit batches and pipeline volumes given in Tables S7 and S8, and technical parameters used in the MILP model (Tables S9 to S11). Table 1 summarizes the computational results from both: the proposed Base model (sec-


Page 34 of 63

Page 35 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tion 4.1) and the Minimization of Interfaces (MI) model (section 4.2). This table indicates: total number of variables; number of integer variables; number of constraints; number of batches for the MILP model; objective function values (respectively, expressions 1 and 32); integrality gap (%); computational time, in seconds; and, finally, the number of necessary plugs within the proposed sequencing. It is possible to observe in Table 1 that the sequencing of 75 batches was obtained in a relative small computational time (≈ 1 second for the Base model and ≈ 34 seconds for the MI model). Table 1: Results for 30-day scenario - Case study 1 Variables Int.Var. Constraints Batches Obj.Value Gap Plugs Time

(#) (#) (#) (#) (104 ) (%) (#) (s)

Base Model 4477 3133 13021 75 4.64 0 6 1

MI Model 5384 3987 15933 75 5.06 0 4 34

The proposed Assignment/Sequencing is embedded into the hierarchical decomposition approach of Figure 2. Figures 11 to 14 are used to illustrate the main outputs obtained for case study 1. For these outputs, the Base model was used within the Assignment/Sequencing Block. A computational time of less than ten seconds was spent in the entire running procedure. Figure 11 presents the obtained Gantt chart of the final scheduling (the output of the Timing Block 1 of Figure 2), which is influenced by the proposed assignment/sequencing methodology. The pipelines used during the scheduling horizon are vertically presented in Figure 11. The scenario starts on January/06th, and 30 days are necessary to complete the scheduling. Different colors are used to represent each involved product. In Figure 11, pipeline D4, product P1 is incompatible with product P2 (the incompatibility table is given in the Supporting Information). A plug of P7 is inserted between them. Inclusions of P7 as a plug occurred four times during the scheduling horizon at this pipeline. In addition, in pipeline D5, product P6 is incompatible with product P4. A plug of P2 is in35 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 11: Gantt chart influenced by the Base model - Case study 1 serted between them. Following a similar reasoning, a plug of P8 is inserted between batches of P9 and P1 in D5. According to Table 1, 6 plugs were necessary in all pipelines during the entire scheduling horizon. On the other hand, the MI model indicated a solution with a total of 4 plugs (in fact, just the plugs indicated in D4). Thus, the MI model eliminates two plugs compared with the Base model sequencing. It is worth emphasizing that the Base model did not have the explicit goal of minimizing plugs; so that, the obtained value for the number of necessary plugs is just a consequence from a sequencing answer that seeks to minimize violation issues, expression (1). On the other hand, the MI model is an extension of the Base model, which takes into account the additional goal of minimizing the number of necessary plugs, expression (32). In Figure 1, pipelines D4 and D6 are reverted. With the aid of Table S3, we can notice that product P2 is sent from N3 to N1 by a route that involves pipeline D4 (N3→4→N2→1→N1) in, for instance, batches 36, 40, 44, 49, 53, 58, and 61. On the other hand, batches 33 and 34 of product P1 (8000 vu each) have to be sent from N4 to N3 by a route that also involves D4 (N4→5→N2→4→N3). Thus, pipeline D4 has to be reverted. Within the Gantt chart of


Page 36 of 63

Page 37 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 11 we can identify in D4 batches of P1 (in brown) and P2 (in red); thus, this pipeline was reverted to attend the operational purposes. We can also notice that product P7 (in green) was used to generate auxiliary reversion batches. A similar analysis can be made to pipeline D6, which is also reverted during the scheduling horizon. Figures 12 to 14 present some storage profiles along the 30-day scheduling horizon. In fact, we choose three different areas, including a depot (N1), a refinery (N4), and a harbor (N5), and two products of each chosen area, to exemplify the profiles. We highlight, however, that capacity level violations do not occurred for any product/area during the entire horizon. If we take, for instance, Figure 12(a), we can notice 7 positive slopes and 6 negative ones. During the refereed positive slopes, the refinery N4 is producing P4 at a constant rate and the storage level is increasing. During the 6 negative ones, a sending procedure of P4 from N4 can be observed. This information can be also checked in the Gantt chart of Figure 11, pipeline D5, where 6 batches of P4 can be identified. The storage level decreases as the output flow rate for the sending procedure is bigger than the production rate. Similar analysis can be made on Figure 13, but care should be taken to notice that N1 is a depot: positive slopes are generated by the receiving of batches; negative ones are created by sending procedures from the depot and also by attending the “local” demand. In exception of node N5, the harbor, a constant rate was considered for demand/production during the scheduling horizon, as indicated in the Supporting Information, Table S12. For node N5, the analysis of Figure 14 is dependent on stepwise demand or supply from vessels, as indicated in Table S12. The planned arrival/departure of vessels is a parameters given by the strategical planning. For instance, in Table S12 we can notice that in N5, product P3, there is the “production” of 45000 vu from hours 0 to 48 (relatively to a 720-hour horizon). This supply of product P3 is from a vessel. Thus, if we take Figure 14(a), we can visualize this product arrival at N5 in the horizon beginning, namely, January 6th. The proposed approach is able to deal with constant demand/production rates and also stepwise ones. Still in Figure 14(a) we can notice that around January 21nd, the product P3 was used to provide auxiliary batches for



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 63

reversion flow procedures in pipeline D6.

(a)

(b)

Figure 12: CS1 - Storage on N4: (a) Product 4; (b) Product 8

(a)

(b)

Figure 13: CS1 - Storage on N1: (a) Product 3; (b) Product 5

(a)

(b)

Figure 14: CS1 - Storage on N5: (a) Product 3; (b) Product 4 This first case study illustrates to the reader by means of a Gantt chart (Figure 11), inventory profiles (Figures 12 to 14) and the detailed list of batches to be ordered (Table S2 in Supporting Information) that pipeline-scheduling problems can be difficult to be solved,


Page 39 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


even in a reduced pipeline network scenario, as illustrated in Figure 1. In fact, the literature indicates that short-term pipeline-scheduling problems are NP-hard combinatorial optimization problems. 19 The proposed scheduling approach based on a collaborative methodology involving heuristic algorithms and MILP models could, however, solve the proposed problem in very few seconds. In section 6, the methodology is applied in more a complex case study involving a real-world pipeline network, which is tested under a set of different input scenarios.

5.2

Case study 1: Monolithic assignment-sequencing approach 2

This subsection presents the computational results obtained in case study 1 by the monolithic MILP model proposed by ref 2, which is a unified approach intended to find the optimal answer for the assignment-sequencing task in the studied problem. Thus, by the observation of sections 5.1 and 5.2, one can compare the solution quality and computational load from both approaches: the proposed decomposition-based and the monolithic model. 2 In a similar reasoning applied to the previous subsection, the used data for the MILP assignment-sequencing monolithic model is available in the Supporting Information. Table 2 summarizes the computational results from the monolithic MILP assignmentsequencing approach for a 30-day scenario of the pipeline network represented in Figure 1. This table indicates: total number of variables; number of integer variables; number of constraints; number of batches suggested by the MILP model; integrality gap (%); and computational time, in seconds. It is possible to observe in Table 2 that the computational load was a concern: the solution have not evolved from a high integrality gap (≈ 99.9%) within the 10-hour allowed computational time. In fact, a solver execution failure due to out-of-memory occurred in less than five hours. Within the proposed approach (Assignment module + MILP Base sequencing model) it was possible to obtain solutions in relatively small computational times (few seconds) for the same 30-day horizon data, as indicated in Table 1. Therefore, comparisons of the decomposition approach solution with the “optimal” 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

solution obtained by the monolithic approach were, unfortunately, not possible due to the computational load evidenced in Table 2. As an alternative, Table 3 compares the generated batches after running the monolithic approach for one hour. The number of batches and the associated volume (lot-sizing) of generated batches from the monolithic and the proposed decomposition approaches may differ. Table 3 exemplifies these features for the first two involved products. For product 1, route N4→5→N2→4→N3, for instance, 3 batches of 5000 volumetric units (vu) each were suggested by the monolithic model and 2 batches of 8000 vu each were suggested by the proposed decomposition. Thus, the sum of volumes suggested were, respectively, 15000 vu and 16000 vu. As previously stated in section 5.1, relatively small deviations from the total planned volume (given in Table S1: 15943 vu for product 1 and 112175 for product 2) can occur. An analogous analysis can be conducted for product 2, route N3→4→N2→1→N1, indicating the sum of volumes suggested as, respectively, 107175 vu (monolithic) and 116000 vu (decomposition). The solution suggested by the decomposition approach uses time-windows propagation for lot-sizing purposes, as explained in section 3.2. It is worth emphasizing that the computational results of section 5.1 indicate that capacity level violations do not occurred for the decomposition approach suggested solution in any product/area during the entire horizon. The suggested batches/volumes were, therefore, considered adequate for operation purposes. Table 2: Assignment-Sequencing monolithic MILP model: 2 Results for 30-day scenario Case study 1

Variables Int.Var. Constraints Batches Gap Time

(#) (#) (#) (#) (%) (s)

Results from MILP monolithic model 2 19383 2184 24962 73 99.994 16440*

*Solver failure (out-of-memory status).

As the computational load presented in Table 2 was relatively high, we also investigated the monolithic MILP assignment-sequencing approach 2 for smaller time horizons, namely 40 ACS Paragon Plus Environment

Page 40 of 63

Page 41 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Assignment-Sequencing monolithic MILP model 2 versus Proposed Decomposition Approach: Example of generated batches - 30-day scenario - Case study 1 Product

Route

1

N4→5→N2→4→N3

2

† After

N3→4→N2→1→N1

one hour of running.

Results from the MILP† monolithic model 2 Batch label Volume (vu) 2 5000 5 5000 10 5000 32 5000 38 5000 48 22175 50 25000 57 25000 58 25000

? Execution

Results from the Proposed Decomposition (with Base Model)? Batch label Volume (vu) 33 8000 34 8000 36 40 44 49 53 58 61

18000 18000 18000 18000 18000 18000 8000

in a second.

3-day, 4-day, and 5-day, instead of a 30-day horizon. The same set of data available in the Supporting Information for an one month period was used: we just considered the production/demand patterns for the first involved days. Table 4 indicates the obtained results. As it can be observed, the model size (number of variables and constraints) reduced significantly from Table 2. The prove of solution optimality, however, was still a difficult task. Just the 3-day and 4-day horizons converged to optimality within an one-hour running. The integrality gap was ≈32% for the 5-day horizon. The number of batches generated were, respectively, 10, 12, and 17. The proposed Decomposition Approach (with the Base Model) generates answers with, respectively, 10, 11, and 18 batches in a fraction of a second. As also indicated in Table 3, the number of generated batches from the monolithic and the proposed decomposition may differ, as the approaches are essentially different: the MILP monolithic model works directly on the of goal of maintaining under control the storage levels approximated by the model; the proposed decomposition uses time-windows propagation for lot-sizing and sequencing purposes. Thus, the evaluation function of each approach is different. However, both approaches conducted to solutions with no violations in capacity levels for the 3-day to 5-day tested scenarios. As a complementary information, we emphasize that capacity level violations do not occurred for any product/area during the entire 30-day 41 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 63

horizon solved by the proposed decomposition approach, as reported in section 5.1. Table 4: Results from the Monolithic MILP Model 2 and the Proposed Decomposition Approach: 3-day to 5-day scenarios Horizon (days) 3 4 5

Variables (#) 791 1002 1640

Int.Var. (#) 77 98 168

Monolithic MILP Model Constraints Time Obj.Value (#) (s) (#) 1084 41.3 217418.1 1395 3504 406358.3 2180 3600 878971.6

Gap (%) 0.0 0.0 32.0

Number of Batches 10 12 17

Decomposition Time Number of (s) Batches 0.1 10 0.1 11 0.2 18

By the evidenced computational load (Table 2 and Table 4), the use of a decomposition approach (Assignment module + MILP Base sequencing model) was of great value to obtain viable solutions in a computational time of few seconds, as indicated in, for instance, Table 1. The computational gain obtained with the proposed approach is also evidenced in this reduced size pipeline-scheduling scenario (Table 4). Further investigation about computational scalability of the proposed decomposition is still conducted for case study 1 in subsection 5.3. In addition, for the tested real-world pipeline-scheduling scenarios (section 6) the Assignment/Sequencing decomposition strategy played a fundamental role: without it, the attainment of practical solutions was not viable in computational times of few minutes (e.g., Table 12).

5.3

Case study 1: Proposed Decomposition Approach Scalability

Subsection 5.1 evidenced that the proposed Assignment/Sequencing decomposition was able to obtain viable operational answers in few seconds (e.g., Table 1, Figures 11 to 14). On the other hand, subsection 5.2 indicated that for the same 30-day scenario a monolithic MILP model 2 was not able to prove the solution optimality when considering the assignment and the sequencing tasks in a unified framework (e.g., Table 2 and Table 4). Thus, the decomposition strategy is evidenced as an important approach to avoid the computational load. However, the question of how far we can circumvent the computational burden with the proposed decomposition still remains. This section aims at addressing this question. 42 ACS Paragon Plus Environment

Page 43 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For doing this, we tested horizons of, respectively, 40 days, 50 days, and 60 days. The new scenarios were based on the same set of data available in the Supporting Information for an one month period. We just repeated the production/demand patterns for the following days as in the original first 30 days, producing scenarios from 40 to 60 days. Thus, the dimensions of generated models (Base model and MI model) increased significantly. Table 5 summarizes the obtained results. The used labels are similar to the ones explained for Table 1. We can observe in Table 5 that the Base model converged for the three tested scenarios (integrality gap = 0) in less than 17 minutes. On the other hand, the MI model was not able to converge to optimality in one hour of run: integrality gaps of, respectively, 14.08%, 16.87%, and 84.81%. The results for a 30-day horizon (Table 1) indicate that the Base model and the MI model converged to optimality in few seconds; according to Table 5, however, a fast convergence was not evidenced as the horizon increased. Thus, the computational load is indeed a significant issue, even in this decomposition scheme. Section 6 further details the computational performance of the proposed decomposition approach in a real-world pipeline network scheduling scenario. Table 5: Results for 40-day, 50-day, and 60-day scenarios - Case study 1

Variables Int.Var. Constraints Obj.Value Gap Plugs Time

6

(#) (#) (#) (107 ) (%) (#) (s)

40 Days Base Model MI Model 8336 10989 6279 8841 24829 33557 1.07 1.21 0 14.08 8 8 10 3600



Results and discussion: Case study 2

The second case study considered a multiproduct pipeline network that was firstly studied by ref 20 and is schematically depicted in Figure 15. This network presents, in fact, a real case study for a Brazilian pipeline network, which involves 14 areas (nodes). It is composed 43 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of 4 refineries (nodes N3, N4, N5, and N6), 2 harbors (N7 and N10), 2 final clients (N2 and N14), and 6 depots (N1, N8, N9, N11, N12, and N13), which receive or send a set of oil derivatives and ethanol. In particular, area N1 does not have tanks to store products. The areas are interconnected by 30 bidirectional pipelines, each one with a particular volume. Nowadays, more than 35 oil derivatives and ethanol can be transported in this network.

Figure 15: Case study 2: Pipeline Network 20 This section is organized in a way to evidence five main features: (i) Assignment module results; (ii) Comparing the final scheduling derived from heuristic ordering and from the MILP base model; (iii) MILP sequencing model general results; (iv) Comparison between results of MILP base model and MI model; and, finally, (v) Comparison of the assignment/sequencing block results with results obtained by ref 2. The tests were conducted in a computer with Intel core i7, 2.4GHz of clock (octa-core), 4GB of RAM, operational system Windows 7 – 64 bits, and Java version 1.6. Eight scenarios (S1 to S8) involving the data for one month period were used in order to perform the computational tests. These scenarios were also used in ref 2.


Page 44 of 63

Page 45 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6.1

Results for assignment module

Table 6 presents the number of batches generated by the assignment algorithm for each scenario. The batches generated by the post-assignment/sequencing process (executed after the MILP sequencing) are also presented: plug batches and batches for reversion flow procedure. The total time spent for this generation methodology, which involved three main steps (pre-assignment, batches’ generation, and post-assignment/sequencing), is presented in the bottom table line. We can observe that, considering a month horizon and a number of generated batches ranging from 270 (in S5) to 320 (in S2), the execution time remained of few seconds, which was considered adequate for operational purposes. Table 6: Some results of assignment module Scenario # of generated batches # of plug batches # of batches for reversion flow Time (seconds)

S1 288 4 4 2.81

S2 320 11 11 5.22

S3 281 11 7 2.94

S4 301 5 14 3.10

S5 270 7 14 2.23

S6 288 8 16 4.21

S7 290 7 9 4.06

S8 297 7 6 2.99

Table 7 exemplifies details of some batches obtained for scenario S2. Information about batch number, product, route, volume, and three levels of time-windows (in hours), associated for each batch, are indicated. The route allows to identify nodes and pipelines to be considered within the batch transportation process. For instance, batch 95 is pumped from node N6 to node N11 by pipelines 24, 28, and 4, respectively. Table 8 presents auxiliary batches to aid the reversion flow procedure. In these cases, the origin and destination nodes for the suggested routes are the same. In fact, the batch is pumped in only one pipeline to allow delivering other batch (batches) that is (are) inside the considered pipeline. Afterwards, the flow direction can be reverted and the auxiliary batch returns to the route origin node. More details about reversion flow procedures can be found, for instance, in ref 25. Table 9 presents plug batches inserted for compatibility purposes. 25 The operational volume of these batches is 5000 vu, for the considered pipeline network scenarios.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 63

Table 7: Examples of generated batches for scenario S2 Batch Number ...

Product Number ...

Obtained Route ...

Volume (vu) ...

95

31

N6-24-28-4-N11

15000

96

26

N3-12-15-N10

8000

97

14

N3-10-8-N5

2700

...

...

...

...

201

25

N6-24-28-4-N11

8000

202

34

N5-20-19-N13

12000

203

18

N6-25-2-5-18-N14

23000

...

...

...

...

TimeWindow ... Target MnMx Cap Target MnMx Cap Target MnMx Cap ... Target MnMx Cap Target MnMx Cap Target MnMx Cap ...

TED (h) ... 179 156 100 323 277 277 251 230 230 ... 685 651 651 216 216 216 240 89 0 ...

TEC (h) ... 462 462 462 539 604 697 324 370 425 ... 720 720 720 716 720 720 720 720 720 ...

TRD (h) ... 472 472 472 289 294 0 331 331 331 ... 660 653 648 216 216 216 0 0 0 ...

TRC (h) ... 525 562 618 333 465 465 477 535 583 ... 714 720 720 720 720 720 719 720 720 ...

Table 8: Reversion flow batches for scenario S2 Batch 215 216 217 218 221 222 223 224 225 226 227

6.2

Product 35 31 35 31 31 31 21 31 31 31 31

Route N13-22-N13 N4-22-N4 N13-22-N13 N4-22-N4 N13-22-N13 N4-22-N4 N8-3-N8 N7-3-N7 N8-3-N8 N7-3-N7 N8-3-N8

Volume (vu) 9400 9400 9400 9400 9400 9400 22100 22100 22100 22100 22100

Comparing the final scheduling derived from heuristic ordering and from the MILP base model

At the end of assignment block, a set of batches with time-windows, routes, volumes, and a preliminary suggested order is obtained. So, what happens if we take this set of batches and use, according to the decomposition strategy in Figure 2, this previously ordered result


Page 47 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 9: Plug batches for scenario S2 Batch 232 233 234 235 236 237 238 239 240 241 242

Product 21 21 21 21 34 31 21 35 21 35 35

Route N4-1-N8 N4-1-26-N6 N4-1-26-N6 N4-1-26-N6 N4-1-26-N6 N4-1-N8 N4-1-N8 N8-26-N6 N8-26-N6 N7-3-N8 N7-3-N8

Volume (vu) 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000

as an input to simulation/timing blocks without the MILP sequencing phase? Reference 20 used, in essence, only a heuristic module to generate the batches and their sequence. In order to develop this study, we will use the output generated by the timing block, 1 that is the final scheduling solution, in order to compare the scheduling derived from the heuristic ordering and from the MILP base model. The generated Gantt charts and inventory profiles of Scenario 2 (S2) were chosen to illustrate this comparison. Figure 16 presents the obtained Gantt chart of the final scheduling influenced by the heuristic ordering suggested by the assignment module. The number of pipelines used during the scheduling horizon are vertically presented. The scenario starts on the beginning of month 12 (12/01), and more than 38 days are necessary to complete the scheduling. Different colors are used to represent each involved product. The black vertical line marks a 30-day period since the scheduling beginning. Pipelines 1, 3, 21, and 26 are the ones in which the expected 30-day scenario was more exceeded. As a consequence, a tendency of storage lack/overflow will be observed in Figures 18 to 20. In a similar reasoning, Figure 17 presents the final Gantt chart to the same scenario (S2), but now using the MILP base sequencing model to aid the scheduling solution. The MILP base model provided an optimized ordering at an early sequencing phase (Figure 2). Figure 17 indicates that a series of ordering changes were suggested. The scheduling horizon was reduced to 32 days, that is, less six days in comparison to the heuristic solution



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 16: Gantt chart influenced by the heuristic ordering - Scenario S2

Figure 17: Gantt chart influenced by the MILP base sequencing model - Scenario S2


Page 48 of 63

Page 49 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(Figure 16). In Figures 18 to 20 we can observe that the scheduling proposed in Figure 17 caused a tendency of better inventory management, which is a desirable condition. In turn, the MILP sequencing made a “simplified timing” at an early phase of the decomposition strategy of Figure 2. The MILP base model was able to represent the main scheduling features involved in the pipeline network and suggested a better batch ordering than the one proposed in the assignment block. Thus, the use of a MILP sequencing model can be of great value for aiding the timing block to determine better final scheduling answers.

(a)

(b)

Figure 18: Storage on N11, product 21: (a) Heuristic ordering; (b) MILP sequencing

(a)

(b)

Figure 19: Storage on N11, product 25: (a) Heuristic ordering; (b) MILP sequencing Figures 18, 19, and 20 illustrate the storage profile tendency of area N11 for products 21, 25, and 31, respectively. A comparison between the inventory influenced by the heuristic ordering (a) and by the MILP sequencing (b) is evidence in these figures. These pairs of area/product were chosen because they represent situations with the more critical inventory 49 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Page 50 of 63

(b)

Figure 20: Storage on N11, product 31: (a) Heuristic ordering; (b) MILP sequencing tendencies observed within the scenario. The inventory of other pairs were relatively wellconducted. It is possible to verify that the MILP sequencing aided to fix, or improve, the inventory management by changing batch orderings. We can observe, for instance that Figure 18(b) tends to maintain the storage inside operational levels, while Figure 18(a) presents a lack tendency in product supply.

6.3 6.3.1

Results for MILP sequencing MILP Sequencing: Base model and MI model results

Table 10 indicates the computational results from the proposed MILP base sequencing model, and Table 11 summarizes the computational results from the proposed Minimization of Interfaces (MI) model. The same eight scenarios (S1 to S8) were tested in both models and a 30-day horizon was considered in all scenarios. The tables indicate: total number of variables; number of integer variables; number of constraints; number of batches for the MILP model (in comparison with Table 6, except the batches generated in post-assignment/sequencing step); objective function (expression 1 value); integrality gap (%); computational time, in seconds; and, finally, the number of necessary plugs within the proposed sequencing. Table 10 indicates that the generated large-scale MILP formulations presented more than 11 thousand variables, with a great majority of integer ones, and more than 33 thousand


Page 51 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 10: MILP base model: Results for 30-day scenarios Scenario Variables Integer var. Constraints Batches Obj. Function Gap (%) Time (s) Plugs

S1 11473 7035 33236 280 8.19E5 0 1.95 4

S2 14243 9189 44954 298 2.02E6 0 39.23 11

S3 12399 7950 38608 263 2.36E5 0 416 11

S4 13771 8809 43487 282 5.51E5 0 9.83 5

S5 11670 7294 37167 249 1.94E6 0 11.54 7

S6 13724 8942 45680 264 1.37E6 0 507 8

S7 13777 9039 46745 274 2.70E6 0 11.61 7

S8 14332 9392 47204 284 1.80E6 0 20.25 7

Table 11: MILP MI model: Results for 30-day scenarios Scenario Variables Integer var. Constraints Batches Obj. Function Gap (%) Time (s) Plugs

S1 11885 7419 34584 280 8.22E5 0 39 3

S2 15155 10047 47946 298 2.03E6 0 423.15 6

S3 12831 8344 40036 263 2.43E5 0 982 7

S4 14760 9747 46742 282 5.54E5 0 24.66 1

S5 12399 7974 39537 249 1.95E6 0 532.45 6

S6 14611 9766 48588 264 1.38E6 0.03 1000 5

S7 14134 9361 47889 274 2.71E6 0 18.81 3

S8 14986 10000 49294 284 1.81E6 0 288.5 6

constraints. For scenarios S3 and S6 the model obtained the optimal solution within a running time of, approximately, 7 and 9 minutes. For other scenarios, the optimal solution was reached in less than 40 seconds. Thus, optimal solutions were obtained in a reasonable computational time, less than 10 minutes, for all tested cases. Thus, the MILP base sequencing model computational performance, as part of a tool to aid the decision-making process for 30-day scheduling scenarios (Figure 2), was considered adequate for operational purposes. Table 11 indicates that the generated large-scale MILP formulations presented more than 11 thousand variables, with a great majority of integer ones, and more than 34 thousand constraints. For scenarios S1, S4, and S7 the optimal solution was reached in less than 40 seconds. For scenarios S2, S3, S5, and S8, respectively, the optimal solutions were obtained in, approximately, 7, 16, 9, and 5 minutes. Scenario S6 not converged to optimality in a time limit of 1000 seconds (≈ 17 minutes), but the integrality gap was small (0.03%). We verified that for all testes scenarios with both, Base model and MI model, no violations occurred within the capacity storage level. 51 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

6.3.2

MILP Sequencing: Comparing Base model and MI model

The direct comparison of Tables 10 and 11 indicates that the MI model presented more variables and constraints than the Base model, as expected by the formulation structure itself. The MI model computational time also increased in relation to the Base model. However, the number of necessary plugs for the sequencing proposed by the MI model was smaller than the number indicated by the Base model. The observed reductions were the following: S1 from 4 to 3; S2 from 11 to 6; S3 from 11 to 7; S4 from 5 to 1; S5 from 7 to 6; S6 from 8 to 5; S7 from 7 to 3; S8 from 7 to 6. As already mentioned, the Base model did not have the explicit goal of minimizing plugs; so that, the obtained value for the number of necessary plugs is just a consequence from a sequencing answer that seeks to minimize violation issues, expression (1). On the other hand, the MI model is an extension of the Base model which takes into account the additional goal of minimizing the number of necessary plugs, expression (32). The MI sequencing model computational performance was inferior in comparison to the Base model, but the obtained computational times of few minutes for 30-day scheduling scenarios were still considered adequate for operational purposes. The benefit of reducing the number of plugs tends to be worth the few additional spent minutes.

6.4

Comparing results with a former proposed approach 2

Table 12 presents results obtained by the monolithic MILP approach proposed by ref 2 and the results from the proposed collaborative approach based on the Assignment module followed by the MILP Base sequencing model. Alternatively, the MI model can be also used instead of the Base model. Hardware and software configurations were maintained in both approaches for computational tests. For each scenario (S1 to S8), 2 the considered time horizon (H) in days, the total number of batches (Batches), the execution time (Time) in seconds, and the relative gap (Gap) in % are presented. The scheduling horizon (H) in each approach is a feature to be observed. Within the 52 ACS Paragon Plus Environment

Page 52 of 63

Page 53 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 12: Comparing results from the Monolithic MILP Model 2 and the proposed Collaborative Approach: Assignment Module + MILP Base sequencing model

Approach Monolithic MILP Model 2 Assignment module + MILP sequencing model Approach Monolithic MILP Model 2 Assignment module + MILP sequencing model Approach Monolithic MILP Model 2 Assignment module + MILP sequencing model Approach Monolithic MILP Model 2 Assignment module + MILP sequencing model

H (days) 7 30

H (days) 7 30

H (days) 7 30

H (days) 7 30

Scenario 1 (S1) Batches Time (#) (s) 60 36000 2.81 288 1.95 Scenario 3 (S3) Batches Time (#) (s) * 2.94 281 416 Scenario 5 (S5) Batches Time (#) (s) 60 36000 2.23 270 11.54 Scenario 7 (S7) Batches Time (#) (s) 60 36000 4.06 290 11.61

Gap (%) 76.33 -

H (days) 7 30

0 Gap (%)

H (days) 7

30 0 Gap (%) 48.41 -

H (days) 7 30

0 Gap (%) 52.01 -

H (days) 7 30

0

* Solver execution failure (out-of-memory status)


Scenario 2 (S2) Batches Time (#) (s) 60 36000 5.22 320 39.23 Scenario 4 (S4) Batches Time (#) (s) 60 36000 3.1 301 9.83 Scenario 6 (S6) Batches Time (#) (s) * 4.21 288 507 Scenario 8 (S8) Batches Time (#) (s) * 2.99 297 20.25

Gap (%) 33.99 0 Gap (%) 89.40 0 Gap (%) 0 Gap (%) 0


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

proposed approach (Assignment module + MILP Base sequencing model) the considered H value is 30 days. In the former approach, 2 due the computational burden, scheduling horizons of just 7 days were considered. Even with smaller horizons, the computational time for each scenario was much bigger in comparison with the herein proposed approach. The use of assignment module followed by the MILP model returned solutions with relative gap equal to 0% in all tested scenarios, in less than 600 seconds. For the assignment and sequencing monolithic model presented in ref 2, the execution time was limited to 36000 seconds and all solutions have relative gap greater than 33% (best case, Scenario S2). Additionally, in some tested scenarios (S3, S6, and S8), the solver returned execution failure. Thus, the computational gain obtained with the proposed approach is evidenced in Table 12. The novel procedure is embedded into a decision making tool, allowing the suggestion of operational decisions in a computational time adequate to specialists’ analysis.

7

Final Considerations

This paper presents a new collaborative approach to solve the assignment and the sequencing tasks for pipeline-scheduling networks, a hard combinatorial problem. The novel proposed approach, based on the cooperation of an assignment algorithm with a MILP sequencing model, was able to solve hypothetical and real-world instances that were not solved by a former literature approach based on a monolithic MILP model. 2 In a more general standpoint, the proposed assignment/sequencing collaborative approach can be used to define operational batches with their volumes (lot-sizing) in pipeline networks, as indicated in case studies 1 and 2. The approach also provides an optimized ordering of batches from source nodes to destination nodes based on the solution of a MILP sequencing model. In order to obtain the final short term scheduling of activities by an integrated tool, the approach herein proposed is used in cooperation with a planning model 2 and a MILP timing


Page 54 of 63

Page 55 of 63

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


model, 1 as indicated in Figure 2. The developed tool has been applied to case studies that represent the operational conditions of real multiproduct pipeline networks that transport oil derivatives and ethanol. The considered problems can be particularly complex and involve many nodes and pipelines, as indicated in Figure 15. The difficulty of obtaining a feasible scheduling of operational activities is a day-to-day problem faced by the company schedulers, and the proposed solution approach can work as a decision tool to assist the decision makers. The proposed approach is embedded in a customized computational tool with a series of graphical facilities that aid specialists to visualize scheduling details for generated solutions involving a horizon of, approximately, one month, in a reduced computational time (few minutes). In a first step, decision-makers analyze a short/middle-scale operational scheduling to be implemented. Afterwards, tendencies related to the lack of products in customers and supply excesses in the refineries can be analyzed through violations on time-windows. If necessary, operational decisions are taken in order to avoid such tendencies.

Acknowledgement The authors would like to thank financial support from Brazilian Oil Company PETROBRAS (grant 0050.0066666.11.9) and CNPq (grants 305405/2012-8, 305816/2014-4, and 309119/2015-4).

Supporting Information Available This material is available free of charge via the Internet at http://pubs.acs.org/.

8

Nomenclature



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 13: Sets and indexes Set N D P R B F

Index {n, n0 , nx, nx0 } {d, d0 , dx, dx0 } {p, p0 , px, px0 , pr} {r, r0 } {b, b0 , bx} {f x}

Domain {1..N N } {1..N D} {1..N P } {1..N R} {1..N B} {“CAP ”, “M inM ax”, “T arget”}

Description Areas Pipelines Products Routes Batches Storage levels

Table 14: Parameters Parameter H NN ND NP NR NB Ktof x Ktdf x KtdDf x Kplug T EDb,f x T ECb,f x T RDb,f x T RCb,f x V OLb ROT Ab P RODb V Bombb,d T DBatb,d T BBatb,d N BatDestn0 ,pr N Batd BatEstDutob,d Incompatibilidadep,p0 U >> 0 L

Collaborative Approach Based on Heuristic Algorithm and MILP Model

Recommend Documents