Rolling Horizon Approach for ProductionâDistribution Coordination of

Article pubs.acs.org/IECR

Rolling Horizon Approach for Production−Distribution Coordination of Industrial Gases Supply Chains Miguel Zamarripa,† Pablo A. Marchetti,‡ Ignacio E. Grossmann,*,† Tejinder Singh,§ Irene Lotero,§ Ajit Gopalakrishnan,§ Brian Besancon,§ and Jean André∥ †

Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States INTEC (UNL-CONICET), Güemes 3450, Santa Fe, Argentina § American Air Liquide Inc., Delaware Research and Technology Center, Newark, Delaware 19702, United States ∥ Air Liquide, Claude & Delorme R&D Center, 78350 Jouy-en-Josas, France ‡

ABSTRACT: This paper addresses industrial gases supply chains involving multiple products at multiple plants that must be coordinated with multiple depot-truck-routes in order to satisfy customer demands. The full-space optimization problem corresponds to a large-scale mixed-integer linear programming model (MILP). To solve large-scale industrial problems, this paper proposes a rolling horizon approach with two aggregation strategies for solving the smaller subproblems. The first one relies on the linear programming (LP) relaxation for which the binary variables (complicating variables) of the distribution problem are treated as continuous, while the second one uses a novel tailored model for the distribution side constraints that leads to improved solutions. A real case study of an industrial gases supply chain has been addressed obtaining good results in both objective value and with lower computational effort compared with the full-space solution. The extension to longer time horizons through a receding horizon is also considered.

1. INTRODUCTION A production−distribution network for industrial gases consists of production sites, air separation units with attached cryogenic storage tanks, and depots from which trucks (and trailers) depart to satisfy several customer demands. The common issues in supply chain (SC) planning (production, inventory management, and simplified distribution management) are complicated by the need to allocate multiple detailed distribution tasks (allocate trucks to distribute the product, select the routing problem, consider truck availabilities, etc.). Additionally, when dealing with industrial size problems this complexity greatly increases and multiple challenges must be faced: (i) the combinatorial nature of the problem from multiple site availabilities, trucks to be allocated to each site, possible routes (or customer sets) to deliver the products to the customers that leads to an exponential growth of the search space;1 (ii) symmetric solutions that are equivalent which make the branch and bound search inefficient; (iii) the extension for addressing uncertainty and rescheduling problems becomes computationally very expensive. Previous work by Marchetti et al.2 has shown the importance of the production−distribution coordination in industrial gases supply chains. An effective coordination provides optimal solutions with lower total cost. Nevertheless, the computational expense becomes a major limitation when dealing with industrial size problems. Decomposition techniques are mainly applied for solving large scale problems, such as MILP problems with a large number of discrete decision (binary variables); with complicating variables and constraints. Decomposition techniques have been proposed depending on the complexity and characteristics of the problem, such as Lagrangean decomposition,3,4 Benders decomposition,5 bilevel decomposition,6 and rolling horizon.7 © XXXX American Chemical Society

The rolling horizon (RH) approach has been shown to be an effective decomposition technique for planning and scheduling problems.8,9 The RH method decomposes sequentially the problem by dividing the time horizon into a “detailed time block” and an “aggregate time block”, and solving a set of subproblems by fixing the discrete variables of previous time periods. The “detailed time block” uses the representation of the original model, while a simplified model is provided for the “aggregate time block”.8 In Process Systems Engineering, RH has been applied to supply chain optimization and production planning and scheduling, including considerations of uncertainty. Examples include strategic planning,10 capacity expansion problems,11,12 scheduling of multiproduct batch plants,13 inventory management of multiproduct processes,14 planning and scheduling,15,16 scheduling of multiproduct batch plants under demand uncertainty,17 and electricity procurement under uncertainty.18 Algorithmic improvements of the RH method include, as an example, a heuristic procedure to improve the solution quality when applied to integrated planning and scheduling optimization problems.16 The heuristic method is based on the development of a production capacity model through a parametric programming approximation. Another interesting approach to improve the RH methodology is presented by Kostin et al.,10 who developed a heuristic method where some of the integer variables are considered as continuous in the aggregate model. The rolling horizon approach has also been found to be a very Received: January 19, 2016 Revised: January 30, 2016 Accepted: February 3, 2016

A

DOI: 10.1021/acs.iecr.6b00271 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Industrial gases supply chain example.

(iv) the energy consumption per unit of product usppmi for each plant p producing in mode m product i, (v) the maximum and minimum production rate limits max (W min pmi , W pmi ), (vi) the maximum storage capacity Lmax for product i at pi plant p, and the initial inventory Lini pi of each product, at each plant p, (vii) a set of customers c ∈ Ci for each product type i, (viii) a set of depots d ∈ D, and the trucks−depot−product association k ∈ Kdi, (ix) trucks are owned by the company that runs the supply chain, provided are the maximum capacity Uktruck, and the travel cost per unit distance ck of each truck k, (x) the set of time periods t ∈ T, each one with duration Δt , where typically two time periods per day are considered (half day peak and off-peak time periods), (xi) the electricity price forecast upt for each plant p at each time t ∈ T, (xii) the demand forecasts are considered as planned delivdeliv eries Pdc,t1,t2 of the product to deliver at customer c for each time period [t1, t2]. The main goal is to minimize production and distribution costs to satisfy customer demands with known planned deliveries, production capacity, electricity prices, and distances. Three main decisions have to be made at each time period: (1) production (levels and modes) for each plant; (2) inventory levels at each storage tank for every plant; (3) detailed distribution planning, that is, how much product to deliver to each customer, and deciding which route and which trucks to use. An industrial gases supply chain coordinates production/ distribution decisions at multiple plants/depots, while satisfying several customer demands (see Figure 1). Plants operate air separation units, which feature high electricity consumption, producing gaseous and liquid products. Distribution is carried out by pipeline for gaseous products and by truck for liquid products. For truck distribution, depot-truck and truck-product assignments are considered. Moreover, trucks at a given depot can deliver product from multiple sources. To ensure customer storage replenishments, products can be purchased from alternative sources. As the number of sources and customers increases, the selection of the possible routes to be included in the formulation becomes a critical issue. Production decisions include production modes and rates at each plant, with varying efficiencies and power consumptions. To handle electricity price fluctuations during the day, a peak/ off-peak pricing scheme is considered. The main distribution decisions include the assignment of vehicles to routes at each time period, the truck load (availability and truck capacity). The distribution capacity depends on the inventory levels of liquid products at each plant, truck availability and maximum product

efficient solution method for scheduling problems with long time horizons.19,20 To our knowledge the RH approach has not been applied before to production with detailed distribution problems, which is the specific application of this work. This work describes a rolling horizon decomposition strategy presenting two approaches for the aggregate time block. The first one relies on a relaxation of the integrality condition of a subset of the binary variables. The second one consists of a novel model for aggregating the distribution by considering direct shipment with unlimited truck availability to satisfy customer demands. While both approaches provide rigorous lower bounds to the distribution cost, the latter provides tighter lower bounds. The application of the proposed approaches reduce the computational requirements by up to a factor of 6 times, while small deviations from optimality are obtained when compared to full-space solutions, which cannot be solved to optimality. The solutions obtained provide higher coordination among plants/depots to satisfy a common set of shared customer demands. This paper is organized as follows. The problem statement is given in section 2, followed by the description of the mathematical models proposed in this paper and an overview of the rolling horizon approach in section 3. The proposed industrial case study is described in section 4. Two case studies have been solved and the results are shown in section 5. Section 6 describes a receding horizon application to face multiple campaigns using an embedded rolling horizon approach. Finally, the main conclusions are presented in section 7. Appendix A and Appendix B, present a summary of the original PDC model and a proof that the simplified model proposed in section 3.2.2 provides a valid underestimation of the distribution cost, respectively.

2. PROBLEM STATEMENT Given a set of plants, product, depots, customers, time periods, and trucks, the main goal is to determine the optimal production, inventory management, and distribution that minimizes the overall supply chain costs that include production, distribution, shutdown, and subcontracting costs. A detailed production−distribution model has been presented in Marchetti et al.,2, and a summary is presented below. The specific problem of production−distribution coordination (PDC) of industrial gases supply chains can be stated as follows. Given are the following items: (i) a set of production plants p ∈ P, (ii) a set of production modes m ∈ Mp in which plant p can operate, (iii) a set of liquid products i ∈ Ip,m, to be produced when plant p operates in mode m, B


Article


approach (MILP). A summary of the model used for this study is described in Appendix A. The resulting production− distribution model (referred in the future as original PDC model) can become very large. To reduce the computational burden a rolling horizon approach is proposed in this paper. We propose two aggregate models for the distribution side constrains of the problem (LP relaxation and a novel simplified distribution model), while the production side is optimized in detail for all the iterations in the RH approach. 3.1. Rolling Horizon Overview. The rolling horizon method divides the time horizon into a “detailed time block” and an “aggregate time block” and successively solves a set of subproblems in which discrete decisions are fixed in previous time periods. The “detailed time block” uses the representation of the original PDC model and a simplified model has been developed to be used as the “aggregate time block”.8 Figure 2 shows an illustrative example in which we describe a 3 day problem, with 6 time periods of 12 h each (2 time periods

purchase at alternative sources. The detailed distribution management considers availability of trucks and trailers, truck capacities, distance calculations (depot-plant-customer-depot), truck constraints (truck-product, truck-depot, truck-customer assignments), and possible customer sets/routes used to distribute the product to the customer. A customer set S is generated by a route generation algorithm.2 The route generation algorithm considers (i) plant−customer correspondence, (ii) filling ratio (truck load/total truck capacity), (iii) plant−customer distances, and (iv) maximum number of deliveries by time period. At the supply chain level, production and distribution decisions must be coordinated to fulfill several customer demands, internal plant demands, and guarantee a safety stock in the inventory levels, while minimizing the overall cost. When dealing with detailed distribution models like the one in this paper, the inventory routing problem (IRP) must be addressed. The IRP determines the policy used to characterize the customer demands (i.e., detailed inventory management, direct deliveries or planned deliveries) and the way replenishment decisions are considered (i.e., order-up-to level and zeroinventory ordering). This work considers fixed planned deliveries for each customer. On the basis of demand forecast data, the inventory management problem is simplified for all the customers, assuming planned deliveries (fixed demands) and order-up-to level (full replenishment of the customer tanks). The inventory management problem is then simplified by considering a fixed amount of product to be delivered to the customer in a time window (for each day, specifically between two time periods). The main model assumptions are as follows: (i) two time periods per day (peak/peak off) are considered, and as previously mentioned planned deliveries occur per day, so the problem considers two time period windows to satisfy a given fixed demand; (ii) each time period represents 12 h., given the time limitation that trucks do not visit more than three customers in a single route; (iii) customer demands are fixed and each customer has at least one demand during the week and no more than one per day. This work further elaborates on the mixed-integer linear programming (MILP) formulation by Marchetti et al.2 This model minimizes the overall cost of production and distribution of industrial gases supply chains. While the previous work can obtain feasible/optimal solutions for modest sized problems, the computational effort required may be extremely high for large-scale problems. This is mainly due to the combinatorial nature of the problem at the distribution side, where multiple trucks are used to pick up products and visit multiple locations at each time period. In this paper we propose a rolling horizon (RH) approach to reduce the computational effort and to improve the solution quality of production−distribution coordination for large-scale industrial gases supply chain problems.

Figure 2. Illustrative example (3 days problem).

per day). This means that we have 3 subproblems (or iterations) that are solved sequentially: Subproblem 1: time periods 1 and 2 are considered as the detailed time block, while time periods 3 to 6 are considered as the aggregate time block. Subproblem 2: in this subproblem the solution of time period 1 and 2 is known, so the algorithm uses the solutions of subproblem 1 to fix the discrete decisions for time periods 1 and 2. Then the detailed time block corresponds to time periods 1 to 4 and the aggregate time block involves time periods 5 and 6. Note that the size of the detailed time block grows, while the size of the aggregate model is reduced. Subproblem 3: fixes the binary variables for time period 1 to 4 and solves the full-space problem (original PDC model from time periods 1 to 6 with fixed information). It is important to note that once all the iterations are performed, the decisions based on the original PDC model replace all the decisions made at the aggregate time block. The detailed time block model corresponds to the original PDC model with original variables and detailed production (eq A1−A9) and distribution (eq A10−A22) constraints. We refer the detailed time block to the original production− distribution coordinated (PDC) model, which is given in Appendix A by eqs A1−A22 and objective function A23−A25. 3.2. Aggregate Time Block Models. 3.2.1. Relaxed Model. This work considers two different aggregate models. The first one corresponds to the simplest aggregate model, which is represented by the LP relaxation where a subset of

3. MATHEMATICAL MODEL In industrial gases supply chain problems, the production distribution coordination is usually solved in a sequential approach.21,22 The production side is usually optimized through traditional mathematical programming techniques, while the distribution side problem usually solves a vehicle routing problem (VRP) using either mathematical programming23 or heuristics.24 Marchetti et al.2 demonstrate that the overall total cost can be improved from the production−distribution coordination using a simultaneous mathematical programming C


Article

Industrial & Engineering Chemistry Research TNp , c , t = PLp , c , t /max[Uktruck ]

the binary variables are treated as continuous, namely, Ykpt (truck-plant selection), and ykst (route-customer selection). We refer to this model as the Relaxed Model (RM) and it corresponds to the same original PDC model given by eqs A1−A20 with objective function A23−A25. The second one corresponds to a simplification of the distribution side constraints, and it is presented in the next section. 3.2.2. Simplified Distribution Model (SD Model). We propose a novel simplified distribution model (SD model) providing a valid underestimation of the distribution cost of the full-space model. Furthermore, it obtains a tighter underestimation compared with the solution of the LP relaxation of the full-space problem. The SD model uses the original production side model and simplifies the distribution side constraints by disregarding the detailed distribution problem previously defined. The new model disregards the trucks k and customer sets s, considering direct shipment from the plant to the customers. PLpct determines the amount of product delivered from plant p to customer c at each time period t. Therefore, the model estimates the product to be delivered by each plant by a tailored model without the combinatorial complexity of the truck and route assignment problem. It is important to emphasize that the proposed simplified distribution model provides a valid underestimation to the total distribution cost of the full-space model, which is demonstrated by a rigorous proof in Appendix B. New Material Balances. Equation 1 keeps track of the inventory level Lpi,t at plant p of product i at time period t. The inventory level at the time period t is equal to the product stored in a previous time period, the summation of the production over the time period t, minus the total product distributed to the network (amount of product supplied on-site Dsite p,i,t and the total product distributed in an aggregated manner PLp,c,t). Lpi , t = Lpi , t − 1 + Δt

∑ m ∈ M pt , i

Wpmi , t − Dpsite ,i,t −

∑

(3)

Cdt =

c ∈ Jc

∀t∈T (4)

c

t2

∑ ( ∑ PLpct ) = Pdcdeliv ,t ,t

1 2

t = t1

p ∈ Ip

∀ c ∈ Ci , i ∈ I , {t1 , t 2} ⊂ T : Pdcdeliv , t1, t 2 > 0

(5)

The new objective function considers the fixed cost (start-up), no run cost (plants are in a shut-down mode at the end of the time horizon), aggregate distribution, and subcontracting costs. Notice that the production cost remains unchanged as in eq A20, while the distribution changes are as in eq 7:

∑

APCost t =

p∈P

+

(F ptstart ·bptstart + PWpt ·Δt ·upt )

own

∑ p∈P

p∈P

(F ptstart ·(1 − Bpmt ))

own

∑

ADCost t =

(6)

Cdt +

own

(C pdist , c · PLp , c , t )

∑ ∑ c

alt

(7)

p∈P ,

minimize ∑ (AP Cost t + ADCost t )

(8)

t∈T

Finally, the simplified distribution model (SD) is given by the original production constraints (eq A1−A8) and the new aggregate material balances and distribution constraints (1−5), with objective function (eqs 6−8). Table 1 presents a brief Table 1. Decision Making Comparison

PLp , c , t (1)

Additionally, as in the original model, the purchase at the alternative sources is constrained. Equation 2 considers the summation of product purchased at the alternative sources and then distributed to the customers to be lower than the purchase capacity, which is also constrained by the capacity of the trucks. Once the truck and route assignment problem has been simplified, we need to estimate the distribution cost (from the transportation of products directly from plant to customer). TNpct estimates the number of trips made to satisfy the demand from plant p to customer c. The estimated number of trucks is determined in eq 3 by considering the product distributed over a maximum truck capacity. Since the new simplified distribution model considers direct shipments from the plant to the customers, a new distribution cost Cdt must be determined (eq 4). The first term is the total distance covered by the shipments from each plant: determined as 2 times the number of trips from plant p to customer c, and then multiplied by the distance between plant and customer. Once we obtain the total distance covered from each plant, the cost is given by all the shipments from each plant p by the distribution cost (dcp). The new demand satisfaction is given by eq 30, in which the time window is also considered.

∑ PLp,c ,t ≤ Q ppurchase ,i,t

∑ ((∑ (2·TNp,c ,t ·UPc ,p))·dcp) p

c ∈ UIu , i

∀ i ∈ Ip , p ∈ P , t ∈ T

∀ c ∈ d, p ∈ P, t ∈ T

k∈K

decision

original PDC model

RM model

SD model

production scheduling plant inventory estimate truck assignment problem routing problem truck load subcontracting

yes yes yes yes yes yes

yes yes yes yes yes yes

yes yes no no yes yes

comparison of the decisions to be taken by the models presented above. As we described above, the simplified distribution model represents an approximation of the distribution side of the problem, with enough detail to provide an underestimation of the production load at the plants, and then replaced by the original decisions in the rolling horizon approach. Appendix B provides a proof about the validity of the underestimation. Given the detailed time block and the aggregate time block, the proposed rolling horizon algorithm (for the SD model) can be described as follows: Rolling Horizon Algorithm Input time period set T,r, 1. for σ = 1...T 2. use eq A1−A24, ∀ t ≤ σ·r 3. use eq A1−A8 and eq 1−7, ∀ T ≥ t ≥ σ·r 4. solve the model to minimize: ⎡ Tcost = ⎢ ∑ (PCost t + DCost t ) + ⎢⎣ t ≤ σ ·r

∀ i ∈ Ji , p ∈ P alt , t ∈ T (2) D

⎤ (APCost t + ADCost t )⎥ ⎥⎦ t ∈ (T ≥ t ≥ σ ·r )

∑


Article


Figure 3. x and y coordinates: p-plants (depot), As-alternative sources and customers (LIN and LOX).

5. f ix Yk, p,t; yk, s,t, ∀ t ≤ σ·r 6. end for where r represents the number of time periods considered by each subproblem (iteration). In a general application of the RH approach r is equal to 1. For the specific application presented for this paper r = 2, since the model considers 2 time periods per day (peak and off peak electricity prices), the customer demands can be satisfied in a time window (2 time periods), and the product can be distributed any time during this time window (i.e., 6 time periods, r = 2, the RH solves three optimization problems).

than 10 000 possible routes. However, using a route generation algorithm2 the number of customer sets (routes) is reduced to 1882, in which the maximum distance, filling ratio, and maximum distribution time are the main setting parameters used to determine the best set of routes (to be used by the model). To justify the number of routes used in this paper, we performed a sensitivity analysis for the number of routes generated by the algorithm (changing the settings of the algorithm). Different sets (from 1300 to 2200 routes) of routes were tested using the original PDC model. The proposed customer set includes 1882 routes and exhibits the best convergence properties (best solution and CPU time).

4. INDUSTRIAL CASE STUDY To illustrate the capabilities of the proposed rolling horizon method, a case study based on an industrial size test case has been considered. Given is a fixed network superstructure of 4 plants-depots, 15 alternative sources (with fixed locations and purchase limits), 46 trucks with their corresponding location and truck load availabilities to satisfy 229 customer demands of 2 products (liquefied nitrogen and oxygen, LIN and LOX, respectively). The production plants operate with 2 or 3 different production modes (depending on the plant and product specification). The main parameters considered are the peak and off-peak electricity prices, minimum and maximum inventories, plant-depot-customer-alternative source distances, and planned deliveries for customers and the geographical locations of plants, depots, and customers, allowing to compute the route distances. The main goal is to provide an optimal production and distribution plan considering production, storage, and distribution capabilities of the given network. The problem is solved for a time horizon of 7 days (divided into two time periods per day and corresponding to a time horizon H = 14 time periods). Figure 3 shows the two-dimensional network in terms of the plants, alternative sources, and customer locations. If we consider all the possible combinations of plants-customers links, the problem under consideration naturally involves more

5. RESULTS All the models are implemented in GAMS 24.1.3 using CPLEX 12.5.1 as the solver. Solutions are obtained on an Intel Core i7-870 (2.93 GHz, 4 cores) machine with 12 GB of RAM. All instances were solved using parallel processing with eight threads. Given the characteristics of the problem under consideration, we tested the two aggregated models RM and SD described in section 3. The first model RM consists of the relaxation of the integrality constraint of the set of binary variables used in the distribution side constraints. The second model consists of using the simplified distribution model (SD model) instead of the original distribution side constraints. Note that both aggregating strategies (RM model and SD model) only affect the distribution side constraints. Therefore, the production side constraints in the aggregate time block use the original production model for all the iterations. 5.1. Illustrative Example. A reduced version of the industrial case study is first considered, which corresponds to a medium size industrial problem. With a smaller number of customers and alternative sources. This reduced problem (see Figure 4) corresponds to solving the production−distribution coordination of 4 plants-depots (Figure 4a), 11 alternative sources (As, Figure 4a) to satisfy 41 customer demands (Figure 4b) of two products LIN (liquefied nitrogen) and LOX (liquefied E


Article


Figure 4. x and y coordinates: p-plants (depot), As-alternative sources and customers (LIN and LOX).

number of routes also changes at each time period. Given the structure of the problem and the characterization of the demand, the size of the model does not decrease monotonically after each iteration (as it should be expected). The most important parameter is the demand and it has been characterized based on customer needs. Since demands periodicity and load vary among clients, the size of the model is different in each iteration. For example let us analyze three types of customers: customers type A have 3 demands per week with full truck-load capacity needs; customers type B have 1 demand per week with 20% of truck-load capacity needs; and customers type C have 5 demands per week with 50% or more of truck-load capacity needs. Hence, a different number of routes are assigned to each customer in the route set. Therefore, a different number of routes activate a different number of equations/variables when solving the problem. 5.2. Results for the Industrial Case Study. As described in previous references,8,10 all iterations of the RH algorithm correspond to a relaxation to the full-space problem, and the solutions obtained in the first subproblem provide a valid lower bound of the full-space problem. Figure 6, which is for the large industrial example, shows the optimal objective in both approaches (RM model; SD model), as it can be observed the objective value increases from one iteration to the next. The original PDC model corresponds to the exact approximation of the problem while the first iteration of the aggregated models SD and RM yield lower bounds (original PDC model =100; SD model iter 1 = 97.2; RM model iter 1 = 54; and RMIP solution =50). The original PDC model objective value corresponds to the best solution after more than 12 h and is used as a reference to compare the RM and SD solutions. Although both approaches provide valid lower bounds for the first subproblem, the SD model provides a tighter relaxation of the full space problem. However, as it can be observed in Figure 6, at iterations 5 and 6 the SD model is not solved to optimality. Notice that the original PDC model has been included as a reference, but this value corresponds to the best solution of the original PDC model after 41000 CPU seconds.

Figure 5. Cost analysis (medium size example), subcontracting cost is equal to 0 in all models.

oxygen). A reduced number of possible routes (288) are proposed using the route generation algorithm. Figure 5 shows the statistics and the summary of costs for the small case study. It can be observed that the solution of the original PDC model and the SD model are exactly the same, whereas the solution of the relaxed model (RM) has a higher total cost. This demonstrates the capability of the RH approach to provide similar quality solutions to the full-space of the original PDC model with lower computational effort. In summary, this medium size example clarifies how the proposed rolling horizon approach works. From Table 2 we observe that the RH underestimates the feasible region reducing the number of discrete decisions, but significantly increasing the number continuous variables and constraints corresponding to the simplified distribution SD model described in section 3.2.2 (column “original PDC model” vs iter 1 rolling horizon). In this line, as mentioned above the last iteration considers the original problem (full-space problem) with fixed information allowing the problem to reduce the search space (298 nodes). The results show that the decisions fixed during the algorithm are feasible for the full-space problem (original PDC model or detailed time block). In general, the size of the model changes as seen in Table 2 because the F


Article

Industrial & Engineering Chemistry Research Table 2. Summary of Results for the Medium Size Example Rolling Horizon (SD model/RM model) original PDC model

iter 1

iter 2

iter 3

iter 4

iter 5

iter 6

iter 7

equations

10113

11436/9119

10396/9242

10383/9420

10591/9611

10097/9771

10200/9940

10114/10114

cont. variables

16270

17451/17681

16244/17776

16559/17357

17014/16898

16321/16811

16083/16562

16270/16270

5266

718

962 1108 2070

1020 2070 3090

592 3090 3682

766 3682 4448

818 4448 5266

summary

5266

718

390 718 1108

normalized cost

100

69.80/30.82

69.83/36.94

71.22/45.20

72.13/56.75

81.24/71.77

81.62/80.73

100/100.78

time (s)

114

4.1/1.2

1.1/0.4

1.5/1.5

0.7/2.9

1.0/1.0

0.2/1.2

0.8/2.4

nodes

10488

6003/0

17088/0

869/0

140/233

1334/146

56/295

298/1878

optimality gap (%)

0.1

0.1/0.4

0.1/0

0.1/0.8

0.1/0.9

0.1/0.6

0.1/0.3

0.1/0.9

discrete decisions

solved fixed total

total cost and it is solved 6.1 times faster than the original PDC model. global gap = |best estimate (original PDC model) − best integer| max(|best estimate (original PDC model)|, |best integer|)

After reporting the computational efficiency of the proposed rolling horizon approach, we provide a quantitative comparison of the solutions. Figures 7, 8, and 9 compare the solutions of the full-space from the original PDC model versus the rolling horizon approach SD model (power consumption and production rate, total power consumption, and inventory levels, respectively). Figure 7 presents the comparison of the power consumption and production rate for each plant for each of the 14 time periods. Figure 7 also shows the comparison of results for both models: original PDC model (or full-space) versus rolling horizon (SD model). As it can be observed in Figure 7b,c, plants 3 and 4 are the most affected as a result of the coordination. Additionally, plant 1 (Figure 7a) and plant 2 exhibit the same solutions for both models (plant 2 is in shutdown mode in all the time periods). Therefore, when the power consumption changes, the production level at the specific plant also changes (i.e., the results for plant 3 in the full-space problem show an increase in the power consumption at time periods 7 to 9). Consequently the production of product 1 is increased in the same time period, while in the RH using the SD model as the aggregate time block this behavior is presented until time period 7. As we mentioned before, the demand for gaseous product is only considered if a plant is in shut-down mode (only applicable to plants 1 and 4). As is the case for plant 4 in time periods t4−t7 (see Figure 7c), liquid product must be vaporized to meet the gaseous customer pipeline demand. Figure 8 summarizes the power consumption corresponding to each plant: observe that the total power consumption is slightly lower for the rolling horizon approach (∼1%). In most cases the savings obtained by the SD model are due to a lower production result, obtained by redistributing the production load among the plants (mainly plants 3 and 4, see Figure 8). The production cost per unit of volume (PPD: $/scf distributed) is also slightly different. The resulting PPD in the original PDC model (full-space problem) is 0.001185, 0.0, 0.000825, and 0.000515 for plants 1, 2, 3, and 4, respectively. In contrast the PPD obtained by the rolling horizon using the SD

Figure 6. Relaxation of the objective value for a large industrial case study.

Table 3. Statistics (Original PDC Model, RM Model, and SD Model) rolling horizon (last iteration)

model size

cost, $

CPU results

binary variables continuous variables constraints production distribution subcontracting shut down normalized total cost time (s) global gap (%)

original PDC model

RM model

SD model

26 486 69 902 35 233 42.65 54.71 1.45 1.20 100 41 511 4.7%

26 486 69 902 35 233 41.88 56.95 3.21 1.20 103.22 1 297 7.6%

26 486 69 902 35 233 42.37 55.16 1.09 1.20 99.81 6 749 4.51%

(9)

Table 3 summarizes the comparison of the statistics for all the proposed models, both RM and SD model solutions correspond to the last iteration of the RH algorithm. The termination criterion for the original PDC model was 5% optimality gap or 41500 CPUs, while for each iteration in the rolling horizon approach the criteria were less than 1% optimality gap or 1500 CPUs for the relaxed model, and less than 1% optimality gap or 1500 CPUs for the SD model (Table 3). To make a fair comparison of the solutions a global gap for the SD and RM models have been calculated using the lower bound (or best estimate) of the original PDC model (see eq 9; https://support.gams.com/solver:what_is_optca_optcr). Table 3 shows that the rolling horizon approach using the SD model as the aggregate time block yields to the lowest G


Article


Figure 7. Power consumption (PW) and production rate (W).

model is 0.001107, 0.0, 0.000826, and 0.000494 for plants 1, 2, 3 and 4, respectively. Figure 9 presents the inventory levels for each product in each plant and at each time period. The use of this model is critical not only to guarantee the production−distribution satisfaction, but also to propose an efficient decision making tool to provide inventory management at the plants. As it can be observed, the results predicted by the rolling horizon model for plants 1 to 3 are very similar (Figure 9a−c), while there are some differences in plant 4 (Figure 9d). Therefore, the use of the proposed rolling horizon approach improves the solutions and does not lead to infeasibilities, but changes in the quality of the solution can be encountered (differences in the production profiles, inventory levels, and distribution profiles).

Figure 8. Total power consumption. H


Article


It is worth mentioning that the inventory at the end of the time horizon plays a critical role in this kind of problems, and it is an important issue to be considered. The results presented above consider a redline policy as the lower bound for the inventory level (policy 1). An alternative to this operation is to specify the final inventory levels at the plants to be equal to the initial inventory (policy 2). This can be easily incorporated, and the proposed models can be solved under this consideration. The results for this operating policy can be observed in Table 4. Table 4. Summary of Results for the Large Industrial Case Study (Policy 2)

production cost distribution cost shut down cost subcontracting cost normalized total cost ($) CPU time (sec) optimality gap (%)

original PDC model

original PDC model

rolling horizon (SD model)

Policy 1

Policy 2

Policy 2

42.65 54.71 1.2 1.45 100 41 511 4.7

54.5 51.45 0 0.70 106.65 18 000 4.7

48.25 58.05 1.2 3.92 114.42 6 700 1%

As it can be observed, suboptimal solutions are obtained in the distribution side of the model used by the rolling horizon approach. Additionally, the solution of the original PDC model predicts that all the plants are working and lower subcontracting actions are needed, while the SD model proposes a shutdown mode for plant 2 and subcontract products from the alternative sources. For policy 2 the RH approach require more CPU-time when final inventories are specified. Presumably this is because, except for the last subproblem, all others treat the final inventory as part of the approximation given by the relaxation. As presented above, two policies have been solved: policy 1 (inventory lower bound = redline of the tanks) and policy 2 (inventory lower bound = initial inventory). The second alternative is more conservative, while the first alternative is more flexible. Since inventories play a very important role, the next section proposes a receding horizon approach to select the optimal final inventories using information from further time periods.

6. EXTENSION TO RECEDING HORIZON We can exploit the flexibility of the rolling horizon approach as an optimization tool by using this framework under a receding horizon approach (Wang et al.25 define receding horizon optimization as the first layer of model predictive control in which real-time steady state economic optimization of the plant variables is performed). The procedure consists of optimizing the problem by considering information for future time periods, in which this information is updated (similar to a model predictive control strategy). Instead of solving a 7 day scheduling problem, we extend the time horizon to 14 days (28 time periods). This scenario corresponds to use the information on 2 weeks of operation to obtain a solution for the first week, which is very expensive to be solved using the full-space model. Figure 10 illustrates the receding horizon framework. As it is observed in Figure 10a the proposed approach aims to minimize the production−distribution costs of the first week of operation by considering additional information for further

Figure 9. Inventory levels (original PDC model vs rolling horizon SD model).

Additionally, this problem is greatly affected by the symmetry of the solutions, since for given power consumption, production rate and inventory levels, several equivalent solutions are obtained with the same objective value. I


Article


can be observed in Table 5, an important improvement can be found when future operations are included in today’s production plan (14% production cost using the original PDC model into the receding horizon approach, and 18% when using the SD rolling horizon model under the receding horizon framework). The embedded receding horizon with the rolling horizon approach provides flexibility to the receding strategy and the computational issue is not a problem. The proposed receding horizon approach reduces the computational burden (up to ∼25 times faster with respect to the receding horizon using the original PDC model) while improving the objective value (∼15% with respect to the original solution, and ∼2% with respect to the receding horizon using the original PDC model).

time periods (which is useful in order to provide feasible initial inventories for the next week, production plants could be started for the next problem, simulate disruptions, etc.). Furthermore, Figure 10b shows the main steps of the receding horizon

7. CONCLUSIONS This paper has proposed a rolling horizon approach for coordinated production−distribution for industrial gases supply chains. Two aggregate approaches have been investigated: the first approach relies in the LP relaxation of the binary variables in the original PDC model, while in the second approach a novel simplified distribution model has been developed, in which the detailed distribution is simplified but provides tighter lower bounds. The proposed rolling horizon approaches are applied to small and large-scale industrial problems. The solutions of the illustrative example show that the RH approach finds similar solutions as when the original problem is solved to optimality (less than 1% for the large scale problem). Furthermore, when the full-space problem cannot be solved to optimality, which is the case for the large scale industrial problem, the results show that the rolling horizon using the special aggregate model SD, yields the best performance compared with both the original PDC model (full-space problem) and the RM model. There is a clear trade-off between optimality and computational time, which can be translated to the trade-off between computational effort and coordination level, since, in some cases a small degree of suboptimality can be obtained with less CPU time. Additionally, the RH approach was successfully applied within a receding horizon approach in order to select the optimal final inventory levels, tackling the production distribution coordination using information from further time periods.

Figure 10. Receding horizon framework for 14 days.

methodology in which each iteration considers information on the next day, and at the end of the iterations we provide feasible solutions for the desired time horizon problem (week 1 in this case). The receding horizon has been tested under both the simultaneous model (original PDC model, see ReH_PDC column in Table 5) and the rolling horizon approach using the SD model Table 5. Receding Horizon (Summary of Results for Week 1) original PDC model (week 1) Re H_PDC Re H_RHSD model size binary variables continuous variables constraints costs production distribution shutdown subcontracting total cost normalized CPU time (s) results optimality gap (%)

20 625 38 090 20 625 54.43 18.89 26.68 0 100 114 0.1%

40.08 19.74 26.68 0 86.50 153 0.1%

36.77 20.96 26.68 0 84.41 6 0.1%

■

APPENDIX A. SIMULTANEOUS PRODUCTION−DISTRIBUTION MODEL In industrial gases supply chain problems, the production distribution coordination is usually solved in a sequential approach.21,22 The production side is usually optimized through traditional mathematical programming techniques, while the distribution side problem usually solves a vehicle routing problem (VRP) using either mathematical programming (You et al.23) or heuristics (Campbell et al.24). Marchetti et al.2 demonstrate that the overall total cost can be improved from the production−distribution coordination using a simultaneous mathematical programming approach. A detailed production− distribution model has been presented in Marchetti et al.,2 and a summary is presented below.

(Re_H_RHSD column in Table 5). First, the receding horizon approach uses the original PDC model in seven subproblems over the 2 weeks. Second, the receding horizon approach uses the rolling horizon (SD model) in seven problems over the 2 weeks, which means it solves 49 subproblems in total. We use the medium size case study presented in section 5.1 with a reduced number of customer demands to be satisfied during the time horizon. The information (electricity prices, inventory limits, production capacities, electricity prices, planned deliveries, production parameters, etc.) has been extended replicating the first week equal to the second week under study (14 days28 time periods). Table 5 shows results for the receding horizon using the original PDC model (Re H_PDC) and of the rolling horizon with the SD model (Re H_RHSD). It should be mentioned that no integer solution is found for the 28 period problem. Hence, the results for the full space problem column are obtained by solving the first 7 days. As it

A.1. Production Side Constraints

The production model optimizes the operating modes and production rates at each plant considering the power consumptions. Production plants are described by simplified linear models relating production rates with power consumptions for J


Article

Industrial & Engineering Chemistry Research each available production mode. While the production of gaseous product is not considered in detail in the model, the gas demand has to be included since it can affect the liquid inventory. In this model, the production of gaseous product takes place when a plant is in shut-down mode and the liquid product from the storage tanks must be vaporized and sent directly to the pipeline. The model assumes that in such scenario the volume of liquid to be gasified and sent by pipeline , is known. Additionally, the following assumptions are Rpipeline pi,t considered: (i) initial operation state of each plant p (binit p = 1, when plant p is running at the beginning of the time horizon); start determines if plant p needs start(ii) the fixed start-up cost Fp,t up at time t; and (iii) the maximum and minimum inventory max level (Q min pit , Q pit ) or tank capacity and redline, respectively for product i at plant p at any given time period t. Mode Selection. Equation A1 is used to enforce the selection of a valid operation mode at each time period with the binary variable Bpmt. Each plant can operate at most in a single mode during time t; otherwise, if Bpmt = 0 for all m then the plant p is in shut-down mode during period t. Note that set P has been divided into two subsets: the owned plants operated by the supply chain and the alternative plants in which additional product can be bought if production is not enough to satisfy the customer demand (P = Pown ∪ Palt).

∑

Bpmt ≤ 1

∑

∀λ ∈ Lim m, m ∈ M pt , p ∈ P own , t ∈ T

PWp , t =

∑ m ∈ M pt

Bpmt ≤

pipeline Dpsite , i , t = R p , i , t (1 −

(A7)

Plant Inventory. Minimum and maximum inventory levels are given by the red-line and the maximum capacity of the storage tanks, respectively. The red-line (Q pimin) provides flexibility to the operation when there is an excess demand from gaseous customers. Note that the red-line may vary over the planning horizon and is based on the gaseous customer demand profile. The maximum storage limits Q max pi are related to the physical limitations of the storage facilities (a constant value) for product i. min Q pit ≤ Lpit ≤ Q pimax

∀i ∈ Ip , p ∈ P own , t ∈ T

(A8)

Material Balances. Equation A9 keeps track of the inventory level Lpi,t of product i at plant p and time period t. Inventory at the time period t is equal to the product stored in previous time period, plus the production over the time period t (where Δt is the length of time period t and Wpmi,t is the production rate of product i for mode m such that Bpmt = 1), minus both the total amount of product supplied on-site (Dsite p,i,t) and the total product distributed by trucks (Dtruck p,i,t ) at time t.

m ∈ M p(t − 1)

(A3)

Production Rate Limits. Equation A4 establishes the production limits for each plant; that is, the production capacity is limited by the maximum and minimum production rate. Wpmi,t represents the production rate of product i at plant p while operating at mode m in time period t.

∑

Lpi , t = Lpi , t − 1 + Δt

truck Wpmi , t − Dpsite , i , t − Dp , i , t

m ∈ M pt , i

∀i ∈ I p , p ∈ P

own

,t∈T

(A9)

Notice that for t = 1 the initial inventory is considered (Lpi,0). Material balances are the link with the distribution side constrains, where the total volume distributed by trucks Dtruck p,i,t is defined.

min max Bpmt W pmi ≤ Wpmi , t ≤ Bpmt W pmi

∀i ∈ Ipm , m ∈ M pt , p ∈ P own , t ∈ T

Bpmt )

∀i ∈ Ip , p ∈ P own , t ∈ T

Bpm(t − 1) + bptstart

∀ p ∈ P own , t ∈ T : t > t0

∑ m ∈ M pt

(A2)

∑

∀p ∈ P own , t ∈ T (A6)

∀ p ∈ P own: (bpinit = 0)

m ∈ M p , t0

(usppmi ·Wpmi , t )

Gaseous Product Supply during Shut-down Mode. The production of gaseous products is not considered in the model. However, the gas demand has to be included since it can affect liquid inventory: liquid product from the storage tanks should be vaporized when gas cannot be supplied directly from the separation column if the plant is in shut-down mode. Hence, eq A7 ensures that liquid inventory will be used to satisfy the gas demand for all over-the-fence customers when the plant is in shut-down.

Shut-Down Mode. Equations A2 and A3 define the transition from shut-down mode to any other mode m. We start as a binary variable that can be treated as define bp,t continuous in the interval [0,1] ; bstart p,t is strictly 1 if plant p is in shut down mode at a previous time period and is turned on when at time period t begins and 0 otherwise because the transition cost (from shut-down to operation) is considered in the objective function. Bp , m , t0 ≤ bpstart , t0

∑ ∑ m ∈ M pt i ∈ Ipm

(A1)

∑

(A5)

The power consumption of plant p at time period t is given by eq A6 considering the energy requirement per unit of product i (unit specific power) when the plant p operates in production mode m.

∀ p ∈ P own , t ∈ T

m ∈ M pt

αpmi , λWpmi , t ≤ Bpmt πpm , λ

i ∈ Ipm

(A4)

Additionally, constraints for the total liquid production are given for any production mode at plant p. Equation A5 uses parameters αpmi,λ and πpm,λ to describe the coefficients and upper bound, respectively, of a linear combination of the production rates of every product i. The subset Limm stands for the limits of the feasible region of production mode m, where each λ is associated with a limiting hyper plane.

A.2. Distribution Side Constraints

Distribution decisions are driven by the amount of product to be delivered to satisfy single customer replenishment demands. Product can be delivered from any source (plant or alternative source available, sets Pown or Palt, respectively); the model decides the truck being used and the customers to visit within a given time period. The trucks can perform a single round trip K


Article

Industrial & Engineering Chemistry Research purchase Dpitruck , t ≤ Q pi , t

(depot-plant-customers-depot) on each time period t. Trucks loading and unloading tasks may occur anytime within each time period (a maximum number of possible visits per trip is assumed). The combinatorial problem arises from the consideration of multiple possible trucks available to distribute the product to any customer inside a possible route (although the customer sets have been reduced by a route generation algorithm, it still considers many routes to explore). Route Selection. Each truck k is assigned to a fixed depot d and product i by defining the subset Kdi. Truck k ∈ Kdi travels from its depot d to a valid source location (i.e., a related plant p) and picks up the product i (considering a truck capacity), then delivers the product through customer set (routes) s (one or more customers previously defined), and finally the truck returns to the depot d. Equation A10 ensures that each truck k is only assigned to a single customer set on each time interval t.

∑ ykst

∑

Ykpt =

p ∈ Pdi

∀s ∈ Si , i ∈ I , t ∈ T

c∈S

(A17)

∑ ykst

∑ ( ∑ dsct ) = Pdcdeliv ,t ,t

1 2

t = t1

∀ k ∈ Kdi , i ∈ I , d ∈ D , t ∈ T

DISkt =

∀k ∈ Kdi , p ∈ Pdi , i ∈ I , d ∈ D , t ∈ T

where is defined as the minimum distance required to deliver product i to customer set s from the plant p that is more conveniently located, and βkt represents the additional distance needed if a different plant is selected. βkt is constrained by its lower bound in eq A20 based on the source and customer-set when plant p is selected. The additional distance (δd,p,s) is defined in eq A21 considering the difference between the minimum and the complete distance (dismin ds , disdps, respectively), and the maximum distance δmax dp,i can be then calculated in eq A22. Note that eqs A21 and A22 only include parameters that can be calculated off line.

∀k ∈ Kdi , s ∈ Sdi , i ∈ I , d ∈ D , t ∈ T

Ekpt

∀k ∈ Kdi , i ∈ I , d ∈ D , t ∈ T

p ∈ Pdi

(A14)

∑ δdpsykst

Plant Pick-up Amount. Having defined the volume loaded in truck k at plant p during a given time period, eq A15 defines the total amount of product delivered by truck from plant p at each time t (Dtruck pi,t ). Deliveries are constrained by the depot− product-plant correspondence in p ∈ Pdi. Dpitruck ,t =

∑

∑

Ekpt

≤ δdpmax , i (1 − Ypkt ) + βkt

s ∈ Sdi

∀p ∈ Pdi , k ∈ Kdi , i ∈ I , d ∈ D , t ∈ T δdps = disdps − disdsmin

(A20)

∀s ∈ Sdi , p ∈ Pdi , i ∈ I , d ∈ D (A21)

∀i ∈ I p , p ∈ P , t ∈ T

δdpmax ,i

d ∈ D :(p ∈ Pdi) k ∈ Kdi

(A15)

= max[δdps]

∀p ∈ Pdi , i ∈ I , d ∈ D

(A22)

A.3. Objective Function

Note that eq A15 is valid for all P (P = P ∪ P ), including the owned plants and the alternative sources. Alternative Sources. Equation A16 constrains the maximum amount of product i that can be purchased at any alternative source at time t at the alternative sources (p ∈ Palt) own

(A19)

dsmin ds

Finally, eq A14 ensures that the amount of product picked up at a given plant is the same one being delivered to the selected customers, for each truck k and time period t.

∑ ekst = ∑

+ βkt

∀k ∈ Kdi , i ∈ I , d ∈ D , t ∈ T

(A13)

s ∈ Sdi

∑ disdsminykst s ∈ Sdi

(A12)

ykst Uktruck

(A18)

At this point it is important to remark that the route generation algorithm generates at least one customer set for each customer (can be seen as single delivery), and then it generates more customer sets involving the same customer. Distance Calculations. The distance traveled by truck k in time period t is given by the continuous variable DISkt, which is defined in eq A19.

Truck Load. Continuous positive variables Ekpt and ekst are introduced to handle the quantity of product delivered by truck k, with upper bounds given by eqs A12 and A13, respectively. The variable Ekpt represents the product distributed by truck k from plant p in time period t, while ekst is the amount delivered by k to the customer set s in the same time period t. YkptUktruck

s ∈ Si

∀c ∈ Ci , i ∈ I , {t1 , t 2} ⊂ T : Pdcdeliv , t1, t 2 > 0

(A11)

ekst ≤

∑ dsct

t2

s ∈ Sdi

Ekpt ≤

ekst =

(A10)

Notice that since each truck is associated with a specific depot, there is no need to decide the depot to be used. Since, the summation of ykst represents the number of visits to the customers and Ykpt denotes if the truck k loads product at plant p during time period t, eq A11 establishes that a sourcing plant is required if and only if truck k is delivering product at time t.

∑

∑

d ∈ D :(s ∈ Sdi) k ∈ Kdi

s ∈ Sdi

(A16)

Demand Satisfaction. The product load to be provided for each customer set is calculated in eq A17, which corresponds to summing up all the demands of the customers c visited by the customer set s equal to the product carried/provided by truck k to the customer s in time period t (ekst), which considers the truck used to satisfy the customer demands in the customer set. Therefore, eq A18 quantifies the product shipments needed to satisfy the customer demand forecast (planned deliveries) during a time window. Note that the demand satisfaction (dsct) does not yield the truck−individual customer assignment information. But, the truck and load used to satisfy each individual customer demand can be characterized from the results.

∀k ∈ Kdi , i ∈ I , d ∈ D , t ∈ T

≤1

∀i ∈ Ip , p ∈ P alt , t ∈ T

alt

The goal is to minimize the overall total cost, which involves the production and distribution costs over the planning horizon considered (t ∈ T). Equation A23 considers the start-up and production costs of each plant, and shut-down mode is penalized at the end of the time horizon, while the distribution L


Article

Industrial & Engineering Chemistry Research To this end, we first operate on eq A19 to obtain

cost is given by eq A24 and determined by two terms. The first one corresponds to the total traveled distance multiplied by the transportation cost. The last term considers the subcontracting cost, that is, the cost resulting from the acquisition of product from alternative sources. Finally, the objective function is given by eq A25.

∑

PCost t =

p∈P

DCost t =

∑ (disdsminykst ) + ∑ ( ∑ s ∈ Sdi

=

s ∈ Sdi

∑ ykst ·(disdsmin + ∑ s ∈ Sdi

(F ptstart ·(1 − Bpmt ))

own

d∈D i∈I

∑ p∈P

=

(A23)

∑ ∑(∑ +

≥

(F ptstart ·bptstart + PWpt ·Δt ·upt )

p∈P

∑ (disdsminykst ) + βkt s ∈ Sdi

own

∑

+

DISkt =

= =

truck (Cppurch , i , t · Dp , i , t )

alt

∑ ykst ·(disdsmin( ∑ ∑ ykst ·( ∑

δdpsYkpt )

p ∈ Pdi

disdpsYkpt )

p ∈ Pdi

∑ ( ∑ disdpsykst Ykpt ) p∈P

(A24)

∑

Ykpt ) +

p ∈ Pdi

s ∈ Sdi

k ∈ Kdi

δdpsYkpt )

p ∈ Pdi

s ∈ Sdi

cdisk ·DISkt )

δdpsykst Ykpt )

p ∈ Pdi

s ∈ Spi

∀ k ∈ Kdi ,I ∈ I , d ∈ D , t ∈ T minimize ∑ (P Cost t + DCost t ) t∈T

where βkt ≥ Σs ∈ Sdi(Σp ∈ Pdiδdps ykst Ykpt) can be derived from eq A20 because Ykpt only has binary values. Therefore, using eq B4 and starting from the first term of the RHS of eq A24 we obtain

(A25)

■

APPENDIX B. PROOF THAT THE SIMPLIFIED DISTRIBUTION MODEL PROVIDES A VALID UNDERESTIMATION OF THE DISTRIBUTION COST For each time t, we prove that the SD model (SDM) underestimates the objective function value of the detailed fullspace model (M). To this end, we show that for each integer solution of the detailed distribution model there is a feasible solution of the aggregate model with a lower objective function value. Let us assume that an integer solution of model (M) with variables Ykpt, ykst, ekst dsct, and βkt (where Ykpt and Ykst have binary values) is given. Using this solution, we define the variables PLpct for the aggregate model as follows: PLpct =

∑

(ηpst dsct )

∑ ∑(∑ d∈D i∈I

cdisk ·DISkt ) ≥

k ∈ Kdi

∑ ∑(∑ d∈D i∈I

disdps ·ykst Ykpt )

k ∈ Kdi (B)

  (A) ≥

 ∑ (min(ck)· ∑ [ ∑ (disps· (∑ (ykst Ykpt )) )]) i∈I

k ∈ Ki

p∈P

s ∈ Spi

k ∈ K ps

(B5)

Then, in eq B5 we operate on the term denoted as (A) by using eq A13:

∑

∀p ∈ Pc , t , c ∈ C , t ∈ T

(ykst Ykpt ) ≥

k ∈ K p,s

⎛ e ⎞ kst ⎜⎜ truck Ykpt ⎟⎟ ⎠ k ∈ K p , s ⎝ Uk

∑

(B1)

s ∈ S :(c ∈ s)

(B4)

(C)

1 ≥ max(Uktruck )

where ηpst is the fraction of product delivered to customer set s that has been sourced from plant p:

∑

(ekst Ykpt )

k ∈ K p,s

(B6)

k ∈ Ki

ηpst =

(∑k ∈ K ekst Ykpt ) p,s

(∑k ∈ K ekst )

Using definition B2 and eq A17, the term (C) is equiva lent to

∀p ∈ P , s ∈ Sp , t ∈ T

s

(B2)

∑

In eq B2, Ks = {k | ∃d(k ∈ Kdi ∧ s ∈ Sdi)} is the set of trucks that can be used to visit customer set s, while Kp,s is the subset that can do so from plant p. Notice that from eq B2 it is easy to prove that

∑ ηpst = 1 p ∈ Ps

(ekst Ykpt ) = ηpst · ∑ (ekst ) = ηpst ·∑ (dsct ) c∈s

k ∈ Ks

k ∈ K p,s

(B7)

By replacing the above results in the term denoted as (B) of eq B5 we get

∀s ∈ S , t ∈ T

∑ (disps· ∑

(B3)

s ∈ Spi

We prove that the above definition of PLpct provides, by using definitions A25 and eq 1 to obtain TNpct and Cdt, a feasible solution of model (SD) with a lower distribution cost than the one of the integer solution given for the model (M).

(ykst Ykpt ))

k ∈ K p,s (D)

≥

1 max(Uktruck) k ∈ Ki

M

∑ (disps·ηpst ·∑ (dsct )) s ∈ Spi

c∈s

(B8) DOI: 10.1021/acs.iecr.6b00271 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Limm = limits on the feasible region of production mode m M = production modes Mpt = production modes available at plant p in time period t P = all plants Palt = alternative sources Pown = plants owned by the company Pdi = plants associated with depot d and product i S = groups of customers visited in a given route. Thus, s ⊆ C Si = customer sets for product i Sdi = customer sets available for product i and depot d Spi = alternative customer sets to source product i from plant p T = time periods

where in turn, (D) verifies

∑ (disps·ηpst ·∑ (dsct )) = ∑ (∑ dispsηpst dsct ) s ∈ Spi

c∈s

s ∈ Spi

=

∑( ∑ c ∈ Ci

≥

dispsηpst dsct )

s ∈ S :(c ∈ s)

∑ (2·σpc· ∑ c ∈ Ci

=

c∈s

(ηpst dsct ))

s ∈ S :(c ∈ s)

∑ (2·σpc·PLpct ) c ∈ Ci

(B9)

Finally, if we replace the result of (B9) in eq B5 we obtain

∑ ∑(∑ d∈D i∈I

αpmi = linear combination coefficient for the production rate of product i at plant p running mode m δdps = difference between actual distance disdps and minimum distance dismin ds δmax dp,i = maximum δdps for all possible sets s ∈ Sdi Δt = duration of time period t πmin pm = lower bound of linear combination of production rates for plant p running mode m πmax pm = upper bound of linear combination of production rates for plant p running mode m binit p = whether plant p is running (1) or shut down (0) at time t cdisk = traveling cost per distance unit for truck k Cpurchase = cost of product i if purchased at alternative source pi,t p in time t disdps = shortest traveling distance of route (d, p, s) obtained by application of a TSP method (customers of set s are visited using the shortest path starting at plant p and finishing at depot d) dismin ds = minimum distance required for a truck of depot d to deliver product from any valid source to customer set s Fstart p,t = start-up cost of plant p at time t H = time horizon Lini pi = initial inventory of product i at plant p Pdc,t1,t2 = Demand forecast for customer c in time window starting at t1 ending t2 Q min pi,t = safety stock in time period t for product i at plant p Q max = storage capacity of customer c c Q max pi = storage capacity of product i at plant p Q purchase = maximum volume of product i available at pi,t alternative source p in time t Rpipeline = forecast of gaseous customer pipeline demand for pi,t product i at plant p in time period t Utruck = trailer capacity for vehicle k k Q withdrawal = fixed truck withdrawal volume of product i from pi,t plant p at time period t upt = electricity price forecast of plant p during time period t usppmi = unit specific power r = number of time periods to be solved at each iteration (rolling horizon algorithm) Wmax pmi = maximum production rate of product i at plant p running production mode m Wmin pmi = minimum production rate of product i at plant p running production mode m

cdisk · DISkt )

k ∈ Kdi

⎛ ⎜ ≥ ∑ ⎜min(cdisk)· ∑ k∈ Ki i∈I ⎜ p∈P ⎝ ≥

Parameters

⎛ ⎞⎞ Dist pct ⎟⎟ ⎜ ∑ ⎜2·σpc· ⎟⎟ max(Uktruck ) ⎟⎟ c ∈ Ci ⎜ ⎝ ⎠⎠ k∈ Ki

∑ (min(cdisk)· ∑ ∑ (2·σpc·TNpct )) i∈I

k∈ Ki

p ∈ P c ∈ Ci

(B10)

Since the right hand side of eq B10 is equivalent to eq 1, and eq 4 defines the aggregate distribution cost at time t, the proposed feasible solution of the aggregate model underestimates the objective function value of the full-space formulation yielding a lower bound. Thus, the proof is complete.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +01 + 412-268-3642. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The financial support received from the Center for Advanced Process Decision-making (CAPD) at Carnegie Mellon University and American Air Liquide is greatly appreciated.

■

NOMENCLATURE

Subscripts

c = customer d = depot i = product k = truck m = production mode p = plant s = customer set (subset of customers visited in a given route) Sets

C = customers Ci = customers for product i D = depots I = products Ip = products of plant p Ipm = products produced by plant p while running in mode m Jc = product grades that can be delivered to customer c Ji = product grades of product i K = trucks Kdi = trucks for product i available at depot d

Binary Variables

Bpmt = denotes that plant p operates in mode m during time period t N


Article


planning of supply chains Application to the sugar cane industry of Argentina. Comput. Chem. Eng. 2011, 35, 2540−2563. (11) Ryan, S. M. Forecast frequency in rolling horizon hedging heuristics for capacity expansion. Eur. J. Oper. Res. 1998, 109, 550− 558. (12) Garcia-Ayala, G.; Rios-Mercado, R. Z.; Chacon-Mondragon, O. L. A disjunctive programming model and a Rolling horizon algorithm for optimal multiperiod capacity expansion in a multiproduct batch plant. Comput. Chem. Eng. 2012, 46, 29−38. (13) Beraldi, P.; Ghiani, G.; Grieco, A.; Guerreiro, E. Rolling horizon and fix and relax heuristics for the parallel machine lot-sizing and scheduling problem with sequence-dependent set-up costs. Comput. Oper. Res. 2008, 35, 3644−3656. (14) Marques, C. N.; Matos, H. A.; Relvas, S. Inventory Management for multiproduct tank farm systems using a milp model with Rolling horizon. Comput.-Aided Chem. Eng. 2012, 30, 472−476. (15) Zondervan, E.; Kaland, M.; van Elzakker, M. A. H.; Jan, C. F.; Jan, M. Integrating Planning and Scheduling in an Oil Refinery with a Rolling Horizon Approach. Comput.-Aided Chem. Eng. 2014, 33, 439− 444. (16) Li, Z.; Ierapetritou, M. G. Rolling horizon based planning and scheduling integration with production capacity consideration. Chem. Eng. Sci. 2010, 65, 5887−5900. (17) Cui, J.; Engell, S. Scheduling of Multiproduct Batch Plant under Multiperiod Demand Uncertainties by Means of a Rolling Horizon Strategy. Comput.-Aided Chem. Eng. 2009, 26, 423−428. (18) Beraldi, P.; Violi, A.; Scordino, N.; Sorrentino, N. Short-term electricity procurement: A rolling horizon stochastic programming approach. Appl. Math. Model. 2011, 35, 3980−3990. (19) van den Heever, S. A.; Grossmann, I. E. A strategy for the integration of production planning and reactive scheduling in the optimization of a hydrogen supply network. Comput. Chem. Eng. 2003, 27, 1813−1839. (20) Kopanos, G. M.; Mendez, C. A.; Puigjaner, L. MIP-based decomposition strategies for large-scale scheduling problems in multiproduct multistage batch plants: A benchmark scheduling problem of the pharmaceutical industry. Eur. J. Oper. Res. 2010, 207, 644−655. (21) Boudia, M.; Louly, M. A. O.; Prins, C. A reactive GRASP and path relinking for a combined production-distribution problem. Comput. Oper. Res. 2007, 34, 3402−3419. (22) You, F.; Grossmann, I. E.; Wassick, J. M. Multisite Capacity, Production, and Distribution Planning with Reactor Modification: MILP Model, Bilevel Decomposition Algorithm versus Lagrangean Decomposition Scheme. Ind. Eng. Chem. Res. 2011, 50, 4831−4849. (23) You, F.; Pinto, J. M.; Capon, E.; Grossmann, I. E.; Arora, N.; Megan, L. Optimal Distribution-Inventory Planning of Industrial Gases 1. Fast Computational Strategies for Large-Scale Problems. Ind. Eng. Chem. Res. 2011, 50, 2910−2927. (24) Campbell, A. M.; Clarke, L.; Kleywegt, A. J.; Savelsbergh, M. W. P. The inventory routing problem. In Fleet Management and Logistics; Crainic, T. G., Laporte, G., Ed.; Springer: US, 1998; pp 95−112. (25) Wang, X.; Teichgraeber, H.; Palazoglu, A.; El-Rarra, N. H. An economic receding horizon optimization approach for energy management in the chlor-alkali process with hybrid renewable energy generation. J. Process Control 2014, 24, 1318−1327.

Ykpt = denotes that truck k loads product at plant p in time period t ykst = denotes that truck k visits the customers in set s during time period t Continuous Variables

APCostt = aggregate production cost at time period t ADCostt = aggregate distribution cost bptstart = denotes that plant p starts operation at time period t βkt = additional distance traveled by truck k to load product from a given plant at time t Cdt = estimated total cost at each time period (aggregate model) Dsite pi,t = volume of product i to be gasified and sent by pipeline at plant p in time period t Dtruck pi,t = volume of product i withdrawn for truck delivery from plant p at time period t DCostt = total distribution cost at time t DISkt = distance traveled by truck k at time t dsct = volume of product distributed from the customer set s to each customer c in time period t Ekpt = volume of product withdrawn from plant p and loaded into truck k at time t ekst = volume of product delivered by truck k to the customers s in time t Lpit = inventory of product i available at plant p at the end of time period t PCostt = total production cost at time t PLp,c,t = product distributed directly from plant p to customer c at time t (aggregate model) PWp,t = power consumption of plant p at time period t TNp,c,t = number of trips used to satisfy the demand form plant p to customer c at time period t (aggregate model) Wpmi,t = production rate of product i at plant p, when p is running in mode m at time t (zero otherwise)

■

REFERENCES

(1) Engell, S.; Harjunkoski, I. Optimal operation: Scheduling, advanced control and their integration. Comput. Chem. Eng. 2012, 47, 121−133. (2) Marchetti, P. A.; Gupta, V.; Grossmann, I. E.; Cook, L.; Valton, P. M.; Singh, T.; Li, T.; André, J. Simultaneous Production and Distribution of Industrial Gas Supply-Chains. Comput. Chem. Eng. 2014, 69, 39−58. (3) Geoffrion, A. M. Lagrangean relaxation for integer programming. Math. Program. Stud. 1974, 2, 82−114. (4) Guignard, M.; Kim, S. Lagrangean decomposition: A model yielding stronger Lagrangean bounds. Math. Program. 1987, 39, 215− 228. (5) Benders, J. F. Partitioning procedures for solving mixed-variables programming problems. Numersche Math. 1962, 4, 238−252. (6) Iyer, R. R.; Grossmann, I. E. Synthesis and operational planning of utility systems for multiperiod operation. Comput. Chem. Eng. 1998, 22, 979−993. (7) Bassett, M. H.; Pekny, J. F.; Reklaitis, G. V. Decomposition techniques for the solution of large-scale scheduling problems. AIChE J. 1996, 42, 3373−3384. (8) Dimitriadis, A. D.; Shah, N.; Pantelides, C. C. RTN-based Rolling Horizon Algorithms for Medium Term Scheduling of Multipurpose Plants. Comput. Chem. Eng. 1997, 21, S1061−S1066. (9) Castro, P. M.; Harjunkoski, I.; Grossman, I. E. Rolling-Horizon Algorithm for Scheduling under Time-Dependent Utility Pricing and Availability. Comput.-Aided Chem. Eng. 2010, 28, 1171−1176. (10) Kostin, A. M.; Guillén-Gosálbez, G.; Mele, F. D.; Bagajewicz, M. J.; Jiménez, L. A novel rolling horizon strategy for the strategic O


Rolling Horizon Approach for ProductionâDistribution Coordination of

Recommend Documents

Rolling Horizon Approach for ProductionâDistribution Coordination of