Optimization of Detailed Schedule for a Multiproduct Pipeline Using a

Article pubs.acs.org/IECR

Optimization of Detailed Schedule for a Multiproduct Pipeline Using a Simulated Annealing Algorithm and Heuristic Rules Haihong Chen,† Changchun Wu,*,† Lili Zuo,† Feng Diao,† Li Wang,† Dapeng Wang,‡ and Baoqiang Song‡ †

National Engineering Laboratory for Pipeline Safety/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum, Beijing, China ‡ PetroChina Oil & Gas Control Center, Beijing, China ABSTRACT: This paper addresses the optimization of delivery/injection plans of intermediate stations along a unidirectional multiproduct pipeline. A mixed-integer linear programming (MILP) model is developed to minimize pump rate variations in every pipeline segment along the pipeline during a specific scheduling horizon. Given the problem complexity, a computational framework is proposed consisting of three main stages: initial solution finding, solution improvement based on a simulated annealing (SA) metaheuristic, and solution refinement based on tailor-made heuristics. The algorithm for making an initial solution is based on a so-called space recursive method that the delivery/ injection schedule of every intermediate station is gradually generated from upstream to downstream segments. The algorithm for generating a new solution by SA consists of three steps: selecting an intermediate station, choosing a batch passing by the selected station, and adjusting the selected operation. For improving the quality of solutions given by SA, the solutions are further fine-tuned based on a heuristic rule on the proper relay of delivery/injection operations carried out by different intermediate stations. Finally, the model and its algorithms are tested and analyzed using historical monthly plans of two multiproduct pipelines.

1. INTRODUCTION With the development of optimization algorithms and computer technology, optimization models (including mathematical models, simulation models, and heuristic models) for studying the hypothetical scheduling of multiproduct pipelines are being gradually applied to real-world multiproduct pipelines. 1.1. Mathematical Models.1−14 The mathematical models that are used for multiproduct pipeline scheduling mainly include discrete models1,2 and continuous models.3−14 The discrete models indicate that the delivery/injection operations occurring at each station along the pipeline are regarded as a series of discrete events. Furthermore, the implementation processes of these discrete events are described by dividing the scheduling horizon into evenly distributed time intervals and dividing each pipeline segment into packs. The discrete mixed-integer linear programming (MILP) model proposed by Rejowski and Pinto1 is a typical discrete model. In this model, the scheduling period was discretized using equal time steps and the volume of a pipeline was discretized using either equal or unequal volumetric steps. The similarity of the two methods for discretizing pipeline segments was that every pipeline segment was segregated into single-product packs. However, the differences were that, in the method based on equal volumetric steps, all pipeline segments had the same volumetric step and only a unique delivery station can carry out a receiving operation at every batch injection, © 2017 American Chemical Society

whereas, in the method based on unequal volumetric steps, different pipeline segments could have different volumetric steps and several delivery stations can simultaneously receive products from in-transit batches at every injection. In contrast to discrete models, continuous models indicate that the delivery/injection operations are regarded as a series of continuous events, such that the scheduling period can be divided into variable time steps. The continuous MILP model proposed by Cafaro and Cerdá4 is a typical continuous model. In this model, the start and end moments of each batch input by the source of a pipeline are selected as a time node to divide the scheduling period, and only a unique delivery station can carry out a receiving operation at every batch injection. Based on the continuous models, Cafaro et al.5−8 and Cafaro and Cerdá9 established many scheduling models over time for different multiproduct pipeline systems, such as a single multiproduct pipeline with a single source and multiple terminals, a single multiproduct pipeline with multiple-sources and multiple terminals, tree-structure multiproduct pipeline networks with multiple input and exit terminals, and mesh-structure multiReceived: Revised: Accepted: Published: 5092

December 8, 2016 March 22, 2017 April 13, 2017 April 13, 2017 DOI: 10.1021/acs.iecr.6b04745 Ind. Eng. Chem. Res. 2017, 56, 5092−5106

Article

Industrial & Engineering Chemistry Research

Figure 1. Schematic diagram of a multiproduct pipeline.

downloaded by a delivery station at a receiving rate less than the flow rate of the batch passing by the station. Compared with fullstream delivery operations, side-stream delivery operations are more economic, because a full-stream delivery operation inevitably makes the downstream adjacent pipeline segment shutdown, which additionally increases the cost for flow restarts in pipeline segments. Unlike the method published by Cafaro et al.,5,6 the method used in this paper applies the simulated annealing (SA) algorithm and heuristic rules to optimize monthly batch plans of a multiproduct pipeline, which allows intermediate stations to carry out side-stream delivery operations. In the literature, pumping costs, costs of reprocessing contaminations, and backorder costs are usually taken as objective functions. The operational stability of a multiproduct pipeline in which the frequency and magnitude of pump rates vary in every pipeline segment along the pipeline has never been considered. With the increasing number of stations along the pipelines, the delivery/injection operations will dramatically increase. If the delivery/injection operations occurring at adjacent stations are not well matched (for example, if one intermediate station terminates a delivery/injection operation and all the other stations do not immediately start a new operation), the magnitude and frequency of pump rate variations of each pipeline segment will dramatically increase. As a result, the costs of operating pipelines will also increase, and even worse, the safety of pipelines also may be badly affected. Therefore, the operational stability of a multiproduct pipeline is selected as the optimization target in this paper.

product pipeline networks with several input and receiving terminals. Although discrete models are usually intuitive and wellunderstood, the scale and computational burden of discrete models are more onerous than those of continuous models if the time and volumetric steps are too small. Generally, the discrete models can produce only suboptimal solutions, whereas the continuous models can produce optimal solutions. However, establishing continuous models is usually difficult. The emphasis of establishing a continuous model is how to track the moving processes of oil batches in the pipeline during the scheduling horizon. 1.2. Simulation Models. In simulation models, the batch plans of multiproduct pipelines are drafted with the help of simulation software. Given an objective function, the draft batch plan is continuously adjusted by algorithms such as the greedy algorithm and contraindication search algorithm15 until a feasible and better batch plan has been achieved. Methods for simulating batch plans of multiproduct pipelines include the special-order method, the elemental-clock method (which is the most flexible method), the master-clock method, the average-rate method and the displacement method.16,17 1.3. Heuristic Models.18−20 In heuristic models, dispatchers make draft batch plans of multiproduct pipelines based on common rules and practical experience. Generally, several improved batch plans can be simultaneously generated instead of the best batch plan. Although the best batch plan is not achieved, the improved batch plans can meet the flexibility and economy requirements of real-world practical applications. Considering the computational times and burdens, most discrete mathematical models that appear in the literature for scheduling of multiproduct pipelines are more applicable for short scheduling periods of approximately a week. Different from the discrete models, Cafaro et al.11 proposed an MILP continuous-time framework for the dynamic scheduling of pipelines over a rolling horizon of four, 1-week periods. However, the dynamic optimization approach of Cafaro et al. can generate only the overall input schedule and the set of aggregate product deliveries to be performed from in-transit lots to depots at every batch injection. Thus, given the aggregate schedule, Cafaro et al.5 proposed a continuous-time planning approach and three different heuristic rules to make the detailed schedule of a single-source pipeline with multiple distribution terminals, which seeks to minimize flow restart and on−off pump switching costs over the planning horizon. However, the MILP model established by Cafaro et al.5 can allow only a unique delivery station to carry out full-stream operations at every batch injection of the initial station. Then, a later work published by Cafaro et al.6 presents a new mathematical programming formulation for the detailed scheduling of a single-source trunk line with multiple delivery stations that can simultaneously carry out side-stream delivery operations. Here, full-stream delivery means that a delivery station receives product from an in-transit batch by the flow rate of the batch passing by the station, whereas side-stream delivery indicates that part of an in-transit batch is

2. PROBLEM DESCRIPTION 2.1. Assumptions. The following assumptions can be made, regarding the work in this paper: (1) A unidirectional multiproduct pipeline is addressed, as depicted in Figure 1. The source of the pipeline has a unique initial station that injects products into the pipeline. The ending of the pipeline has a unique terminal station that receives products from the line. The intermediate stations can be a delivery station or an injection station. (2) All the pipeline segments are always full with incompressible products. (3) When batches move in the pipeline, contamination will occur inevitably between the two adjacent products. The mixed product section is considered to be an interface. (4) Whenever the initial station is injecting a batch to the pipeline, several delivery stations can simultaneously receive products from the trunk pipeline. (5) Full-stream delivery operations carried out by an intermediate station inevitably make the downstream adjacent pipeline segment shut down. Considering the fact that the cost for flow restarts in pipeline segments is high, full-stream delivery operations are forbidden at any intermediate delivery station. 5093

DOI: 10.1021/acs.iecr.6b04745 Ind. Eng. Chem. Res. 2017, 56, 5092−5106

Article


(1) Constraints of Demands. The total volume that any product p ∈ P is delivered/injected by any intermediate station j (for j < Jmax) must be equal to the demand for product p at station j:

(6) For any intermediate injection station, when the station carries out a full-stream injection operation, the upstream adjacent pipeline segment must be under the idle state. Considering the cost for flow restarts in pipeline segments is high, full-stream injection operations are forbidden at any intermediate injection station. (7) During the scheduling horizon, the overall delivery volumes of products that intermediate delivery stations need download from the pipeline and the overall injection volumes of products that intermediate injection stations need inject into the pipeline are prior given. (8) The inventory constraints of station tanks on batch plans are neglected. Whenever delivery stations receive products from the pipeline, the station tanks are assumed to have enough available storage capacity to store these products that are from the pipeline. Whenever injection stations inject products into the pipeline, the station tanks are assumed to have enough available storage products that are supplied to the pipeline. 2.2. Known Conditions. The following conditions are described for the work performed in this paper: (1) The initial line-fill is known, including the sequence of products and the volumes of batches. (2) The input plan of the initial station is known, including the sequences of products and start−end moments and pump rates of all the batches input by the initial station during the scheduling horizon. (3) The overall demand for every product at each intermediate station along the pipeline is known. (4) The location of each station and the inner diameter of each pipeline segment along the pipeline are known. (5) The minimum and maximum pump rates of all the pipeline segments are known. (6) The minimum and maximum flow rates at which each intermediate station delivers/injects each product are known. (7) The minimum and maximum volumes at which each intermediate station delivers/injects each product in every single operation are known. According to the aforementioned conditions, an MILP model is developed to determine the start moment, end moment, flow rate and volume of every delivery/injection operation carried out at each intermediate station along the pipeline. The model is solved using the SA algorithm and heuristic rules. The objective function is to minimize the magnitude of pump rate variations in each pipeline segment along the pipeline during the scheduling horizon.

VP jdemand = ,p

Imax k max − 1

∑ ∑ i=1

yi , p ·Vj , i , k

2 ≤ j ≤ Jmax − 1, p ∈ P

k=1

(1)

(2) Delivering/Injecting Volume Constraints. (a) Intermediate Stations. For intermediate stations, Vj , i , k = Q j , i , k ·Δt

2 ≤ j ≤ Jmax − 1, i ∈ I , k ≤ k max − 1 (2)

To reduce the number of single delivery/injection operations carried out by intermediate station j (for j < Jmax), the volume of any in-transit batch containing product p that station j delivers/ injects must not be less than the minimum delivery/injection volume. yi , p ·Vj , i , k ≥ oj , i , k ·yi , p ·VI jmin ,p 2 ≤ j ≤ Jmax − 1, i ∈ I , p ∈ P , k ≤ k max − 1

(3)

Although the influences of tank farms on batch plans are neglected, the maximum inventory limitations of tank farms connected to intermediate stations should be considered. Therefore, the total volume of any in-transit batch containing product p that station j delivers/injects must not be greater than the maximum delivery/injection volume. yi , p ·Vj , i , k ≤ oj , i , k ·yi , p ·VI jmax ,p 2 ≤ j ≤ Jmax − 1, i ∈ I , p ∈ P , k ≤ k max − 1

(4)

(b) Terminal Station. The terminal station can be desecribed using VJ

max

,i,k

≤ FJ ·o J max

max

,i,k

i ∈ I , k ≤ k max − 1

(5)

(3) Delivery/Injection Flow Rate Constraints. The delivery/ injection flow rates of operations are limited by pipeline delivery/ injection infrastructures, such as pumps, valves, flow meters, etc. The flow rate of any intermediate station j (for 2 ≤ j < Jmax) delivering/injecting any product p ∈ P must be controlled within the allowable range. max oj , i , k ·yi , p ·QI jmin , p ≤ yi , p · Q j , i , k ≤ oj , i , k · yi , p · QI j , p

2 ≤ j ≤ Jmax − 1, i ∈ I , p ∈ P , k ≤ k max −1

(6)

(4) Delivering/Injecting Conditions. The conditions for an intermediate station j (2 ≤ j < Jmax) delivering/injecting a batch i ∈ I in the time interval (tk,tk+1) can be described as follows. The upper coordinate of batch i at time tk must be greater than the coordinate of station j, and the lower coordinate of batch i at time tk+1 must be less than the coordinate of station j.

3. SCHEDULING MODEL 3.1. Nomenclature. The nomenclature used in this work is detailed in the Nomenclature section near the end of this paper. 3.2. Time Representation. In this paper, the discrete time formulation is used to describe the transportation of batches in the pipeline during the scheduling horizon, in which the scheduling horizon is divided into evenly distributed time steps. 3.3. Constraints. The following constraints are observed in this work: the constraints of demands, the delivering/injecting volume constraints, the delivery/injection flow rate constraints, the delivering/injecting conditions, the constraints of material balance, the constraints of batch tracking, and the pump rate constraints of pipeline segments.

Fj ≤ Hi , k + FJ ·(1 − oj , i , k ) max

2 ≤ j ≤ Jmax − 1, i ∈ I , k ≤ k max − 1

(7)

Hi , k + 1 − Wi , k + 1 ≤ Fj + FJ ·(1 − oj , i , k ) max

2 ≤ j ≤ Jmax − 1, i ∈ I , k ≤ k max − 1 5094

(8)


Article

Industrial & Engineering Chemistry Research In any time interval (tk,tk+1), any intermediate station cannot simultaneously deliver/inject more than one batch.

QAj =

Imax

∑ oj , i ,k ≤ 1

2 ≤ j ≤ Jmax − 1, k ≤ k max − 1

i=1

max

max

+ Hi , k + 1 ≥ FJ

,i,k)

max

k max − 1

∑

QDj , k

j ≤ Jmax − 1

(21)

k=1

If intermediate station j is a delivery station, the pump rate in the upstream adjacent pipeline segment j − 1 must be not less than that in the downstream adjacent pipeline segment j. If intermediate station j is an injection station, the pump rate in the upstream adjacent pipeline segment j − 1 must be not greater than that in the downstream adjacent pipeline segment j.

(9)

The conditions for the terminal station Jmax delivering a batch i ∈ I in time interval (tk,tk+1) can be described as follows. The upper coordinate of batch i at tk+1 must be equal to the coordinate of station Jmax. FJ ·(1 − o J

1 · k max − 1

bj ·QDj − 1, k ≤ bj ·QDj , k

i ∈ I , k ≤ k max − 1

2 ≤ j ≤ Jmax − 1, k ≤ k max − 1

(10)

(22)

(5) Constraints of Material Balance. At any time interval, the sum of (a) the overall volume of all the batches injected by the initial solution and (b) the total volume of all the batches delivered/injected by the intermediate and terminal stations must be zero (0).

To ensure safe operations of a pipeline, the pump rates of any pipeline segment along the pipeline must be not greater than the maximum pump rate. Moreover, considering economic operations of a pipeline, the pump rates of any pipeline segment along the pipeline must be not less than the minimum pump rate (>0). Since the cost for flow restarts in pipeline segments is high, all the pipeline segments are always kept active during the scheduling horizon. For an intermediate delivery station, when the station carries out a full-stream delivery operation, the downstream adjacent pipeline segment will be shut down and its pump rate will become zero, which violates constraint (23). For an intermediate injection station, when the station carries out a full-stream injection operation, the upstream adjacent pipeline segment must be under the idle state and its pump rate must be zero, which violates violates constraint (23). Here, assumptions (5) and (6) in Section 2.1 indicate that the full-stream delivery and full-stream injection operations that are forbidden at any intermediate station can be satisfied.

Imax Jmax

Imax

∑ Vini ,k+ ∑ ∑ bj ·Vj ,i ,k = 0 i=1

k ≤ k max − 1 (11)

i=1 j=2

(6) Constraints of Batch Tracking. The volume of any batch i in the pipeline at any moment tk must equal the sum of (a) the volume of the batch i at tk−1 and (b) the volume of batch i input by the initial station during period (tk−1, tk), and (c) the total volume of the batch i delivered/injected by the intermediate and terminal stations during period (tk−1, tk). Jmax

WBi , k = WBi , k − 1 + Vini , k − 1 +

∑ bj ·Vj ,i ,k− 1 j=2

i ∈ I , 2 ≤ k ≤ k max

QSjmin ≤ QDj , k ≤ QSjmax

(12)

j ≤ Jmax − 1, k ≤ k max − 1 (23)

WBi , k = WBi0 WBi , k = 0

i ∈ I old , k = 1

i∈I

new

,k=1

(13)

3.4. Objective Function. To ensure that a multiproduct pipeline is smoothly operated, the magnitude of pump rate variations in each pipeline segment along the pipeline during the scheduling horizon should be minimized. Achieving this target can reduce the frequency of pressure variations of each pipeline segment, and the times that pump stations along the pipelines are restarted and shut down also can be reduced. Obviously, the cost of operating a pipeline also will be decreased.

(14)

The upper volumetric coordinate of any batch i at any moment tk must equal the sum of the upper volumetric coordinate of batch i + 1 at time tk and the volume of batch i at time tk. Hi , k = Hi + 1, k + WBi , k Hi , k = WBi , k

1 ≤ i < Imax , k ≤ k max

i = Imax , k ≤ k max

Hi , k = Fj

i = 1, k ≤ k max , j = Jmax

Hi , k ≤ Fj

2 ≤ i ≤ Imax , k ≤ k max , j = Jmax

(15)

Jmax − 1 k max −1

(16)

min f =

j=1

(17) (18)

j = 1, k ≤ k max − 1

i=1

Imax

QDj , k =

j

J

(19)

Imax

j ′= 2 i = 1

2 ≤ j ≤ Jmax − 1, k ≤ k max −1

cpj ·|QDj , k − QDj , k + 1|

(24)

k=1 −1

k

−1

max max The intent of using the expression ∑j=1 ∑k=1 cpj·|QDj,k − QA j| is to reduce the summation of overall deviations between the real-time pump rates and the mean pump rate in each pipeline max−1 max−2 segment. The expression ∑Jj=1 ∑kk=1 cpj·|QDj,k − QDj,k+1| is used to reduce the summation of overall deviations of the realtime pump rates between adjacent periods in each pipeline segment, which can avoid wild variations in the pump rates of the pipeline segments in a short period. The weight factor cpj is used to avoid the pump rate fluctuations being concentrated in one pipeline segment, although, overall, they are minimized. Since the objective function is nonlinear, two auxiliary variablesf 1j,k and f 2j,kare introduced to linearize the objective function.

∑ Qini ,k+ ∑ ∑ bj ′·Q j ′ ,i ,k i=1

∑ ∑ j=1

Imax

∑ Qini ,k

cpj ·|QDj , k − QAj|

k=1


+

(7) Pump Rate Constraints of Pipeline Segments. At any time interval, the pump rate of pipeline segment j must be equal the sum of the injecting rate of the initial station and the overall delivering/injecting flow rate of intermediate stations 2−j. QDj , k =

∑ ∑

(20) 5095


Article

Industrial & Engineering Chemistry Research Jmax − 1 k max −1

min f =

∑ ∑ j=1

cpj ·f 1j , k +


k=1

∑ ∑ j=1

−f 1j , k ≤ QDj , k − QAj ≤ f 1j , k

cpj ·f j2, k

k=1

station from upstream to downstream. The processes involved in making an initial solution are shown in Figure 2. (25)

j ≤ Jmax − 1, k ≤ k max − 1 (26)

− f j2, k ≤ QDj , k − QDj , k + 1 ≤ f j2, k

j ≤ Jmax − 1, k ≤ k max − 2 (27)

4. SOLVING ALGORITHM 4.1. General Approach. The developed model is a mixedinteger linear programming (MILP) model, which can be solved using CPLEX. However, in this paper, the SA algorithm is also used to solve the same problem that is addressed by the MILP formulation. In a robust SA algorithm, eight essential aspects must be considered, including a fitness function, how to make an initial solution, how to generate new solutions, the acceptance criteria for new solutions, the initial temperature, the final temperature, the cooling schedule, and the length of the Markov chain. Each of these criteria is discussed in the following sections. 4.2. Fitness Function. The penalty function method is used to transform constrained optimization problems to nonconstrained optimization problems, because the SA algorithm cannot be used directly to solve constrained optimization problems. Although the initial solution and new solutions can easily satisfy constraints (1)−(22), the pump rate constraint (constraint (23)) of pipeline segments are difficult to satisfy. Therefore, the pump rate constraints of pipeline segments are transformed to the fitness function. According to the objective function (constraint (24)) and constraint (23), the fitness function is expressed as described in eqs 28 and 29.

Figure 2. Processes involved in making an initial solution of the proposed model.

To reduce the complexity of the solving algorithm, any intermediate station is assumed to deliver/inject any batch at most once during the scheduling horizon, while the flow rate of any delivery/injection operation must be constant. Moreover, to better illustrate the solving algorithm, the variables qstation , vstation , j,i j,i station end head tail pass dtstation , t , q , tt , tt , and W are used. In the processes j,i j,i j,i j,i j,i j,i involved in making an initial solution, the methods used to determine the initial values of some decision variables are described as follows. (1) xj,i: Determine whether the batch i (yi,p = 1) is suitable to be delivered/injected at intermediate station j, according to the following conditions: (a) station j has demand for product p; (b) batch i can arrive at station j during the batch injections; and (c) the volume of station j delivering/injecting batch i must not be less than the minimum delivery/injection volume. (2) Variable vterminal (yi,p = 1) is determined by the rule that the j,i demand for product p at station j is positively allocated to each batch containing product p, based on the ratio of volume Wpass j,i of each batch containing the product p that passes by station j during the scheduling horizon. (3) Variable qterminal (yi,p = 1) is preliminarily set as (QImin j,i j,p + QImax )/2. If the volume that station j delivers/injects into j,p max terminal product p by (QImin j,p + QIj,p )/2 cannot meet the demand, qj,i min max is added among (QIj,p ∼ QIj,p ) until the demand can be met. end (4) tstart j,i ,tj,i : All delivery/injection operations are carried out at the middle of the periods when batches pass by stations, which ensures that the adjustable space of operations is wider. For example, if the period of a batch passing by an intermediate station is (100−200 h), the volume of the batch delivered by the

Jmax − 1 k max − 1

min f =

∑ ∑ j=1

cpj ·|QDj , k − QAj|

k=1


+

∑ ∑ j=1

cpj ·|QDj , k − QDj , k + 1| + M ·U

k=1

(28)


U=

∑ ∑ j=1

[|min(QDj , k − QSjmin , 0)|

k=1

+ |max(QDj , k − QSjmax , 0)|]

(29)

The parameter M is a penalty factor, which should be positively infinite. Notably, if the M value is set too small, the final solution may violate the pump rate constraints of pipeline segments, and a feasible solution may not be obtained. According to the results of the case study described subsequently in Section 5, the M value is recommended to be >105. 4.3. Making an Initial Solution. The quality of the initial solution not only affects the computational time of an intelligent algorithm but also affects the quality of the final solution. To obtain a better initial solution, the space recursive method is used. The delivery/injection operations that occur at the upstream stations affect the timing of batches arriving at downstream stations, but the delivery/injection operations that occur at downstream stations have no effect on the delivery/ injection operations that occur at upstream stations. In other words, the batch plan has no aftereffects on space. The space recursion method takes advantage of this aftereffects feature and is used to gradually generate the batch plan of each intermediate 5096


Article


Very Fast Simulated Re-Annealing algorithm.21 Giving the same chance to adjust the time, the volume and the flow rate variables ensure that each variable can be fully adjusted and make the final solution closer to the optimal solution.

station is 2000 m3, and the delivery rate is 100 m3/h, then the station will deliver the batch in the period of (140−160 h). The specific steps for making an initial solution are described as follows. Step 1. Read the batch plan of the initial station. Then, set j = 2. Step 2. According to the batch plans of the upstream stations, calculate the detailed information on each batch passing by tail pass station j, including tthead j,i , ttj,i , and Wj,i . Step 3. Set p = 1. Step 4. If Vdemand = 0, go to Step 10. j,p Step 5. Determine whether station j delivers/injects each head batch i (∀i ∈ I, yi,p = 1): If tttail j,i − ttj,i > 0, set xj,i = 1; otherwise, set xj,i = 0. Step 6. Set vj,istation to vjstation = xj , i· ,i

W jpass ,i I

max ∑i ′= y · x · W jpass ,i′ 1 i′,p j,i′

·V jdemand ,p

⎧ α′ = α + kk·(B − A) ⎪ ⎪ ⎡⎛ ⎤ ⎞|2u ′− 1| ⎨ 1 ⎢ − 1⎥⎥ ⎪ kk = sgn(u′ − 0.5) ·Tf ·⎢⎜⎜1 + ⎟⎟ Tf ⎠ ⎪ ⎝ ⎣ ⎦ ⎩

In eq 30, α is the decision variable (α ∈ [A,B]). The variable α′ represents the neighborhood solution; Tf is the temperature parameter in the SA algorithm; u′ is a random number within the range of 0−1, and sgn is the sign function. The specific steps for generating a new solution are described as follows. Step 1. Randomly select an intermediate station sj(sj ∈ J, sj < Jmax). Step 2. Randomly select a batch si ∈ I (ysi,sp = 1) passing by station sj. Step 3. If xsj,si = 0, go to Step 4; otherwise, go to Step 5. Step 4. Allow station sj to deliver/inject batch si and set the volume, flow rate, and time of the operation using eqs 31.

∀ i ∈ I , yi , p = 1

station Step 7. If vstation < VImin j,i j,p (∃i ∈ I, xj,i = 1, yi,p = 1), set xj,i = 0, vj,i = 0, and go to Step 6. Step 8. For ∀i ∈ I, yi,p = 1, set max ⎧ x ·(QI jmax , p + QI j , p ) ⎪ q station = j , i 2 ⎪ j ,i ⎪ station ⎪ station xj , i·vj , i = ⎪ dt j , i qjstation ,i ⎪ ⎪ ⎨ ⎡ tt tail + tt head ⎤ dt jstation ,i j ,i ⎪ start head ⎥ ⎢ j ,i = − t x max , tt , , , j i j i j i ⎪ ⎢⎣ ⎥⎦ 2 2 ⎪ ⎪ ⎡ tt tail + tt head ⎤ ⎪ end dt jstation ,i j,i j,i ⎥ + , tt jtail ⎪ t j , i = xj , i max⎢ ,i ⎢⎣ ⎥⎦ 2 2 ⎪ ⎩

⎧ xsj , si = 1 ⎪ Pmax ⎪ ⎪ vsjstation = xsj , si· ∑ ysi , p ·VIsjmin , si ,p ⎪ p=1 ⎪ Pmax ⎪ ⎪ q station = xsj , si· ∑ y ·QIsjmin ,p si , p ⎪ sj , si p=1 ⎪ ⎪ station ⎨ station xsj , si·vsj , si d t = sj , si ⎪ qsjstation , si ⎪ ⎪ ⎤ ⎡ (tt tail + tt head) ⎪ start dtsjstation j,i j,i , si ⎥ , tt jtail − ⎪ tsj , si = xsj , si·min⎢ ,i ⎥⎦ ⎢⎣ 2 2 ⎪ ⎪ ⎤ ⎡ (tt tail + tt head) ⎪ dtsjstation j,i , si tail ⎥ ⎢ j,i ⎪ tsjend x min , tt = · + sj , si j,i ⎥⎦ ⎢⎣ ⎪ , si 2 2 ⎩

Step 9. For ∀i ∈ I, yi,p = 1, set qjstation ,i

(30)

⎡ ⎛ ⎞⎤ V jdemand ,p ⎢ station ⎜ max ⎟⎥ = xj , i max⎢qj , i , min⎜ i ′= I , QI j , p ⎟⎥ station max ⎢⎣ ⎝ ∑i ′= 1 yi ′ , p ·xj , i ′·dt j , i ′ ⎠⎥⎦

Step 10. Set p = p + 1. If p ≤ Pmax, go to Step 4. Step 11. Set j = j + 1. If j ≤ Jmax − 1, go to Step 2; otherwise, exit the loop. As seen in the processes involved in making an initial solution, the reason why a batch cannot be delivered/injected by an intermediate station is either that the volume of the batch passing by the station is too low or the station has no demand for the product that the batch contains. 4.4. Generating New Solutions by the Simulated Annealing Algorithm. 4.4.1. Construction Process. The framework for generating a new solution consists of three steps: (i) randomly select an intermediate station; (ii) randomly select a batch passing by the selected station; and (iii) determine whether to add a new delivery/injection operation or adjust the selected original delivery/injection operation. With regard to step (iii): (a) If the selected batch has not been delivered/injected by the selected station, then add a new delivery/injection operation. (b) Otherwise, adjust the originally selected operation. Randomly select a variable from the continuous variables tstart j,i , vstation , qstation by equal probability and adjust the selected variable j,i j,i by the formulations expressed by eq 30, which comes from the

(31)

Step 5. Generate a random number u within the range 0−1. If 1 ≥ u > 2/3, go to Step 6. If 2/3 ≥ u > 1/3, go to Step 7. Otherwise, go to Step 8. Step 6. Keeping the variables vstation and qstation constant, sj,si sj,si randomly adjust the start and end moments of the selected tail operation xsj,si = 1 within the range [tthead sj,si , ttsj,si], using eq 32. ⎧ t ′start = t start + kk ·(tt tail − tt head) sj , si sj , si sj , si ⎪ sj , si ⎪ ⎡⎛ ⎤ ⎞|2u ′− 1| ⎪ 1 ⎢ − 1⎥⎥ ⎪ kk = sgn(u′ − 0.5)·Tf ·⎢⎜⎜1 + ⎟⎟ Tf ⎠ ⎨ ⎣⎝ ⎦ ⎪ ⎪ station ⎪ t ′end = t ′start + vsj , si sj , si ⎪ sj , si qsjstation , si ⎩

(32)

Then, go to Step 9. end Step 7. Keeping the variables tstart sj,si and ttsj,si constant, randomly adjust the volume of the selected operation xsj,si = 1, using eq 33. 5097


Article

Industrial & Engineering Chemistry Research ⎧ v′station = v station + 3·q station ·sgn(u′ − 0.5) sj , si sj , si ⎪ sj , si ⎪ ⎨ v′station sj , si ⎪ q′station = sj , si end ⎪ (tsj , si − tsjstart , si ) ⎩

moment of the operation performed by station j, when the start moment of the operation performed by station j + 1 is being randomly adjusted. To avoid pump rates of pipeline segments varying in a very short period, a heuristic rule derived from dispatchers is introduced into the process of generating a new solution. The rule requires that once an intermediate station terminates a delivery/injection operation, another intermediate station should immediately start a new operation. The algorithm for implementing the proper relay of operations is described as follows. (1) After adjusting the selected operation xsj,si in the processes of generating a new solution, attempt to adjust the selected operation by a relay to delivery/injection operations performed by other stations. (2) Randomly choose the relay direction. If u > 0.5, an upstream station will be randomly selected and the operation xsj,si will be adjusted by a relay to an operation performed by the selected upstream station; otherwise, a downstream station will be randomly selected and the operation xsj,si will be adjusted by a relay to an operation performed by the selected downstream station. (3) The situations that the operation xsj,si must be adjusted by a relay to the operation xrj,ri performed by station rj are shown in Figure 3. The red block indicates the selected operation xsj,si, the

(33)

Then, go to Step 9. end Step 8. Keeping the variables tstart sj,si and tsj,si constant, randomly adjust the flow rate of the selected operation xsj,si = 1 within the max range [QImin sj,sp, QIsj,sp ] using eq 34. ⎧ q′station = q station + +kk·(QI max − QI min) sj , sp sj , sp sj , si ⎪ sj , si ⎪ ⎡⎛ ⎞|2u ′− 1| ⎪ 1 ⎢ ⎨ kk = sgn(u′ − 0.5)·Tf · ⎜⎜1 + ⎟⎟ − ⎢ Tf ⎠ ⎪ ⎝ ⎣ ⎪ ⎪ station station end start ⎩ v′sj , si = q′sj , si ·(tsj , si − tsj , si )

⎤ 1⎥⎥ ⎦ (34)

Then, go to Step 9. Step 9. After adding a new operation or adjusting the originally selected operation, the completion of the demand for product sp at station sj must be evaluated again. If the demand is exceeded (or unsatisfied), the delivery/injection flow rate of each operation (∀i′ ∈ I, yi,sp = 1, xsj,i′ = 1) must be gradually reduced max (or increased) within the range [QImin sj,sp, QIsj,sp ] until the demand for product sp has been just met. Step 10. After adding a new operation or adjusting the selected operation, the moments that some batches reach and leave the downstream stations and the pump rates of downstream pipeline segments will change. Thus, a new batch plan must be simulated tail to calculate the detailed information (tthead j,i , ttj,i ) of each batch passing by each station and the change law of pump rate varying with time in each pipeline segment. 4.4.2. Heuristic Rule. During the above process of generating a new solution, an unsatisfactory solution may arise in which the pump rate of a certain pipeline segment varies in a very short period. As a result, the operational scheme of the pump station located at the end of the pipeline segment may have to change frequently in a short period, which is obviously not economical. In reality, once an intermediate station terminates a delivery/ injection operation, another station should immediately start a new operation. By implementing the relay of delivery/injection operations performed by different stations, the pump rates of downstream pipeline segments can be stabilized. For example, assume that the inlet flow rate of station j during the period (100 h, 150 h) is 1000 m3/h and there are two different batch plans: • Case A, where two sequential delivery operations are performed by station j at (100 h, 125.66 h) and (127.88 h, 150 h), and the flow rates are 200 m3/h. • Case B, where two sequential delivery operations are performed by station j at (100 h, 125.66 h) and (125.66 h, 147.78 h), and their flow rates are 200 m3/h. The pump rate of the pipeline segment between station j + 1 and station j + 2 undergoes two changes in Case A (at 125.66 and 127.88 h), but only one change in Case B (at 125.66 h). The operations performed by stations j and j + 1 are better matched in Case B than in Case A, and Case B is more reasonable. However, the relay of operations in Case B is very difficult to be implemented directly using the SA algorithm. According to the formulations for generating a neighborhood solution of a continuous variable in the SA algorithm, it is almost impossible to randomly generate a number that is just equal to the end

Figure 3. Situations for which operations must be synchronized.

blue block indicates the operation xrj,ri, and the left and right sides of a filled block represent the start and end moments, respectively, of an operation. I. Situation a−b: Set tsj,si ′start = tstart ′end = tsj,si ′start + dtstation rj,ri , tsj,si sj,si . start end end II. Situation c−d: Set tsj,si ′ = trj,ri , tsj,si ′ = tsj,si ′start + dtstation sj,si . station end start end III. Situation e−f: Set t′sj,si = tend rj,ri , t′sj,si = t′sj,si − dtsj,si . start end start end IV. Situation g−h: Set t′sj,si = trj,ri , t′sj,si = t′sj,si − dtstation sj,si . V. Situation i: Set tsj,si ′start = tstart ′end = tsj,si ′start + dtstation rj,ri , tsj,si sj,si . When the time deviation of two delivery/injection operations performed by two different intermediate stations is less than τ, the relay of the two operations will be implemented according to the above formulas for different situations. In this paper, the parameter τ is set to be 8 h. The improved algorithm for generating a new solution is shown in Figure 4. 4.5. The Acceptance Criterion for New Solutions. After adding a new operation or adjusting the selected operation, the moments that some batches reach and leave the downstream stations will change. As a result, the new solution may become unfeasible in that the start and end moments of the operations performed by the downstream stations may no longer satisfy the time constraints of delivery/injection operations. Thus, before determining whether to accept a new solution, the feasibility of the new solution must be judged. If a new solution does not satisfy constraints (1)−(21), then refuse it; otherwise, determine 5098


Article


Figure 4. Improved simulated annealing algorithm for generating a new solution.

whether to accept the new solution by using the Metropolis rule, which is described as follows. (I) If the new solution is better than the old solution, accept the new solution. (II) Otherwise, accept the new solution by acceptance probability. The formula for calculating the acceptance probability is pro = exp[−(f ′ − f )/Tf ], wherein f′ and f respectively indicate the fitness function values of the new solution and the old solution. Tf represents the current temperature of the SA algorithm. 4.6. Iteration Termination Conditions. According to the SA algorithm theory, the outer loop and the inner loop are used to control the number of iterations. The outer loop controls the number of iterations by setting the initial temperature, final temperature, and cooling schedule. The inner loop controls the number of iterations by setting the length of the Markov chain. In this paper, the parameters in the SA algorithm are set as follows. (I) Initial temperature: (Tf)max = 1000. (II) Final temperature: (Tf)min = 1. (III) Cooling schedule: Tf ′ = r·Tf. Parameter r represents a cooling factor with a recommended range of 0.75−0.9, which will be discussed by the case study in Section 5. (IV) The length of Markov chain: num = 1000.

Figure 5. Schematic diagram of the multiproduct pipeline XA. D0−D3 are stations; B1−B3 are batches of products; P1 is #0 diesel, and P2 is #93 gasoline.

Table 1. Demand for Each Product at Every Intermediate Station along Pipeline XA (m3) Demand (m3) station D1

station D2

11 700 7800

7800 5200

P1 P2

Table 2. Input Plan of the Initial Station of Pipeline XA during 0−260 h batch product B4 B5 B6 B7

5. CASE STUDY In this section, the MILP formulation and the SA algorithm are respectively applied to determine detailed schedules of two realworld multiproduct pipelines, called XA and XB. Pipeline XA is presented as Case 1 in Section 5.1; pipeline XB is presented as Case 2 in Section 5.2. In Section 5.3, the influences of the cooling factor on the converging processes, the final objective function values, and the computational times are analyzed. Section 5.4 presents statistical results about the computational times for different instances based on the improved SA algorithm. In Cases

P1 P2 P1 P2

start moment (h)

end moment (h)

flow rate (m3/h)

volume (m3)

0 70 150 220

70 150 220 260

1000 1000 1000 1000

70 000 80 000 70 000 40 000

1 and 2, the weight factors cpj of all intermediate station are set to be 1. 5.1. Case 1. The model was first tested on the multiproduct pipeline XA illustrated in Figure 5; the pipeline consists of four stations (D0−D3) and can simultaneously transport the oil products #0 diesel (P1) and #93 gasoline (P2). The total length of the pipeline is 417 km, and the inner diameter is 642 mm. The pump rate range of each pipeline segment along the pipeline is 400−1200 m3/h. The delivery rate ranges of stations D1 and D2 are 100−200 m3/h, and the volume ranges at which stations D1 5099


Article


Figure 6. Transportation of batches in pipeline XA for the solution given by the improved simulated annealing algorithm.

Figure 7. Pump rate profiles of all pipeline segments along pipeline XA for the solution given by the improved simulated annealing algorithm.

Table 3. Solutions for Cases 1 and 2 Given by CPLEX and the Improved SA Algorithm example

methodology

No. of constraints

No. of binary variable

No. of continuous variables

objective function (m3)

CPU time (s)

gap

Case 1

MILP SA

3733

546

1522

2525 2560

1.08 8.39

0%

Case 2

MILP SA

38 764

5904

14 292

36 368 39 364

360.13 60.19

0%

Figure 8. Schematic diagram of the multiproduct pipeline XB. D0−D6 are stations; B1−B5 are batches of products; P1 is #0 diesel, P2 is #93 gasoline, and P3 is #97 gasoline.

moments, the flow rates, and the delivery volumes. The initial line-fill consists of three batches: (B1, P2, 16 504 m3), (B2, P1, 53 709 m3), and (B3, P2, 64 710 m3); these are depicted in Figure 5. Given the input schedule of the initial station, the demands of all of the intermediate stations, and the initial line-fill of the

and D2 deliver each product in every single operation are 1000− 20 000 m3. Table 1 indicates the demand for each product at every intermediate station along pipeline XA. Table 2 indicates the input plan of the initial station (D0) during the scheduling horizon (0−260 h), including the sequences, the start and end 5100


Article


operations in turn. For example, station D2 terminates a receiving operation at 20 h and station D1 immediately starts a new receiving operation. Then, station D1 ends a receiving operation at 40 h, which is followed by a new operation performed by station D2. Stations D1 and D2 continue to take turns performing receiving operations until the end of the scheduling horizon. As seen in Figure 7, the good bridging of receiving operations keeps the pump rate of the pipeline segment D2−D3 at a constant value (870 m3/h) during the scheduling horizon. As a comparison to the improved SA algorithm developed in this paper, the proposed MILP model is used to solve Case 1 with the help of CPLEX 12.6.3. Before solving the MILP model using CPLEX, the scheduling horizon is divided into evenly distributed time steps (Δt = 10). The optimal solution is also given by a MILP solved using CPLEX, as seen in Table 3. The computational time to determine the final solution given by the improved SA algorithm is 7.31 s longer than the time required to determine the optimal solution given by CPLEX, while the objective function value for the final solution given by the improved SA algorithm is just 1.4% higher than that of CPLEX. In the solution given by the improved SA algorithm, the pump rate variations are fully concentrated in the pipeline segment D1−D2, whereas those distributed in two pipeline segments D1−D2 and D2−D3 for the optimal solution are given by CPLEX. 5.2. Case 2. In this section, a historical monthly plan of pipeline XB is tested, which consists of seven stations (D0−D6; see Figure 8) and can simultaneously transport three oil products: #0 diesel (P1), #93 gasoline (P2), and #97 gasoline (P3). The total length of the pipeline is 1840 km, and the inner diameter is 540 mm. The pump rate of each pipeline segment along the pipeline ranges from 550 m3/h to 1450 m3/h. The delivery/injection rate ranges of all of the intermediate stations are 110−250 m3/h, and the volume ranges that all of the intermediate stations deliver/inject each product in every single operation are 1000−40 000 m3. Table 4 indicates the demand for each product at every intermediate station along pipeline XB. Table 5 presents the

Table 4. Demand for Each Product (P1−P3) at Every Intermediate Station (D1−D5) along Pipeline XB Demand (m3) P1 P2 P3

D1

D2

D3

D4

D5

51 900 6900 0

21 500 2700 3000

29 900 23 500 0

3000 3960 0

29 600 8400 0

Table 5. Input Plan of the Initial Station of Pipeline XB during 0−830 h batch

product

start moment (h)

end moment (h)

flow rate (m3/h)

volume (m3)

B5 B6 B7 B8 B9 B10 B11 B12

P1 P2 P1 P2 P1 P2 P3 P1

0 20 70 250 310 540 590 630

20 70 250 310 540 590 630 820

1100 1100 1100 1100 1100 1100 1100 1100

22 000 55 000 198 000 66 000 253 000 55 000 44 000 209 000

pipeline, the detailed delivery plans of stations D1 and D2 are obtained using the improved SA algorithm and are shown in Figure 6, which visually describes the transportation of batches in the pipeline during the scheduling horizon. In Figure 6, the abscissa indicates time and the vertical axis represents distance. The names of all the stations are shown on the left of the diagram and the positions of these stations are on the right. The black oblique lines represent the trajectory of the contamination segment between two adjacent batches. The short bold bars represent delivery/injection operations. The labels under short bold bars represent the start moments, end moments, and flow rates of operations. For example, the red bold bar (0/20/130) in the upper left corner of the diagram denotes that station D2 receives batch B2 in the period (0−20 h) at a flow rate of 130 m3/ h. As seen in Figure 6, stations D1 and D2 perform receiving

Figure 9. Transportation of batches in pipeline XB for the solution given by the improved SA algorithm. 5101


Article


Figure 10. Transportation of batches in pipeline XB for the initial solution.

Given the input schedule of the initial station (D0), the demands of all of the intermediate stations and the initial line-fill of the pipeline, the final solution given by the improved SA algorithm, and the initial solution, are shown in Figures 9 and 10, respectively. As seen in Figure 10, the delivery/injection rates of all operations performed by intermediate stations are equal to the average values of the maximum and minimum delivery/injection rates. The execution times of all operations performed by intermediate stations are arranged at the middle periods of batches passing by stations. For example, station D1 delivers batch B6 in the period (160−171.67 h), while the period of batch B6 passing by station D1 is (140.73−190.94 h). Compared with those in Figure 9, the bold bars in Figure 10 located at station D6 are more uneven, which indicates that the pump rates of the pipeline segment D5−D6 frequently change. By using the improved SA algorithm to adjust the initial solution, the pump rates of the pipeline segment D5−D6 become much smoother. The SA algorithm and the SA algorithm that was improved by introducing a heuristic rule are compared in terms of the objective function value and pump rates in pipeline segments. Figure 11 presents the objective function value changes in iterations for the solutions given by the SA algorithm and the improved SA algorithm. Comparisons of pump rates for all the

Figure 11. Comparison of objective function values given by the simulated annealing (SA) algorithm and the improved SA algorithm.

input plan of the initial station of pipeline XB during the scheduling horizon (0−830 h), including the sequences, the start and end moments, the flow rates, and volumes. The initial line-fill consists of five batches: (B1, P2, 20 970 m3), (B2, P1, 186 395 m3), (B3, P2, 46 599 m3), (B4, P3, 25 630 m3) and (B5, P1, 149 116 m3); these are shown in Figure 8.

Figure 12. Pump rate profile of the pipeline segment D1−D2 given by the simulated annealing (SA) algorithm and the improved SA algorithm. 5102


Article


Figure 13. Pump rate profile of the pipeline segment D2−D3 given by the simulated annealing (SA) algorithm and the improved SA algorithm.




pipeline segments are shown in Figures 12−16. These comparisons are discussed below. (1) In Figure 11, when iteration terminates, the final objective function value given by the improved SA algorithm is ∼8% lower than that given by the SA algorithm. As seen in Figures 12−16, compared to the solution given by the SA algorithm, the magnitudes of pump rate variations in all of the pipeline

segments for the solution given by the improved SA algorithm are almost always lower. (2) For the solution given by the SA algorithm, the pump rates of pipeline segments always drastically change in a short period of several hours. For example, the pump rate of the pipeline segment D3−D4 in the period (64.33−133.64 h) is 990 m3/h, and becomes 1100 m3/h at 133.64 h. After 4.36 h, the pump rates of the pipeline segment D3−D4 revert to 990 m3/h. The pump 5103


Article

Industrial & Engineering Chemistry Research Table 6. Summary of Pump Rate Variations in Every Pipeline Segment for the Three Different Batch Plans: The Initial Solution, the Simulated Annealing (SA) Algorithm, and the Improved SA Algorithm segment

initial solution

solution given by the SA algorithm

solution given by the improved SA algorithm

D0−D1 D1−D2 D2−D3 D3−D4 D4−D5 D5−D6 total

0 14 22 36 40 54 166

0 13 20 32 34 44 143

0 12 14 20 22 16 84

Figure 18. Computational processes for r ranging from 0.05 to 0.35.

rates changing at D3−D4 are unacceptable variations. This and similar situations are avoided using the improved SA algorithm. The number of pump rate variations in every pipeline segment for the initial batch plan, the solution given by the SA algorithm, and the solution given by the improved SA algorithm are compared in Table 6. Compared with those in the initial solution, 14% fewer pump rate variations occur in all pipeline segments for the solution given by the SA algorithm and 49% fewer occur in the solution given by the improved SA algorithm. Generally, introducing the heuristic rule into the SA algorithm reduces both the magnitudes of pump rate variations of pipeline segments and also the number of pump rate variations. As a comparison to the improved SA algorithm developed in this paper, the proposed MILP model is used to solve Case 2 with the help of CPLEX 12.6.3. Before solving the MILP model using CPLEX, the scheduling horizon is divided into evenly distributed time steps (Δt = 10). The results given by CPLEX and the improved SA algorithm are shown in Table 3 and Figure 17. The computational time to determine the optimal solution given by CPLEX is ∼6 times longer than the time required to determine the final solution given by the improved SA algorithm, while the objective function value for the final solution given by the improved SA algorithm is just 7.6% higher than that of CPLEX. In the detail schedule given by the MILP formulation, the total number of single operations performed by all of the intermediate stations is 52, whereas that number is 28 for the detailed schedule given by the improved SA algorithm. Here, if the flow rates that a delivery station receiving product from an in-transit batch or an injection station injecting product to an in-transit batch in two consecutive periods are the same, the two operations can be seen as a single operation. In the improved SA algorithm, to reduce the complexity of the solving algorithm, any intermediate station is assumed to deliver/inject any batch by, at most, one single operation during the scheduling horizon, which reduces the



quality of solutions. However, the MILP formulation has no such a limitation. In reality, pipeline operators sometimes prefer the detailed schedule with lower number of single operations, because their workloads for changing the opening of valves are lighter. The comparison indicates that the improved SA

Figure 17. Evolution of results given by CPLEX and the improved simulated annealing (SA) algorithm. 5104


Article


Through analyzing the influences of different cooling factors, the range of the cooling factor r is recommended to be 0.75−0.9. 5.4. Statistics for Different Instances. The computational burdens of the algorithms proposed in this paper are mainly affected by the number of stations along a pipeline and the number of batches to be delivered. The variations in computational times for variant numbers of stations and batches are shown in Figure 22. When the numbers of stations and the batches increase, the computational time increases linearly rather than exponentially. For the case with 15 stations and 15 batches, the computational time is relatively short (∼270 s).

6. CONCLUSIONS This paper develops a mixed-integer linear programming (MILP) model for optimizing the delivery/injection plans of intermediate stations along a multiproduct pipeline. The model seeks to minimize the magnitude of pump rate variations in each pipeline segment along the pipeline and is solved using the Simulated Annealing (SA) algorithm and heuristic rules. The solution procedure incorporates the following features: (1) A novel algorithm called the space recursive method is used to make an initial solution. (2) To overcome the shortcoming of the SA algorithm and improve the quality of solutions, a heuristic rule is introduced into the SA algorithm. The rule implements the proper relay of any two delivery/injection operations performed by different stations. The proposed model can generate improved batch plans of a multiproduct pipeline in a short time, as demonstrated using data from two actual multiproduct pipelines. The heuristic rule, which is derived from the experience of dispatchers, makes the computational results more suitable for real-world applications than purely hypothetical algorithms. The model and algorithms proposed in this paper can be incorporated into optimization software that dramatically reduces the workloads of dispatchers. In the proposed model, the input plan of the initial station of a multiproduct pipeline is a known condition. However, in the real world, the input plan of the initial station of a multiproduct pipeline is difficult to draft, because of the storage limitation of the tank farm that is linked to the initial station. Therefore, the problem of how to optimize the input plan of the initial station of a multiproduct pipeline deserves further research and will be addressed in the future.

Figure 21. Influences of the cooling factor r on the final objective function value and the computational time.

Figure 22. Variation in computational times as a function of number of stations (Jmax) and number of batches.

algorithm can provide solutions of excellent quality in efficient computation times. 5.3. The Influence of the Cooling Factor. In Cases 1 and 2 described in Sections 5.1 and 5.2, the value of the cooling factor r was set to be 0.8. Next, the method for setting the value of r in the improved SA algorithm will be discussed using Case 2 as an example. The cooling factor r not only affects the quality of solutions, but also affects computational times. Thus, choosing an appropriate value of r is crucial. Different values of r from 0.05 to 0.95 are used to solve Case 2. The influences of the cooling factor r on the converging processes, the final objective function values, and the computational times are analyzed, as shown in Figures 18−21, respectively. These influences are briefly discussed as follows. (1) When r increases from 0.05 to 0.35, the objective function value is continuing an obvious downward trend, even though the iterations have been terminated. (2) When r ranges from 0.4 to 0.7, the computational processes are almost convergent. However, as seen in Figure 21, the final objective function value can be improved further. (3) When r increases from 0.75 to 0.95, the final objective function value barely changes. (4) With the increase of r, the computational time exponentially grows. As seen in Figure 19, the computational time for r = 0.95 is approximately twice the computational time for r = 0.9.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Changchun Wu: 0000-0002-6735-006X Notes

The authors declare no competing financial interest.

■

NOMENCLATURE

Indexes and Sets

j ∈ J = set of stations, which ranges from the initial station (j = 1) to the terminal station (j = Jmax) p ∈ P = set of products, p ≤ Pmax i ∈ I = set of batches, I = Iold ∪ Inew, i ≤ Imax, where Iold is the set of old batches in the pipeline at the beginning of the time horizon and Inew is the set of new batches to be injected during the scheduling horizon

5105


Article

Industrial & Engineering Chemistry Research k ∈ K = set of time nodes, k ≤ kmax Parameters

Fj = volumetric coordinate of station j cpj = weight factor that is used to avoid the pump rate fluctuations being concentrated in one pipeline segment bj = type of station j; the value assigned is −1 if station j is a delivery station or 1 if the station j is an injection station QImin j,p = minimum flow rate at which intermediate station j delivers/injects product p QImax j,p = maximum flow rate at which intermediate station j delivers/injects product p VImin j,p = minimum volume at which intermediate station j delivers/injects product p in every single operation VImax j,p = maximum volume at which intermediate station j delivers/injects product p in every single operation QSmin = minimum pump rate of pipeline segment j between j station j and station j + 1 QSmax = maximum pump rate of pipeline segment j between j station j and station j + 1 demand VPj,p = delivery/injection demand for product p at intermediate station j tk = time required for time node k Qini,k = flow rate at which that the initial station injects batch i in time interval (tk, tk+1) Vini,k = volume at which that the initial station injects batch i in time interval (tk, tk+1) WB0i = volume of old batch i in the pipeline at the start moment of the scheduling horizon yi,p = binary parameter that is assigned a value of 1 if batch i contains product p Δt = discrete time step τ = time interval length that a delivery/injection operation needs to synchronize with another delivery/injection operation carried out by another station Tf = temperature parameter of the SA algorithm (Tf)max = initial temperature of the SA algorithm (Tf)min = iteration terminating temperature of the SA algorithm r = cooling factor of the SA algorithm num = Markov Chain’s length of the SA algorithm pro = acceptance probability for a new solution of the SA algorithm u, u′ = random numbers between 0 and 1

■

dtstation = duration over which station j delivers/injects batch i j,i tstart = start moment at which station j delivers/injects batch i j,i tend j,i = end moment at which station j delivers/injects batch i tthead = moment that batch i reaches station j j,i tttail j,i = moment that batch i leaves station j Wpass j,i = total volume of batch i that passes by the station j during the scheduling horizon

REFERENCES

(1) Rejowski, R.; Pinto, J. M. Scheduling of a multiproduct pipeline system. Comput. Chem. Eng. 2003, 27, 1229−1246. (2) Rejowski, R.; Pinto, J. M. Efficient MILP formulations and valid cuts for multiproduct pipeline scheduling. Comput. Chem. Eng. 2004, 28, 1511−1528. (3) Rejowski, R.; Pinto, J. M. A novel continuous time representation for the scheduling of pipeline systems with pumping yield rate constraints. Comput. Chem. Eng. 2008, 32, 1042−1066. (4) Cafaro, D. C.; Cerdá, J. Optimal scheduling of multiproduct pipeline systems using a non-discrete MILP formulation. Comput. Chem. Eng. 2004, 28, 2053−2068. (5) Cafaro, V. G.; Cafaro, D. C.; Mendez, C. A.; Cerdá, J. Detailed scheduling of operations in single-source refined products pipelines. Ind. Eng. Chem. Res. 2011, 50, 6240−6259. (6) Cafaro, V. G.; Cafaro, D. C.; Méndez, C. A.; Cerdá, J. Detailed Scheduling of Single-Source Pipelines with Simultaneous Deliveries to Multiple Offtake Stations. Ind. Eng. Chem. Res. 2012, 51, 6145−6165. (7) Cafaro, V. G.; Cafaro, D. C.; Méndez, C. A.; Cerdá, J. MINLP model for the detailed scheduling of refined products pipelines with flow rate dependent pumping costs. Comput. Chem. Eng. 2015, 72, 210−221. (8) Cafaro, D. C.; Cerdá, J. Operational scheduling of refined products pipeline networks with simultaneous batch injections. Comput. Chem. Eng. 2010, 34, 1687−1704. (9) Cafaro, D. C.; Cerdá, J. Efficient tool for the scheduling of multiproduct pipelines and terminal operations. Ind. Eng. Chem. Res. 2008, 47, 9941−9956. (10) Relvas, S.; Matos, H. A.; Barbosa-Povoa, A. P. F. D.; Fialho, J.; Pinheiro, A. S. Pipeline scheduling and inventory management of a multiproduct distribution oil system. Ind. Eng. Chem. Res. 2006, 45, 7841−7855. (11) Cafaro, D. C.; Cerdá, J. Dynamic scheduling of multiproduct pipelines with multiple delivery due dates. Comput. Chem. Eng. 2008, 32, 728−753. (12) Cafaro, D. C.; Cerdá, J. Optimal Scheduling of Refined Products Pipelines with Multiple Sources. Ind. Eng. Chem. Res. 2009, 48, 6675− 6689. (13) Shah, N. Mathematical programming techniques for crude oil pipeline. Comput. Chem. Eng. 1996, 20, S1227−S1232. (14) Herrán, A.; de la Cruz, J. M.; de Andrés, B. A mathematical model for planning transportation of multiple petroleum products in a multipipeline system. Comput. Chem. Eng. 2010, 34, 401−413. (15) Garcia-Sanchez, A.; Arreche, L. M.; Ortega-Mier, M. Combining simulation and Tabu search for oil-derivatives pipeline Scheduling. Comput. Intell. 2008, 128, 301−325. (16) Techo, R. Product-line computer scheduler. Oil Gas J. 1973, 71, 12−41. (17) Techo, R.; Taylor, R. F.; Armstrong, F. T. Digital computer schedules Colonial operations from startup. Oil Gas J. 1964, 4, 94−100. (18) Sasikumar, M.; Ravi Prakash, P.; Patil, S. M.; Ramani, S.; et al. Pipes, a heuristic search model for pipeline schedule generation. Knowledge-Based Syst. 1997, 10, 169−175. (19) Hane, C. A.; Ratliff, H. D. Sequencing inputs to multi-commodity pipelines. Ann. Oper. Res. 1995, 57, 73−101. (20) Chen, H.; Wu, C.; Zuo, L.; Li, W.; Yu, Y. Applying the Simulated Annealing Algorithm to Optimize the Scheduling of Products Pipelines. Presented at the PSIG Annual Meeting; Pipeline Simulation Interest Group (PSIG): Vancouver, British Columbia, 2016; Paper No. 1605. (21) Ingber, L. Very fast simulated re-annealing. Math. Comput. Modelling 1989, 12, 967−973.

Variables

xj,i = binary variable that is assigned 1 if the station j delivers/ injects batch i during the scheduling horizon oj,i,k = binary variable that is assigned 1 if the station j delivers/ injects batch i in time interval (tk,tk+1) Vj,i,k = volume at which that the station j delivers/injects batch i in time interval (tk,tk+1) Qj,i,k = flow rate at which the station j delivers/injects batch i in time interval (tk,tk+1) Hi,k = upper volumetric coordinate of the batch i at the moment tk WBi,k = volume of batch i in the pipeline at moment tk QDj,k = pump rate of pipeline segment j in time interval (tk,tk+1) QAj = average pump rate of the pipeline segment j during the scheduling horizon f 1j,k, f 2j,k = auxiliary variables qstation = flow rate at which station j delivers/injects batch i j,i vstation = volume at which station j delivers/injects batch i j,i 5106


Optimization of Detailed Schedule for a Multiproduct Pipeline Using a

Recommend Documents