4948
Ind. Eng. Chem. Res. 2010, 49, 4948–4965
Integration of Operational Planning and Medium-Term Scheduling for Large-Scale Industrial Batch Plants under Demand and Processing Time Uncertainty Peter M. Verderame and Christodoulos A. Floudas* Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544
Operational planning, medium-term scheduling, and short-term scheduling are closely related activities. Despite prior advances, the effective integration of planning and scheduling under demand and processing time uncertainty remains a challenging problem, especially for large-scale industrial processes. To address this issue, a novel framework for the integration of planning and scheduling under demand due date and amount uncertainty, as well as batch processing time uncertainty is proposed. Both the robust operational planning with production disaggregation model and the operational conditional value-at-risk planning with production disaggregation model have been successfully interfaced with an industrially validated medium-term scheduling model by means of a rolling horizon framework used in conjunction with a novel feedback loop. Demand uncertainty is taken into account at the planning level, processing time uncertainty is addressed within the selected scheduling model, and a novel feedback loop allows for the two-way interaction between the planning and scheduling levels. An industrial case study has been conducted for a large-scale multipurpose and multiproduct batch plant capable of producing hundreds of products over a time horizon of 3 months in order to demonstrate the viability of the proposed approach. 1. Introduction Planning and scheduling deal with the allocation of plant resources; however, the effective integration of planning and scheduling has proven to be difficult because of the disparate time scales of operation. For example, operational planning operates over a time horizon of several months and has the objective of determining the production targets and raw material requirements for the plant in question. Scheduling occurs over a shorter time horizon and deals in greater detail with the dayto-day operations of the plant such as the allocation of tasks to the appropriate plant resources like process units. Kallrath1 discussed in detail the various aspects of planning and scheduling, Pochet and Wolsey2 provided an in-depth discussion regarding the modeling and solving of planning and scheduling problems, and Grossmann3 and Maravelias and Sung4 presented an overview of the various frameworks for the integration of planning and scheduling. The scheduling level utilizes the production profile supplied by the planning level in order to schedule the plant, and unrealistic production targets can cause the scheduling level to allocate resources in a suboptimal manner, which can lead to lost profits. The ability of the planning model to provide a tight upper bound on the production capacity of the plant and the effectiveness of the framework used to integrate the planning and scheduling levels can greatly affect a company’s bottom-line as noted by Shobrys and White.5 There exist several integration frameworks which consider all system parameters to be deterministic in nature. For example, Papageorgiou and Pantelides6,7 proposed a decomposition scheme having an upper level planning model and a lower level scheduling model, which addresses the planning and scheduling of batch and semicontinuous plants. Subrahmanyam et al.8 developed an iterative framework for the planning and scheduling of a chemical plant based in part upon the premise that the aggregate planning model is a relaxation of the scheduling * To whom correspondence should be addressed. Tel: 1-609-2584595. E-mail:
[email protected].
model. Within a review paper by Bassett et al.,9 an integration strategy composed of an aggregate planning model, which determines the upper level production targets and the time horizon for each scheduling problem, was presented. Similarity analysis was used subsequently in order to group scheduling subproblems into families where representative problems were solved in order to determine the feasibility of the overall scheduling problem. Dimitriadis et al.10 developed two frameworks for the integration of planning and scheduling by means of a rolling horizon approach. Zhu and Majozi11 proposed a formulation for the planning and scheduling of a multipurpose plant. Munawar and Gudi12 developed a multilevel planning, scheduling, and rescheduling framework which utilized control theory in order to minimize the revision of aggregate production targets and thus avoid the violation of customer commitments. Dogan and Grossmann13,14 presented a decomposition based method for the simultaneous planning and scheduling of a multiproduct, single-stage plants. A fundamentally important problem within the field of integrated planning and scheduling is the consideration of system uncertainty. For example, demand and processing time uncertainty should be explicitly taken into account within a planning and scheduling framework in order to ensure the appropriate allocation of plant resources. Petkov and Maranas15 applied chance constraint programming within a multiperiod planning and scheduling framework for multiproduct batch plants with uncertain demand profiles due at the end of each period. Wu and Ierapetritou16 presented a hierarchical approach to planning and scheduling under uncertainty. The proposed approach modeled uncertainty by means of multistage programming where the overall horizon was discretized into time periods based upon the given system parameters’ degree of certainty which varied inversely with the time horizon. As the relative uncertainty increased, the number of scenarios, which were generated by sampling the underlying probability distribution of each uncertain parameter, increased. Despite the recent and continued dedication to the field of planning and scheduling, the integration
10.1021/ie901973e 2010 American Chemical Society Published on Web 04/12/2010
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
of planning and scheduling under uncertainty remains a challenging problem within the literature as noted in the reviews of Sahinidis17 and Li and Ierapetritou.18 The proposed approach toward the integration of planning and scheduling under demand due date and amount uncertainty, as well as processing time uncertainty entails extending the framework of Verderame and Floudas.19 Demand due date and demand amount uncertainty will be taken into account at the operational planning level by replacing the nominal operational planning with production disaggregation model (PPDM) developed by Verderame and Floudas19 with either the robust operational PPDM (see Verderame and Floudas),20 which was developed by means of the robust optimization framework originally presented by Lin et al.21 and Janak et al.22 and later extended by Verderame and Floudas,20 or the operational CVAR-PPDM (see Verderame and Floudas),23 which utilized in its formulation the conditional value-at-risk theory originally presented by Artzner et al.,24 Ogryczak,25 and Rockafellar and Uryasev.26,27 Processing time uncertainty will be taken into account at the scheduling level by applying the robust optimization framework developed by Lin et al.21 and Janak et al.22 to the medium-term scheduling model developed by Janak et al.28,29 The rolling horizon approach presented within the nominal integrated planning and scheduling framework of Verderame and Floudas19 will be adapted to the planning and scheduling uncertainty problem by implementing a modified feedback loop. It will be demonstrated that the proposed approach explicitly takes into account the aforementioned forms of system uncertainty and allows the planning level to be gradually refined so that the selected operational planning model provides a tight upper bound on the production capacity of the plant in question and avoids the misallocation of plant resources. The remainder of the paper will take the following form. First, a description of the problem will be presented followed by an exposition of how the various forms of system uncertainty will be taken into account at the planning and scheduling levels. The rolling horizon approach along with a feedback loop, which allows the scheduling level results to refine the given operational planning model, will be presented, and the proposed framework will be implemented within an industrial case study in order to demonstrate the viability of the proposed approach toward dealing with the integration of operational planning and scheduling under demand and processing time uncertainty. Finally, concluding remarks highlighting the fundamental contributions of the paper will be made. 2. Problem Description The industrial plant under investigation is a multipurpose and multiproduct batch plant having the capability of producing hundreds of different products. The plant under consideration has previously been studied in the work of Janak et al.28,29 and Verderame and Floudas.19 The given batch plant produces both made-to-order (bulk) and made-to-stock (packed) goods. For the made-to-order products, the customer supplies the plant with intermediate demand due dates, which specify on a given day the amount of a product which needs to be sent to market. The supplied intermediate demand due dates have uncertainty associated with both the requested amount of product, as well as the day of demand realization, and both of these uncertain characteristics need to be explicitly modeled at the operational planning level. For the made-to-stock goods, the customer supplies the plant with weekly demand totals which should be satisfied by the end of the given week, and there is uncertainty associated with the amount of demand requested at the end of
4949
Figure 1. Plant state-task-network (STN).
each week. The plant in question has many different types of operations, and the state-task-network (STN) for the plant can be seen in Figure 1 where the rectangles denote the various unit operations occurring within the plant and the circles denote the different states with F, I, and P standing for feed, intermediate, and product, respectively. An in-depth description of the various production recipes for the given classes of final states can be found in the work of Janak et al.28,29 The batch processing times for the various unit operations capable of being undertaken within the plant are considered to be stochastic in nature, and these uncertain processing times need to be taken into account at the scheduling level. Having been supplied the plant and product characteristics as outlined by Janak et al.,28,29 the customer-supplied demand distributions and the processing time distributions, the proposed integration framework should effectively interface the planning and scheduling levels while concurrently modeling in an explicit manner the aforementioned forms of system uncertainty. 3. Planning Level Uncertainty At the operational planning level, demand due date and demand amount uncertainty need to be taken into account explicitly so that the production profile supplied to the scheduling level is immune to uncertainty. To be concise and maintain the convention set forth within the work of Verderame and Floudas,20,23 it is assumed that the amount associated with each demand due date parameter for both the bulk and packed product classes follows a normal distribution. The normal distribution case has been chosen because it is perhaps the most likely distribution for the demand due date parameters’ amount to follow, since demand is typically affected by several stochastic elements and the sum of these stochastic elements should in the limit follow a normal distribution due to the Central Limit Theorem.15 It should be noted, however, that both the conditional value-at-risk and robust optimization frameworks can take into account several other types of probability distributions related to amount uncertainty as demonstrated within the work of Verderame and Floudas.20,23 The demand due date uncertainty related to each bulk product customer order is assumed to follow a uniform, discrete distribution with bounds of plus or minus one day from the nominal day while the due dates for packed product customer orders are considered to be deterministic in nature. Processing time uncertainty is not modeled at the planning level because the unit aggregation approach used in part to formulate the PPDM19 outlined within the Appendix A.1 utilizes the maximum historically observed processing time for each task playing a role in the production recipe for a final product when calculating the production time for the product’s aggregate
4950
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
task. As a result, the aggregate processing time already represents a worst-case scenario so as to prevent the nominal PPDM from overestimating the production capacity of the plant, and additional processing time uncertainty analysis does not need to be taken into account at the planning level. 3.1. Nomenclature. The following indices, sets, parameters, and variables are required when modeling demand uncertainty within the PPDM by means of robust optimization or conditional value-at-risk frameworks: Indices. s ) states d ) days w ) weeks m ) months p ) planning periods or ) orders sc ) scenarios Sets. D ) days in the overall planning horizon W ) weeks in the overall planning horizon M ) months in the overall planning horizon P ) planning periods Or ) customer orders Sc ) demand due date parameters scenarios Sp ) states which are final products Sbp ) states which are bulk final products Spp ) states which are packed final products Parameters. rs,d ) demand for state s on day d (intermediate demand due date parameter) r_ors,d,or ) demand for state s on day d related to order or (intermediate demand due date order parameter) σd ) time horizon dependent standard deviation probs,d,or ) probability that the order or for state s will be realized on day d r_scs,d,sc ) demand for state s on day for scenario sc (intermediate demand due date scenario parameter) lasts,d,or ) takes on a value of 1 if day d is the last day that the customer order or for state s could be realized prices ) price for state s initial_ww ) first day of week w term_ww ) last day of week w initial_mm ) first day of month m term_mm ) last day of month m intial_pp ) first day of planning period p term_pp ) last day of planning period p ε ) degree to which uncertain parameter is known κ1 ) packed product underproduction probabilistic level of constraint violation (robust optimization parameter) κ2 ) packed product overproduction probabilistic level of constraint violation (robust optimization parameter) δr ) allowable degree of constraint violation (robust optimization parameter) ω ) user-specified probability level (CVAR parameter) δc ) user-specified risk aversion parameter (CVAR parameter) Continuous Variables. tot(s, d) ) amount of state s produced on day d sl1a(s, d) ) daily underproduction demand slack variable for bulk products sl1b(s, d) ) daily overproduction demand slack variable for bulk products sl2a(s, w) ) weekly underproduction demand slack variable for packed products
ξbulk_under(w) ) acceptable bulk underproduction loss level for week w ξbulk_over(m) ) acceptable bulk overproduction loss level for month m ξpack_under(m) ) acceptable packed underproduction loss level for month m ξpack_over(m) ) acceptable packed overproduction loss level for month m zbulk_under(w, sc) ) bulk underproduction loss function variable for week w and scenario sc zbulk_over(m, sc) ) bulk overproduction loss function variable for month m and scenario sc zpack_under(m, sc) ) packed underproduction loss function variable for month m and scenario sc zpack_over(m, sc) ) packed overproduction loss function variable for month m and scenario sc. 3.2. Nominal Demand Constraints. 3.2.1. Bulk Products. For the nominal PPDM presented within Verderame and Floudas,19 the demand restrictions for the bulk products are modeled as daily underproduction and overproduction demand constraints as shown in eqs 1 and 2. The aforementioned constraints enforce that the amount of bulk state s produced on day d should be equivalent to the demand for the given state while still maintaining feasibility with the presence of demand slack variables which are penalized within the objective function. The bulk demand constraints allow for the carryover of outstanding demand, as well as inventory from one day to the next. d')d
∑
d')d
∑r
tot(s, d') + sl1a(s, d) g
d')1
∀s|s ∈ Sp ∩ Sbp, ∀d
s,d'
d')1
(1) d')d
∑
d')d
tot(s, d') - sl1b(s, d) e
d')1
∑r
∀s|s ∈ Sp ∩ Sbp, ∀d
s,d'
d')1
(2) 3.2.2. Packed Products. The demand requirements for the packed products are enforced as weekly underproduction and monthly overproduction demand constraints within the nominal PPDM. Equation 3 represents the weekly underproduction demand constraint enforcing that the amount of a packed product produced by the end of week w should be greater than or equal to the weekly demand for the given product while still allowing for the carryover of outstanding demand, as well as inventory. The weekly packed product underproduction demand constraint is not strictly enforced, and a weekly underproduction slack variable is present within eq 3 and penalized within the PPDM’s objective function. The monthly overproduction demand constraint for each packed product is strictly enforced, however, and within any month, the amount produced has to be less than or equal to 120% of the total demand for the given packed product within the month m (eq 4).
∑
tot(s, d) + sl2a(s, w) g
determ_ww
∑
rs,d
determ_ww
∀s|s ∈ Sp ∩ Spp, ∀w (3) where term_ww is the last day of week w. 5 · 6
∑
initial_mmedeterm_mm
tot(s, d) e
∑
rs,d
initial_mmedeterm_mm
∀s|s ∈ Sp ∩ Spp, ∀m (4)
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
where initial_mm and term_mm are the first and last days of month m, respectively. 3.3. Robust Optimization Framework. The robust optimization framework consists of formulating probabilistic constraints for those nominal demand constraints which contain uncertain parameters, and then those probabilistic constraints are converted into their deterministic equivalents, which are subsequently used to augment the proposed model ensuring that any obtained solution is not only feasible for the nominal set of system conditions but also robust in terms of the uncertainty associated with the system’s stochastic elements. At the operational planning level, the demand due date and demand amount associated with each bulk product demand due date parameter is considered to be uncertain while the demand amount associated with each packed product demand due date parameter is considered to be uncertain. For an in-depth explanation of how the robust optimization framework can be applied to the problem of demand due date and amount uncertainty, consult the work of Verderame and Floudas.20 A summary of the theoretical underpinnings of robust optimization also can be found in the Appendix A.2. 3.3.1. Bulk Products. To model concurrently the uncertainty associated with both the due date and amount for each bulk product demand due date parameter (rs,d), each demand due parameter is given a unique order designation, or. The order has a 33% chance of being realized by the previous day (r_ors,d-1,or), 66% chance of being realized by the nominal day (r_ors,d,or), and a 100% chance of being realized by the following day (r_ors,d+1,or) in accordance with the assumption that the due date follows a discrete, uniform distribution with bounds of plus or minus one day from the nominal. Each possible realization of the customer order has a realized amount which follows the user-specified probability distribution (eq 5) where ε stands for an uncertainty level and ξs,d represents the random variable having a distribution which might be state and/or time specific. r_ors,d-1,or ) (1 + εξs,d-1)rs,d r_ors,d,or ) (1 + εξs,d)rs,d r_ors,d+1,or ) (1 + εξs,d+1)rs,d
(5)
To take into account the possibility of a customer order being realized on the previous, nominal, or following day, the daily demand constraints for the bulk products (eqs 1 and 2) need to be changed into the form shown in eqs 6 and 7. The total demand for a bulk state s up to a given day d has now been modified to take into account the possibility of alternative order due date realization without double counting the same customer order. d')d
{
Pr
d')d
d')d
∑ tot(s, d') + sl1a(s, d) < ∑ ∑ [r˜_or
d')1
or
or
[ (∑ ∑
δr max 1,
or
r_ors,d'-1,or(1 - lasts,d'-1,or)] d')d
∀s|s ∈ Sp ∩
(6)
d')d
∑ tot(s, d') - sl1b(s, d) e ∑ ∑ [r_or
s,d',or
or
-
d')1
r_ors,d'-1,or(1 - lasts,d'-1,or)]
∀s|s ∈ Sp ∩ Sbp, ∀d (7)
where the parameter lasts,d,or takes on a value of 1 if day d is the last day that the customer order or for state s could be realized. Equations 8 and 9 represent the probabilistic underproduction and overproduction demand constraints, respectively. It is important to note that the term κ which represents the
)]}
e
(
)]}
e
r_ors,d',or - r_ors,d'-1,or(1 - lasts,d'-1,or)
d')1
∏
1-
)
probs,d,or
or|probs,d,or>0
{∑ [ (∑ ∑ d')d
d')1
∀s|s ∈ Sp ∩ Spb, ∀d (8)
d')d
tot(s, d') - sl1b(s, d) >
Pr
∑ ∑ [r˜_or or
s,d',or
-
d')1
r˜_ors,d'-1,or(1 - lasts,d'-1,or)] +
d')d
δrmax 1,
or
(
r_ors,d',or - r_ors,d'-1,or(1 - lasts,d'-1,or)
d')1
∏
1-
)
probs,d,or
or|probs,d,or>0
∀s|s ∈ Sp ∩ Spb, ∀d (9)
Using the robust optimization framework originally developed by Lin et al.21 and Janak et al.,22 the generic robust bulk product underproduction and overproduction demand constraints can be written as seen in eqs 10 and 11, respectively. The function f depends upon the nominal values of the uncertain right-hand side parameters, as well as the underlying probability distribution that these aforementioned stochastic elements follow. d')d
d')d
∑ tot(s, d') + sl1a(s, d) g ∑ ∑ [r_or
s,d',or
d')1
or
-
d')1
r_ors,d'-1,or(1 - lasts,d'-1,or)] + εf(r˜_ors,d',or) -
[ (∑ ∑
d')d
δrmax 1,
or
r_ors,d',or - r_ors,d'-1,or(1 -
d')1
d')d
Sbp, ∀d
-
r˜_ors,d'-1,or(1 - lasts,d'-1,or)] -
-
d')1
s,d',or
d')1
d')d
lasts,d'-1,or) s,d',or
d')1
probabilistic level of constraint violation in the general formulation described in the work of Janak et al.22 has been replaced with the term (1 - ∏or|probs,d,or>0probs, d,or), and in doing so, the level of constraint reliability is now directly related to the likelihood of the sum of the pertinent orders being realized by the day in question, where probs,d,or is the probability of an individual order or for state s being realized by day d. The probabilistic constraint reliability level now explicitly takes into account the discrete, uniform distribution of the demand due date, and hence, the robust optimization framework allows for the coupling of the uncertain demand due date with the uncertain demand amount within standard probabilistic constraints which can in turn be reformulated into their deterministic robust counterparts.
d')d
∑ tot(s, d') + sl1a(s, d) g ∑ ∑ [r_or
d')1
4951
)]
∀s|s ∈ Sp ∩ Sbp, ∀d (10)
d')d
∑ tot(s, d') - sl1b(s, d) e ∑ ∑ [r_or
s,d',or
d')1
or
-
d')1
r_ors,d'-1,or(1 - lasts,d'-1,or)] - εf(r˜_ors,d',or) +
[ (∑ ∑
d')d
δr max 1,
or
r_ors,d',or - r_ors,d'-1,or(1 -
d')1
lasts,d'-1,or)
)]
∀s|s ∈ Sp ∩ Sbp, ∀d (11)
If each bulk demand due date parameter’s amount follows a normal distribution, then the sum of the aforementioned demand due date parameters also follows a normal distribution. Therefore, the generic function f can be replaced as shown in eqs 12 and 13 which represent the deterministic robust counterpart
4952
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
constraints for the bulk overproduction and underproduction demand constraints, respectively. The linearity of the robust demand constraints is preserved because the terms within the square root operator are deterministic parameters which can be ascertained before model execution. d')d
due parameter, which can be expressed as shown in eq 17 in terms of a random variable ξs,d with a given distribution, the user-specified parameter ε denoting the degree of parameter uncertainty, and the nominal demand due parameter rs,d as shown in eq 17. r˜s,d ) (1 + εξs,d)rs,d
d')d
∑ tot(s, d') + sl1a(s, d) g ∑ ∑ [r_or
s,d',or
d')1
or
(∑ ∑
-
d')1
r_ors,d'-1,or(1 - lasts,d'-1,or)] +
d')d
ελ
or
[(σd')(r_ors,d',or)]2 - [(σd'-1)(r_ors,d'-1,or)(1 -
)
d')1
[( )]
1/2
lasts,d'-1,or)]2
- δrmax 1,
r_ors,d'-1,or(1 - lasts,d'-1,or) d')d
d')d
∑ ∑ r_or or
s,d',or
-
d')1
∀s|s ∈ Sp ∩ Sbp, ∀d (12)
Given the probabilistic constraints shown in eqs 15 and 16, the generic packed product robust demand constraints can be written as seen in eqs 18 and 19 where the function f depends upon the nominal values of the uncertain demand parameters, as well as the underlying distribution of these uncertain parameters.
∑
determ_ww
tot(s, d) + sl2a(s, w) g
[
s,d',or
d')1
or
ελ
5 · 6
d')1
r_ors,d'-1,or(1 - lasts,d'-1,or)] -
(
d')d
∑ ∑ [(σ )(r_or d')1
2
s,d',or)]
d'
or
)
lasts,d'-1,or)]2
or
r_ors,d'-1,or(1 - lasts,d'-1,or)
r_ors,d',or -
∀s|s ∈ Sp ∩ Sbp, ∀d (13)
where σd is the time horizon dependent standard deviation defined by eq 14 of the random variable ξs,d used to represent the uncertain nature of the parameter r_ors,d,or and λ ) -1 F-1 n (∏or|probs,d,or>0probs,d,or) with Fn being the inverse standard normal cumulative distribution function.
{
Pr
{
Pr
∑
tot(s, d) + sl2a(s, w)
initial_mmedeterm_mm
δr max 1,
∑
determ_ww
rs,d
]}
r˜s,d +
initial_mmedeterm_mm
e κ2
∀s|s ∈ Sp ∩ Spp, ∀m
(16)
where κ1 and κ2 represent the probabilistic levels of constraint violation for the packed product weekly underproduction and monthly overproduction demand constraints, respectively, and r˜s,d represents the true realization of the packed product demand
∀s|s ∈ Sp ∩ Spp, ∀w (18)
∑
rs,d -
initial_mmedeterm_mm r
∑
]
rs,d
initial_mmedeterm_mm
∀s|s ∈ Sp ∩ Spp, ∀m (19)
If each packed product demand due date parameter’s amount follows a normal distribution, then eqs 20 and 21 represent the weekly and monthly deterministic robust counterpart demand constraints, respectively.
∑
determ_ww
rs,d +
determ_ww
ελ1
[
δr max 1, 5 · 6
∑
tot(s, d) + sl2a(s, w) g
σd ) 0.05 + 0.025(p - 1)∀d, ∀p|initial_pp e d e term_pp (14) where p stands for planning period and initial_pp and term_pp are the first and last days within the considered planning period, respectively. 3.3.2. Packed Products. The packed products have weekly demand requirements which are represented by weekly underproduction and monthly overproduction demand constraints. Since the due date is considered to be deterministic for all packed products, the probabilistic weekly and monthly demand constraints can be written in a form as shown in eqs 15 and 16, respectively.
]
tot(s, d) e
1,
d')1
rs,d + εf(r˜s,d) -
εf(r˜s,d) + δ max
[
d')d
+ δr max 1,
rs,d
determ_ww
∑
∑
determ_ww
initial_mmedeterm_mm
- [(σd'-1)(r_ors,d'-1,or)(1 -
[ (∑ ∑ )]
1/2
∑
δr max 1,
d')d
∑ tot(s, d') - sl1b(s, d) e ∑ ∑ [r_or
(17)
∑
∑
determ_ww
∑
]
∀s|s ∈ Sp ∩ Spp, ∀w (20)
∑
tot(s, d) e
ελ2
[
rs,d
determ_ww
initial_mmedeterm_mm
δr max 1,
(σdrs,d)2 -
rs,d -
initial_mmedeterm_mm
∑
(σdrs,d)2 +
initial_mmedeterm_mm
∑
initial_mmedeterm_mm
rs,d
]
∀s|s ∈ Sp ∩ Spp, ∀m (21)
where λ1 and λ2 can be determined through F-1 n (1 - κ1) and -1 F-1 n (1 - κ2), respectively, with Fn being the inverse standard normal cumulative distribution function. 3.4. Conditional Value-at-Risk Framework. Following the recent work of Verderame and Floudas,23 a scenario-based conditional value-at-risk approach also has been applied to the problem of demand due date and amount uncertainty at the operational planning level. Loss functions representing the underproduction and overproduction demand constraints are formulated and then constrained by means of encapsulated conditional value-at-risk constraints. The suggested approach has the ability to adapt to any demand due date and demand amount distributions since the forms of the distributions only affect the generation of the demand due date parameter scenarios. When generating demand due date parameter scenarios for a given bulk product, both the demand due date and demand amount distributions are sampled concurrently. For packed products, the due date is considered to be deterministic,
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
and so, demand due date parameter scenarios are generated by means of sampling only the related amount distribution. Having described how the demand due date parameter scenarios should be generated, the following sections now will present the development of both the bulk and packed CVAR demand constraints which replace the nominal PPDM demand constraints in order to formulate the operational CVAR-PPDM presented previously by Verderame and Floudas.23 For an outline of scenario-based conditional value-at-risk theory, consult the Appendix A.3, as well as the work of Verderame and Floudas.23 3.4.1. Bulk Products. To maintain computational tractability, both the underproduction and overproduction loss functions represent an aggregation of the nominal daily underproduction and overproduction demand constraints present within the deterministic PPDM (see Appendix A.1) and shown in eqs 1 and 2, respectively. As opposed to writing underproduction and overproduction loss functions for each bulk product on a daily basis, the underproduction and overproduction loss functions shown in eqs 22 and 23, respectively, are written on a weekly and monthly basis, respectively, for the entire bulk product class. The state specific terms found within the bulk product underproduction loss function (eq 22) are weighted by the respective price of the bulk product in order to emphasize the importance of meeting customer demand for profitable products. The state specific terms found within the bulk product overproduction loss function (eq 23) are not weighted by their respective prices because the inventory cost for each bulk product is the same regardless of how expensive the product is on the open market. The underproduction loss function for a given week w and scenario sc takes on a positive value when the total weighted production of bulk product is less than the total weighted demand for the scenario, and the overproduction loss function for a month m and scenario sc takes on a positive value when the total production is greater than the total demand for the scenario.
∑
fbulk_under(w, sc) )
∑
prices · (r_scs,d,sc -
s∈Sp∩Spb initial_wwedeterm_ww
tot(s, d))
∑
fbulk_over(m, sc) )
∑
∀w, ∀sc (22)
(tot(s, d) -
s∈Sp∩Spb initial_mmedeterm_mm
r_scs,d,sc)
∀m, ∀sc (23)
where fbulk_under(w, sc) and fbulk_over(m, sc) are the bulk product underproduction and overproduction loss functions for a given scenario sc, respectively. The bulk underproduction and overproduction loss functions are used in conjunction with the bulk underproduction and overproduction acceptable loss variables (ξbulk_under(w) and ξbulk_over(m), respectively) to define the loss function underproduction and overproduction variables (zbulk_under(w, sc) and zbulk_over(m, sc), respectively) as shown in eqs 24 and 25, respectively. If zbulk_under(w, sc) takes on a positive value, then the plant underproduced beyond the acceptable threshold for a given week w and scenario sc, and similarly, if zbulk_over(m, sc) is positive for a given month m and scenario sc, then the plant overproduced beyond the acceptable threshold for the given time period and scenario. zbulk_under(w, sc) g
∑
∑
prices · (r_scs,d,sc - tot(s, d)) -
s∈Sp∩Spb initial_wwedeterm_ww
ξbulk_under(w)
∀w, ∀sc (24)
∑
zbulk_over(m, sc) g
∑
4953
(tot(s, d) -
s∈Sp∩Spb initial_mmedeterm_mm
r_scs,d,sc) - ξbulk_over(m)
∀m, ∀sc (25)
Having defined the requisite loss function variables (zbulk_under(w, sc) and zbulk_over(m, sc)), the encapsulated CVAR constraints for the underproduction and overproduction demand constraints can be written in a form as shown in eqs 26 and 27, respectively.
∑ξ
bulk_under(w)
+
w
1 · (1 - ω)|Sc| δc
(w, sc) e ∑ ∑z · ∑ ∑ price · r (26) bulk_under
w
sc
s
s∈Sp∩Spb
∑ξ
bulk_over(m)
+
m
1 · (1 - ω)|Sc|
s,d
d
∑ ∑z
bulk_over(m, sc)
m
sc
δc ·
∑ ∑r
e
s,d
s∈Sp∩Spb
(27)
d
The left-hand sides of the above CVAR constraints represent the scenario-based approximation of the conditional value-atrisk for the particular loss function, and the right-hand sides are the acceptable loss thresholds as specified by the user which are partially dependent upon the historical demand distribution, the price of each of the products under investigation, and the user-specified parameter δc. This user specified parameter can take on a value between 0 and 1 and is used to express the user’s level of risk aversion. The parameter ω represents a probabilistic confidence level with bounds between 0 and 1. Equations 24-27 replace the nominal daily overproduction and underproduction demand constraints (eqs 1 and 2) within the operational PPDM so as to take into account bulk product demand uncertainty by means of conditional value-at-risk theory. 3.4.2. Packed Products. Having presented the bulk product underproduction and overproduction CVAR demand constraints, the packed product demand constraints (eqs 3 and 4) will be reformulated into their respective overproduction and underproduction CVAR demand constraints, and the packed product monthly underproduction and overproduction loss functions are defined by eqs 28 and 29, respectively. fpack_under(m, sc) )
∑
∑
prices · (r_scs,d,sc -
s∈Sp∩Spp initial_mmedeterm_mm
tot(s, d)) fpack_over(m, sc) )
∑
∑
∀m, ∀sc (28)
(tot(s, d) -
s∈Sp∩Spp initial_mmedeterm_mm
r_scs,d,sc)
∀m, ∀sc (29)
where fpack_under(m, sc) and fpack_over(m, sc) are the packed product underproduction and overproduction loss functions for a given scenario sc, respectively. The packed product loss functions represent an aggregation of the original packed product demand constraints on both the time and state levels in order to maintain computational tractability, and the terms comprising the two loss functions are weighted by their respective packed product price in accordance with the same standards used to generate the bulk product loss functions. The packed product underproduction and overproduction loss functions (fpack_under(m, sc) and fpack_over(m, sc), respectively) can be used in tandem with the packed product underproduction and overproduction acceptable
4954
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
loss variables (ξpack_under(m) and ξpack_over(m), respectively) in order to provide a lower bound on the packed product underproduction and overproduction loss function variables (zpack_under(m, sc) and zpack_over(m, sc), respectively) as shown in eqs 30 and 31. zpack_under(m, sc) g
∑
∑
prices · (r_scs,d,sc - tot(s, d)) -
s∈Sp∩Spp initial_mmedeterm_mm
ξpack_under(m) zpack_over(m, sc) g
∑
∑
∀m, ∀sc (30)
(tot(s, d) -
s∈Sp∩Spp initial_mmedeterm_mm
r_scs,d,sc) - ξpack_over(m)
∀m, ∀sc (31)
The underproduction and overproduction loss function variables are further constrained by the encapsulated underproduction and overproduction CVAR demand constraints shown in eqs 32 and 33, respectively.
∑ξ
pack_under(m)
+
m
∑ ∑z
1 · (1 - ω)|Sc|
pack_under(m, sc)
m
∑ ∑ price
δc ·
pack_over(m)
m
+
1 · (1 - ω)|Sc|
· rs,d (32)
s
s∈Sp∩Spp
∑ξ
d
∑ ∑z
pack_over(m, sc)
m
e
sc
δ ·
s,d,sc∀s|s
∈ Sp
(34)
∑ tot(s, d) g ∑ min r_sc
s,d,sc∀s|s
∈ Sp
(35)
d
d
d
d
sc
sc
∑ tot(s, d) g 54 · ∑ r d
s,d∀s|s
∈ Sp
(36)
d
3.4.4. Sample Average Approximation. One proposed approach toward modeling explicitly demand uncertainty at the operational planning level utilizes scenario-based conditional value-at-risk theory in order to maintain linearity. The issue of adequately sampling the given probability space and ensuring that the obtained solution is feasible when considering the complete uncertainty space for the bulk and packed customer orders arises as a result. To maintain computational tractability, only a finite number of scenarios can be considered within the operational CVAR-PPDM, and therefore, an iterative approach needs to be adopted in order to guarantee that the final production profile satisfies the conditional value-at-risk demand constraints when considering the entire uncertainty space. The sample average approximation (SAA) developed by Wang and Ahmed30 provides such an iterative framework for determining a valid, feasible upper bound on the original minimimization problem when considering the entire uncertainty space, and the reader is directed to the work of Verderame and Floudas23 for a detailed explanation of how the operational CVAR-PPDM can be implemented within the SAA algorithm.
e 4. Scheduling Level Uncertainty
sc
c
∑ tot(s, d) e ∑ max r_sc
∑ ∑r
s∈Sp∩Spp
s,d
(33)
d
The scenario-based approximation of the conditional valueat-risk represented by the left-hand side of the underproduction and overproduction CVAR demand constraints (eqs 32 and 33, respectively) should be less than or equal to the maximum allowable risk the user is willing to be exposed to which is expressed as a fraction of the potential packed profit for the underproduction case and a fraction of the total packed demand for the overproduction case. Equations 30-33 replace eqs 3-4 so that the operational PPDM explicitly takes into account packed product demand uncertainty by means of conditional value-at-risk theory. Overall, eqs 24-27 and 30-33 have replaced the nominal demand constraints represented by eqs 1-2 and 3-4 in order to formulate the operational CVAR-PPDM. Also, the objective function for the operational CVAR-PPDM now only contains the reactor underutilization penalization and gross profit maximization terms due to the elimination of all the demand slack variables present within the nominal operational PPDM outlined within the Appendix A.1. 3.4.3. Additional Production Constraints. Because of computational tractability issues, the CVAR demand constraints have been written in an aggregated time basis and in terms of bulk and packed product classes as opposed to individual states. As a result, additional production constraints need to be added to the operational CVAR-PPDM in order to ensure an even production of all products having customer orders within the time horizon. Equations 34 and 35 bound the production of each final product based upon the minimum and maximum demand due date parameter scenarios for the state s on day d. Equation 36 mandates that at least 80% of the total demand for a given final product s be satisfied within the time horizon.
The medium-term scheduling model developed by Janak et al.28 decomposes the scheduling horizon into subhorizons which are then solved within a rolling horizon framework using a continuous-time unit-specific event-based model based upon the work of Floudas and co-workers.31-41 The rolling horizon framework involves successively solving each scheduling subhorizon and carrying over any unsatisfied demand to the following subhorizon. When provided the plant’s characteristics and the daily production profile supplied by the planning level, the medium-term scheduling model can determine the appropriate allocation of resources. However, the various forms of uncertainty need to be taken into account in order to ensure that the final schedule is consistent. Demand due date and demand amount uncertainty is explicitly taken into account at the operational planning level by means of the proposed robust optimization framework or the conditional value-at-risk theory ensuring that the production profile supplied to the scheduling level is immune to the aforementioned forms of demand uncertainty. As a result, only processing time uncertainty needs to be taken into account at the scheduling level. The robust optimization technique originally developed by Lin et al.21 and Janak et al.22 will be applied to task duration constraints present within the medium-term scheduling model developed by Janak et al.,28 and in doing so, the selected medium-term scheduling model will generate a schedule that is not only feasible for the nominal system conditions but also robust with regard to processing time uncertainty. 4.1. Nomenclature. The following indices, sets, parameters, and variables are required when modeling processing time uncertainty at the medium-term scheduling level. Indices. i ) tasks j ) units
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
n ) event points Sets. I ) tasks capable of being undertaken within the plant J ) units present within the plant N ) events points Parameters. Ri,j ) constant batch processing time for task i in unit j βi,j ) batch size proportional processing time for task i in unit j ε ) degree of parameter uncertainty κ ) probabilistic level of constraint violation δs ) allowable degree of robust task duration constraint violation Continuous Variables. B(i, j, n) ) batch size for task i within unit j at event point n Ts(i, j, n) ) start time of task i within unit j at event point n Tf(i, j, n) ) finish time of task i within unit j at event point n Binary Variables. wV(i, j, n) ) unit allocation variable denoting activation of task i within unit j at event point n. 4.2. Processing Time. Equation 37 represents the nominal task duration constraint present within the formulation of Janak et al.,28 which indicates that the finishing time of a task i taking place within unit j at event point n should be equal to the start time plus the requisite processing time for the given activated batch task. Following the work of Janak et al.,22 eq 37 is relaxed into a form as shown in eq 38 so that the robust optimization framework can be applied. T f(i, j, n) ) T s(i, j, n) + Ri,j · wV(i, j, n) + βi,j,n · B(i, j, n) ∀i, ∀j|(Ri,j + βi,j) > 0, ∀n (37) T f(i, j, n) g T s(i, j, n) + Ri,j · wV(i, j, n) + βi,j,n · B(i, j, n) ∀i, ∀j|(Ri,j + βi,j) > 0, ∀n (38) It should be noted that every considered batch processing task falls into one of two possible categories where either the batch processing time is constant regardless of the respective batch size, which is denoted by a nonzero Ri,j, or is proportional to the batch size of the given task, which is denoted by a nonzero βi,j. As a result, two probabilistic constraints can be written in order to capture the uncertainty related to the processing time parameters Ri,j and βi,j as seen in eqs 39 and 40. Pr{T f(i, j, n) < T s(i, j, n) + R˜ i,j · wV(i, j, n) δs max[1, Ri,j]} e κ ∀i, ∀j|Ri,j > 0, ∀n (39) Pr{T f(i, j, n) < T s(i, j, n) + β˜ i,j · B(i, j, n) δs max[1, βi,j]} e κ ∀i, ∀j|βi,j > 0, ∀n (40) where the entities R˜ i, j and β˜ i, j, can be expressed in terms of their nominal values, the parameter ε representing a parameter’s degree of uncertainty, and the random variables ξi,R j and ξi,β j as shown in eqs 41 and 42. R˜ i,j ) (1 + εξRi,j)Ri,j
(41)
β˜ i,j ) (1 + εξβi,j)βi,j
(42)
The generic robust deterministic counterparts of eqs 39 and 40 can be seen in eqs 43 and 44, respectively, with the function
4955
f being dependent upon the distribution that the given batch processing time parameter obeys and the probabilistic level of constraint violation, κ. T f(i, j, n) g T s(i, j, n) + (Ri,j + εf(R˜ i,j)) · wV(i, j, n) δs max[1, Ri,j] ∀i, ∀j|Ri,j > 0, ∀n (43) T f(i, j, n) g T s(i, j, n) + (βi,j + εf(β˜ i,j)) · B(i, j, n) δs max[1, βi,j] ∀i, ∀j|βi,j > 0, ∀n (44) Making the assumption that all pertinent batch processing time parameters each follow a normal distribution, then the generic deterministic robust counterpart constraints can be rewritten into the form as shown in eqs 45 and 46. Tf(i, j, n) g Ts(i, j, n) + Ri,j · (1 + ελ) · wV(i, j, n) δs max[1, Ri,j] ∀i, ∀j|Ri,j > 0, ∀n (45) T f(i, j, n) g T s(i, j, n) + βi,j · (1 + ελ) · B(i, j, n) δs max[1, βi,j] ∀i, ∀j|βi,j > 0, ∀n (46) -1 where λ ) F-1 n (1 - κ) with Fn being the inverse standard normal cumulative distribution function.
5. Rolling Horizon Framework The integration of operational planning and scheduling under processing time, demand due date, and demand amount uncertainty requires the interaction of models which operate on disparate time scales. Both the robust operational PPDM and the operational CVAR-PPDM operate over a longer time horizon when compared to the medium-term scheduling model developed by Janak et al.,28 and because of model complexity limits, the operational planning model cannot represent the plant’s production capacity as rigorously as the selected scheduling model. An effective framework for integrating operational planning and scheduling under uncertainty should allow for the feedback of scheduling level production results in the form of additional operational planning model constraints. These additional constraints should reduce the feasible space of the operational planning problem and lead to the operational planning model providing a tighter upper bound on the plant’s production capacity. As the scheduling of the time horizon progresses, the operational planning model should provide the scheduling model with daily production targets which are not only theoretically immune to demand uncertainty but also more accurate with regard to reflecting the plant’s true production capacity. The proposed forward rolling horizon approach addresses this requirement by allowing for the two-way interaction between the operational planning and scheduling levels. As stated by Verderame and Floudas,19 the given rolling horizon approach minimizes the number of required iterations in order to schedule the entire time horizon when compared to traditional iterative frameworks, because, once a portion of the time horizon has been scheduled, no further modifications are undertaken for the given period of the time horizon. As a result, the final schedule and production profile obtained from the rolling horizon approach may not be the global optimum; however, it does represent a fair compromise between global optimality and computational effort. The integration of the robust operational PPDM20 or the operational CVAR-PPDM23 with the medium-term scheduling
4956
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
formulation of Janak et al.28 using the forward rolling horizon approach requires the addition of the following parameters and variables to the robust operational PPDM or the operational CVAR-PPDM: Parameters count_p ) order of the planning period following the last scheduled planning period Ups, p ) aggregate scheduling level production up to planning period p for state s Continuous Variables slm1_bulk ) medium-term scheduling feedback bulk overproduction slack variable slm2_bulk ) medium-term scheduling feedback bulk underproduction slack variable slm1_pack ) medium-term scheduling feedback packed overproduction slack variable slm2_pack ) medium-term scheduling feedback packed underproduction slack variable. When applying the rolling horizon approach, the overall time horizon is discretized into planning periods having end points which represent the moments in time when the planning and scheduling levels directly interact with one another. The length of the planning periods should reflect the inherent time horizon restrictions of the planning and scheduling models, as well as the computational and optimality tolerances of the user. A finer level of planning period discretization could theoretically lead to tighter production bounds at the expense of increased computational time. The first iteration of the forward rolling horizon approach entails the application of either the robust operational PPDM or operational CVAR-PPDM to the entire operational planning horizon. At this stage, the aggregate production totals for the overall operational planning time horizon can be used to determine an upper bound on the required amount of raw material, as well as the daily production profile for the first planning period. The forms of demand uncertainty are explicitly taken into account at the operational planning level so that the supplied production profile is immune to demand due date and demand amount uncertainty. The second iteration involves the medium-term scheduling model determining the appropriate allocation of resources along with a feasible daily production profile for the first planning period while simultaneously taking into account batch processing time uncertainty. Then, the selected operational PPDM is solved again for all those planning periods after the first planning period. The time horizon considered by the operational PPDM is still equal to the overall operational planning horizon; however, the production totals to be determined by the planning model for the second planning period are constrained by the production totals realized at the scheduling level for the first planning period. For example, eqs 47 and 48 provide upper and lower bounds, respectively, on the bulk production totals for the planning period following the last scheduled planning period, and these bounds are based upon the medium-term scheduling level’s aggregate production total for the bulk product class up to the given planning period, as well as the nominal demand profile. Equations 49 and 50 provide analogous upper and lower bounds, respectively, for the packed product class. The ratio of the mediumterm scheduling level’s aggregate production total up to the planning period of interest to the total nominal demand for the same period of time helps to restrict bulk and packed production in the next planning period to be scheduled. It should be noted that the underproduction and overproduction
feedback slack variables are penalized within the selected operational planning model’s objective function in order to constrain the production totals to be within the specified limits when feasible. Also, the planning period count parameter (countp) is set equal to 2 for the second application of the given operational planning model within the rolling horizon framework.
(∑
∑
) (∑
slm1_bulk e
(∑
)
∑
rs,d
s∈Sp∩Spb initial_ppedeterm_pp
∀p|p > 1, p ) count_p (47)
rs,d ·
s∈Sp∩Spb d 1
d)1
Parameters.
subject to
(54)
d)d'-1
Tu(u, d') e Tud' +
s∈Sp∩Spb
∀u, ∀d (53)
s∈Sp∩Su
Sets.
β·(
4961
∀s|s ∈ Sp ∩ Spp, ∀w (62) 5 · 6
∑
initial_mmedeterm_mm
∑
tot(s, d) e
rs,d
initial_mmedeterm_mm
∀s|s ∈ Sp ∩ Spp, ∀m (63) A.2. Robust Optimization Framework. The discussion to follow presents the general robust optimization approach, which has the goal of explicitly taking into account the various forms of parameter uncertainty found within model constraints and/ or objective function in a robust, deterministic manner. The first step in the robust optimization framework is to represent the true realization of the uncertain system parameters in the following manner: a˜ ) (1 + εξ)a
(64)
where a is an uncertain parameter with a˜ being its true realization, ε representing the relative degree of uncertainty, and ξ standing for a random variable with a known, but currently, undefined distribution. The true realization of the uncertain parameters can then be inserted into the relevant constraint(s). For example, eq 65 represents a generic model constraint in matrix form containing a vector of continuous variables x, a vector of binary variables y, matrices A and B which both contain uncertain parameters, and a vector p containing uncertain parameters. The given constraint is a part of a generic mixedinteger linear programming problem (eq 66). Ax + By e p
(65)
4962
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
min/max cTx + dTy
(66)
x,y
transform the mixed-integer linear programming problem into a mixed-integer linear or nonlinear programming problem.
subject to
min/max cTx + dTy subject to
Now, we observe in eqs 67 and 68 the replacement of the uncertain parameters with their respective true realizations. It should be noted that a small degree of constraint violation is allowed within the given constraints as evidenced by the presence of the δ term, and the subsets Ml and Kl encompass those parameters which possess a degree of uncertainty for a given constraint (l).
∑ a˜
lmxm
+
m
∑a
lmxm
∑ b˜ y
lk k
- p˜l e δ · max[1, |pl |]
∀l (67)
k
∑b y
+
lk k
m
- pl + ε(
k
∑ξ
lmalmxm
lk lk k
lmxm
+
m
∑b y
- ξlpl) e δ · max[1, |pl |]
lk k
- pl + ε(
k
∑ξ b y
lk lk k
∑ξ
lmalmxm
+
m∈Ml
- ξlpl) > δ · max[1, |pl |]} e κ (69)
lmalmxm
+
∑ξ b y
lk lk k
- ξlpl
- κ) ) f(λ1, |alm |xm, |blk |yk, |pl |)
lmxm
m
+
∑b y
lk k
k
∀l
∀l (73)
The final step is to augment the generic mixed-integer linear programming problem (eq 66) with the additional generic robust deterministic counterpart constraints (eq 74) so that any solution obtained is feasible for the nominal system conditions, as well as immune to the various forms of uncertainty present within the problem. Depending upon the form of uncertainty, the addition of the robust deterministic counterpart constraints may
2
lk
σlk2yk + pl2σl2 (75)
where σlm, σlk, and σl are the standard deviations for the independent normal random variables used to represent the uncertainty associated with the parameters alm, blk, and pl, respectively. Replacing the function f with its equivalent form for random variables following independent normal distributions, the deterministic robust counterpart constraint can be written as seen in eq 76. When considering only uncertainty associated with the right-hand side parameters, the robust formulation maintains its linearity.
∑a
+
∑b y
lk k
- pl +
k
ελ1
∑
alm2σlm2xm2 +
m∈Ml
∑b
2
lk
σlk2yk + pl2σl2 e
k∈Kl
δ · max[1, |pl |]∀l (76)
where λ1 can be determined from the following relationships:
κ)1-
(72)
∑b
k∈Kl
m∈Ml
(71)
- pl + εf(λ1, |alm |xm, |blk |yk, |pl |) e δ · max[1, |pl |]
σlm2xm2 +
2 lm
κ ) 1 - Pr{ξ e λ1} κ ) 1 - Fn(λ1)
Equation 73 represents the generic robust deterministic counterpart of the probabilistic constraint shown in eq 69 for any distribution of random variables.
∑a
∑a
(70)
k∈Kl
Fξ1(λ1) ) Pr{ξ1 e λ1} ) 1 - κ Fξ-1 (1 1
- pl + εf(λ1, |alm |xm, |blk |yk, |pl |) e
Making the assumption that all the relevant random variables are independent and follow a normal distribution, the function f can be represented in the following manner:
lmxm
The respective random variables can then be grouped into an aggregate random variable (eq 70), and the aggregate random variable’s cumulative distribution Fξ1 and inverse cumulative distribution F-1 ξ1 functions (eqs 71 and 72, respectively) can be used in order to reformulate the probabilistic constraint into its deterministic robust counterpart.
∑ξ
lk k
δ · max[1, |pl |] _x e x e jx y ) 0, 1
m
m∈Ml
∑b y k
∀l (68)
k∈Kl
ξ1 )
+
λ1
Equation 68 can be rewritten into a probabilistic form (eq 69) where the degree of constraint violation is captured by the user-specified parameter κ, which has bounds between 0 and 1. If the tolerance for constraint violation is low, then the value of κ should be similarly small.
∑a
lmxm
m
Ex + Fy ) e Ax + By e p
f(λ1, |alm |xm, |blk |yk, |pl |) )
+
k∈Kl
Pr{
∑a
m∈Ml
∑ξ b y
(74)
x,y
Ex + Fy ) e Ax + By e p _x e x e jx y ) 0, 1
∫
λ1
1
-∞
√2π
( )
exp
-x2 dx 2
(77)
For completeness and to demonstrate the flexibility of the robust optimization framework, the deterministic robust counterpart constraint for eq 78 will be developed. It should be noted that eq 78 is analogous to eq 65 with the exception that the inequality sign is inverted. Ax + By g p
(78)
To transform eq 78 into a form more amenable to the robust optimization framework, the inequality sign will be flipped by negating both sides of the given inequality constraint (eq 79). -Ax - By e -p
(79)
Similar to the previous development of the generic deterministic robust counterpart, we observe in eqs 80 and 81 the replacement of the uncertain parameters with their respective true realizations. It should be noted that a small degree of
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
constraint violation is allowed within the given constraints as evidenced by the presence of the δ term, and the subsets Ml and Kl encompass those parameters which possess a degree of uncertainty for a given constraint l. -
∑ a˜
∑ b˜ y
-
lmxm
lk k
m
-
∑a
lmxm
+ p˜l e δ · max[1, |pl |]
∀l (80)
-
∑b y
lk k
+ pl + ε(-
k
∑ξ
lmalmxm
-
+ ξlpl) e δ · max[1, |pl |]
lk lk k
∑a
lmxm
m
∑b y
-
lk k
+ pl + ε(-
k
∑ξ b y
lk lk k
∀l (81)
∑ξ
lmalmxm
-
m∈Ml
+ ξlpl) > δ · max[1, |pl |]} e κ (82)
k∈Kl
The respective random variables can then be grouped into an aggregate random variable (eq 83). It should be noted that the aggregate random variable for the given situation, ξ2, is the negative of the original grouped random variable, ξ1, defined in eq 70, and it is the negation which facilitates the evaluation of the left tail of the known distribution. As before, the given aggregate random variable’s cumulative distribution Fξ2 and inverse cumulative distribution Fξ2-1 functions (eqs 84 and 85, respectively) can be used in order to reformulate the probabilistic constraint into its deterministic robust counterpart. ξ2 ) -
∑ξ
lmalmxm
-
m∈Ml
∑ξ b y
lk lk k
+ ξlpl
(83)
k∈Kl
Fξ2(λ2) ) Pr{ξ2 e λ2} ) 1 - κ
(84)
Fξ2-1(1 - κ) ) f(λ2, |alm |xm, |blk |yk, |pl |)
(85)
Equation 86 represents the generic robust deterministic robust counterpart of the probabilistic constraint shown in eq 82 for any distribution of random variables. lmxm
+
∑b y
lk k
- pl - εf(λ2, |alm |xm, |blk |yk, |pl |) g
k
-δ · max[1, |pl |]
∀l (86)
Making the assumption that all the relevant random variables are independent and follow a normal distribution, the function f can be be expressed through the use of eq 87. Since a normal distribution is symmetric about the mean, the expression of function f is the same as in eq 75. Note that this would not be the case if the known distribution was not symmetric. f(λ2, |alm |xm, |blk |yk, |pl |) ) λ2
∑
m∈Ml
+
∑b y
- pl -
lk k
k
∑a
σlm2xm2 +
2 lm
m∈Ml
Equation 81 can be rewritten into a probabilistic form (eq 82) where the degree of constraint violation is captured by the userspecified parameter κ.
m
lmxm
m
alm2σlm2xm2 +
∑b
2 lk
2
lk
σlk2yk + pl2σl2 g
-δ · max[1, |pl |]
∀l (88)
The reader is directed to the work of Lin et al.,21 Janak et al.,22 and Verderame and Floudas20 for a complete exposition of how to formulate the deterministic robust counterpart for several other types of distributions besides the normal distribution. A.3. Conditional Value-at-Risk. Both value-at-risk and conditional value-at-risk seek to guard against unfavorable realization of uncertain parameters by going beyond expectation evaluation when expressing the uncertainty of system parameters. For example, in the portfolio management problem the investor is attempting to select those commodities which will maximize the portfolio’s profit; however, the values of these commodities are stochastic. In order to guard against potential losses, the investor defines a loss function, f(x, u), where x represents a decision variable (i.e., the vector of obtained commodities) and u is a stochastic vector (i.e., the vector of random commodity values) representing uncertainty and having a probability distribution, p(u). Both value-at-risk and conditional value-at-risk theory can be used in order to constrain the level of the loss function depending upon the investor’s degree of risk aversion. Value-at-risk can be defined as the minimum loss that is expected to be exceeded with a probability of (1 ω) while conditional value-at-risk is the expected loss given that the loss is greater than the value-at-risk for a given confidence level, ω. Conditional value-at-risk is the average of the probabilistic loss function’s tail defined by the lower and upper bounds of the value-at-risk and maximum loss, respectively. Conditional value-at-risk represents a more pessimistic loss level when compared to value-at-risk, and so, it is often the preferable representation of loss. Conditional value-at-risk also has the modeling advantage of maintaining convexity regardless of what probability distribution is being utilized, while value-at-risk is difficult to implement within an optimization model when the probability distribution is not normal or lognormal, because the resulting problem becomes nonconvex.44 Given the enhanced modeling capabilities, as well as the more pessimistic level of risk evaluation, CVAR has been applied to the demand uncertainty problem at the operational planning level, and the text to follow will provide an overview of the relevant conditional value-at-risk theory. When applying both conditional value-at-risk, as well as value-at-risk, a loss function, f(x, u), needs to be first defined where x represents a decision vector and u is a stochastic vector representing uncertainty and having a probability distribution, p(u). Having generically defined the loss function, the valueat-risk problem can be formulated as follows:
σlk2yk + pl2σl2 (87)
ξω ) minξ ξ
k∈Kl
where σlm, σlk, and σl are the standard deviations for the independent normal random variables used to represent the uncertainty associated with the parameters alm, blk, and pl, respectively, λ2 ) λ1, and λ1 can be determined from eq 77. Replacing the function f with its equivalent form for random variables following independent normal distributions, the de-
∑b
k∈Kl
m∈Ml
∑ξ b y
∑a
∑a
ελ2
k∈Kl
Pr{-
terministic robust counterpart constraint can be written as seen in eq 88, which is linear when there is only uncertainty associated with the right-hand side parameter, pl.
k
m
4963
(89)
subject to
∫
f(x,u)eξ
p(u) du g ω
where ω is a user-specified confidence level taking on a value between 0 and 1.
4964
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010
Having defined the ω -value-at-risk, ξω, by means of eq 89, the conditional expectation of the loss function being greater than ξω is expressed in eq 90.
ξ+
∫
Literature Cited
Φω(x) )
1 · (1 - ω)
f(x,u)>ξω
f(x, u) · p(u) du
(90)
Following the work of Rockafellar and Uryasev,26,27 a special function, Fω(x, ξ), where ξ represents a given loss level can be defined (eq 91). Fω(x, ξ) has the properties illustrated within eq 92, and due to these properties, it suffices to use Fω(x, ξ) instead of Φω(x) when applying conditional value-at-risk theory. Fω(x, ξ) ) ξ +
1 · (1 - ω)
∫
u∈Rn
max(0, f(x, u) - ξ) · p(u) du
(91) Φω(x) ) min Fω(x, ξ) ξ
ξω ) left endpoint of arg min Fω(x, ξ) ξ
Φω(x) ) Fω(x, ξω(x)) min Φω(x) ) min Fω(x, ξ) x
(92)
x,ξ
Fω(x, ξ) can be approximated by having the integral replaced by a summation over a number of scenarios where sc stands for a given scenario and πsc is the probability that a given scenario will occur (eq 93). ∼
Fω(x, ξ) ) ξ +
1 · (1 - ω)
∑π
· max(0, f(x, sc) - ξ)
sc
sc
(93) The loss function scenarios as denoted by f(x, sc) can be determined in part by sampling the underlying probability distribution of the stochastic vector u; consequently, assuming that the probability of each scenario occurring is equal to one another, then eq 93 can be rewritten into the form as shown in eq 94. ∼
Fω(x, ξ) ) ξ +
1 · (1 - ω) · |Sc|
∑ max(0, f(x, s) - ξ) sc
(94) where |Sc| is the number of scenarios. To include eq 94 representing the scenario-based conditional value-at-risk in a generic mixed-integer linear programming problem, the inner maximization term must be reformulated as seen in eq 95 where z(sc) is a positive variable representing the adjusted loss function for a given scenario sc. z(sc) g f(x, sc) - ξ
∀sc
(95)
With the aforementioned reformulation (eq 95), the scenariobased conditional value-at-risk (F˜ω(x, ξ)) can be restricted by means of eq 96 to be below some threshold defined by the parameter risk_aversion. Overall, eqs 95 and 96 are used to represent and bound the scenario-based conditional value-atrisk within an optimization model. The inclusion of the two given constraints within an optimization model helps to restrict the evaluation of the system’s variables in accordance with the level of risk aversion expressed by the parameter risk_aversion, which is a user-specified entity representing the user’s tolerance of risk. If the user has a low risk tolerance, then the value of risk_aversion should be similarly low when compared to a user with a high risk tolerance.
1 · (1 - ω) · |Sc|
∑ z(sc) e risk_aversion
(96)
sc
(1) Kallrath, J. Operational Plannning and Scheduling in the Process Industry. OR Spectrum 2002, 24, 219. (2) Pochet, Y.; Wolsey, L. A. Production Operational Planning by Mixed Integer Programming; Springer: New York, 2006. (3) Grossmann, I. E. Enterprise-wide Optimization: A New Frontier in Process Systems Engineering. AIChE J. 2005, 51, 1846. (4) Maravelias, C. T.; Sung, C. Integration of Production Planning and Scheduling: Overview, Challenges and Opportunities. Comput. Chem. Eng. 2009, 33, 1919. (5) Shobrys, D. E.; White, D. C. Operational Planning, Scheduling and Control Systems: Why They Cannot Work Together. Comput. Chem. Eng. 2002, 26, 149. (6) Papageorgiou, L. G.; Pantelides, C. C. Optimal Campaign Operational Planning/Scheduling of Multipurpose Batch/Semicontinuous Plants. 1. Mathematical Formulation. Ind. Eng. Chem. Res. 1996, 35, 488. (7) Papageorgiou, L. G.; Pantelides, C. C. Optimal Campaign Operational Planning/Scheduling of Multipurpose Batch/Semicontinuous Plants. 2. A Mathematical Decomposition Approach. Ind. Eng. Chem. Res. 1996, 35, 510. (8) Subrahmanyam, S.; Pekny, J. F.; Reklaitis, G. V. Decomposition Approaches to Batch Plant Design and Operational Planning. Ind. Eng. Chem. Res. 1996, 35, 1866. (9) Bassett, M. H.; Dave, P., III; 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, 821. (10) Dimitriadis, A. D.; Shah, N.; Pantelides, C. C. RTN-Based Rolling Horizon Algorithms for Medium Term Scheduling of Multipurpose Plants. Comput. Chem. Eng. 1997, 35, S1061. (11) Zhu, X. X.; Majozi, T. Novel Continuous-Time MILP Formulation for Multipurpose Batch Plants. 2. Integrated Operational Planning and Scheduling. Ind. Eng. Chem. Res. 2001, 40, 5621. (12) Munawar, S. A.; Gudi, R. D. A Multilevel, Control-Theoretic Framework for Integration of Planning, Scheduling, and Rescheduling. Ind. Eng. Chem. Res. 2005, 44, 4001. (13) Dogan, M. E.; Grossmann, I. E. A Decomposition Method for the Simultaneous Planning and Scheduling of Single-Stage Continuous Multiproduct Plants. Comput. Chem. Eng. 2006, 45, 299. (14) Dogan, M. E.; Grossmann, I. E. Simultaneous Planning and Scheduling of Single-Stage Multiproduct Continuous Plants with Parallel Lines. Comput. Chem. Eng. 2008, 32, 2664. (15) Petkov, S. B.; Maranas, C. D. Multiperiod Planning and Scheduling of Multiproduct Batch Plants under Demand Uncertainty. Ind. Eng. Chem. Res. 1997, 36, 4864. (16) Wu, D.; Ierapetritou, M. Heirarchical Approach for Production Planning and Scheduling under Uncertainty. Comput. Chem. Eng. 2007, 46, 1129. (17) Sahinidis, N. V. Optimization under Uncertainty: State-of-the-Art and Opportunities. Comput. Chem. Eng. 2004, 28, 971. (18) Li, Z.; Ierapetritou, M. Process Scheduling under Uncertainty: Review and Challenges. Comput. Chem. Eng. 2008, 32, 715. (19) 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, 4845. (20) Verderame, P. M.; Floudas, C. A. Operational Planning of LargeScale Industrial Batch Plants under Demand Due Date and Amount Uncertainty: I. Robust Optimization Framework. Ind. Eng. Chem. Res. 2009, 48, 7214. (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. (22) Janak, S. L.; Lin, X.; Floudas, C. A. A New Robust Optimization Approach for Scheduling under Uncertainty: II. Uncertainty with Known Probability Distribution. Comput. Chem. Eng. 2007, 31, 171. (23) Verderame, P. M.; Floudas, C. A. Operational Planning of LargeScale Industrial Batch Plants under Demand Due Date and Amount Uncertainty: II. Conditional Value-at-Risk Framework. Ind. Eng. Chem. Res. 2009, 49, 260. (24) Artzner, P.; Delbaen, F.; Eber, J. M.; Heath, D. Coherent Measure of Risk. Math. Finance 1999, 9, 203. (25) Ogryczak, W. Stochastic Dominance Relation and Linear Risk Measures; Progress and Business Publication: Cracow, Poland, 1999.
Ind. Eng. Chem. Res., Vol. 49, No. 10, 2010 (26) Rockafellar, R.; Uryasev, S. Optimization of Conditional Valueat-Risk. J. Risk 2000, 2, 21. (27) Rockafellar, R.; Uryasev, S. Conditional Value-at-Risk for General Loss Distributions. J. Bank. Finance 2002, 26, 1443. (28) 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. (29) Janak, S. L.; Floudas, C. A.; Kallrath, J.; Vormbrock, N. Production Scheduling of a Large-Scale Industrial Batch Plant. II. Reactive Scheduling. Ind. Eng. Chem. Res. 2006, 45, 8253. (30) Wang, W.; Ahmed, S. Sample Average Approximation of Expected Value Constrained Stochastic Programs. Operat. Res. Lett. 2008, 36, 515. (31) 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. (32) Ierapetritou, M. G.; Floudas, C. A. Effective Continuous-Time Formulation for Short-Term Scheduling. 2. Continuous and Semicontinuous Processes. Ind. Eng. Chem. Res. 1998, 37, 4360. (33) Ierapetritou, M. G.; Hene’, T. S.; Floudas, C. A. Effective Continuous-Time Formulation for Short-term Scheduling. 3. Multiple Intermediate Due Dates. Ind. Eng. Chem. Res. 1999, 38, 3446. (34) Lin, X.; Floudas, C. Design, Synthesis and Scheduling of Multipurpose Batch Plants via an Effective Continuous-Time Formulation. Comput. Chem. Eng. 2001, 25, 665. (35) Lin, X.; Floudas, C. A.; Modi, S.; Juhasz, N. M. Continuous-Time Optimization Approach for Medium-Range Production Scheduling of a Multiproduct Batch Plant. Ind. Eng. Chem. Res. 2002, 41, 3884. (36) Janak, S. L.; Lin, X.; Floudas, C. A. Enhanced Continuous-Time Unit-Specific Event-Based Formulation for Short-Term Scheduling of
4965
Multipurpose Batch Processes: Resource Constraints and Mixed Storage Policies. Ind. Eng. Chem. Res. 2004, 43, 2516. (37) Shaik, M. A.; L., J. S.; Floudas, C. A. Continuous-Time Models for Short-Term Scheduling of Multipurpose Batch Plants: A Comparative Study. Ind. Eng. Chem. Res. 2006, 45, 6190. (38) Shaik, M. A.; Floudas, C. A. Unit-Specific Event-Based ContinuousTime Approach for Short-Term Scheduling of Batch Plants Using RTN Framework. Comput. Chem. Eng. 2008, 32, 260. (39) 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. (40) Shaik, M. A.; Floudas, C. A.; Kallrath, J.; Pitz, H. Production Scheduling of a Large-Scale Industrial Continuous Plant: Short-Term and Medium-Term Scheduling. Comput. Chem. Eng. 2009, 33, 670. (41) Shaik, M. A.; Floudas, C. A. Novel Unified Modeling Approach for Short-Term Scheduling. Ind. Eng. Chem. Res. 2009, 48, 2947. (42) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMS: A User’s Guide; GAMS: San Francisco, CA, 2003. (43) CPLEX, ILOG CPLEX 9.0 User’s Mannual; ILOG, Inc.: Mountain View, CA, 2005. (44) Krokhmal, P.; Uryasev, S.; Zrazhevsky, G. Risk Management for Hedge Fund Portfolios: A Comparative Analysis of Linear Rebalancing Strategies. J. Altern. InVest. 2002, 5, 10.
ReceiVed for reView December 11, 2009 ReVised manuscript receiVed March 5, 2010 Accepted March 25, 2010 IE901973E