Deterministic Planning Model and Robust Optimization for Demand

Feb 24, 2012 - Operational Planning of Large-Scale Continuous Processes: Deterministic Planning Model and Robust Optimization for Demand Amount and ...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/IECR

Operational Planning of Large-Scale Continuous Processes: Deterministic Planning Model and Robust Optimization for Demand Amount and Due Date Uncertainty Jie Li,†,‡ Peter M. Verderame,† and Christodoulos A. Floudas†,* †

Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544-5263, United States State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, P. R. China



S Supporting Information *

ABSTRACT: The operational planning problem typically determines both the plant’s aggregate and daily production targets to meet product demands over a long time horizon varying from one to three months. The daily production profile generated from the operational planning level is used to schedule the day-to-day operations of the plant. In this paper, we first extend the production disaggregation model framework proposed by Verderame and Floudas, Ind. Eng. Chem. Res. 2008, 47, 4845−4860 to develop a novel mixed-integer linear programming operational planning model based on discrete-time representation for a largescale multiproduct continuous plant. We allow unused processing time to carry over from one day to the next and thus capture the continuous time nature of the plant. The production totals are disaggregated into a feasible distribution of daily production requirements, and daily production profiles are provided to meet the requirements of the medium-term scheduling model developed by Shaik and Floudas, Comput. Chem. Eng. 2009, 33, 670−686. Then, we extend the work of Lin et al. Ind. Eng. Chem. Res. 2004, 28, 1069−1085 and Janak et al. Ind. Eng. Chem. Res. 2007, 31, 171−195 to address various forms of demand amount and due date uncertainty. The computational results show that the proposed deterministic and robust planning models successfully provide tight upper bounds on the plant’s true production capacity.

1. INTRODUCTION The planning decisions can be divided into strategic planning, tactical planning (operational planning) and short-term planning.1 While the strategic (long-term) planning determines the structure of the supply chain, for instance, facility location over one or several years, short-term planning, which is referred to as scheduling at production levels, determines the allocation of tasks to units and their sequencing and timing in each unit on a daily or weekly basis. The operational planning problem typically determines both the plant’s aggregate and daily production targets for the plant to meet product demands over a long time horizon from one to three months. In this paper, we mainly concentrate on operational planning. The operational planning problems have received much attention in the literature. Many planning models have been developed for various industrial problems. Wilkinson et al.2 proposed an aggregate planning model for multipurpose batch plants. McDonald and Karimi3 developed a multiperiod midterm planning models for parallel semicontinuous processes. Kallrath and Timpe,4 Kallrath,5 and Kallrath6 developed a multisite, discretetime planning model, which divides the planning horizon into commercial periods and then further discretized those commercial periods into production periods. Erdirik-Dogan and Grossmann7 proposed a planning model for batch processes which addressed the sequencing of tasks accurately. Recently, Verderame and Floudas8,9 developed novel discrete-time mixedinteger linear programming (MILP) operational planning with production disaggregation models (PPDM) for a large-scale multipurpose and multiproduct batch plant. They disaggregated © 2012 American Chemical Society

the production totals into a feasible distribution of daily production requirements and hence provided the daily production profile. They also allowed any unutilized reactor time to carry over from one day to the next and hence captured the continuous time nature of the plant. The computational results show that the proposed framework provides a tighter upper bound on the true capacity of the plant. An overview of various aspects on planning can be referred to Kallrath,10 Pochet and Wolsey,11 and the recent review of Verderame et al.12 All parameters in the aforementioned models are assumed to be deterministic in nature. In real practice, frequent uncertainties are unavoidable such as demand amount fluctuations, due date variations, and uncertainty on processing rate. Neglecting the uncertain nature may cause a detrimental effect on company operations.13 Factoring the uncertain parameters into the analysis would greatly enhance the robustness of the planning decision. Verderame et al.12 presented a detailed review on planning and scheduling under uncertainty. Verderame and Floudas14 applied the robust optimization framework to address demand due date and demand uncertainty in the problem of operational planning of large-scale industrial batch plants.8 Note that the robust optimizaiton framework is originally proposed for linear and quadratic programming by Ben-Tal, Nemirovsky, and co-workers.15−18 Received: Revised: Accepted: Published: 4347

November 19, 2011 February 4, 2012 February 24, 2012 February 24, 2012 dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

Figure 1. Schematic of an industrial continuous plant.

a tight upper bound on the plant’s production capacity. Then, we extend the work of Lin et al.,21 Janak et al.,22 and Li et al.23 to address various forms of demand uncertainty including demand due date and amount uncertainty.

A similar robust optimization framework was also independently proposed by El Ghaoui and co-workers.19,20 Later it is extended for general mixed-integer linear programming problems by Lin et al.,21 Janak et al.,22 and Li et al.23 and for linear and discrete programming by Bertsimas and co-workers.24−26 Verderame and Floudas27 also applied the robust optimization framework to address demand and processing time uncertainty in the problem of integration of operational planning and medium-term scheduling for large-scale industrial batch plants,8 and demand and transportation uncertainty in a multisite planning problem.28 In addition, an alternative framework based on conditional valueat-risk theory was also introduced to address demand due date and demand amount uncertainty29 in the operational planning of batch processes,8 and demand and transportation uncertainty28 in the operational planning of multisite batch plants.9 To the best of our knowledge, no work exists in the literature that develops deterministic and robust operational planning models for large-scale continuous processes to generate daily production profiles. In this paper, we first extend the PPDM framework developed by Verderame and Floudas8 to develop a novel mixed-integer linear programming (MILP) operational planning model based on discrete-time representation for a large-scale industrial semicontinuous/continuous multiproduct plant. In the model, we allow unused processing time to carry over from one day to the next and capture the continuous time nature of the plant. The production totals are disaggregated into a feasible distribution of daily production requirements and hence daily production profiles are provided to meet the requirements of the medium-term scheduling model developed by Shaik and Floudas,30 which is a continuous-time unit-specific event-based short-term scheduling model.31−46 The computational results show that the proposed planning model provides

2. PROBLEM STATEMENT Figure 1 shows a schematic of an industrial continuous plant. It involves a total of U units (u = 1, 2, 3, ..., U), i tasks (i = 1, 2, 3, ..., I) and s (s = 1, 2, 3, ..., S) states. All tasks can be classified into five basic operations including continuous feed supply (FT), feed storage (FS), main production (M), product storage (PS), and final production packing/filling (PF). The units up to 100 in total comprise feed supply units, feed storage units, main production units, product storage units, and final packing/fill units. Each unit u processes Iu tasks. The main units are represented by Um. The states include raw material states, intermediate product states, and final product states. There are also some operational restrictions and rules that are described in detail in Shaik et al.30 The plant produces more than 100 articles that belong up to 100 different bulk products (denoted as Sp). These bulk products that are produced in the main units are packed/filled into different types of either bags, big bags, or truck articles, which are used to meet two types of demands: make-to-order (MTO) and make-to-stock (MTS). The MTO and MTS products are included into sets SpMTO and SpMTS, respectively. All truck articles (included in sets SpTR) belong to MTS. Each article has its own demand denoted as dem(s,d), which includes requested amount and intermediate due date. Compared to MTS articles, severe penalties are imposed for lateness of MTO articles. For MTS articles, severe penalties are also imposed for lateness of truck 4348

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

articles relative to nontruck articles (included in sets SpNTR). The entire problem can be stated as follows, Given the production recipe of the plant, suitable units and their capacities, the processing rates, sequence-dependent cleanup times, storage policy for the intermediates, the planning horizon, and demands and due dates of different articles that can be produced, the deterministic operational planning model is to determine the plant’s daily production profiles. During the planning stage, the most common uncertainty that needs to be considered arises from demand, processing times, and price parameters. These uncertain parameters can be described using discrete or continuous distributions. In some cases, only limited knowledge about the distribution is available, for example, the uncertainty is bounded, or the uncertainty is symmetrically distributed in a certain range. In the best situation, the distribution function for the uncertain parameter is given, for instance, as a normal distribution with known mean and standard deviation. The robust operational planning model is to explicitly address the various forms of demand uncertainty and determine the plant’s daily production profiles simultaneously.

The operational planning horizon is discretized into daily production periods, but the discretization of time does not lead to unnecessary reactor downtime due to the formulation of the timing constraints. Unsatisfied demand for bulk products is carried over and penalized on a daily basis while unsatisfied demand for packed goods is penalized on a weekly basis. Overproduction of bulk products is penalized daily to minimize the amount of bulk product inventory for any given day. The operational planning model allows for the overproduction of packed goods to provide an upper bound on the plant’s capacity, but a strict upper bound on the amount of overproduction is present within the set of production constraints.

4. DETERMINISTIC OPERATIONAL PLANNING WITH PPDM We extend the PPDM schematic proposed by Verderame and Floudas8 to this large-scale multiproduct continuous plant. Since the processing rates of the main units (M1-M10) are much smaller than those of other units, they are identified to be bottleneck units. The entire deterministic model is presented as follows. 4.1. Duration Constraints. To assign product states to units in each day, similar binary variables Y(u,s,d) to Verderame and Floudas8 are defined as follows:

3. OVERVIEW OF OPERATIONAL PLANNING WITH PRODUCTION DISAGGREGATION MODEL (PPDM) Verderame and Floudas8 developed an operational PPDM schematic (illustrated in Figure 2) for a multiproduct and

⎧1 if a state s is produced in unit u on day d Y (u , s , d ) = ⎨ ⎩ 0 otherwise ∀ u ∈ Um, s ∈ Sp , d

Similar to Verderame and Floudas,8 the total time required in each main unit u on day d can be calculated using eq 1 and it is bounded using eqs 2 and 3. Equation 2 states that the total processing time of a main unit u in the first day (i.e., d = 1) cannot exceed the total available processing time for this main unit u in day d. Equation 3 means the total processing time of a main unit u in other days (d ≠ 1) cannot exceed the total accumulated available processing times for this main unit u in day d including the accumulated unutilized processing times of previous days (1

(3)

u

where T (u,d) denotes the total time that main unit u is utilized on day d, RL(u,s,d) is defined as the processing time of main unit u to produce state s in day d, and T(u,d) is a parameter to denote the total available processing time in main unit u on day d. To ensure that the operational planning model provides an upper bound on production capacity of the plant and that all reactors are efficiently utilized, Vederame and Floudas8 4349

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

enforced that the total processing time of a unit u must exceed a user-supplied aggregate time. However, we demand the total processing time of a main unit u up to week w must exceed some percentage of total available time up to this week.

since these articles are required to be met more strictly than the nontruck MTS articles. d



T u(u , d) + slu(u , w) ≥ η· term _w(w)

∑ d ≤ term _w(w) m

∀u∈U ,w

d ′= 1



T(u , d′) (5)

tot1(u , s , d)





d ′= 1





cost = α·[

∀ s ∈ SMTO ,d p

(13)



sl1a _b(s , d) ·price(s) D



+



sl1a _p(s , d) ·price(s)]

TR s ∈ SMTS p , Sp d = 1

D

+ β·[

dem(s , d′)





sl1b _b(s , d) ·price(s)

d=1 s ∈ SMTO p D



+



sl1b _p(s , d) ·price(s)]

TR s ∈ SMTS p , Sp d = 1

d

d ′= 1

dem(s , d)

d=1 s ∈ SMTO p

(8)



∑ init _w(w) ≤ d ≤ term _w(w)

D

(7)

d ′= 1

tot(s , d′) − sl1b _b(s , d) ≤

(12)

tot(s , d) ≤ 1.2·

∀ s ∈ Sp , d

∀ s ∈ SMTO ,d p d

dem(s , d′)

where, init_w(w) is the beginning time of week w. 4.4. Objective Function. Similar to Verderame and Floudas,8 our objective is also to penalize all the slack variables in the model and maximize the gross profit from the sale of the products scheduled to be produced in the planning horizon.

(6)



∑ d ′= 1

init _w(w) ≤ d ≤ term _w(w) ∀ s ∈ SMTS ∩ SNTR ,d p p

d

tot(s , d′) + sl1a _b(s , d) ≥

d

tot(s , d′) + sl1a _p(s , d) ≥

∀ s ∈ SMTS ∩ SNTR ,d p p

4.3. Customer Demand Constraints. The daily underproduction and overproduction demand constraints for MTO products are represented by eqs 8−9, which are similar to eqs 13−14 from Verderame and Floudas8 and eqs 85−86 of Verderame and Floudas.14 d

(11)

d ′= 1

The total amount of state s produced on day d [denoted as tot(s,d)] is computed as u ∈ Um

dem(s , d′)

d ′= 1

d



tot1(u , s , d) = Rate(u , s) ·RL(u , s , d)





where, sl1a_p(s,d) and sl1b_p(s,d) are underproduction and overproduction slack variables. For the nontruck MTS articles, we enforce underproduction demand constraints and allow for some overproduction, which provides an upper bound on the production of these articles.

where, RLmin is the minimum processing time, which is used to ensure that if binary variable Y(u,s,d) is active, then the state s must be processed in unit u for some minimum time. 4.2. Production and Capacity Constraints. We use Rate(u,s) to denote the processing rate of product state s in main unit u. The total amount of product state s processed in main unit u on day d [tot1(u,s,d)] can be calculated by

tot(s , d) =

d

tot(s , d′) − sl1b _p(s , d) ≤

∀ s ∈ SMTS ∩ STR p p ,d

d ′= 1

∀ u ∈ Um, d , s ∈ Sp

(10)

d ′= 1

d

∀ u ∈ Um, d , s ∈ Sp

dem(s , d′)

d ′= 1

d

where, slack variable slu(u,w) is used to indicate that strict satisfaction of eq 4 is not always required. Parameters η is userdefined utilization percentage and term_w(w) is the total available processing time up to week w. The processing time of a main unit u to produce state s in day d is bounded by





∀ s ∈ SMTS ∩ STR p p ,d

(4)

Y (u , s , d) ·RLmin ≤ RL(u , s , d) ≤ Y (u , s , d) ·

d

tot(s , d′) + sl1a _p(s , d) ≥

D

dem(s , d′)



+ β·[

d ′= 1



sl1a _p(s , d) ·price(s)]

NTR d = 1 s ∈ SMTS p ,Sp

(9)

+ γ·

∑ ∑ slu(u , w) u ∈ Um w

where, sl1a_b(s,d) and sl1b_b(s,d) are underproduction and overproduction slack variables, respectively. We also enforce daily underproduction and overproduction demand constraints for truck MTS articles using eqs 10 and 11,



∑ ∑ tot(s , d)·price(s) s ∈ Sp d

4350

(14)

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

through scheduling represent the true production capacity of the plant. Table 1 also illustrates the aggregate production level from the medium-term scheduling model. While the operational model provides production totals of 18818.513 mu, it is 18107.800 mu from the medium-term scheduling model. Therefore, the proposed operational planning model provides an upper bound on the production capacity of the plant. We divide the entire planning horizon into 14 periods with 3 days for the first period, 7 days for the next 12 periods, and 4 days for the last period. The aggregate production totals of all articles at both planning and scheduling levels in each period are illustrated in Figure 3. From Figure 3, the largest difference of the aggregate production totals from the planning and medium-term scheduling models are around 9.9% in period 3. It demonstrates that the proposed operational planning model provides a tight upper bound on the production capacity of the plant. Figure 4 depicts the true production schedule for those bottleneck units (i.e., main units M1-M10) from the mediumterm scheduling model. Table 3 presents the utilization of bottleneck units. From Table 3, the bottleneck unit M6 reaches the highest utilization of 80.92%, and M1 has the lowest utilization of 35.61%. From Figure 4, the idle times on M1-M10 include transition times. It can also be observed that M1, M3, M8 and M9 have more idle times than other bottleneck units. This is because many articles are not allowed to be processed in M1, M3 and M8-M9, as shown in Table 4.

In eq 14, the first term is to penalize the underproduction for MTO and truck MTS articles. The second term is to penalize the overproduction for MTO and truck MTS articles. The third term is to penalize the underproduction for nontruck MTS articles. The fourth term is to penalize the time utilization and the last term is to maximize gross profit from products that are to be produced. Up to this point, we complete the deterministic operational planning model (denoted as M) comprising eqs 1−14 with cost minimization in eq 14.

5. COMPUTATIONAL STUDY: DETERMINISTIC PLANNING We use an example from a large-scale multiproduct continuous plant to evaluate the proposed operational planning model. The planning horizon is about three months. This example is solved using GAMS 22.647/CPLEX 11.0.0 on Dell OPTIPLEX 960 of Intel Xeon CPU 3.0 GHz with 2 GB RAM running Linux. We impose an upper bound of 14400 CPU seconds for the model M. The model M involves 31 155 binary variables, 89 054 continuous variables, and 122 315 constraints. Table 1 presents the aggregate Table 1. Aggregate Production Levels from Model M and Medium-Term Scheduling Model with Three-Month Planning Horizon model the proposed deterministic operational planning model the medium-term scheduling

aggregate production level (mu) 18818.513

6. ROBUST OPTIMIZATION FRAMEWORK FOR OPERATIONAL PLANNING WITH PPDM

18107.800

6.1. Overview of Robust Optimization Framework. The robust optimization framework15−18,21−23 is used to address the various forms of demand uncertainty present in the planning model. In the following, we present in brief the general robust optimization approach, which explicitly takes into account the various forms of parameter uncertainty within constraints and/or objective function.

production levels for a three-month period from the planning model. We also use the daily production profile generated from the operational planning model as an input parameter to solve the medium-term scheduling developed by Shaik et al.30 Since the medium-term scheduling model is a continuous-time formulation and assumes no unit aggregation, the production results generated

Figure 3. Aggregate production totals of all articles from model M and the medium-term scheduling model. 4351

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

Figure 4. Operational schedule for bottleneck main units from model M with three-month time horizon.

min/max c Tx + dTy

Table 2. Aggregate Production Levels from Model ROM and Medium-Term Scheduling Model with Three-Month Planning Horizon

x ,y

Ex + Fy = e

s.t.

aggregate production level (mu)

model the proposed robust operational planning model the medium-term scheduling

Ax + By ≤ p x ≤ x ≤ x̅ y = 0, 1

16877.907 16571.278

(15)

Assume that the left-hand-side coefficients A and B, and the right-hand-side parameters p of the inequality constraints are uncertain parameters. The true realization of an uncertain parameter is represented by

Table 3. Utilization of Each Bottleneck Unit for Deterministic and Robust Cases utilization (%)

a ̃ = (1 + ε·ξ) ·a

unit

deterministic

robust

M1 M2 M3 M4 M5 M6 M7 M8 M9 M10

35.61 58.28 46.94 60.09 79.40 80.92 50.67 43.26 39.56 58.99

27.22 53.18 44.35 76.35 78.39 73.95 48.32 38.58 34.68 51.77

(16)

where a is an uncertain parameter with ã being its true realization, ε represents a given (relative) uncertain level, and ξ stands for a random variable. In the robust optimization framework,15−18,21−23 a solution (x, y) is called robust if (i) (x, y) is feasible for the nominal problem, (ii) whatever are the true values of the coefficients and right-hand-side parameters, (x, y) must satisfy the lth inequality constraint with an error of at most δmax[1, |pl|], where δ is a given infeasibility tolerance. Then, any inequality constraint l in eq 15 becomes



al mx m +

m ∉ Ml

Consider a generic deterministic mixed-integer linear programming (MILP) problem:



al̃ mx m +

blk yk +

k∉Kl

m ∈ Ml

≤ pl̃ + δmax[1, |pl |]





̃ y blk k

k∈Kl

∀l

(17)

Table 4. Units, Suitable Articles, and Amounts MTO articles

truck articles

nontruck articles

unit

numbers

amount (mu)

numbers

amount (mu)

numbers

amount (mu)

total number of articles

total amount (mu)

F41 F42 F43 F44 F45 G81 G82 G83 G84 G85

4 1 1 4 0 0 0 2 0 0

90.00 6.00 30.00 45.00 0 0 0 35.40 0 0

1 3 1 6 9 1 1 2 2 3

60.00 828.00 180.00 1137.00 1844.70 546.00 594.00 810.00 990.00 945.00

10 6 9 6 18 13 12 8 9 12

1565.16 561.27 1987.36 522.42 2894.22 2783.22 1657.49 2355.83 2729.70 3177.09

15 10 11 16 27 14 13 12 11 15

1715.17 1395.27 2197.36 1704.43 4738.92 3329.22 2251.89 3201.22 3719.69 4122.09

4352

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

where, Ml and Kl are the subsets that encompass those uncertain parameters for the given constraint l. Known Probability Distribution. The true values of the uncertain parameters are obtained from their nominal values by random variables:

The deterministic robust optimization model for the original uncertain MILP problem can be derived as follows: min/max c Tx + dTy s.t. Ex + Fy = e

al̃ m = (1 + εξl m)al m ̃ = (1 + εξ )blk blk

Ax + By ≤ p

∑ al mxm + ∑ blkyk

lk

pl̃ = (1 + εξl)pl

≤ δmax[1, |pl |] x ≤ x ≤ x̅

Assume that the probability distributions of the random variables ξlm, ξlk, and ξl are known. Then, eq 17 becomes



al mx m +

m ∉ Ml





+



al m(1 + εξl m)x m +

yk = 0, 1

blk(1 + εξlk)yk ≤ pl (1 + εξl) + δmax[1, |pl |]

∀l

∀k

Unknown Probability Distribution. In this case, the righthand side parameters in eq 15 are uncertain and follow an unknown probability distribution, that is,

blkyk

k∉Kl

m ∈ Ml

k

m

(18)

− pl + ε·f (λ , |al m|x m , |blk |yk , pl )

∀l

∑ al mxm + ∑ blkyk

k∈Kl

(19)

∑ plĩ



k

m

+

i ∈ Il

∑ pli

+ δmax[1, |pl |]

∀l

i ∉ Il

(25)

Equation 19 is assumed to be violated by the user-specified parameter κ (0 ≤ κ ≤ 1). In other words, Pr{∑ al mx m +

∑ blkyk

+

− pl + ε(

k

m





where, pl = ∑ipli. Its probabilistic form with constraint violation of the userspecified parameter κ is represented by Pr{∑ al mx m +

al m ξl mx m

m ∈ Ml

blk ξlk yk − pl ξl) > δmax[1, |pl |]} ≤ κ

∀l

i ∈ Il

al m ξl mx m +



Fξ−1(1 − κ) = f (λ , |al m|x m , |blk |yk , |pl |) l



m

∑ blkyk

+

k

∑ pli

+ δmax[1, |pl |]}

i ∉ Il

−E[∑i ∈ I plĩ ] l

−∑m al mx m − ∑k blk yk + ∑i ∉ I pli + δmax[1, |pl |]

i ∈ Il

(27)

i ∈ Il

∑ pli ] i ∈ Il

(28)

Combine eqs 27 and 28, Pr{− ∑ plĩ > −∑ al mx m −

(23)

i ∈ Il



m

∑ blkyk k

+

∑ pli

+ δmax[1, |pl |]}

i ∉ Il

−max[0, 2κ·∑i ∈ I pli − ∑i ∈ I pli ] l

l

−∑m al mx m − ∑k blk yk + ∑i ∉ I pli + δmax[1, |pl |] l

If let

− pl + ε·f (λ , |al m|x m , |blk |yk , pl )

∀l

the left-hand side of eq 26 is

E[ ∑ plĩ ] ≥ max[0, 2κ· ∑ pli −

− max[0, 2κ·∑i ∈ I pli − ∑i ∈ I pli ] l l ≤κ −∑m al mx m − ∑k blk yk + ∑i ∉ I pli + δmax[1, |pl |] l

k

≤ δmax[1, |pl |]

∑ pli i ∉ Il

If the unknown distribution is approximately symmetric and the extremes are equidistant from the mean, for a given value of κ,

Combing eqs 21−23, the generic deterministic robust counterpart of the probabilistic constraint can be formulated as follows for random variables following any distribution.

∑ al mxm + ∑ blkyk

+

l

(22)

∀l

14

m

∀l

We define Fξl as the cumulative distribution for the aggregate random variable ξl and Fξl−1 to be inverse cumulative distribution. ∀l

∑ blkyk

∀l

Pr{− ∑ plĩ > −∑ al mx m −

(21)

Fξ (λ) = Pr{ξl ≤ λ} = 1 − κ l

i ∉ Il

k

m

From Markov Inequality,

k∈Kl

m ∈ Ml

∑ pli

(26)

i ∈ Il

blk ξlk yk − pl ξl

i ∈ Il

+ δmax[1, |pl |]} ≤ κ

Note that if the tolerance for the constraint violation is low, then the value of κ should be similarly small. We aggregate the respective random variables into one random variable.

+

∀l

⇒ Pr{− ∑ plĩ > −∑ al mx m −

(20)



k

m

∑ plĩ

>

+ δmax[1, |pl |]} ≤ κ

k∈Kl

ξl =

∑ blkyk

(24) 4353

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

For simplicity, universal values for ε, κ, and δ are adopted. However, Janak et al.22 highlighted that the robust optimization framework can easily be extended for cases where ε can vary from parameter to parameter and κ, and δ can vary from constraint to constraint. Next, we extend and apply the robust optimization framework to address the various forms of demand uncertainty present during the operations. 6.2. Robust Counterpart for Demand Uncertainty: Amount and Due Date. Verderame and Floudas 14 considered two forms of demand uncertainty for bulk products, that is, the realized date and the required amount associated with each bulk product demand due date parameter. While the actual date of demand realization were treated to follow a discrete, uniform distribution with bounds of plus or minus one day from the nominal day, the realized amount of demand followed a user-specified probability distribution. To model concurrently both the amount and day realization, they assigned a unique order designation, or, to each demand due parameter, which is illustrated in Figure 5. Therefore, the order

then the deterministic robust counterpart for eq 25 is represented by

∑ al mxm + ∑ blkyk



k

m

∑ pli

∑ pli i ∉ Il

+ δ·max[1, |pl |]

⎡ 2κ − 1 ⎤ + max⎢0, · ⎣ κ ⎦⎥

∀l

i ∈ Il

(29)

The deterministic robust counterpart optimization model for the original uncertain MILP problem can be derived as follows, min/max c Tx + dTy s.t. Ex + Fy = e Ax + By ≤ p

∑ al mxm + ∑ blkyk k

m

∑ pli



∑ pli i ∉ Il

+ δmax[1, |pl |]

⎡ 2κ − 1 ⎤ + max⎢0, · ⎣ κ ⎦⎥

∀l

i ∈ Il

x ≤ x ≤ x̅ yk = 0, 1 ∀ k

Bounded Uncertainty. Assume that the uncertain parameters vary in a bounded interval. They are represented by |al̃ m − al m| ≤ ε|al m| ̃ − blk | ≤ ε|blk | |blk

Figure 5. Bulk product order designation diagram from Verderame and Floudas (2010)14.

|pl̃ − pl | ≤ ε|pl |

(30)

has a 33% chance of being realized by the previous day [dem_or(s, d − 1, or)], a 66% chance of being realized by the nominal day [dem_or(s, d, or)], and a 100% chance of being realized by the following day [dem_or(s, d + 1, or)]. Each possible realization of the customer order has a realized amount, which follows the user-specified probability distribution given as follows,

The deterministic robust counterpart for eq 17 is derived as follows,

∑ al mxm + ∑ blkyk k

m



+ ε(

|al m|u m +

≤ pl − ε|pl | + δmax[1, |pl |]



|blk |yk )

k∈Kl

m ∈ Ml

∀l

(31)

where, −um ≤ xm ≤ um The deterministic robust counterpart of the original uncertain MILP problem can be derived as follows:

͠ _or(s , d − 1, or ) = [1 + ε·ξ(s , d − 1)]·dem(s , d) dem ͠ _or(s , d , or ) = [1 + ε·ξ(s , d)]·dem(s , d) dem

min/max c Tx + dTy x ,y,u

͠ _or(s , d + 1, or ) = [1 + ε·ξ(s , d + 1)]·dem(s , d) dem

s.t. Ex + Fy = e

where ε denotes an uncertainty level, and ξ(s,d) represents the random variable following a distribution. 6.2.1. Robust Demand Counterparts for MTO Articles. Following the approach of Verderame and Floudas,14 the daily demand constraints for the MTO articles (i.e., eqs 8−9) need to be changed into the following forms:

Ax + By ≤ p

∑ al mxm + ∑ blkyk k

m

+



+ ε(



|al m|u m

m ∈ Ml

|blk |yk ) ≤ pl − ε|pl | + δmax[1, |pl |]

∀l

k∈Kl

−um ≤ xm ≤ um x ≤ x ≤ x̅ yk = 0, 1

d

∀ m ∈ Ml

∑ d ′= 1

d

tot(s , d′) + sl1a _b(s , d) ≥

∑ ∑

͠ _or(s , d′, or ) (dem

or d ′= 1

͠ _or(s , d′ − 1, or )·[1 − L(s , d′ − 1, or )]) − dem

∀k

∀ s ∈ SMTO ,d p

(33)

(32) 4354

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

d

d



Pr{

tot(s , d′) − sl1b _b(s , d)

d

d

∑ ∑

>

͠ _or(s , d′, or ) − dem ͠ _or(s , d′ − 1, or ) · (dem

[1 − L(s , d′ − 1, or )])

,d ∀ s ∈ SMTO p

d

(34)

[1 − L(s , d′ − 1, or )])]}

The generic robust deterministic counterparts for eqs 33 and 34 can be formulated as follows for random variables following any distribution. d



d

͠ _or(s , d′, or ) (dem

tot(s , d′) + sl1a _b(s , d)

d ′= 1 d

or d ′= 1



͠ _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )]) − dem

∑ ∑

[1 − L(s , d′ − 1, or )]) + ε·f [λ1, dem _or(s , d′, or ), or ]

(dem _or(s , d′, or )

d

or d ′= 1

− δ·max[1,

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])] ,d ∀ s ∈ SMTO p

d ′= 1

,d ∀ s ∈ SMTO p

͠ _or(s , d′, or ) (dem

or d ′= 1



d

∑ ∑

(dem _or(s , d′, or )

d



,d ∀ s ∈ SMTO p

(36)

d

+ δ·max[1,

tot(s , d′) + sl1a _b(s , d)

∑ ∑

,d ∀ s ∈ SMTO p

͠ _or(s , d′, or ) − dem ͠ _or(s , d′ − 1, or ) · (dem

or d ′= 1 d

∑ ∑

f [λ1, dem _or(s , d′, or ), or ]

(dem _or(s , d′, or )

= Fn−1[

or d ′= 1

prob(s , d , or )]

prob(s , d , or )]

f [λ2, dem _or(s , d′, or ), or ]

or : prob(s , d , or ) > 0

∀ s ∈ SMTO ,d p

∏ or : prob(s , d , or ) > 0

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])]}



(40)

Normal Distribution. If the demand due parameters follow an independent normal distribution, then

[1 − L(s , d′ − 1, or )])

≤ [1 −

(dem _or(s , d′, or )

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])]

d

−δ·max[1,

∑ ∑ or d ′= 1

d ′= 1


0

tot(s , d′) − sl1a _b(s , d)

∑ ∑



≤ [1 −

d

+ δ·max[1,

(dem _or(s , d′, or ) − dem _or(s , d′ − 1, or )·

or d ′= 1

Now the robust optimization framework can be applied to eqs 33−34.



∑ ∑

+δ· max[1,

⎧1 if day d is the last day that the customer ⎪ L(s , d , or ) = ⎨ order or for state s could be realized ⎪ ⎩ 0 otherwise

≤ −∑

͠ _or(s , d′, or ) − dem ͠ _or(s , d′ − 1, or )· (dem

[1 − L(s , d′ − 1, or )])

where, the parameter L(s, d, or) is defined as follows,



∑ ∑ or d ′= 1

or d ′= 1



tot(s , d′) − sl1b _b(s , d)

d ′= 1

d ′= 1





= Fn−1[



prob(s , d , or )]

or : prob(s , d , or ) > 0

(37) 4355

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

Unknown Probability Distribution. If the summation of the demand due date parameters for a given MTO article is assumed to follow an unknown distribution, then the robust counterparts for eqs 33 and 34 are given by

where Fn−1[



prob(s , d , or )]

or : prob(s , d , or ) > 0

d



is the inverse standard normal cumulative distribution function and σ(d) is the time horizon dependent standard deviation of the random variable ζ(s, d) used to represent the uncertain parameter dem_or(s, d, or). The robust counterparts are represented by

tot(s , d′) + sl1a _b(s , d)

d ′= 1



1 1 − ∏or : prob(s , d , or ) > 0 prob(s , d , or )

∑ ∑ d



(dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

or d ′≤ d tot(s , d′) + sl1a _b(s , d)

[1 − L(s , d′ − 1, or )]) − δ·

d ′= 1

d

d



∑ ∑

{dem _or(s , d′, or ) − dem _or(s , d′ − 1, or )·

max[1,

or d ′= 1

∑ ∑

(dem _or(s , d′, or ) − dem _or

or d ′= 1

[1 − L(s , d′ − 1, or )]} + ε·λ(s , d)· d ⎛ ⎜ ∑ ∑ [σ(d′)·dem _or(s , d′, or )]2 ⎜ ⎝ or d ′= 1

(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])] ,d ∀ s ∈ SMTO p

{

(43)

⎞1/2

2 ⎟

}

d

− (σ(d′ − 1)·dem _or(s , d′ − 1, or )[1 − L(s , d′ − 1, or )]) ⎟ ⎠



d

− δmax[1,

∑ ∑

d ′= 1

(dem _or(s , d′, or ) − dem _or(s , d′ − 1, or )·

or d ′= 1

[1 − L(s , d′ − 1, or )])]

,d ∀ s ∈ SMTO p

tot(s , d′) + sl1b _b(s , d)

≤ (41)

2·(1 − ∏or : prob(s , d , or ) > 0 prob(s , d , or )) − 1 1 − ∏or : prob(s , d , or ) > 0 prob(s , d , or ) d

∑ ∑ d



tot(s , d′) + sl1b _b(s , d)

[1 − L(s , d′ − 1, or )]) + δ·

d ′= 1

d

d



(dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

or d ′= 1

∑ ∑

max[1,

{dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])]

[1 − L(s , d′ − 1, or )]} − ε·λ(s , d) · d ⎛ ⎜ ∑ ∑ [σ(d′) ·dem _or(s , d′, or )]2 ⎜ ⎝ or d ′= 1

,d ∀ s ∈ SMTO p

{

⎞1/2 [1 − L(s , d′ − 1, or )])2 ⎟ ⎟ ⎠

}

ξ(s , d′ − 1, or′) = α(s , d − 1, d′ − 1, or , or′) ·ξ(s , d − 1, or ) + β(s , d − 1, d′ − 1, or , or′)

d

(dem _or(s , d′, or )

ξ(s , d′, or′) = α(s , d , d′, or , or′) ·ξ(s , d , or ) + β(s , d , d′, or , or′)

or d ′= 1

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])] ,d ∀ s ∈ SMTO p

(44)

Correlated Demand Distribution. If two separate orders for the same state are linearly correlated with one another, then these two correlated orders can be written:

− (σ(d′ − 1) ·dem _or(s , d′ − 1, or )

∑ ∑

(dem _or(s , d′, or ) − dem _or

or d ′= 1

or d ′= 1

+ δmax[1,

∑ ∑

ξ(s , d′ + 1, or′) = α(s , d + 1, d′ + 1, or , or′) ·ξ(s , d + 1, or ) + β(s , d + 1, d′ + 1, or , or′)

(42)

The deterministic robust counterparts for other known probability distributions are presented in the Supporting Information section S1.

∀ or ∈ ORI , or′ ∈ ORD, match(or , or′) > 0, d , nominal(s , d , or ) > 0, d′, nominal(s , d′, or′) > 0 4356

(45)

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

d

where, sets ORI denotes independent orders and ORD denotes dependent orders, parameters α(s,d,d′,or,or′) is order correlation coefficient and β(s,d,d′,or,or′) is order correlation constant, parameter nominal(s,d,or) is used to denote whether order or for state s is satisfied on the nominal day d; Match(or,or′) identifies the presence of correlation between independent order or and dependent order or′. From Kenney and Keeping,48 each α and β can be determined by means of finding the minimum of a related least-squares function:



tot(s , d′) + sl1b _b(s , d)

d ′= 1 d



∑ ∑

(dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

or d ′= 1

[1 − L(s , d′ − 1, or )]) − ε·Fn−1[

prob(s , d , or )]· σ2

∏ or : prob(s , d , or ) > 0 d

{ξ(s , d′, or′) − [α(s , d , d′, or , or′) ·ξ(s , d , or ) + β(s , d , d′, or , or′)]}2

+ δ·max[1,

∑ ∑

(dem _or(s , d′, or )

or d ′= 1

(46)

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])] ,d ∀ s ∈ SMTO p

The solution of one of these three systems of two linear algebraic equations obtained through separate first-order partial differentiation of the least-squares function with respect to the given α coefficient and β constant is

α(s , d , d′, or , or′) =

where, σ2 = ⎧⎡ ⎪⎢ d ⎪⎢ ⎪ ∑ ∑ ⎨⎢dem_or(s , d′, or) + ∑ ⎢ or ′∈ ORD, or ∈ ORI d ′= 1 ⎪⎢ ⎪ ⎪⎢⎣ or ′ :match(or , or ′) > 0 ⎩

cov(s , d , d′, or , or′) σd2

β(s , d , d′, or , or′) = μ(s , d′, or′) − α(s , d , d′, or , or′) ·μ(s , d , or )

⎫2 ⎤ ⎪ ⎥ d ⎪ ⎪ cov(s , d′, d″ , or , or′)·dem _or(s , d″ , or′) ⎥⎥ × ∑ ·σ(d′)⎬ 2 ⎥ ⎪ σ ′ ( ) d d ′= 1 ⎥ ⎪ ⎪ ⎥⎦ ⎭

(47)

where, cov(s,d,d′,or,or′) is the covariance of the two random variables ζ(s,d′,or) and ζ(s,d,or), σ(d) is the standard deviation of the random variable ζ(s,d,or), μ(s,d,or) is the average of the random variable ζ(s,d,or). Note that the given β constant can be set to zero if the means of the random variables are equal to zero. If each of the independent and dependent random variables follows a normal distribution, and the independent and dependent random variables are linearly correlated, then the robust demand counterparts for underproduction and overproduction constraints are represented by,

⎧ ⎪ ⎪ ∑ − ⎨dem _or(s , d′ − 1, or ) + ⎪ or ′∈ ORD, ⎪ or ′ :match(or , or ′) > 0 ⎩ d

(dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

d

[1 − L(s , d′ − 1, or )]) + ε·



prob(s , d , or )]· σ2



(dem _or(s , d′, or )



dem(s , d′)

d ′= 1

d

or d ′= 1



− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )])] ,d ∀ s ∈ SMTO p

tot(s , d′) + sl1a _p(s , d) ≥

∀ s ∈ SMTS ∩ STR p p ,d

d

∑ ∑

d

d ′= 1

or : prob(s , d , or ) > 0

− δ·max[1,

·σ(d′ − 1)·

6.2.2. Robust Counterparts for Truck MTS Articles. In section 3, the following underproduction and overproduction constraints (i.e., eqs 50 and 51) are enforced for truck MTS articles, which are similar to those of MTO articles,

or d ′≤ d

Fn−1[

σ2(d′ − 1)

⎫2 ⎪ ⎪ ⎪ [1 − L(s , d′ − 1, or )]⎬ ⎪ ⎪ ⎪ ⎭

tot(s , d′) + sl1a _b(s , d)

∑ ∑

cov(s , d′ − 1, d″ , or , or′)·dem _or(s , d″ , or′)

d ″= d ′− 1

d ′= 1





×

d



(49)

d

tot(s , d′) − sl1b _p(s , d) ≤

d ′= 1

4357



dem(s , d′)

d ′= 1

∀ s ∈ SMTS ∩ STR p p ,d

(48)

(50)

(51)

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

Following the approach of Vederame and Floudas,14 eqs 50 and 51 change to

6.2.3. Robust Counterparts for Nontruck MTS Articles. For nontruck MTS articles, underproduction constraints are also

d



imposed, which are similar to that of MTO articles: tot(s , d′) + sl1a _p(s , d)

d ′= 1

d d



∑ ∑



͠ _or(s , d′, or ) − dem ͠ _or(s , d′ − 1, or ) · {dem

d ′= 1

or d ′= 1

[1 − L(s , d′ − 1, or )]}

∀ s ∈ SMTS ∩ STR p p ,d



dem(s , d′)

d ′= 1

∀ s ∈ SMTS ∩ SNTR ,d p p

(56)

(52)

Using the approach of Vederame and Floudas (2010),14 eq 56

d



d

tot(s , d′) + sl1a _p(s , d) ≥

tot(s , d′) − sl1b _p(s , d)

becomes

d ′= 1 d



∑ ∑

d

͠ _or(s , d′, or ) − dem ͠ _or(s , d′ − 1, or ) · {dem



or d ′= 1

[1 − L(s , d′ − 1, or )]}

∀ s ∈ SMTS ∩ STR p p ,d

tot(s , d′) + sl1a _p(s , d)

d ′= 1 d

(53)



Applying the robust optimization framework, the generic robust deterministic counterparts for eqs 52 and 53 can be generated as follows for random variables following any distribution.

͠ _or(s , d′, or ) − dem ͠ _or(s , d′ − 1, or ) · (dem

∑ ∑ or d ′= 1

[1 − L(s , d′ − 1, or )])

,d ∀ s ∈ SMTS ∩ SNTR p p (57)

d



tot(s , d′) + sl1a _p(s , d)

The generic robust deterministic counterpart of eq 57 from the robust optimization framework can be formulated as follows for random variables following any distribution.

d ′= 1 d



∑ ∑

{dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

or d ′= 1

[1 − L(s , d′ − 1, or )]} + ε·f [λ1, dem _or(s , d′, or ), or ]

d

d

− δ·max[1,

∑ ∑

∑ {dem _or(s , d′, or )

tot(s , d′) + sl1a _p(s , d) ≥

d ′= 1

or d ′= 1

∑ ∑

{dem _or(s , d′, or )

or d ′= 1

− dem _or(s , d′ − 1, or )· [1 − L(s , d′ − 1, or )]}+

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )]}] ∀ s ∈ SMTS ∩ STR p p ,d

ε·

f [λ1, dem _or(s , d′, or ), or ] − δ· d

(54)

max[1,

∑ ∑

{dem _or(s , d′, or ) − dem _or(s , d′ − 1, or )·

or d ′= 1

d



d

[1 − L(s , d′ − 1, or )]}]

tot(s , d′) + sl1b _p(s , d)

,d ∀ s ∈ SMTS ∩ SNTR p p

(58)

d ′= 1 d



∑ ∑

Supporting Information, section S3, presents the deterministic robust counterparts for some known probability distributions such as uniform and normal distributions. Although similar overproduction constraints are not proposed for nontruck MTS articles, an upper bound constraint as follows is imposed to allow for some overproduction of these articles.

{dem _or(s , d′, or ) − dem _or(s , d′ − 1, or ) ·

or d ′= 1

[1 − L(s , d′ − 1, or )]} − ε·f [λ2, dem _or(s , d′, or ), or ] d

+ δ·max[1,

∑ ∑

{dem _or(s , d′, or )

or d ′= 1

− dem _or(s , d′ − 1, or ) ·[1 − L(s , d′ − 1, or )]}] ∀ s ∈ SMTS ∩ STR p p ,d



(55)

tot(s , d) ≤

initial _w(w) ≤ d ≤ ter m_w(w)



The deterministic robust counterparts for some known probability distributions such as uniform and normal distributions are presented in Supporting Information, section S2.

initial _w(w) ≤ d ≤ ter m_w(w) ∀ s ∈ SMTS ∩ SNTR ,d p p 4358

6 · 5

dem(s , d)

(59)

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

Following the approach of Verderame and Floudas (2010),14 eq 59 becomes 5 6



The deterministic robust counterparts for some known probability distributions are presented in Supporting Information, section S3. The robust planning model (denoted as ROM) is defined below,

tot(s , d)

initial _w(w) ≤ d ≤ ter m_w(w)

(ROM) Min COST

͠ _or(s , d , or ) {dem





s.t. eqs 1−13

initial _w(w) ≤ d ≤ ter m_w(w)

͠ _or(s , d − 1, or ) ·[1 − L(s , d − 1, or )]· − dem ,d ∀ s ∈ SMTS ∩ SNTR p p

YW (s , w , d − 1)}

eqs 39−40, 43−44, and 48−49 for corresponding robust counterparts

(60)

where, YW(s,w,d) is defined by

eqs 54−55, and S23−S26 for corresponding robust counterparts

⎧1 if day d is in week w YW (s , w , d) = ⎨ ⎩ 0 otherwise

eqs 58, 61, S33−S34, and S37−S38 for corresponding robust counterparts

Using the robust optimization framework, the generic robust deterministic counterpart of the probabilistic constraint eq 60 can be formulated as follows for random variables following any distribution. 5 6



7. COMPUTATIONAL STUDY: PLANNING UNDER UNCERTAINTY The proposed deterministic robust counterpart optimization formulation ROM is applied to solving the large-scale industrial example in section 5. The uncertainty on demand amount is assumed to follow a normal distribution. The standard deviation [σ(d)] is calculated using the following equation:

tot(s , d)

initial _w(w) ≤ d ≤ ter m_w(w)



∑( or



{dem _or(s , d , or )

initial _w(w) ≤ d ≤ ter m_w(w)

σ(d) = 0.05 + 0.025·(p − 1)

− dem _or(s , d − 1, or ) ·[1 − L(s , d − 1, or )]·

∀ d , p , initial _p(p) ≤ d ≤ term _p(p)

YW (s , w , d − 1)}) − ε·f [λ , dem _or(s , d , or ), or ] + δ·max{1,

∑( or



where, p stands for planning period and initial_p(p) and term_p(p) are the first and last days within the given planning period, respectively. Each planning period has a time duration of approximately two weeks. The average for each utilized random variable is set to zero. The uncertainty level is ε = 100% which expresses the high degree of demand due date

dem _or

initial _w(w) ≤ d ≤ ter m_w(w)

(s , d , or ) − dem _or(s , d − 1, or ) ·[1 − L(s , d − 1, or )] ·YW (s , w , d − 1))}

∀ s ∈ SMTS ∩ SNTR ,d p p

(62)

(61)

Figure 6. Aggregate production totals of all articles from model ROM and the medium-term scheduling model. 4359

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

Figure 7. Operational schedule for bottleneck main units from model ROM with three-month time horizon.

parameter uncertainty; the infeasibility tolerance is δ = 0.01 to maintain a stringent level of constraint feasibility. We impose an upper bound of 14400 CPU seconds for the model ROM. The model ROM involves 31 155 binary variables, 89 054 continuous variables, and 137 999 constraints. The aggregate production levels from ROM for a three-month planning horizon are presented in Table 2. The production profiles generated from model ROM are supplied to the medium-term scheduling framework proposed by Shaik et al. (2009)30 and the aggregate production levels are obtained which are also presented in Table 2. From Table 2, the robust operational model produces a total of 16877.907 mu articles, whereas the total articles from the medium-term scheduling framework are 16571.278 mu. This demonstrates that the supplied production profile from the robust operational planning model provides a tight upper bound on the true production capacity of the plant. Compared to the aggregate production levels in Table 1, the aggregate production level from the planning model is reduced from 18818.513 to 16877.907 mu because of demand amount and due date uncertainty. Figure 6 shows the aggregate production totals of all articles at both planning and scheduling levels in each period. The aggregate production totals from the scheduling model deviate at most 11.7% from those of the robust planning model, which is a bit larger than the deviation of 9.9% between the deterministic planning model and medium-term scheduling model. Similar to the deterministic planning model, the robust planning model also provides a tight upper bound on the production capacity of the plant mainly because some amounts of overproduction of nontruck MTS articles are allowed. Figure 7 depicts the true operational production schedule from the medium-term scheduling framework. The utilization of each bottleneck unit is presented in Table 4. From Table 4, M5 reaches the highest utilization of 78.39% and M1 has the lowest utilization, which is about 27.22%. Compared to the deterministic case, the highest utilization is reduced from 80.92% to 78.39% and the lowest utilization is reduced from 35.61% to 27.22%. This is possibly because the total production profiles provided from the robust planning model is reduced. From Figure 7, M1−M10 have some idle times, which represent the transition times in those units. M1, M3, M7, M8, and M9 have more idle times than other units because many articles are not allowed to be processed in M1, M3, M7, M8, and M9, as shown in Table 4.

8. CONCLUSIONS In this paper, we successfully extended the operational PPDM framework proposed by Verderame and Floudas8 to develop deterministic operational planning model for a large-scale multiproduct continuous plant. The robust optimization framework developed by Lin et al.,21 Janak et al.,22 and Li et al.23 was also successfully utilized and applied to develop robust optimization models for demand due date parameters uncertainty. The computational studies show that the production profiles supplied from the deterministic and robust planning model provided tight upper bounds on the true production capacity of the plant. In the future, we will extend the robust optimization framework to address other forms of uncertainty such as processing rates variation and develop integration framework for planning and scheduling to reduce the gap between the planning and scheduling.



ASSOCIATED CONTENT

S Supporting Information *

Additional calculations as described in the text. This material is available free of charge via the Internet at http://pubs.acs.org.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +1-609-258-4595. Fax: +1-609-258-0211. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors gratefully acknowledge support from the National Science Foundation (CMMI-08856021). NOMENCLATURE Indices u = extruder s = states d, d′ = days w = weeks m = months p = periods Sets

Su = states that can be produced in a given main unit u Sp = states that are final articles SpMTO = states that are MTO final articles 4360

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

SpMTS = states that are MTS final articles SpTR = final articles that are delivered by trucks SpNTR = final articles that are not delivered by trucks Us = main units that can produce a given state s Um = main units

(10) Kallrath, J. Operational Planning and Scheduling in the Process Industry. OR Spectrum 2002, 24, 219−250. (11) Pochet, Y.; Wolsey, L. A. Production Operational Planning by Mixed Integer Programming; Springer: New York, 2006. (12) Verderame, P. M.; Elia, J. A.; Li, J.; Floudas, C. A. Planning and Scheduling under Uncertainty: A Review across Multiple Sections. Ind. Eng. Chem. Res. 2010, 49, 3993−4017. (13) Joung, J. Y.; Blau, G.; Pekny, J. F.; Reklaitis, G. V.; Eversdyk, D. A Simulation Based Optimization Approach to Supply Chain Management under Demand Uncertainty. Comput. Chem. Eng. 2004, 28, 2087−2106. (14) Verderame, P. M.; Floudas, C. A. Operational Planning of Large-Scale Industrial Batch Plants under Demand Due Date and Amount Uncertainty: I. Robust Optimization Framework. Ind. Eng. Chem. Res. 2010, 48, 7214−7231. (15) Ben-Tal, A.; Nemirovski, A. Robust Convex Optimization. Math. Oper. Res. 1998, 23, 769−805. (16) Ben-Tal, A.; Nemirovski, A. Robust Solutions of Uncertain Linear Programs. Oper. Res. Lett. 1999, 25, 1−13. (17) Ben-Tal, A.; Nemirovski, A. Robust Solutions of Linear Programming Problems Contaminated with Uncertain Data. Math. Prog. A. 2000, 88, 411−424. (18) Ben-Tal, A.; Goryashko, A.; Guslitzer, E.; Nemirovski, A. Adjustable Robust Solutions for Uncertain Linear Programs. Math. Prog. A. 2004, 99, 351−376. (19) Ghaoui, L. H. E. Robust Solutions to Least-Square Problems with Uncertain Data. SIAM J. Matrix Anal. Appl. 1997, 18, 1035−1064. (20) Ghaoui, L. H. E.; Oustry, F.; Lebret, H. Robust Solutions to Uncertain Semidefinite Programs. SIAM J. Opt. 1998, 9, 33−52. (21) Lin, X.; Janak, S. L.; Floudas, C. A. A New Robust Optimization Approach for Scheduling under Uncertainty: I. Bounded Uncertainty. Comput. Chem. Eng. 2004, 28, 1069−1085. (22) Janak, S. L.; Lin, X.; Floudas, C. A. A New Robust Optimization Approach for Scheduling II. Uncertainty with Known Probability Distribution. Comput. Chem. Eng. 2007, 31, 171−195. (23) Li, Z. K.; Ding, R.; Floudas, C. A. A Comparative Theoretical and Computational Study on Robust Counterpart Optimization: I. Robust Linear Optimization and Robust Mixed-Integer Linear Optimization. Ind. Eng. Chem. Res. 2011, 50, 10567−10603. (24) Bertsimas, D.; Sim, M. Robust Discrete Optimization and Network Flows. Math. Prog. B. 2003, 98, 49−71. (25) Bertsimas, D.; Sim, M. The Price of Robustness. Oper. Res. 2004, 52, 35−53. (26) Bertsimas, D.; Pachamanova, D.; Sim, M. Robust Linear Optimization under General Norms. Oper. Res. Lett. 2004, 32, 510− 516. (27) Verderame, P. M.; Floudas, C. A. Integration of Operational Planning and Medium-Term Scheduling for Large-Scale Industrial Batch Plants under Demand and Processing Time Uncertainty. Ind. Eng. Chem. Res. 2010, 49, 4948−4965. (28) Verderame, P. M.; Floudas, C. A. Multisite Planning under Demand and Transportation Uncertainty: Robust Optimization and Contitional Value at Risk Framework. Ind. Eng. Chem. Res. 2011, 50, 4959−4982. (29) Verderame, P. M.; Floudas, C. A. Operational Planning of Large-Scale Industrial Batch Plants under Demand Due Date and Amount Uncertainty: II. Contional Value-at-Risk Framework. Ind. Eng. Chem. Res. 2010, 49, 260−275. (30) Shaik, M. A.; Floudas, C. A.; Kallrath, J.; Pitz, H.-J. Production Scheduling of a Large-Scale Industrial Continuous Plant: Short-Term and Medium-Term Scheduling. Comput. Chem. Eng. 2009, 33, 670− 686. (31) Floudas, C. A.; Lin, X. Continuous-Time versus Discrete-Time Approaches for Scheduling of Chemical Processes: A Review. Comput. Chem. Eng. 2004, 8, 2109−2129. (32) Floudas, C. A.; Lin, X. Mixed-Integer Linear Programming in Processing Scheduling: Modeling, Algorithms, and Applications. Ann. Oper. Res. 2005, 139, 131−162.

Parameters

CleanTime(u) = maximum changeover time in main unit u T(u,d) = total available processing time in main unit u on day d price(s) = price for state s Rate(u,s) = processing rate of state s in main unit u dem(s,d) = demand for state s on day d (intermediate due date parameter) H = planning horizon Binary Variables

Y(u,s,d) = one if the production of state s is assigned to main unit u in day d Continuous Variables

Tu(u,d) = total time that extruder u is utilized in day d RL(u,s,d) = processing time of state s in main unit u during day d Slt(u) = processing time underutilization slack variable for main unit u tot1(u,s,d) = amount of state s that is produced on day d in main unit u tot(s,d) = amount of state s that is produced on day d sl1a_b(s,d) = daily underproduction demand slack variable for MTO articles sl1b_b(s,d) = daily overproduction demand slack variable for MTO articles sl1a_p(s,d) = daily underproduction demand slack variable for MTS articles sl1b_p(s,d) = daily overproduction demand slack variable for MTS articles z = total cost



REFERENCES

(1) Tempelmeier, H. In Supply Chain Planning with Advanced Planning Systems, The 3rd Aegean International Conference on Design and Analysis of Manufacturing Systems, Tinos Island, Greece, 19−22 May, 2001. (2) Wilkinson, S. J.; Shah, N.; Pantelides, C. C. Aggregate Modeling of Multipurpose Plant Operation. Comput. Chem. Eng. 1995, 19, S583−S588. (3) McDonald, C. M.; Karimi, I. A. Planning and Scheduling of Parallel Semicontinuous Processes. 1. Production Planning. Ind. Eng. Chem. Res. 1997, 36, 2691−2700. (4) Kallrath, J.; Timpe, C. H. Optimal Operational Planning in Large Multi-site Production Networks. Eur. J. Oper. Res. 2000, 126, 422−435. (5) Kallrath, J. Combined Strategic and Operational PlanningAn MILP Success Story in Chemical Industry. OR Spectrum 2002, 24, 315−341. (6) Kallrath, J. Solving Planning and Design Problems in the Process Industry Using Mixed Integer and Global Optimization. Ann. Oper. Res. 2005, 140, 339−373. (7) Erdirik-Dogan, M.; Grossmann, I. E. Operational Planning Models for Parallel Batch Reactors with Sequence-Dependent Changeovers. AIChE J. 2007, 53, 1298−2300. (8) 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. (9) Verderame, P. M.; Floudas, C. A. Operational Planning Framework for Multisite Production and Distribution Networks. Comput. Chem. Eng. 2009, 33 (5), 1036−1050. 4361

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362

Industrial & Engineering Chemistry Research

Article

(33) Ierapetritou, M. G.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling: 1. Multipurpose Batch Processes. Ind. Eng. Chem. Res. 1998, 37, 4341−4359. (34) Ierapetritou, M. G.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 2. Continuous and SemiContinuous Process. Ind. Eng. Chem. Res. 1998, 37, 4360−4374. (35) Ierapetritou, M. G.; Hene, T. S.; Floudas, C. A. Effective Continuous-Time Formulation for Shor-Term Scheduling. 3. Multiple Intermediate Due Dates. Ind. Eng. Chem. Res. 1998, 38, 3446−3461. (36) Lin, X.; Floudas, C. A. Design, Synthesis and Scheduling of Multipurpose Batch Plants via an Effective Continuous-Time Formulation. Comput. Chem. Eng. 2001, 25, 665−674. (37) Lin, X.; Floudas, C. A.; Modi, S.; Juhasz, N. M. ContinuousTime Optimization Approach for Medium-Range Production Scheduling of a Multiproduct Batch Plant. Ind. Eng. Chem. Res. 2002, 41, 3884−3906. (38) Lin, X.; Chajakis, E. D.; Floudas, C. A. Scheduling of Tanker Lightering via a Novel Continuous-Time Optimization Framework. Ind. Eng. Chem. Res. 2003, 42, 4441−4451. (39) Janak, S. L.; Lin, X.; Floudas, C. A. Enhanced Continuous-Time Unit-Specific Event-Based Formulation for Short-Term Scheduling of Multipurpose Batch Processes: Resource Constraints and Mixed Storage Policies. Ind. Eng. Chem. Res. 2004, 43, 2516−2533. (40) Janak, S. L.; Floudas, C. A.; Kallrath, J.; Vormbrock, N. Production Scheduling of a Large-Scale Industrial Batch Plant. I. Short-Term and Medium-Term Scheduling. Ind. Eng. Chem. Res. 2006, 45, 8234−8252. (41) Janak, S. L.; Floudas, C. A. Improving Unit-Specific Event-Based Continuous-Time Approaches for Batch Processes: Integrality Gap and Task Splitting. Comput. Chem. Eng. 2008, 32, 913−955. (42) Shaik, M. A.; Janak, S. L.; Floudas, C. A. Continuous-Time Models for Short-Term Scheduling of Multipurpose Batch Plants: A Comparative Study. Ind. Eng. Chem. Res. 2006, 45, 6190−6209. (43) Shaik, M. A.; Floudas, C. A. Improved Unit-Specific EventBased Model Continuous-Time Model for Short-Term Scheduling of Continuous Processes: Rigorous Treatment of Storage Requirement. Ind. Eng. Chem. Res. 2007, 46, 1764−1779. (44) Shaik, M. A.; Floudas, C. A. Unit-Specific Event-Based Continuous-Time Approach for Short-Term Scheduling of Batch Plants using RTN Framework. Comput. Chem. Eng. 2008, 32, 260− 274. (45) Shaik, M. A.; Floudas, C. A. Novel Unified Modeling Approach for Short-Term Scheduling. Ind. Eng. Chem. Res. 2009, 48, 2947−2964. (46) Li, J.; Floudas, C. A. Optimal Event Point Determination for Short-Term Scheduling of Multipurpose Batch Plants via Unit-Specific Event-Based Continuous-Time Approaches. Ind. Eng. Chem. Res. 2010, 49, 7446−7469. (47) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R., GAMS: A User’s Guide. GAMS: San Francisco, CA, 2003. (48) Kenney, J. F.; Keeping, E. S. Linear Regression and Correlation. In Mathematics of Statistics, 3rd ed.; Van Nostrand: Princeton, NJ, 1962; Part 1.

4362

dx.doi.org/10.1021/ie202670a | Ind. Eng. Chem. Res. 2012, 51, 4347−4362