A Progressive Hedging-Based Solution Approach for Integrated

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

A Progressive Hedging-Based Solution Approach for Integrated Planning and Scheduling Problems under Demand Uncertainty Zedong Peng, Yi Zhang, Yiping Feng,* Gang Rong, and Hongye Su State Key Laboratory of Industrial Control Technology, Department of Control Science and Engineering, Zhejiang University, Hangzhou, 310027, China

Downloaded via NOTTINGHAM TRENT UNIV on August 12, 2019 at 11:48:01 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Progressive hedging (PH) is a classical decomposition algorithm for solving multistage stochastic problems. However, due to the exponentially growing model size of real-world enterprise-wide optimization problems, critical issues arise when implementing PH in practice. In this work, we propose a novel PH-based algorithm to address integrated planning and scheduling problems under demand uncertainty in a general mathematical formulation. Strategies are proposed to accelerate and guarantee the convergence of the algorithm. Through application of the enhanced PH to solve variants of a typical state−task network example and a real-world ethylene plant case, computational results demonstrate that the proposed algorithm outperforms directly invoking commercial solvers and gets a better solution within nearly two-thirds of the direct solution time on a serial computer. The advantage of the multistage stochastic programming method is also demonstrated by comparing the model solution with the counterparts of an expected value-based deterministic model and a two-stage stochastic model.

1. INTRODUCTION In process industry, production planning and scheduling are two important parts of enterprise-wide decision making.1 The quality of these decisions directly determines the profit and performance of a company.2 However, there are many factors that influence the quality of decisions in production planning and scheduling. The key factor lies in the fundamental differences between production planning and scheduling. Production planning determines the optimal allocation of resources and production targets in a medium-range time scale at the company’s management level. On the contrary, production scheduling involves more of the plant level and determines assignment of tasks to units and the sequencing of tasks to meet the production target over a short-term horizon.3 Because of these differences, a traditional strategy in the process industry is to first solve the planning problem to determine production targets and then solve the scheduling problem to achieve the goal. The lack of interaction between these two decision-making processes often results in higher operating costs or even infeasible production plans.4,5 Driven by this motivation, many researchers have explored ways to develop a production planning and scheduling integration framework. The most common and precise way to integrate production planning and scheduling is to build a detailed full-space model that spans the entire planning horizons, in which new constraints or variables are created to link planning model and scheduling model.6 However, the high precision of the full-space model generally leads to a large-scale © XXXX American Chemical Society

mixed-integer linear or nonlinear programming model, and various decomposition methods are designed to solve these full-space problems.7,8 Although the advances in decomposition methods reduce the computational burden, the simplification of a full space model also serves as a solution to the full space problem. An optional simplification method is to remove some constraints or aggregate some decisions of the original scheduling model, such as keeping task assignment variables and constraints, and abandoning sequencing variables and constraints.4,9 Moreover, in order to replace individual simulations, surrogate models can also be adopted to economize on the overall computational effort.10 Another compromise between modeling accuracy and computational burden is to use a hybrid modeling method in a rolling horizon frame. When a planning period covers several scheduling periods, a detailed integration model is implemented in the first several periods and can be replaced by a surrogate model in the rest of the periods. In this way, the production target is updated as the completion of production operations in every period and the approximation of the later periods will not significantly affect the quality of the solution.11−13 Received: Revised: Accepted: Published: A

May 12, 2019 July 17, 2019 July 18, 2019 July 18, 2019 DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 1. Integrated planning and scheduling problem under demand uncertainty.

decomposition22 and the alternating direction method of multipliers (ADMM),23 for stochastic programs. The progressive hedging algorithm was initially proposed by Rockafellar and Wets in 1991 as a decomposition strategy for solving large-scale stochastic linear programs.21 The core of the PH algorithm is to relax the nonanticipativity constraints and solve the scenario subproblems independently. Solving scenario subproblems separately is naturally much easier than solving the extensive form directly. However, the difficulty lies in satisfying nonanticipativity constraints for all subproblems. The common situation is that variables of the optimal solution are not equal for the same stage even though the solution can be easily found. To solve this problem, the PH algorithm penalizes deviations from the average values of decision variables. Then scenario subproblems are iteratively solved to enforce nonanticipativity.24 The progressive hedging algorithm has been proven to possess theoretical convergence properties when all decision variables are continuous and even if there is a linear convergence rate given a convex reference scenario optimization model. Besides the successful applications in stochastic linear programming, the PH algorithm has also been proven to be a very effective heuristic algorithm for solving stochastic mixed-integer programs (SMIP). Løkketangen and Woodruff proposed a hybrid algorithm combined PH algorithm and tabu search to solve multistage stochastic mixed-integer programming in 1996.25 Teodor and Xiaorui proposed a two-stage stochastic program to model the stochastic fixed-charge capacitated multicommodity network design problem with demand uncertainty, and used the PH algorithm to solve it.26 Goncalves and Lamghari respectively applied the PH algorithm to solve the operation planning problem of a hydrothermal system and the production scheduling problem in open-pit mines.27,28 Besides the application of the progressive hedging algorithm on different fields, the computational aspects of the PH algorithm in SMIP have also been studied. In 2010, Watson and Woodruff introduced several enhancements of the basic PH algorithm to ensure and accelerate convergence in a mixedinteger case. The enhancements involve an effective penalty factor computation, convergence accelerator, and termination criteria.29 In 2015, Watson and Munoz proposed a scalable decomposition algorithm to solve stochastic transmission and generation planning problems, including a simple scenario reduction framework based on a clustering algorithm.30 In 2016, Gade et al. proposed a method for computing lower bounds in the PH algorithm for multistage stochastic mixed-

Besides the lack of interaction between different decision levels, uncertainties also affect the quality of decisions. The optimal solution of a deterministic model may result in a suboptimal performance in the actual production environment due to uncertainties.14 Stochastic programming and robust optimization are the two main mathematical programming techniques to handle uncertainties. In general, stochastic programming is more appropriate for the planning level and robust optimization is more appropriate for the scheduling level due to the solution procedure difference of these two techniques. Robust optimization uses worst-case analysis to guarantee the feasibility over a specified uncertainty set, which guarantees meeting the production target in all uncertainty scenarios in the scheduling level.15 Stochastic programming optimizes the expectation of the objective function and allows recourse decisions to be made at future time points to adapt to reveal uncertainties.16 We focus on the demand uncertainty in this paper, so stochastic programming is a more appropriate modeling method. To furtherly improve the quality of planning and scheduling decisions, integrated models should be adjusted or improved to lower the impact of uncertainties. Wu and Ierapetritou proposed a hierarchical framework to integrate production planning and scheduling and considered uncertainty utilizing a three-stage stochastic programming formulation.17 Chu and You proposed a bilevel programming method for integrating planning and scheduling problems under production uncertainties. Production uncertainties, such as uncertain task processing times and sequence-dependent transition times, are mainly considered.18 The above literature review reveals the necessity and challenges of integrated planning and scheduling optimization under uncertainty. In this work, we proposed a general formulation to integrate production planning and scheduling under demand uncertainty based on multistage stochastic programming (MSSP). Since MSSP models are generally difficult to solve due to the considerable model scale, decomposition methods for MSSP have been a popular research topic in recent years. Decomposition methods for stochastic programs generally fall into two groups: stage-based methods and scenario-based methods. Stage-based methods, such as the L-shaped method,19 decompose stochastic programs by stages, while scenario-based methods, such as the dual decomposition algorithm20 and progressive hedging algorithm,21 decompose stochastic programs by scenarios. In a sense, the L-shaped method and PH algorithm are the specialization of two more generalized algorithms, the Benders B

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 2. Scenario tree of demand uncertainty.

scenario d represents the possible combination of demand realizations in each planning periods, which corresponds to a path from the root node to the leaf node in the scenario tree. Apparently, the number of scenarios increases with the planning periods, which corresponds to the fact that it is easier for us to predict the demand tomorrow than a month later. A general MSSP formulation (P1) for integrated planning and scheduling problem under demand uncertainty is represented in the following compact form:

integer programs and theoretically proved that the lower bound can be as tight as the one from a dual decomposition method.31 On the basis of the successful implementation and enhancement of the PH algorithm stated above, we choose it to solve the multistage stochastic mixed-integer linear program. In this paper, a progressive hedging algorithm-based solution methodology is proposed. Several accelerating strategies are designed to improve the convergence of the PH algorithm. The outline of this paper is as follows. We present a general multistage stochastic formulation for the integrated planning and scheduling problem under demand uncertainty in section 2. The detailed implementation of the progressive hedging algorithm and several accelerating strategies are presented in section 3. The efficiency of the proposed method is tested through a state−task network (STN) case and an industrial case in section 4. The study is concluded and several future research directions are outlined in section 5.

min

∑ ∑ pd ((cP)T xt ,d + (cS)T yt ,d ) d∈D

t

(1)

Subject to AP xt , d ≤ bP , d

∀ d ∈ D , t = 1, .., T

DP1xt , d = eP1, d

∀ d ∈ D , t = 1, .., T

DP 2xt − 1, d + DP 3xt , d = eP 2, d

2. PROBLEM STATEMENT AND FORMULATION In this paper, we focus on multiperiod integrated planning and scheduling problems under demand uncertainty as shown in Figure 1. The planning horizon is divided into a set of planning periods t = {1, 2, ···, T} with equal time length. In every period of the planning horizon, the scheduling constraints are incorporated into the planning model. The planning decisions and the scheduling decisions are connected by the coupling constraints. At the beginning, we only know the exact demand of the first planning period. The demands of future periods are stochastic, and the only available information is their probability distribution. Therefore, a multistage stochastic program with full recourse to demand uncertainty is proposed, which means every planning period is taken as a stage in the scenario tree and the production decisions are made sequentially upon the realization of demand uncertainty over time. Here-and-now decisions are made in the first planning period and wait-and-see decisions are determined in later planning periods to take recourse actions.32 There are several common approaches to generate a scenario tree, such as bound-based construction, Monte Carlo-based schemes, and the moment-matching principle.33 In this paper, we assume that the discrete probability distribution of demand uncertainty is known and the scenario tree is directly generated from it. An example of a demand uncertainty scenario tree is provided in Figure 2. Each node refers to a certain realization of demand uncertainty. The probability of each node refers to the conditional probability given that all the previous demand uncertainty up to current stage is observed or revealed. A

∀ d ∈ D , t = 2, .., T

(2) (3) (4)

AS yt , d ≤ bS

∀ d ∈ D , t = 1, .., T

(5)

DS1yt , d = eS

∀ d ∈ D , t = 1, .., T

(6)

DP 4 xt , d + DS2yt , d = e x t , d = x t , d′ yt , d = yt , d′

∀ d ∈ D, ∀ d ∈ D,

∀ d ∈ D , t = 1, .., T ∀ d′ ∈ Dt , d , t = 1, .., T ∀ d′ ∈ Dt , d , t = 1, .., T

(7) (8) (9)

Here, D is the set of scenarios and Dt,d is the sets of scenarios that contain the same node as scenario d in stage t in the scenario tree. pd is the probability of scenario d and ∑d pd = 1. xt,d and yt,d, respectively, represent the planning level variables and scheduling level variables in planning period t in scenario d. Coefficient matrixes (AP, DP1, DP2, DP3, DP4) or vectors (bP,d, eP1,d, eP2,d) with subscript P are related to the planning level and the others (AS, bS, DS1, DS2, eS) with subscript S are related to scheduling level. The objective function (1) of the integrated model is often concerned with financial performance indicators, such as minimizing the total cost and maximizing the profit. These indicators often include sales revenue, production cost, inventory cost, backorder cost, raw materials cost, and operation mode switch cost, etc.34 Equations 2−4 represent the constraints in the planning level, which are discrete-time formulated. These constraints mainly include inventory balance, backorder balance, and capacity constraints. Equation 4 describes the spatial and C

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

ρ i min ∑ ∑ jjjcPxt , d + cSyt , d + ωP , t , dxt , d + xt , d − X̅ t , d 2 d t k

temporal relationships on planning variables xt−1,d and xt,d between consecutive time stages. For example, assuming that backorder is not permitted and products are delivered at the end of each planning period, the inventory level of product at the end of planning period t − 1 equals to the sum of the delivery amount and the inventory level of product at the beginning of planning period t. Equations 5 and 6 represent the constraints in the scheduling level. These constraints often include equipment allocation, production capacity limitations, storage capacity limitations, material balance and time limitations. The scheduling model may have various formulations, such as discrete-time formulation, continuous-time formulation and TSP-based formulation.35−37 Both the detailed scheduling model and the aggregated scheduling model could be included in the general integration formulation and will not affect the solution methodology in section 3. Equation 7 represents the logical connections between the planning variables xt,d and scheduling variables yt,d, which incorporates the planning model and scheduling model in the same time scale. A common logical connection is that the planning model determines the production target according to the demand and the scheduling model determines the assignment and sequence of tasks to meet the production target. Equations 8 and 9 are the nonanticipativity constraints, which guarantee the decisions taken at any stage should only depend on information observed up to the respective stage, without the hindsight of future realizations of uncertainties or future actions. For the sake of understandability, the nonanticipativity constraints are explicitly formulated. Although the explicit formulation is more complex than implicit formulation, the presolvers in commercial solvers are able to quickly identify and eliminate most of the redundant variables and constraints in explicit formulation.38,39 With appropriate adjustment, the formulation can also apply to integrated production planning and the scheduling model under other exogenous uncertainties, such as supply uncertainty.

+ ωS , t , dyt , d +

ρ y − Yt̅ , d 2 t ,d

2

2y z

zz {

(10)

Because we only consider demand uncertainty in this model, there are no constraints that link variables across scenarios and no further relaxation is needed here. The relaxed scenario subproblems with penalizations are denoted as (P2). l o i o o min ∑ jjjcPxt , d + cSyt , d + ωP , t , dxt , d o o o o t k o o ρ o 2 o o o + 2 xt , d − X̅ t , d + ωS , t , dyt , d (P2)m o o ρ 2y o o + yt , d − Yt̅ , d zzz o o o 2 { o o o o o o s.t. eqs 2−7 n

(11)

To formally describe the procedure, we introduce several (k) notations. X(k) t,d and Yt,d denote the vectors of planning and scheduling decisions of scenario d in planning period t at (k) iteration k. ω(k) P,t,d and ωS,t,d denote vectors ωP,t,d and ωS,t,d at (k) (k) iteration k. X̅ t,d and Y̅ t,d denote vectors X̅ t,d and Y̅ t,d at iteration k. Then the PH algorithm for the integrated planning and scheduling problem is designed as follows:

3. SOLUTION METHODOLOGY 3.1. Progressive Hedging Algorithm-Based Scheme. In this section, the progressive hedging algorithm is modified to solve the stochastic mixed-integer linear program (P1). The scenario-based decomposition is accomplished by introducing scenario-specific decision variables and relaxing the nonanticipativity restrictions, that is, eqs 8 and 9. Nonanticipativity is then restored iteratively by penalizing deviations of scenario subproblem solutions from the aggregated solution in each stage. Thus, each scenario subproblem in the PH algorithm corresponds to a deterministic version of the stochastic problem, with a penalty-augmented objective function. The objective function with penalty terms is presented as eq 10. The deterministic objective function augments terms that penalize the lack of implementability using a subgradient estimator for the nonanticipativity constraints and a squared proximal term. We introduce a set of multiplier vectors ωP,t,d, ωS,t,d for the subgradient estimator and a penalty parameter ρ for the squared proximal term. X̅ t,d and Y̅ t,d are the probability average of planning and scheduling variables.

First, the algorithm initializes the penalty factor and the multipliers to zero. Then the relaxed individual scenario problems (P2) are solved to obtain the initial value of decision (0) variables X(0) t,d and Yt,d . At the next step, the algorithm D

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (0) calculates a probability-average solution X̅ (0) t,d , Y̅ t,d for each node over adjacent scenarios according to the probability distribution. The scope of the aggregation includes scenarios that contain the node. Then, the penalty factor ρ is set to 0.01 (0) and the multipliers ω(0) P,t,d, ωS,t,d are computed based on the results from step 2 and step 3. At this point, the iteration begins. The values of ωP,t,d, ωS,t,d, X̅ t,d, Y̅ t,d in eq 10 are set to the corresponding values obtained in the last iteration. For each scenario d, scenario subproblems (P2) are solved again to get (k) scenario solutions X(k) t,d and Yt,d in step 6. Then, the weighted (k) (k) (k) average decisions X̅ t,d , Y̅ t,d and multipliers ω(k) P,t,d, ωS,t,d are updated in steps 7 and 8. If the termdiff convergence criterion is satisfied, then the algorithm is terminated. Otherwise, go back to step 5 and continue the iteration. 3.2. Accelerating Strategies. As we have mentioned above, PH algorithm serves as a heuristic algorithm for solving stochastic mixed-integer programs. Besides the basic implementation of PH algorithm, we also adopt some strategies to accelerate the convergence of PH algorithm. 3.2.1. Effective Parameter ρ Computation. In the progressive hedging algorithm, the selection of the penalty factor ρ of quadratic terms plays an important role in the convergence speed. According to Watson and Woodruff,29 the optimal ρ value need not to be fixed at a constant value and variable-specific ρ may be more appropriate for some problems. Three ρ computation methods are proposed by these authors, including fixed method, cost-proportional method, and heuristic method. In the cost-proportional method, ρ for each variable equals to the product of a multiplier λ and the cost parameter in the objective function as eqs 12 and 13. In the heuristic method, ρ is calculated via eqs 14 and 15 for discrete variables and eqs 16 and 17 for continuous variables in each iteration. There is no conclusion to which one of these three methods is the best for different stochastic programs. It is often the case that there might be a great difference among the cost parameters. If the ρ value is not close in magnitude to the cost parameter, the computation of the initial multipliers (ω(0) P,t,d and ω(0) S,t,d in Step 4) will yield a small fraction of xt,d and yt,d. The per-iteration change in the penalty term (ωP,t,dxt,d and ωS,t,dyt,d) will be comparatively small. Slow changes in the penalty terms necessarily yield little movement in xt,d and yt,d, which in turn significantly delays PH convergence. Therefore, the cost-proportional method and the heuristic method are generally better than the fixed method. In this paper, we try all these three methods for the case study and choose the best one.

ρP = λcP

(12)

ρS = λcS

(13)

ρP =

ρS =

ρP =

ρS =

− mind ∈ Dd ,t (X t(0) , d ) + 1)

− mind ∈ Dd ,t (Y t(0) , d ) + 1)

((∑

) )

(0) p |X t(0) , d − X̅ t , d | , 1

d ∈ Dd , t d

(17)

Figure 3. Linearization of penalty terms.

3.2.3. Binary Variable Fixation. As we know, the computational complexity of mixed-integer programming lies in integer variables. In PH iterations, binary variable fixing is an empirically effective heuristic way to accelerate convergence. When the binary variable has satisfied nonanticipativity for a number of iterations, an alternative way is to fix these variables to reduce the scenario subproblem size. Here, satisfying nonanticipativity means the values of the binary variables in the same node all equal 1 or 0. For example, we consider a multistage stochastic program with the scenario tree structured as Figure 2 and a binary variable x in the first stage. When the value of x in the solution of each subproblem equals 1 or 0 for 10 PH iterations, variable x will be fixed. Because we only check the current value and previous value of each variable, there are no coupling relationships in the checking process, which can be implemented simultaneously in PH iterations. 3.2.4. Mipgap Tolerance of Scenario Subproblem. Considering the fact that high accuracy of the scenario

(14)

(15)

cP max

) )

(0) p |Y t(0) , d − Yt̅ , d | , 1

d ∈ Dd , t d

3.2.2. Proximal Penalty Terms Linearization. The PH algorithm involves the nonanticipativity constraints relaxation, which results in quadratic penalty terms in the objective function. Commercial solvers such as CPLEX and Gurobi, can handle problems with quadratic objectives, but the efficiency is often dramatically worse than the linear one according to the experience of solving the two case studies in section 4. To address this issue, we linearize all quadratic penalty terms. As for binary variable penalty terms, the linearization is simple due to only feasible 0−1 values. Supposing that xt,d is a binary variable, the value of ∥xt,d − X̅ t,d∥2 equals to X̅ 2t,d − xt,d. For continuous variables, a piecewise linear function is adopted to fit the quadratic term. Breakpoints are uniformly sampled between the maximum value maxd∈Dd,t(X(k) t,d ) and minimum value mind∈Dd,t(X(k) ) of each variable in the same node. The t,d segments between lower bound and min value, max value, and upper bound are respectively fitted with a linear function. The goal of the above linearization for continuous variables is to increase the fitting accuracy over an interval that may contain an optimal solution and reduce the fitting accuracy over an interval that may not. An example of a linear interpolation for a continuous variable is represented in Figure 3.

cS (maxd ∈ Dd ,t(Y t(0) ,d )

((∑

max

cP (maxd ∈ Dd ,t(X t(0) ,d )

cS

(16) E

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 4. Flowchart of PH algorithm embedded with accelerating strageties.

linearization only exists in PH iterations and will not affect the final solution precision of the partly fixed extensive form. 3.2.6. Warm Starts. For each subproblem (P2), the difference between the subproblem in iteration k and the subproblem in iteration k + 1 only lies in the penalty terms in the objective function. The feasible region remains the same. Therefore, feasibility is first guaranteed when using the previous computed solution as the warm start point. Because X̅ t,d and Y̅ t,d are probability average value and ωP,t,d and ωS,t,d are updated according to Step 4 and Step 8, the objective function in adjacent iterations does not change significantly. Therefore, the previously computed solution in iteration k can be used as an efficient warm start in iteration k + 1.38 There is a trade-off between the optimality and the computational complexity, which is controlled by parameters in binary variable fixation strategy and the early termination strategy. When strict parameters are chosen, it takes more iterations to fix binary variables and a more optimal solution may be obtained. When loose parameters are chosen, the algorithm may terminate quickly with a relatively suboptimal solution. The flowchart of the progressive hedging algorithm imbedded with accelerating strategies is presented in Figure 4.

subproblem solutions is not necessary in early iterations when the current solution is far from the optimal one, it is more desirable to allow each subproblem solution within an acceptable gap with the optimal solution. This will save computational effort and accelerate the convergence without significant loss in the quality of the final solution. In later iterations, the gap can be tightened with a small value of the mipgap parameter to guarantee the optimality of the final solution. 3.2.5. Early Termination of PH Algorithm. The convergence of the PH algorithm is determined by the convergence criteria, such as the termdiff convergence metric. However, for problems with discrete decision variables, the convergence of the last few variables often requires more computational effort than it is worth. A practical way is to solve the significantly smaller extensive form with fixed variables, after a limited number of iterations of PH algorithm based on the binary variable fixing strategy. The criteria of early termination often depend on the proportion of the fixed binary variable to total binary variables. For example, it is designed that a progressive hedging iteration will be early terminated when 90% of the binary variables are fixed in previous iterations. The solution time of the partly fixed extensive form will be massively reduced due to the smaller scale of integer variables. Because the linearization strategy is only adopted in PH iterations to solve subproblems and the partly fixed extensive form contains no quadratic penalty terms, the possible error brought by the approximation of

4. CASE STUDY To test the efficiency of the proposed algorithm, two examples are introduced in this section. The first example is based on a classic STN proposed by Kondili and Pantelides in 1993. The second example is an industrial case of a real-world ethylene F

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 5. State-task network of the case study.

plant in China. We encode the multistage stochastic mixedinteger model and the PH algorithm using PySP module in Pyomo5.3 Python3.6.38,40,41 All scenario subproblems and extensive forms are solved by commercial solvers, CPLEX 12.7.1 and Gurobi 7.0.2. 4.1. Motivating Example. 4.1.1. Case Study Description and Model Formulation. The STN of this example is presented in Figure 5. Three feedstocks are used to produce two different products through five processing stages: heating, reactions 1, 2, and 3, and separation. Process data, storage capacity, and cost data of the STN case are presented in the Supporting Information. In each planning period, production targets are set for the production process depending on the market demands and production capacity. Task assignments, execution sequences, task starting time, and task ending time are computed according to the production target within a scheduling horizon which is equal to the length of the planning period. Backorders are permitted and penalized through the backorder cost. The actual delivery of products is determined jointly by the inventory and demand. The MSSP model is an extension of the work of Li8 and it is given as follows. The objective function 18 minimizes the expectation of the total cost over the set of scenarios and over a set of periods. The total cost is given by the sum of inventory cost, backorder cost, and production cost. ij j min ∑ pd jjj ∑ hsInvst , d + jj d k s ∈ SP , t +

∑ i ,j,n,t

∑

zyz (fci·wvit,,jd, n + vci·bit,,jd, n)zzz zz {

The planning level constraints and scheduling level constraints are connected by eqs 21 and 22. Equation 21 expresses the connection between the production target Pst,d and the inventory of products in each planning period. Equation 22 represents the connection constraints for the initial product inventory for different adjacent periods. stst,,nd= N − stinst , d = Pst , d

stinst , d = Invst − 1, d

∑ wvit,,jd,n ≤ 1

(22)

∀ j ∈ J , n ∈ N , d ∈ D, t

i ∈ Ij

(23)

t ,d t ,d max t ,d vimin , j , n wvi , j , n ≤ bi , j , n ≤ vi , j , n wvi , j , n

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t

(24)

Equations 25−32 express time constraints and eqs 33−34 represent the material balances for each state s expressing that at each event point n the amount of state s is equal to that at event point n − 1, adjusted by the produced and consumed amounts between event points n − 1 and n. Equation 35 represents the storage capacity constraints. More details about the above scheduling constraints can be found in the paper by Ierapetritou and Floudas (1998).42

usUst , d

Tf it, j, d, n = Tsit,,jd, n + αi , jwvit,,jd, n + βi , jbit,,jd, n ∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t (18)

(25)

Tsit,,jd, n + 1 ≥ Tf it, j, d, n − H(1 − wvit,,jd, n)

According to the formulation in section 2, the model contains planning level constraints 19−20, connection constraints 21− 22, scheduling level constraints 23−35 and nonanticipativity constraints. For a particular scenario, eq 19 represents the inventory balance and eq 20 represents the backorder balance.

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t

(26)

Tsit,,jd, n + 1 ≥ Tf it′,,dj , n − H(1 − wvit′,,dj , n) ∀ i , i′ ∈ Ij , j ∈ Ji , n ∈ N , d ∈ D , t

∀ s ∈ SP , d ∈ D , t

(27)

Tsit,,jd, n + 1 ≥ Tf it′,,dj ′ , n − H(1 − wvit′,,dj ′ , n)

(19)

Ust , d = Ust − 1, d + Demst , d − Delst , d

∀ s ∈ SP , d ∈ D , t

(21)

In the scheduling level, eq 23 is the allocation constraint and eq 24 represents the capacity limitations of production units.

s ∈ SP , t

Invst , d = Invst − 1, d + Pst , d − Delst , d

∀ s ∈ SP , d ∈ D , t

∀ i , i′ ∈ Ij , i ≠ i′, j ∈ Ji , j′ ∈ Ji , n ∈ N , d ∈ D , t

∀ s ∈ SP , d ∈ D , t (20)

(28) G

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Tsit,,jd, n + 1 ≥ Tsit,,jd, n

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t

Tf it, j, d, n + 1 ≥ Tf it, j, d, n

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t

Table 1. Expectation of Demand Data Distribution

(29)

(30)

Tsit,,jd, n ≤ H

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t

(31)

Tf it, j, d, n ≤ H

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , t

(32)

stst,,nd = stst,,nd‐ 1 −

∑ φsC,i ∑ bit,,jd,n + ∑ φsP,i ∑ bit,,jd,n− 1, i ∈ Is

j ∈ Ji

i ∈ Is

stst,,nd= 1 = stinst , d −

∑ φsC,i ∑ bit,,jd,n= 1 i ∈ Is

(33)

∀ s ∈ S , d ∈ D, t

j ∈ Ji

(34)

stst,,nd

≤

stsmax

∀ s ∈ S , n ∈ N , d ∈ D, t

(35)

Equations 36−45 are nonanticipativity constraints. The equations represent that the decision variables at the same node in the scenario tree must have the same value. The t,d t,d t,d t,d t,d decision variables consist of Invt,d s , Ps , Ds , Us , wvi,j,n, bi,j,n, t,d t,d t,d , Ts , st , stin . Tf t,d i,j,n i,j,n s,n s Invst , d = Invst , d′

∀ s ∈ S , d ∈ D , d′ ∈ Dt , d , t

(36)

Pst , d = Pst , d′

∀ s ∈ S , d ∈ D , d′ ∈ Dt , d , t

(37)

Dst , d = Dst , d′

∀ s ∈ S , d ∈ D , d′ ∈ Dt , d , t

(38)

Ust , d = Ust , d′

∀ s ∈ S , d ∈ D , d′ ∈ Dt , d , t

(39)

wvit,,jd, n

=

bit,,jd, n = bit,,jd, n′

(40)

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , d′ ∈ Dt , d , t

Tf it, j, d, n = Tf it, j, d, n′ (42)

Tsit,,jd, n = Tsit,,jd, n′ ∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , d′ ∈ Dt , d , t stst,,nd = stst,,nd′

∀ s ∈ S , n ∈ N , d ∈ D , d′ ∈ Dt , d , t

stinst , d = stinst , d′

∀ s ∈ S , d ∈ D , d′ ∈ Dt , d , t

1 2 3 4 5 6

600 400 2200 400 2200 800

1400 400 800 200 2200 200

no. of periods

scenarios

binary variables

continuous variables

constraints

3 4 5 6

9 27 81 243

5 400 21 600 81 000 291 600

15 007 58 045 211 891 747 625

45 117 186 381 712 233 2 595 969

Table 3. Model Statistics of the Subproblem of (P1) with Different Planning Periods

(41)

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , d′ ∈ Dt , d , t

expectation of demand for product 2

Table 2. Model Statistics of the Extensive Form of (P1) with Different Planning Periods

wvit,,jd, n′

∀ i ∈ I , j ∈ Ji , n ∈ N , d ∈ D , d′ ∈ Dt , d , t

expectation of demand for product 1

We suppose that the demand uncertainty in each planning period fluctuates and follows a multinomial distribution with three possible values. Therefore, every parent node in the scenario tree owns three child nodes, average node, 10% above average node, and 10% below average node. There are two products in the STN case and three demand uncertainties for each product in each planning period. Apparently, there could be nine combinations of child node for each parent node. We suppose that the demands for product 1 and product 2 are perfectly correlated and the two products have the same behavior at each stage of the scenario tree, while both products have its own demand probability distribution. In other words, if the market condition is good at stage t, the demand scenario for products P1 and P2 can be expected to be high and only three combinations are considered in this case study. The statistical data of the extensive form and the subproblem with planning periods ranging from 3 to 6 are shown in Table 2 and Table 3. With the increase of planning

j ∈ Ji

∀ s ∈ S , n ≥ 2, d ∈ D , t

planning period

no. of periods

scenarios

binary variables

continuous variables

constraints

3 4 5 6

9 27 81 243

600 800 1000 1200

1375 1833 2291 2749

3697 4929 6161 7393

(43)

periods, the integrated problem scales up not only from the increase of scenarios but also from the increasing complexity of scenario subproblems. These two factors result in an intractable model size. Obviously, scenario subproblems are much easier to solve compared to the integrated problem due to its significantly reduced scale. The most straightforward method to solve a multistage stochastic program involves generating the deterministic equivalent form and then invoking a standard deterministic mixed-integer programming solver. The direct approach is first implemented to provide a baseline for the proposed method. To test the performance of the enhanced PH algorithm with different commercial solvers, we adopt both CPLEX and Gurobi as the scenario subproblem solver in the PH algorithm.

(44) (45)

4.1.2. Computational Results. In this example, the planning horizon is divided into a number of planning periods with equal one-week lengths, and the number of event points is set to 12 in the continuous time scheduling model. Four cases with 3-weeks, 4-weeks, 5-weeks, and 6-weeks planning horizon are considered. Correspondingly, the integrated planning and scheduling model changes from a three-stage stochastic model into six-stage stochastic model. The expectation of demand data distribution of six planning periods is presented in Table 1. H

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

scenarios. Moreover, the application of the proximal penalty terms linearization strategy will not affect the solution accuracy of the final solution due to the early terminate strategy. Three ρ computation methods are tested and the cost-proportional method performs better than the fixed method and heuristic method in this case study. The coefficient of the costproportional method is set to 0.01. According to the mipgap tolerance strategy proposed in section 3, the mipgap parameter is set to 0.3 for the first iteration to accelerate the convergence and 0.1 for the remainder iterations to guarantee the optimality. The binary variables are fixed when they have satisfied nonanticipativity for three iterations. The iterations are terminated when 90% binary variables are fixed. The detailed computational performance of the enhanced PH algorithm is given in Table 7 and Table 8.

All the computations in this example are performed on an 8core windows workstation with 3.5 GHz CPU and 16 Gb RAM. The mipgap parameter is set to 0.03, which implies that the algorithm will terminate when the dual gap reaches 3%.Table 4 and Table 5 give the computational performance of CPLEX and Gurobi. Table 4. Results of Solving the Extensive Form Directly with CPLEX no. of periods

no. of scenarios

objective value

gap

bound

solution time (s)

3 4 5 6

9 27 81 243

974 507.65 1 223 088.04 2 071 403.91 3 010 031.61

0.42% 0.99% 1.40% 2.87%

970 403.21 1 211 010.70 2 042 271.97 2 923 714.58

15.61 1 198.17 1 638.03 2 498.09

Table 7. Results of the PH Algorithm with Accelerating Strategies Using CPLEX

Table 5. Results of Solving the Extensive Form Directly Using Gurobi no. of periods

no. of scenarios

objective value

gap

bound

solution time (s)

no. of periods

scenario numbers

objective value

PH execution time (s)

variable fixed solving time (s)

total time (s)

3 4 5 6

9 27 81 243

974 438.75 1 220 446.38 2 069 471.29 2 978 391.86

2.38% 0.75% 1.26% 1.88%

951 224.12 1 211 356.51 2 043 483.56 2 922 310.12

12.32 1460.00 1674.00 2321.00

3 4 5 6

9 27 81 243

974 457.81 1 222 907.93 2 066 420.51 2 992 222.56

17.34 241.85 507.36 824.59

10.92 421.25 608.49 448.06

28.26 663.1 1115.85 1272.65

The results in Table 4 and Table 5 show that the time for both CPLEX and Gurobi to solve a multistage stochastic problem grows greatly with the increase of planning periods. It is reasonable according to the analysis of the problem scale in Table 2. This demonstrates the necessity for a more efficient algorithm to solve large-scale mixed-integer linear programs. The standard progressive hedging algorithm is then implemented to solve the integrated problem. In this experiment, a penalty factor ρ is set as 0.1. We choose the termdiff convergence metric as the convergence criteria, which is defined as the unscaled sum of differences between a variable value and mean. The computational results in Table 6 show

Table 8. Results of PH Algorithm with Accelerating Strategies Using Gurobi

no. of periods

scenario numbers

objective value

PH execution time(s)

CPLEX Gurobi

3 3

9 9

975 304.56 975 304.56

22.96 19.55

scenario numbers

objective value

PH execution time (s)

variable fixed solving time (s)

total time (s)

3 4 5 6

9 27 81 243

974 438.75 1 220 343.34 2 067 897.68 2 970 047.74

29.80 84.26 201.76 634.00

9.90 857.00 826.00 908.00

39.70 941.26 1027.76 1542.00

We compare the objective value and computation time of the PH algorithm with accelerating strategies and the direct approach. To help us analyze the algorithm performance, Rate1 and Rate2 are introduced, which are calculated by eqs 46 and 47. Objdirect and ObjPH represent the objective value form direct approach and accelerated PH algorithm. Timedirect and TimePH represent the computation time form direct approach and accelerated PH algorithm. The values of Rate1 and Rate2 are presented in Figure 6.

Table 6. Results of Standard PH Algorithm solver

no. of periods

that the standard PH algorithm is only capable of solving the three-planning-periods case and failed to converge when the number of planning periods grows to 4, 5, and 6. More precisely, it is hard for the aggregated decision variables (X̅ t,d and Y̅ t,d) to satisfy all scenarios when the number of scenarios increases. This is primarily due to the increasing deviation of accumulated demands between different scenarios. Therefore, we may regard that the standard PH algorithm is not suitable for a multistage stochastic MILP with complex scenario tree structure in this case study. The progressive hedging algorithm with accelerating strategies is then applied to accelerate the convergence of the PH algorithm. Among all strategies stated in section 3.2, the variable fixing strategy and early termination of the PH strategy are the key factors that solve the divergence problem of the standard PH algorithm, because the solution of the extension form with partial variables fixed destines to satisfy all

Rate1 =

Rate 2 =

Objdirect − ObjPH Objdirect Timedirect − TimePH Timedirect

(46)

(47)

It is not hard to find that the effect of the PH algorithm is not obvious or even negative when the model scale is small like the 3-planning-periods case. This is due to the overhead associated with communicating with scenario subproblem solvers in each PH iteration. However, this overhead becomes less significant when the total solution time increases. When the model scale expands, the implementation of progressive hedging algorithm shortens about one-third of the computation time compared to the direct approach in this case study. Taking the five-planning-periods case study as an example, the I

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 6. Rate1 and Rate2 with different planning periods.

Figure 7. Flowchart of a typical ethylene plant.

PH algorithm with accelerating strategies finds a 0.24% better solution in only 69% running time using the CPLEX solver. The algorithm finds almost the same solution in only 62% running time using the Gurobi solver. The same situation applies to other case studies in Figure 6. Therefore, we can draw the conclusion that the PH algorithm with accelerating strategies outperforms the directly invoking solvers. The

advantages of the PH algorithm with accelerating strategies do not decline with the growth of scenarios. 4.2. Process Industry Example. 4.2.1. Case Study Description and Model Formulation. The topology of an ethylene plant is shown in Figure 7. The plant can be divided into four parts, the in-plant site, the ethylene unit area, the downstream production area, and the out-plant site. The inplant site imports and stores the raw materials for the ethylene J

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research units, which include the liquefied petroleum gas (LPG), light diesel (LD), naphtha (NA), ethane (ETH), and cracked C4 fraction (CC4). Six major products are produced, including ethylene, propylene, hydrogasoline, glycol, ethylene oxide, and butadiene. The objective function 48 of the optimization problem is to maximize the profit in terms of sales revenue, raw materials cost, backorder cost, operating cost, and operation mode switch cost. ij j max profit = ∑ pd jjjj ∑ ∑ ∑ pricelprodl , n , t jj d k t ∈ T l ∈ Lo n ∈ N − ∑ ∑ ∑ pricer consl , n , t − ∑ ∑ ulUl , t t ∈ T l ∈ Lc , n ∈ N

−

t

l ∈ Lo

∑ ∑ ∑ ∑ ∑ αi ′ , i , j , n, tcoj − ∑ ∑ ∑ ∑ t ∈ T i ′∈ Ij i ∈ Ij j ∈ J n ∈ N

t ∈ T i ∈ Ij j ∈ J n ∈ N

SPlt,,id, j , n = ρl P, i , j

∑

DoKlt,,kd, o , n +

k ∈ Ko

SClt,,id, j , n = ρlc, i

conslt,,nd

=

∑ k ∈ Kc

+

∑

yz z upi , j ti , j , n , t zzzz zz { (48)

∑ SClt′,,di ,j ,n ,

l ′∈ lic

∑ SPlt,,id,j ,n = ∑ +

FJJlt, ,jd, j ′ , n +

j ′∈ Jj ∩ Jlc

i ∈ Ij

∑ o ∈ Ojp

DoJlt, ,jd, o , n

∑ k ∈ Kl ∩ K jp

FiKlt,,jd, k , n

∀ l ∈ Ljp , j ∈ Jl p , n ∈ N , d ∈ D , t

The feeding material of the unit j can come either from the upstream units or tanks or from the in-plant site c.

∑ SClt,,id,j ,n = ∑ i ∈ Ij

+

j ′∈ Jj ∩ Jl

∑

InJlt, ,cd, j , n

p

FJJlt, j, d′ , j , n +

∑ k ∈ Kl ∩ K jc

FoKlt,,jd′ , k , n

∀ l ∈ Lj , j ∈ Jlc , n ∈ N , d ∈ D , t

c ∈ Cj

(59)

InJlt, ,cd, j , n ,

j ∈ Jc

∑ prodlt,,nd ,

(57)

(58)

(49)

∀ l ∈ Lc , c ∈ C , n ∈ N , d ∈ D , t

Ult , d = Ult − 1, d + rlt , d −

(56)

The product l of a unit j can be sent to the next processing unit j′ or the storage tanks, or it can be transferred to the out-plant site o.

j ∈ Jo

InKlt,,cd, k , n

SPlt′,,di , j , n ,

∀ l ∈ Ljc , i ∈ Ij , j ∈ J , n ∈ N , d ∈ D , t

∑ DoJlt,,jd,o,n ,

∀ l ∈ Lo , o ∈ O , n ∈ N , d ∈ D , t

l ′∈ liP

∀ l ∈ Ljp , i ∈ Ij , j ∈ J , n ∈ N , d ∈ D , t

Equation 49 presents the total amount of product deliveries, including those delivered from storage tanks and those delivered directly from the processing unit. Equation 50 represents the similar constraint for raw material. Equation 51 represents a backorder balance for each product l. prodlt,,nd =

∑

Equation 60 represents the balance between the total produced material and the total consumed material.

(50)

∑ SPlt,,id,j ,n = ∑ SClt,,id,j ,n ,

∀ l ∈ Lo , d ∈ D , t

i ∈ Ij

n∈N

(51)

∀ l ∈ Lj , j ∈ Jlc , n ∈ N , d ∈ D , t

Equations 52−54 are used to impose the logical relation between the amount of transported materials and the corresponding time duration.

(52)

DoKlt,,kd, o , n ≥ FLok , o , l (Tof kt ,,od, n − Toskt ,,od, n), ∀ l ∈ Lo , k ∈ Ko , o ∈ O , n ∈ N , d ∈ D , t

t ,d t ,d R imin , j (Tf i , j , n − Tsi , j , n) ≤

(53)

l ′∈ liP

SPlt′,,di , j , n

∀ i ∈ Ij , j ∈ J , n ∈ N , d ∈ D , t (54)

(61)

Equation 62 enforces the logical relation between the time duration of each slot n and the utilization of unit j with operation mode i.

For each raw material, there is an upper bound of the average daily supply capacity. t ,d t ,d InKlt,,cd, k , n ≤ FLincup , k , l (Tif c , k , n − Tisc , k , n),

∀ l ∈ Lc , k ∈ Kc , c ∈ C , n ∈ N , d ∈ D , t

∑

t ,d t ,d ≤ R imax , j (Tf i , j , n − Tsi , j , n),

InJlt, ,cd, j , n ≥ FLinc , j , l(Tif ct,,jd, n − Tisct,,jd, n), ∀ l ∈ Lc , j ∈ Jc , c ∈ C , n ∈ N , d ∈ D , t

(60)

min The upper and lower bounds of the running capacity Rmix i,j , Ri,j for each unit j with operation mode i are given in eq 61. The production amount in each unit should be between the minimum and the maximum running capacity multiplied by the time duration.

DoJlt, j, d, o , n ≥ FLoj , o , l (Tof jt,,od, n − Tosjt,,od, n), ∀ l ∈ Lo , j ∈ Jo , o ∈ O , n ∈ N , d ∈ D , t

i ∈ Ij

Tf it, j, d, n − Tsit,,jd, n ≥ RH minwvit,,jd, n − H(1 − wvit,,jd, n),

(55)

∀ i ∈ Ij , j ∈ J , n ∈ N , d ∈ D , t

For each unit j, the amount of materials consumed/produced is calculated by the running capacity ∑l ′∈ l P SPlt′,,di , j , n and the i P production recipe ρl,i,j , ρCl,i,j in a different mode.

(62)

Equations 63 and 64 are the material balance constraints for the storage tank. K

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research stlt,,kd, n = stinlt,,kd + −

j ∈ JKp

FoKlt,,kd, j , n

∑ j ∈ Jkc

−

c∈C

∑

i

DoKlt,,kd, o , n

−

=

∑

+

j ∈ JKp

FiKlt,,jd, k , n

+

(63)

∑ c ∈ Ci

αit′ ,, di , j , n ≤ wvit′,,dj , n ,

InKlt,,cd, k , n

∀ j ∈ J , i , i′ ∈ Ij ,i ≠ i′, n < N , d ∈ D , t

∑ FoKlt,,kd,j ,n − ∑ DoKlt,,kd,o,n j ∈ Jkc

αit′ ,, di , j , n ≤ 1 + wvit,,jd, n′ −

o ∈ ol

(64)

+

Vmax k

≤

∀ j ∈ J, k ∈

×

αit′ ,, di , j , n ≥ wvit′,,dj , n + wvit,,jd, n′ − 1 −

n ∈ N , d ∈ D, t

(66)

k∈ Kj

(67)

stlt,,kd, n − 1 +

∑ Kolt,,kd, j , n ≤ ∑

(68)

k∈ Kj

i ∈ Ij ∩ Ilc

c∈C

FoKlt,,kd, j , n

Tf it, j, d, n ≥ Tsit,,jd, n ,

∀ i ∈ Ij , j ∈ Ji , n ∈ N , d ∈ D , t

Tf it, j, d, n + 1 ≥ Tf it, j, d, n ,

+

∑

≤

∀ i ∈ Ij , j ∈ Ji , n < N , d ∈ D , t

Tsit,,jd, n + 1 ≥ Tf it′,,dj , n − H(1 − wvit′,,dj , n), (69)

∀ i ∈ Ij , i′ ∈ Ij , i′ ≠ i , j ∈ J , n ∈ N , n < N , d ∈ D , t (80)

Tjf jt,,kd, n ≥ Tjsjt,,kd, n ,

(70)

∀ k ∈ K , j ∈ JKp , n ∈ N , d ∈ D , t (81)

λVkmax ,

Tjsjt,,kd, n + 1 ≥ Tjf jt,,kd, n ,

o ∈ Ol

∀ l ∈ L, k ∈ K , n ∈ N , d ∈ D, t

(78)

(79)

c ∈ Ci

DoKlt,,kd, o , n

(77)

Equations 78−88 express time constraints about the relationship between starting time and ending time for each task and unit.

∑ FiKlt,,jd, k ,n + ∑ InKlt,,cd, k , n ≤ λVkmax , j ∈ JKp

(76)

wvit,,jd, n ,

∀ j ∈ J , l ∈ LCj , n ∈ N , d ∈ D , t

i

∀ l ∈ L , k ∈ Kl , n > 1, d ∈ D , t j ∈ Jkc

wvit,,jd, n ,

∑ FiKlt,,jd, k ,n + ∑ InKlt,,cd, k , n ≤ λVkmax ,

∀ l ∈ L , k ∈ Kl , n = 1, d ∈ D , t

∑

i ∈ Ij ∩ IlP

∀ j ∈ J , l ∈ LjP , n ∈ N , d ∈ D , t

Outkt ,,od, n ,

In practical operations, the amount of material transfer of each tank k should not exceed the maximum inventory capacity. Equations 69−71 are enforce that the production plan is reliable under these restrictions. Considering there are cases in which materials are transferred in and out at the same time slot, it is not reasonable to strictly limit the transferred amount to less than the maximum capacity. The parameter λ is used to control the reliability level of the inventory capacity constraints (e.g., λ = 2). When λ is greater than 1, an extra transferred out amount is considered in the inventory capacity constraints. j ∈ JKp

wvit″,,dj , n″

Equations 76−77 enforce that if the unit j is not processing during time slot n, there is no material flows in or flows out from this unit.

∑ Kilt,,jd, k ,n ≤ ∑

∀ l ∈ L , k ∈ Kl , o ∈ Ol , n ∈ N , d ∈ D , t

Initl , k +

∑

(75)

∀ j ∈ J , k ∈ K jc , n ∈ N , d ∈ D , t

×

∑

∀ j ∈ J , i , i′ ∈ Ij ,i ≠ i′, n < N , n < n′ ≤ N , d ∈ D , t

(65)

Foklt,,kd, j , n ≤ Vkmax × Kolt,,kd, j , n ,

≤

wvit″,,dj , n″

i″∈ Ij , n″∈ n < n″< n′

∀ j ∈ J , k ∈ Kjp , n ∈ N , d ∈ D , t

Vkmax

∑

(74)

FiKlt,,jd, k , n ≤ Vkmax × Kilt,,jd, k , n ,

DoKlt,,kd, o , n

∑

∀ j ∈ J , i , i′ ∈ Ij ,i ≠ i′, n < N , n < n′ ≤ N , d ∈ D , t

Inct,,kd, n ,

Kjp ,

∑ wvit′,,dj ,n′

i″∈ Ij , n″∈ n < n″< n′

Equations 65−68 give the upper bound for the material transfer amount in each storage tank k in time slot n. The binary variables on the right-hand side are used to designate whether there will be a material transfer operation from/to the storage tank k during the time slot n. Vkmax

(73)

i ′∈ Ij

∀ l ∈ L , k ∈ Kl , n > 1, d ∈ D , t

InKlt,,cd, k , n

(72)

Since production mode switches can cause production yield disturbance and generate additional cost, it is penalized in the objective function.

o ∈ ol

stlt,,kd, n − 1

∀ i ∈ Ij , j ∈ J , n ∈ N , d ∈ D , t

i

∀ l ∈ L , k ∈ Kl , n = 1, d ∈ D , t stlt,,kd, n

∑ wvit,,jd,n ≤ 1,

∑ FiKlt,,jd, k , n + ∑ InKlt,,cd, k , n

∀ k ∈ K , j ∈ JKp , n ∈ N , n < N , d ∈ D , t

(71)

Equation 72 is the logical constraint for processing units. It enforces that only one operation mode i is active for each unit j during time slot n.

Tjf kt ,,jd, n ≥ Tjskt ,,jd, n ,

(82)

∀ k ∈ K , j ∈ JKp , n ∈ N , d ∈ D , t (83)

L

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Tf it, j, d, n ≤ Tjf kt ,,jd, n + H(2 − wvit,,jd, n − Kolt,,kd, j , n)

Tjskt ,,jd, n + 1 ≥ Tjf kt ,,jd, n , ∀ k ∈ K , j ∈ JKp , n ∈ N , n < N , d ∈ D , t

Tof kt ,,od, n ≥ Toskt ,,od, n ,

∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

(84)

Tf it, j, d, n ≥ Tjf kt ,,jd, n − H(2 − wvit,,jd, n − Kolt,,kd, j , n)

∀ k ∈ Ko , o ∈ O , n ∈ N , d ∈ D , t (85)

Toskt ,,od, n + 1

≥

∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

Tof kt ,,od, n ,

∀ k ∈ Ko , o ∈ O , n ∈ N , n < N , d ∈ D , t

Tif ct,,kd, n ≥ Tisct,,kd, n ,

(86)

(87)

Tisct,,kd, n + 1 ≥ Tif ct,,kd, n , (88)

For each tank that is connected to more than one unit, the starting time of the next transfer operation from unit j′ in time slot n + 1 should be larger than the ending time of current transfer operation from unit j in time slot n. Tjsjt′,,dk , n + 1 ≥ Tjf jt,,kd, n − H(2 − Kilt,,jd, k , n − Kilt,,jd′ , k , n)

Table 9. Expectation of Demands Major Products

∀ k ∈ K , j , j′ ∈ JKp , j′ ≠ j , n ∈ N , n < N , d ∈ D , t (89)

Tjsk , j ′ , n + 1, t ≥ Tjf kt ,,jd, n − H(2 − Kolt,,kd, j , n − Kolt,,kd, j ′ , n) ∀ k ∈ K , j , j′ ∈ JKp , j′ ≠ j , n ∈ N , n < N , d ∈ D , t (90)

For the intermediate storage tanks which are connected by both upstream and downstream units, the processing operation of the units and the material transfer operations should occur t,d simultaneously. Notice that the timing variables Tst,d i,j,n, Tf i,j,n should be less than the scheduling horizon H. Tsit,,jd, n

≤

Tjsjt,,kd, n

+ H(2 −

wvit,,jd, n

−

(91)

Tsit,,jd, n ≥ Tjsjt,,kd, n − H(2 − wvit,,jd, n − Kilt,,jd, k , n) ∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

(93)

(94)

Tsit,,jd, n ≤ Tjskt ,,jd, n + H(2 − wvit,,jd, n − Kolt,,kd, j , n) ∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

(95)

Tsit,,jd, n ≥ Tjskt ,,jd, n − H(2 − wvit,,jd, n − Kolt,,kd, j , n) ∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

second quarter

third quarter

21 600 17 550 10 800 14 400 14 400 1287 000

25 200 14 850 10 800 7 200 9 000 14 040

21 600 19 170 9 000 10 800 10 800 9 900

raw material

flow rate (ton/day)

LPG naphtha light diesel ethane cracked C4 fraction

400 1460 300 500 810

Table 10. The variation between the low level, median level, and high level of demand is set at 5%. The production yield of units in ethylene plant with different operation modes is presented in the Supporting Information. The extensive form of the nine-months case problem has 19 278 binary variables, 225 643 continuous variables, and 509 760 constraints. The corresponding subproblem has 2142 binary variables, 23 791 continuous variables, and 51 238 constraints. This example is also first solved directly by CPLEX and Gurobi. Then the enhanced PH algorithm is applied to solve it. All the configurations of parameters and accelerating strategies of the enhanced PH algorithm in this example are the same as the configurations in the STN example. All the computations in this example are performed on a 16-core Linux server with 2.10 GHz CPU and 128 Gb RAM. The computational performance of the direct approach and the enhance PH algorithm are shown in Table 11 and Table 12.

Tf it, j, d, n ≥ Tjf jt,,kd, n − H(2 − wvit,,jd, n − Kilt,,jd, k , n) ∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

first quarter

(92)

Tf it, j, d, n ≤ Tjf jt,,kd, n + H(2 − wvit,,jd, n − Kilt,,jd, k , n) ∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

products ethylene propylene hydrogasoline glycol ethylene oxide butadiene

Table 10. Flow Rate Data for Feeding Material

Kilt,,jd, k , n)

∀ l ∈ L , k ∈ Kl , j ∈ JKp , i ∈ Ij , n ∈ N , d ∈ D , t

(98)

The nonanticipativity constraints for this example are similar to eqs 36−45 and they are neglected here. 4.2.2. Computational Results. In this example, the length of each planning period is set to 120 days and a nine-months planning horizon is considered. The data used for this case are based on the realistic production of a Sinopec plant. The parameter λ in eqs 69−71 is set to 2. As in the first example, we suppose that the demands of the six main products in each planning period follow multinomial distributions with three possible values (low, median, and high) and the fluctuations of these demands are perfectly correlated in each planning period. The history data on the product demand orders were collected to fit the distribution of the demand uncertainty. The expectation of demand distribution information is given in Table 9. The flow rate data for feeding materials is given in

∀ k ∈ Kc , c ∈ C , n ∈ N , d ∈ D , t

∀ k ∈ Kc , c ∈ C , n ∈ N , n < N , d ∈ D , t

(97)

(96) M

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

is first solved and then the expected results RP2 are computed by fixing the first-stage variables to the two-stage solution in the multistage stochastic model. To analyze the sensitivity of the multistage stochastic model, we create two new sets of demand scenarios with higher variations, 15% and 20%, between demand uncertainty realizations.44 The expected value-based deterministic, twostage stochastic and multistage stochastic models are solved again with new scenario sets, and the results are reported in Table 13. Compared to the expected value-based deterministic

Table 11. Results of Solving the Extensive Form Directly solver

objective value

gap

bound

solution time(s)

CPLEX Gurobi

1.37821 × 108 1.40702 × 108

8.18% 5.56%

1.49089 × 108 1.43320 × 108

13 736.07 9 657.08

Due to the intractable scale of this example, it takes about ten thousand seconds to get an optimal solution when the problem is directly solved by CPLEX and Gurobi. The enhanced Progressive Hedging algorithm shortens more than 35% solution time and provides a better solution compared to the direct method. For example, Gurobi finds an objective function value of 1.40702 × 108 in 9657.08 s, while the enhanced PH algorithm finds an objective function value of 1.41332 × 1008 only in 6265.08 s. The results of the ethylene plant case reinforce the conclusion that the PH algorithm with accelerating strategies outperforms both the conventional PH algorithm and directly invoking solvers. Besides the strategies stated in section 3, the parallelization of PH is an effective strategy to solve problems with a large number of scenarios. In the above implementations of the PH algorithm, subproblems are sequentially solved by commercial solvers. This process can be parallelized through MPI and the Pyro4 module. However, CPU utilization has closely reached 100% even in the sequential condition, because both Gurobi and CPLEX use all of the available CPU cores on the computer. Therefore, the implementation of parallelization instead results in a slow-down due to the overhead associated with the contention among the threads. 4.3. The Value of the Multistage Stochastic Model. To quantify the value of the multistage stochastic program proposed in section 2, we take the five-planning-periods STN case as an example to compare the solution of the multistage stochastic program with an expected value-based deterministic model and a two-stage stochastic model. The value of the stochastic solution (VSS) is first proposed by Brige and Louneaux to compare stochastic solutions with deterministic solutions.43 The VSS is defined as the difference between the expected result of using the expected value problem solution (EEV) and the solution value of the multistage stochastic program (RPM, where subscript “M” denotes multistage). The expected value (EV) problem is derived by replacing the demand uncertainties in each stage by their expected values. The expected result of using EV solution is obtained through solving the multistage stochastic program with first-stage variables fixed to the EV solution. VSS = EEV − RPM

Table 13. Value of Multistage Stochastic Model cases

10% variation

15% variation

20% variation

EEV RP2 RPM EEV-RPM RP2-RPM

2 071 170.80 2 069 874.83 2 066 420.51 4 750.29 3 454.32

2 065 292.45 2 063 329.35 2 062 017.25 3 275.20 1 312.10

2 070 504.61 2 069 155.32 2 068 034.03 2 470.58 1 121.29

model, the multistage stochastic model always provides a more robust solution under different demand uncertainty variations. The superiority is less significant when it is compared to the two-stage stochastic solution. This is reasonable because the scenario tree of the two-stage stochastic program does not describe the evolvement of demand uncertainty as precisely as that of the multistage stochastic program. The multistage stochastic program also extends superiority by taking multistep recourse actions according to the realization of uncertainties. The results justify the adoption of multistage stochastic programming in integrated production planning and scheduling under demand uncertainty.

5. CONCLUSION In this paper, we mainly demonstrate the advances of the accelerated progressive hedging (PH) algorithm for integrated planning and scheduling problems under demand uncertainty. A general multistage stochastic formulation is first presented to address the problem of integrated production planning and scheduling under demand uncertainty. On the basis of the general structure of multistage stochastic mixed-integer programming, we apply the PH algorithm to decompose the integrated model into scenario subproblems. Several accelerating strategies are designed to guarantee the convergence of the PH algorithm. The results from case studies show that the accelerated PH algorithm outperforms a directly invoking commercial solver by generating an even better solution within about 65% of the computation time. The value of the multistage stochastic model is demonstrated through the comparison with expected value based deterministic model and the two-stage stochastic model. It is possible to extend the study in several ways, which can be suggested as future research areas. The first direction for further research is to deal with the scenario tree generation method. In this paper, we focus on solution efficiency of the integrated model and directly generate a scenario tree from the

(99)

Furtherly, we formulate and solve a two-stage stochastic program with 81 scenarios. In the two-stage stochastic model, the five planning periods are divided into two parts. The first stage contains the decisions in the first planning periods and the second stage contains the remainder of the decisions. The realizations of demand in these 81 scenarios are the same as the multistage stochastic model. Therefore, the scenario tree of the two-stage model has one node in the first stage and 81 leaf nodes in the second stage. Similar to EEV, the two-stage model Table 12. Results of Enhanced PH Algorithm solver

objective value

bound

PH execution time(s)

variable fixed solving time(s)

total time (s)

CPLEX Gurobi

1.40883 × 108 1.41332 × 108

1.47854 × 108 1.42660 × 108

1283.58 3024.22

5469.63 3240.86

6753.21 6265.08

N

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research cS pd T

probability distribution of the demand uncertainty. Efficient scenario tree generation methods, such as moment matching method and clustering method can also predigest model complexity. What is more, only demand uncertainty is considered in this paper. However, there are many other uncertainties that need to be considered in production planning and scheduling, such as the fluctuation of temperature, kinetic constants, and product yields. The modeling of multiuncertainty and the interaction between these uncertainties are worth studying. The parallelization of the PH algorithm has not been fully explored due to the restriction of computing environments, which can also be studied in the future.

■

Variables

xt,d planning level variable yt,d scheduling level variable 2. State-Task-Network Case Indices

i j n s I hs hj J hi N S hp

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.9b02620. Process data of units in motivating example; storage capacity of states in motivating example; cost data in motivating example; production yield of units in ethylene plant with different operation modes (PDF)

set set set set set set set set

hs hs φPs,i φCs,i stsmax s

AUTHOR INFORMATION

*Tel.: 086-0571-87952277. E-mail: [email protected]. ORCID

Zedong Peng: 0000-0001-6001-1738 Yi Zhang: 0000-0002-0249-0547 Gang Rong: 0000-0002-6779-2968

vmin i,j vmax i,j αi,j βi,j fci vci H Demt,d s

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge financial support from the Intelligent Manufacturing Integrated Standardization and New Model Application Project (No.2017-ZJ-003) and National Natural Science Foundation of China (61621002).

tasks tasks which produce or consume state s tasks which can be performed in unit j units units which are suitable for performing task i event points within the time horizon all involved state s products

inventory unit cost of state s backorder unit cost of product s proportion of state s produced in task i proportion of state s consumed in task i available maximum storage capacity for state s when processing task i minimum amount of unit j when processing task i maximum capacity of unit j when processing task i constant term of processing time of task i in unit j variable term of processing time of task i in unit j fixed production unit cost of task i dynamic production unit cost of task i scheduling horizon demand of product s in planning period t

Variables

Invt,d s Pt,d s Delt,d s Ut,d s wvt,d i,j,n

NOMENCLATURE

General Formulation Indices

t planning periods/stages in scenario tree d scenarios

stt,d s,n bt,d i,j,n

Sets

stint,d s Tf t,d i,j,n

D set of scenarios Dt,d set of scenarios that are equivalent (or indistinguishable) to scenario d at planning period t

Tst,d i,j,n

Parameters

AP,DP1,DP2,DP3,DP4 coefficient variables AS,DS1,DS2 coefficient variables bP,d,eP1,d,eP2,d parameter constraints bS,eS parameter constraints e parameter straints cP cost vector

of of of of of of of of

Parameters

Corresponding Author

■

tasks units event points representing the beginning of a task states

Sets

ASSOCIATED CONTENT

S Supporting Information *

■

cost vector for scheduling variables probability of scenerio d the number of planning periods

inventory level of state s at the end of planning period t production target of state s in planning period t delivery of product s in planning period t backorder of product s in planning period t whether or not task i in unit j start at event point n (binary) amount of state s at event point n amount of material undertaking task i in unit j at event point n initial inventory for state s in planning period t time at which task i finishes in unit j while it starts at event point n time at which task i starts in unit j at event point n

3. Ethylene Plant Case Indices

matrixs for planning level

k j l n c o i

matrixs for scheduling level vectors for planning level vectors for scheduling level vector for connection con-

storage tanks processing units in ethylene plant materials event point slots in-plant site ports out-plant site ports operation mode

Sets

N batch time slots

for planning variables O

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research stt,d l,k,n

O J Jj Jc Jo Jpk Kj Kl Kc Ko KPj

out-plant site ports processing units units connected with other unit j units connected with in-plant site port c units connected with out-plant site port o units with produced material connected with tank k tanks connected with other unit j tanks that store material l tanks connected with in-plant site port c tanks connected with out-plant site port o tanks connected with unit j which has produced material in tank Kcj tanks connected with unit j which has raw material in tank Lo material sending to the out-plant site port o LP major products under demand uncertainty Lpj products from unit j Lc feeding material from in-plant site ports c to downstream units Lcj consumed material from unit j Lj material produced by unit j Ij operation mode with unit j Cj in-plant site ports connected with unit j Oj out-plant site ports connected with unit j

Tof t,d k,o,n Tost,d j,o,n Tst,d i,j,n Tf t,d i,j,n Tjst,d j,k,n Tjf t,d j,k,n Tjst,d k,j,n Tjf t,d k,j,n Tjst,d c,k,n Tjft,d c,k,n Tost,d k,o,n Tof t,d k,o,n

Parameters

H FLoj,o,l FLok,o,l FLinc,j,l FLinup c,k,l rt,d l ρPl,i,j ρcl,i,j Rmax i,j Rmin i,j Initl,k Vmax k

scheduling horizon of the integrated model out flow rate for material l from unit j to port o out flow rate for material l from tank k to port o feeding flow rate for material l from port c to unit j upper bound of feeding flow rate for material l from port c to tank k demand orders for material l during planning horizon roduction ratio of material l of operation mode i in unit j consumption ratio of material l of operation mode i in unit j maximum running capacity of unit j in operation mode i minimum running capacity of unit j in operation mode i initial amount of tank k with material l maximum storage capacity for tank k

current amount of material l in storage tank k in time slot n finishing time of material transfer from storage tank k to out-site plant o in time slot n starting time of material transfer from unit j to outsite plant o in time slot n starting time of operation mode i of unit j of time slot n finishing time of operation mode i of unit j of time slot n starting time of material transfer from unit j to storage k in time slot n finishing time of material transfer from unit j to storage k in time slot n starting time of material transfer from storage k to unit j in time slot n finishing time of material transfer from storage k to unit j in time slot n starting time of material transfer from in-plant site c storage k in time slot n finishing time of material transfer from in-plant site c storage k in time slot n starting time of material transfer from storage k to out-plant site o in time slot n finishing time of material transfer from storage k to out-plant site o in time slot n

Variables (Binary)

Int,d c,k,n t,d Kil,j,k,n

Kot,d l,k,j,n Outt,d k,o,n wvt,d i,j,n t,d αi′,i,j,n

■

binary variable which indicates the material transfer from in-plant site c to storage tank k in time slot n binary variable which indicates the material transfer of material l from unit j to storage tank k in time slot n binary variable which indicates the material transfer of material l from storage tank k to unit j in time slot n binary variable which indicates the material transfer from storage tank k to out-plant site o in time slot n binary variable which indicates the utilization of operation mode i in unit j binary variable indicates the mode switching from i to i’ in unit j of time slot n

REFERENCES

(1) Grossmann, I. Enterprise-Wide Optimization: A New Frontier in Process Systems Engineering. AIChE J. 2005, 51 (7), 1846−1857. (2) Grossmann, I. E. Advances in Mathematical Programming Models for Enterprise-Wide Optimization. Comput. Chem. Eng. 2012, 47, 2−18. (3) Dias, L. S.; Ierapetritou, M. G. From Process Control to Supply Chain Management: An Overview of Integrated Decision Making Strategies. Comput. Chem. Eng. 2017, 106, 826−835. (4) Maravelias, C. T.; Sung, C. Integration of Production Planning and Scheduling: Overview, Challenges and Opportunities. Comput. Chem. Eng. 2009, 33 (12), 1919−1930. (5) Bassett, M. H.; Dave, P.; Doyle, F. J.; Kudva, G. K.; Pekny, J. F.; Reklaitis, G. V.; Subrahmanyam, S.; Miller, D. L.; Zentner, M. G. Perspectives on Model Based Integration of Process Operations. Comput. Chem. Eng. 1996, 20 (6), 821−844. (6) Chen, P.; Papageorgiou, L. G.; Pinto, J. M. Medium-Term Planning of Single-Stage Single-Unit Multiproduct Plants Using a Hybrid Discrete/Continuous-Time MILP Model. Ind. Eng. Chem. Res. 2008, 47 (6), 1925−1934. (7) Dogan, M. E.; Grossmann, I. E. A Decomposition Method for the Simultaneous Planning and Scheduling of Single-Stage Continuous Multiproduct Plants. Ind. Eng. Chem. Res. 2006, 45 (1), 299− 315.

Variables (Continuous)

prodt,d total produced material of event point n of material l l,n const,d total consumed material of event point n of material l l,n InKt,d l,c,k,n feeding material from in-plant site c to raw material storage tank k with material l in time slot n InJt,d feeding material from in-plant site c to unit j with l,c,j,n material l in time slot n DoKt,d l,k,o,n produced product l from storage tank k to out-plant site o in time slot n DoJt,d produced product l from unit j to out-plant site o in l,j,o,n time slot n FJJt,d produced product l from unit j to unit j’ in time slot n l,j,j′,n FiKt,d produced product l from unit j to tank k in time slot l,j,k,n n FoKt,d l,k,j,n material l transferred from tank k to unit j in time slot n SPt,d material l produced in unit j of operation mode i in l,i,j,n time slot n SCt,d material l consumed in unit j of operation mode i in l,i,j,n time slot n P

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (8) Li, Z.; Ierapetritou, M. G. Production Planning and Scheduling Integration through Augmented Lagrangian Optimization. Comput. Chem. Eng. 2010, 34 (6), 996−1006. (9) Váncza, J.; Kis, T.; Kovács, A. Aggregation - the Key to Integrating Production Planning and Scheduling. CIRP Ann. 2004, 53 (1), 377−380. (10) Chu, Y.; You, F. Integrated Planning, Scheduling, and Dynamic Optimization for Batch Processes: MINLP Model Formulation and Efficient Solution Methods via Surrogate Modeling. Ind. Eng. Chem. Res. 2014, 53 (34), 13391−13411. (11) Verderame, P. M.; Floudas, C. A. Integrated Operational Planning and Medium-Term Scheduling for Large-Scale Industrial Batch Plants. Ind. Eng. Chem. Res. 2008, 47 (14), 4845−4860. (12) Lima, R. M.; Grossmann, I. E.; Jiao, Y. Long-Term Scheduling of a Single-Unit Multi-Product Continuous Process to Manufacture High Performance Glass. Comput. Chem. Eng. 2011, 35 (3), 554−574. (13) Sung, C.; Maravelias, C. T. An Attainable Region Approach for Production Planning of Multiproduct Processes. AIChE J. 2007, 53 (5), 1298−1315. (14) Apap, R. M.; Grossmann, I. E. Models and Computational Strategies for Multistage Stochastic Programming under Endogenous and Exogenous Uncertainties. Comput. Chem. Eng. 2017, 103, 233− 274. (15) Zhang, Y.; Feng, Y.; Rong, G. Data-Driven Chance Constrained and Robust Optimization under Matrix Uncertainty. Ind. Eng. Chem. Res. 2016, 55 (21), 6145−6160. (16) Grossmann, I. E.; Apap, R. M.; Calfa, B. A.; García-Herreros, P.; Zhang, Q. Recent Advances in Mathematical Programming Techniques for the Optimization of Process Systems under Uncertainty. Comput. Chem. Eng. 2016, 91, 3−14. (17) Wu, D.; Ierapetritou, M. Hierarchical Approach for Production Planning and Scheduling under Uncertainty. Chem. Eng. Process. 2007, 46 (11), 1129−1140. (18) Chu, Y.; You, F.; Wassick, J. M.; Agarwal, A. Integrated Planning and Scheduling under Production Uncertainties: Bi-Level Model Formulation and Hybrid Solution Method. Comput. Chem. Eng. 2015, 72, 255−272. (19) Van Slyke, R.; Wets, R. L-Shaped Linear Programs with Applications to Optimal Control and Stochastic Programming. SIAM J. Appl. Math. 1969, 17 (4), 638−663. (20) Carøe, C. C.; Schultz, R. Dual Decomposition in Stochastic Integer Programming. Operations Research Letters 1999, 24 (1), 37− 45. (21) Rockafellar, R. T.; Wets, R. J.-B. Scenarios and Policy Aggregation in Optimization Under Uncertainty. Mathematics of Operations Research 1991, 16 (1), 119−147. (22) Geoffrion, A. M. Generalized Benders Decomposition. J. Optim Theory Appl. 1972, 10 (4), 237−260. (23) Boyd, S. Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers. FNT in Machine Learning 2010, 3 (1), 1−122. (24) Veliz, F. B.; Watson, J.-P.; Weintraub, A.; Wets, R. J.-B.; Woodruff, D. L. Stochastic Optimization Models in Forest Planning: A Progressive Hedging Solution Approach. Ann. Oper Res. 2014, 232 (1), 259−274. (25) Løkketangen, A.; Woodruff, D. L. Progressive Hedging and Tabu Search Applied to Mixed Integer (0,1) Multistage Stochastic Programming. J. Heuristics 1996, 2 (2), 111−128. (26) Crainic, T. G.; Fu, X.; Gendreau, M.; Rei, W.; Wallace, S. W. Progressive Hedging based Metaheuristics for Stochastic Network Design. Networks 2011, 58 (2), 114−124. (27) Gonçalves, R. E. C.; Finardi, E. C.; Silva, E. L. da. Applying Different Decomposition Schemes Using the Progressive Hedging Algorithm to the Operation Planning Problem of a Hydrothermal System. Electr. Power Syst. Res. 2012, 83 (1), 19−27. (28) Lamghari, A.; Dimitrakopoulos, R. Progressive Hedging Applied as a Metaheuristic to Schedule Production in Open-Pit Mines Accounting for Reserve Uncertainty. European Journal of Operational Research 2016, 253 (3), 843−855.

(29) Watson, J.-P.; Woodruff, D. L. Progressive Hedging Innovations for a Class of Stochastic Mixed-Integer Resource Allocation Problems. Computational Management Science 2011, 8 (4), 355−370. (30) Munoz, F. D.; Watson, J.-P. A Scalable Solution Framework for Stochastic Transmission and Generation Planning Problems. Comput. Manag Sci. 2015, 12 (4), 491−518. (31) Gade, D.; Hackebeil, G.; Ryan, S. M.; Watson, J.-P.; Wets, R. J.B.; Woodruff, D. L. Obtaining Lower Bounds from the Progressive Hedging Algorithm for Stochastic Mixed-Integer Programs. Math. Program. 2016, 157 (1), 47−67. (32) Shapiro, A.; Dentcheva, D.; Ruszczyński, A. Lectures on Stochastic Programming: Modeling and Theory; Society for Industrial and Applied Mathematics, 2009. DOI: 10.1137/1.9780898718751. (33) Heitsch, H.; Römisch, W. Scenario Tree Modeling for Multistage Stochastic Programs. Math. Program. 2009, 118 (2), 371−406. (34) Kallrath, J. Planning and Scheduling in the Process Industry. OR Spectrum 2002, 24 (3), 219−250. (35) Kopanos, G. M.; Puigjaner, L.; Maravelias, C. T. Production Planning and Scheduling of Parallel Continuous Processes with Product Families. Ind. Eng. Chem. Res. 2011, 50 (3), 1369−1378. (36) Liu, S.; Pinto, J. M.; Papageorgiou, L. G. A TSP-Based MILP Model for Medium-Term Planning of Single-Stage Continuous Multiproduct Plants. Ind. Eng. Chem. Res. 2008, 47 (20), 7733−7743. (37) Vieira, M.; Pinto-Varela, T.; Moniz, S.; Barbosa-Póvoa, A. P.; Papageorgiou, L. G. Optimal Planning and Campaign Scheduling of Biopharmaceutical Processes Using a Continuous-Time Formulation. Comput. Chem. Eng. 2016, 91, 422−444. (38) Hart, W. E.; Laird, C. D.; Watson, J.-P.; Woodruff, D. L.; Hackebeil, G. A.; Nicholson, B. L.; Siirola, J. D. PyomoOptimization Modeling in Python; Springer Optimization and Its Applications; Springer International Publishing: Cham, 2017; Vol. 67. DOI: 10.1007/978-3-319-58821-6. (39) Stellato, B.; Banjac, G.; Goulart, P.; Bemporad, A.; Boyd, S. OSQP: An Operator Splitting Solver for Quadratic Programs. In 2018 UKACC 12th International Conference on Control (CONTROL) 2018, 339−339. (40) Watson, J.-P.; Woodruff, D. L.; Hart, W. E. PySP: Modeling and Solving Stochastic Programs in Python. Math. Prog. Comp. 2012, 4 (2), 109−149. (41) Hart, W. E.; Watson, J.-P.; Woodruff, D. L. Pyomo: Modeling and Solving Mathematical Programs in Python. Math. Prog. Comp. 2011, 3 (3), 219. (42) Ierapetritou, M. G.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 1. Multipurpose Batch Processes. Ind. Eng. Chem. Res. 1998, 37 (11), 4341−4359. (43) Birge, J. R.; Louveaux, F. Introduction to Stochastic Programming, 2nd ed.; Springer series in operations research and financial engineering; ; Springer: New York, 2011. (44) Xie, F.; Huang, Y. A Multistage Stochastic Programming Model for a Multi-Period Strategic Expansion of Biofuel Supply Chain under Evolving Uncertainties. Transportation Research Part E: Logistics and Transportation Review 2018, 111, 130−148.

Q

DOI: 10.1021/acs.iecr.9b02620 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX