ARTICLE pubs.acs.org/IECR
Multisite Planning under Demand and Transportation Time Uncertainty: Robust Optimization and Conditional Value-at-Risk Frameworks Peter M. Verderame and Christodoulos A. Floudas* Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States ABSTRACT: The operational planning of a multisite production and distribution network entails determining the daily production and shipment profiles for the supply chain under consideration. The generated profiles should provide a tight upper bound on the true production capacity of the supply chain to ensure the maximization of customer satisfaction along with the minimization of resource misallocation. With an operational planning time horizon of 1-3 months, it is also imperative to take into account pertinent parameter uncertainty, so that the production and shipment profiles are not only a tight upper bound on the supply chain’s production capacity but also immune to the different forms of system uncertainty, such as demand due date, demand amount, and transportation time uncertainty. These types of uncertainty have been explicitly considered by means of the robust optimization framework and conditional value-at-risk theory. An industrial case study has been undertaken to demonstrate the viability of the proposed multisite operational planning under uncertainty approaches.
1. INTRODUCTION The operational planning of a supply chain entails determining the daily production and shipment profiles for each production facility under consideration, to satisfy the demand requirements of the various customer distribution centers within the supply chain. Levi et al.,1 Shah,2,3 Varma et al.,4 and Papageorgiou5 presented some of the key issues and challenges related to supply chain management. To avoid the misallocation of resources, it is paramount that the generated profiles represent a tight upper bound on the true production capacity of the supply chain. One potential option is to apply a scheduling model that also takes into account the transportation considerations associated with shipping final product from the production facilities to the customer distribution centers. However, this option is computationally prohibitive when considering an operational planning time horizon of 1-3 months, along with the inherently complicated supply chain infrastructure. To address this issue, certain aggregation or relaxation schemes have been adopted when formulating operational planning models. Wilkinson et al.6 and Bassett et al.7 applied time and model aggregation approaches, which involve the discretization of the time horizon into time periods and the grouping of related scheduling level constraints and variables into aggregate resource, capacity, and production constraints and variables. The values of variables are only calculated at the end of a discrete time period. Kallrath and co-workers8-13 developed a comprehensive multisite, discrete-time planning model that divides the planning horizon into commercial periods and then further discretizes those commercial periods into production periods. A unit aggregation approach is another modeling alternative where the production of pertinent final states is expressed through a facility’s set of bottleneck reactors, and Sung and Maravelias14 noted that unit aggregation schemes often have been adopted in discrete manufacturing industries. Verderame and Floudas15 developed a discrete-time operational planning model that was r 2010 American Chemical Society
based on the unit aggregation approach and showed, through computational studies, that it provides a production profile, which is a tight upper bound on the production capacity of the supply chain under investigation. The viability of any proposed modeling approach is ultimately determined by its ability to provide a tight upper bound on the true production capacity of the supply chain, to ensure the maximization of customer satisfaction, along with the minimization of resource misallocation. However, because of the length of the operational planning time horizon, parameter uncertainty must be explicitly addressed as well. Li and Ierapetritou16 and Verderame et al.17 reviewed and highlighted the importance of explicitly considering the multiple forms of parameter uncertainty when addressing supply chain management problems. One approach toward modeling parameter uncertainty is two-stage stochastic programming that is defined by an outer deterministic model and an inner recourse model dealing with uncertain parameter realization, and Gupta and Maranas,18,19 Levis and Papageorgiou,20 Wu and Ierapetritou,21 Colvin and Maravelias,22 and Al-Qahtani and Elkamel23 all have applied two-stage stochastic programming to address parameter uncertainty in various planning problems. Clay and Grossmann24 noted that, within the field of planning, it is often neither optimal nor computationally efficient to decouple the problem in the aforementioned manner. An alternative to two-stage stochastic programming is parametric programming, which is based on the theory of sensitivity analysis, where the objective is to define a function that maps the uncertain Special Issue: Puigjaner Issue Received: July 1, 2010 Accepted: September 9, 2010 Revised: September 7, 2010 Published: October 12, 2010 4959
dx.doi.org/10.1021/ie101401k | Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research parameter values to a given optimal solution for the entire uncertain parameter space. Pistikopoulos and Dua25 and Ryu et al.26 demonstrated that parametric programming can be used to address certain variants of the planning under uncertainty problem. The potential downside of parametric programming is that it entails solving a series of optimization problems in an iterative approach, and, therefore, the computational time required may make the application of parametric programming to large-scale industrial problems challenging. Based on the work of Charnes and Cooper,27-29 chance constraint programming does not possess the same potential computational burdens of twostage stochastic programming or parametric programming, and Petkov and Maranas,30 Li et al.,31 and Mitra et al.32 utilized chance constraint programming when solving specific aspects of the supply chain management under uncertainty problem. Yet, chance constraint programming does not guarantee the feasibility of an obtained solution with a specific realization of the uncertain parameters. The robust optimization framework that is based on the work of Ben-Tal et al.,33-35 Ghaoui et al.,36,37 Bertsimas et al.,38-40 Lin et al.,41 and Janak et al.42 overcomes the inherent shortcoming of chance constraint programming without the additional computational burden commonly associated with two-stage stochastic programming and parametric programming. Verderame and Floudas43 applied the robust optimization framework to take into account demand due date and demand amount uncertainty as it pertains to a single production facility. Verderame and Floudas44 also introduced conditional value-at-risk (CVAR) theory, which was developed by Artzner et al.,45 Ogryczak,46 and Rockafellar and Uryasev,47,48 to the same operational planning problem under the aforementioned forms of demand uncertainty. Carneiro et al.49 utilized CVAR when considering the risk management for an oil supply chain at the strategic planning level. Verderame and Floudas50 extended the work of Verderame and Floudas43,44 to apply the robust optimization framework and CVAR theory to the integration of operational planning and medium-term scheduling under demand and processing time uncertainty. Overall, all of the contributions indicated that the reliability of planning decisions is greatly enhanced through the explicit consideration of parameter uncertainty. The primary objective of this work is to apply the robust optimization framework, as well as CVAR theory, to the multisite operational planning problem under multiple forms of system uncertainty, which has not been addressed previously within the literature. The multisite planning with production disaggregation model (PPDM) developed by Verderame and Floudas15 will be used as an initiation point. Both demand due date and demand amount uncertainty, along with transportation time uncertainty, will be explicitly modeled, to ensure that the obtained production and shipment profiles are immune to the perceived forms of parameter uncertainty. The remainder of the paper will take the following form. The supply chain planning problem under consideration will be presented, followed by the proposed robust multisite PPDM. A novel multisite production shifting planning model (PSPM), utilizing both the robust optimization framework and CVAR theory in a hierarchical manner, then will be covered. A computational study will be presented to demonstrate the viability of the proposed approaches toward addressing the supply chain operational planning under uncertainty problem. Finally, some concluding remarks will be made to highlight the key contributions of the given work.
ARTICLE
Table 1. Supply Chain Network Representation
2. PROBLEM DESCRIPTION The supply chain of interest was previously studied within the work of Verderame and Floudas,15 and Table 1 shows the supply chain topology, which is composed of three batch production facilities that are capable of producing hundreds of both made-toorder (bulk) and made-to-store (packed) products, which are then shipped to three customer distribution centers. Each production facility is a multipurpose and multiproduct batch plant characterized by its set of 13 bottleneck reactors, which can be used to model the final state production at the operational planning level. The bulk products have intermediate demand due dates, which have uncertainty associated with both the requested demand due date and demand amount. The packed products have weekly demand totals, which possess demand amount uncertainty. There is also transportation time uncertainty, meaning that there is uncertainty associated with the required time to ship an allotment of final product from a given production facility to a customer distribution center, where the customer receives the order of interest. Due to the fact that the operational time horizon is on the order of 1-3 months, it is important to take into account these aforementioned forms of system uncertainty, which may manifest themselves at the operational planning level. Both the robust optimization framework and conditional valueat-risk theory (CVAR) can be utilized to take into account the supply chain uncertainty, and the following sections will present the robust and CVAR multisite planning models. 3. ROBUST OPTIMIZATION FRAMEWORK Before presenting the robust multisite planning model, an overview of the robust optimization framework is presented within this section. The reader should consult the works of Lin et al.,41 Janak et al.,42 and Verderame and Floudas43 for a complete elucidation of the given robust optimization framework. The goal of the robust optimization framework is to 4960
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
explicitly take into account the various forms of parameter uncertainty found within the 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 þ EξÞa
Ax þ By e p
ð2Þ
min=max cT x þ dT y x, y s:t: Ex þ Fy ¼ e Ax þ By e p _x e x e x
ð3Þ
y ¼ 0, 1 Now, in eqs 4 and 5, we observe the replacement of the uncertain parameters with their respective true realizations. Note 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. ~alm xm þ ~blk yk - ~pl e δ max 1, jpl j "l ð4Þ
∑k
∑m alm xm þ ∑k blk yk - pl þ E m ∑∈ M ξlm alm xm þ k ∑∈ K ξlkblkyk - ξlpl l
!
∑
m ∈ Ml
ξlm alm xm þ
∑k
∑
!
> δ max 1, jpl j
ξlk blk yk - ξl pl
ξlk blk yk - ξl pl
ð7Þ
Fξ1 ðλ1 Þ ¼ Prfξ1 e λ1 g ¼ 1 - k
ð8Þ
Fξ-1 ð1 - kÞ ¼ f ðλ1 , jalm jxm , jblk jyk , jpl jÞ 1
ð9Þ
"l
ð10Þ
The final step is to augment the generic MILP problem (eq 3) with the additional generic robust deterministic counterpart constraints (eq 11), 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 on the form of uncertainty, the addition of the robust deterministic counterpart constraints may transform the MILP problem to a mixed-integer nonlinear programming (MINLP) problem or a different MILP problem.
s:t:
min=max cT x þ dT y x, y Ex þ Fy ¼ e Ax þ By e p
∑m alm xm þ ∑k blkyk - pl þ Ef ðλ1, jalm jxm, jblk jyk , jpl jÞ e δ max½1, jpl j
" l x e x e xy ¼ 0, 1 _
ð11Þ
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: f ðλ1 , jalm jxm , jblk jyk , jpl jÞ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ λ1 alm 2 σlm 2 xm 2 þ blk 2 σ lk 2 yk þ pl 2 σl 2
∑
m ∈ Ml
∑
k ∈ Kl
ð12Þ 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 13. When considering only uncertainties associated with the right-hand side parameters, the robust formulation maintains its linearity.
ek
∑
k ∈ Kl
e δ max½1, jpl j
ð5Þ
"l
k ∈ Kl
ξlm alm xm þ
∑m alm xm þ ∑k blk yk - pl þ Ef ðλ1 , jalm jxm , jblkjyk, jpl jÞ
e
Equation 5 can be rewritten into a probabilistic form (eq 6), where the degree of constraint violation is captured by the userspecified parameter κ, which has bounds between 0 and 1. If the tolerance for constraint violation is low, then the value of κ should be similarly small. alm xm þ blk yk - pl Pr þE
∑
m ∈ Ml
Equation 10 represents the generic robust deterministic counterpart of the probabilistic constraint shown in eq 6 for any distribution of random variables:
l
δ max 1, jpl j
∑m
ξ1 ¼
ð1Þ
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 2 represents a generic model constraint in matrix form that contains 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 that contains uncertain parameters. The given constraint is a part of a generic mixedinteger linear programming (MILP) problem (eq 3).
∑m
distribution Fξ1-1 functions (eqs 8 and 9, respectively) can be used to reformulate the probabilistic constraint into its deterministic robust counterpart.
∑m alm xm þ ∑k blk yk - pl
ð6Þ þ Eλ1
The respective random variables can then be grouped into an aggregate random variable (eq 7), and the aggregate random variable’s cumulative distribution Fξ1 and inverse cumulative
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi alm 2 σlm 2 xm 2 þ blk 2 σlk 2 yk þ pl 2 σl 2
∑
m ∈ Ml
∑
k ∈ Kl
e δ max½1, jpl j 4961
"l
ð13Þ
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
Figure 1. Schematic diagram detailing the robust multisite planning with production disaggregation model (robust multisite PPDM).
Here, λ1 can be determined from the following relationship: k ¼ 1 - Prfξ e λ1 g k ¼ 1 - Fn ðλ1 Þ ! Z λ1 1 -x2 pffiffiffiffiffiffi exp k ¼ 1dx ð14Þ 2 2π -¥ Although the given overview assumes, in part, that the uncertain parameters each followed independent normal distributions, the robust optimization approach can be applied to problems where the uncertain parameters do not follow a normal distribution, and Janak et al.42 have presented the theoretical underpinnings for taking into account uncertain parameters, which follow a variety of continuous, as well as discrete, distributions. For the sake of simplicity of presentation, universal values for ɛ, κ, and δ were adopted; however, as noted by Janak et al.,42 the robust optimization framework can easily be extended to take into account scenarios where ɛ can vary from parameter to parameter and κ and δ can vary from constraint to constraint.
4. ROBUST MULTISITE PLANNING WITH PRODUCTION DISAGGREGATION MODEL (ROBUST MULTISITE PPDM) 4.1. Overview. To address the multiple forms of parameter uncertainty that must be explicitly taken into account at the operational planning level, the multisite planning model developed by Verderame and Floudas15 has been transformed to its robust counterpart by means of the robust optimization framework. To facilitate the given transformation, the original customer distribution center demand constraints have been reformulated into a truncated set of underproduction and overproduction demand constraints that encompass both the demand due date and demand amount uncertainty, as well as the aforementioned transportation time uncertainty associated with the supply chain of interest. As a result, only the given set of demand constraints must be transformed to their robust counterparts and augmented to the nominal multisite planning model to ensure that any obtained solution is feasible for the nominal set of system conditions, as well as being robust with regard to the unfavorable realization of the given uncertain parameters. The proposed discrete-time planning model is based on a unit aggregation approach, where the production of all final states is taken into account through the set of bottleneck units present at each production facility. To capture the continuous-time nature of each production facility within the discrete-time planning model, the task duration constraints allow for the carryover of unused processing time from one day to the next, and the following text will provide a detailed description of all relevant model
ARTICLE
constraints. Figure 1 provides an overview regarding the required supply chain input data and the pertinent results generated through the application of the robust multisite planning with production disaggregation model (robust multisite PPDM). The following indices, sets, parameters, and variables are required for the robust multisite PPDM. 4.2. Nomenclature Indices i = production facilities c = customer distribution centers u = reactors s = states d = days t = weeks m = months p = planning periods or = customer orders ap = delivery summation possibilities l = delivery combinations Sets I = production facilities C = customer distribution centers U = reactors S = material states D = days in the overall planning horizon T = weeks in the overall planning horizon M = months in the overall planning horizon P = planning periods in the overall planning horizon OR = customer order AP = shipment addition possibility; AP = 0, 1, 2, 3 L = delivery combinations Is = facilities that can produce state s Cs = customer distribution centers that can receive state s Us = reactors that can produce state s Ui = reactors that are present at facility i Si = states that can be produced at facility i Sc = states that can be shipped to center c Sp = states that are final products Sbp = states that are bulk final products Spp = states that are packed final products Su = states that can be produced in reactor u Parameters dft = last day of week t dfm = last day of month m dfp = last day of planning period p Fixedtimei,u,s = batch processing time for state s in reactor u at facility i Ti, u, d = total available processing time for reactor u at facility i on day d Hi, u = lower bound on aggregate processing time for reactor u at facility i Capmini, u, s = minimum batch size for state s in reactor u at facility i Capmaxi, u, s = maximum batch size for state s in reactor u at facility i prices = price for state s Cti, c, s = fixed transportation cost from facility i to center c for state s σ_transi, c, s = proportional transportation cost from facility i to center c for state s 4962
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Cleanup = batch changeover downtime Δi, c = transportation time from facility i to center c initial_pp = first day within planning period p term_pp = last day within planning period p rc, s, d = demand for state s on day d at center c (intermediate demand due date parameter) r_orc, s, d, or = demand for state s at center c on day d related to order or (intermediate demand due date order parameter) prob_orc, s, d, or = probability that the order or for state s will be realized on day d at center c prob_lc, d, s, l = probability that delivery combination l for state s will arrive at center c on day d lastc, s, d, or = takes on a value of 1 if day d is the last day that the customer order or for state s can be realized at center c σd = time-horizon-dependent standard deviation ɛ = degree to which an uncertain parameter is known κ1 = packed product underproduction probabilistic level of constraint violation κ2 = packed product overproduction probabilistic level of constraint violation δr = allowable degree of constraint violation Continuous Variables Tu(i, u, d) = total time reactor u is utilized on day d within facility i sla(i, u) = processing time underutilization slack variable for reactor u within facility i tot1(i, u, s, d) = amount of state s produced on day d in reactor u within facility i tot(i, s, d) = amount of state s produced on day d within facility i D(i, c, s, d) = amount of state s sent from facility i to center c on day d Inv_bulk(c, s, d) = inventory of bulk product state s within center c on day d Under_bulk(c, s, d) = unsatisfied customer demand for bulk state s at center c on day d Under_pack(c, s, t) = unsatisfied customer demand for packed state s at center c on day d del_sum(c, s, d, l) = takes on a value of 1 if delivery combination l is realizable for state s at center c on day d Binary Variables y(i, u, s, d) = assigns the production of state s to reactor u on day d within facility i w(i, c, s, d) = assigns an amount of state s to be sent from facility i to center c on day d del_sum_a(c, s, d, ap) = indicates if delivery summation possibility ap has occurred at center c for state s on day d. Using the above notation, the robust multisite PPDM is composed of the following constraints and objective function. 4.3. Model. 4.3.1. Duration Constraints. Equation 15 determines, on a daily basis, the total time that bottleneck unit u located at facility i is occupied with processing tasks or is down because of requisite cleanup operations:
∑
s ∈ Sp ∩ Su ∩ Si
For the first day within the time horizon, eq 16 places an upper bound on the available time a given unit can be in an active state characterized by undertaking processing or cleanup activities as being less than or equal to the user-specified parameter Ti, u, d, which denotes the total operation time for unit u located at facility i on day d. T u ði, u, dÞ e Ti, u, d
After the initial day, eq 17 provides an adjusted bound on the total available processing time for each unit under consideration, which allows for the carryover of unutilized time from one day to the next. In doing so, eq 17 enables the proposed discrete-time model to approach the continuous-time nature of plant operations and avoid forced downtime due to discretization of the time horizon. Ti, u, d0 - T u ði, u, d0 Þ T u ði, u, dÞ e Ti, u, d þ 0 d ed-1 ð17Þ " i, " uju ∈ Ui , " djd > 1
∑
Equation 18 ensures an efficient utilization of each bottleneck unit; however, the given constraint is not strictly enforced, as denoted by the presence of the slack variable sla(i, u), which is penalized within the objective function.
∑d T uði, u, dÞ þ slaði, uÞ g Hi, u
∑
s ∈ Sp ∩ Su ∩ Si
Capmini, u, s 3 yði, u, s, dÞ e tot1ði, u, s, dÞ e Capmax i, u, s 3 yði, u, s, dÞ " i, " uju ∈ Ui , " sjs ∈ Sp ∩ Su ∩ Si , " d
∑
u ∈ U i ∩ Us
ð15Þ
ð19Þ
tot1ði, u, s, dÞ ¼ totði, s, dÞ " i, " sjs ∈ Sp ∩ Si , " d
ð20Þ
The total amount of state s produced at facility i on day d is then shipped to the multiple customer distribution centers, as represented by eq 21. Equations 22 and 23 link the shipment activation variables (w(i, c, s, d)) with the actual shipment of product, to explicitly take into account the fixed cost associated with the shipment of product. totði, s, dÞ ¼
∑
c ∈ Cs
Dði, c, s, dÞ
" i, " sjs ∈ Sp ∩ Si , " d ð21Þ
Dði, c, s, dÞ e
∑
u ∈ Ui ∩ Us
Capmax i, u, s 3 wði, c, s, dÞ
"c, " i, " sjs ∈ Sp ∩ Si ∩ Sc , " d Dði, c, s, dÞ g wði, c, s, dÞ " c, " i, " sjs ∈ Sp ∩ Si ∩ Sc , " d
yði, u, s, dÞ ¼ T u ði, u, dÞ
" i, " uju ∈ Ui , " d
" i, " uju ∈ Ui ð18Þ
4.3.2. Production and Capacity Constraints. The maximum and minimum batch size for task i occurring within unit u producing state s on day d is taken into account through eq 19, and eq 20 provides the balance to determine the total amount of state s produced within facility i on day d:
½Fixedtimei, u, s 3 yði, u, s, dÞ þ Cleanup 3
" i, " uju ∈ Ui , " djd ¼ 1 ð16Þ
ð22Þ
ð23Þ
4.3.3. Nominal Demand Constraints. Equations 24-27 are the reformulated bulk and packed underproduction and 4963
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
overproduction demand constraints, which represent an enhancement over the set of demand constraints originally presented by Verderame and Floudas,15 that required additional variables and constraints to express the given demand satisfaction restrictions. Equations 24 and 25 are the bulk daily underproduction and overproduction demand constraints, respectively, which allow for the carryover of inventory and unsatisfied demand. For the packed product class, eq 26 is the weekly underproduction demand constraint, and eq 27 is the monthly overproduction demand constraint.
∑
∑
i ∈ Is d0 e ðd - Δi, c Þ
Dði, c, s, d0 Þ þ Under_bulkðc, s, dÞ g
∑
d0 e d
" c, " sjs ∈ Sbp ∩ Sc , " d
∑
∑
i ∈ Is d0 e ðd - Δi, c Þ
rc, s, d0
ð24Þ
Dði, c, s, d0 Þ - Inv_bulkðc, s, dÞ e
∑
d0 e d
rc, s, d0
" c, " sjs ∈ Sbp ∩ Sc , " d ð25Þ
∑
∑
i ∈ Is d e ðd f - Δ Þ i, c t
Dði, c, s, dÞ þ Under_packðc, s, tÞ g
∑
i ∈ Is d e ðdf - Δ Þ i, c m
f
rc, s, d
ð26Þ
Dði, c, s, dÞ e
∑
f
d e dm
r c, s , d
ð27Þ
" c, " sjs ∈ Spp ∩ Sc , " m 4.3.4. Robust Demand Constraints. For bulk products, there exist demand due date and demand amount uncertainty, and to model these two forms of demand uncertainty concurrently by means of the robust optimization framework, each demand due date parameter (rc, s, d) is given a unique order designation. The demand due date uncertainty is assumed to follow a uniform discrete distribution with bounds of plus or minus one day from the nominal demand due date, and, therefore, a given order has a 33% chance of being realized on the previous day (r_orc, s, d-1, or), a 67% chance of being realized on the nominal day (r_orc, s, d, or), and a 100% chance of being realized on the following day (r_orc, s, dþ1, or). Each possible realization of the customer order has a realized amount that follows a specified probability distribution (eq 28), where ε denotes an uncertainty level, ξc, s, d represents the random variable having a distribution that might be state-, center-, and/or time-specific, and ~r _orc, s, d, or is the true realization of the r_orc, s, d, or parameter. ~r _orc, s, d - 1, or ¼ ð1 þ Eξc, s, d - 1 Þrc, s, d ~r _orc, s, d, or ¼ ð1 þ Eξc, s, d Þrc, s, d ~r _orc, s, d þ 1, or ¼ ð1 þ Eξc, s, d þ 1 Þrc, s, d
σd ¼ 0:05 þ 0:025ðp - 1Þ
ð29Þ
" d, " pjinitial_pp e d e term_pp
For the packed product class, demand due dates are considered to be deterministic in nature, while the demand amount is considered to be uncertain. As a result, the true realization of each packed product demand due date parameter (~r c, s, d) can be represented as shown in eq 30. As with the bulk product class, each ξc, s, d is assumed to follow a normal distribution with zero as its mean and a standard deviation defined by eq 29. However, note that the robust optimization framework can take into account many different distribution types for both demand due date and amount uncertainty, as demonstrated within the works of Janak et al.42 and Verderame and Floudas.43 ~r s, d ¼ ð1 þ Eξc, s, d Þrc, s, d
d e dt
" c, " sjs ∈ Spp ∩ Sc , " t
∑
∑
should conform to a normal distribution. Therefore, each ξc, s, d is modeled as a normally distributed random variable with a mean of zero and a standard deviation (σd) that varies linearly with the time horizon (eq 29):
To address transportation time uncertainty, it is assumed that each transportation time parameter (Δi, c) follows a uniform discrete distribution with bounds of plus or minus one day from the nominal value. If a given shipment of product is shipped from facility i to center c on day d, then it has a 33% chance of arriving by day (d þ Δi, c - 1), a 67% chance of arriving by (d þ Δi, c), and 100% chance of arriving by day (d þ Δi, c þ 1). The multisite planning model originally proposed by Verderame and Floudas15 allows for the consideration of the aforementioned transportation and demand uncertainty solely within the set of underproduction and overproduction demand constraints (eqs 24-27). However, to take into account both transportation time and demand uncertainty by means of the robust optimization framework, the given demand constraints must be reformulated further. Following the work of Verderame and Floudas,43 the right-hand side of the nominal bulk product demand constraints is adapted to facilitate the modeling of the demand due date uncertainty, while the right-hand side of the packed product demand constraints remain unchanged, because the packed products do not have due date uncertainty. The single summation on the left-hand side of both the bulk and packed nominal demand constraints is broken into three summations, to take into account the possibility of shipments arriving a day ahead of the nominal transportation time or a day after the nominal transportation time. Equations 31 and 32 are the modified bulk product underproduction and overproduction demand constraints, respectively, while eqs 33 and 34 are the modified packed product underproduction and overproduction demand constraints, respectively.
∑
∑
Dði, c, s, d0 Þ þ
∑
∑
Dði, c, s, d0 Þ þ Under_bulkðc, s, dÞ
i ∈ Is d0 e d - ðΔi, c þ 1Þ
ð28Þ
ð30Þ
i ∈ Is d0 ¼ d - ðΔi, c - 1Þ
g
Given the multiple forms of uncertainty that play a role in determining the demand amount for a given order, it follows from the Central Limit Theorem that the demand amount
∑
∑
i ∈ Is d0 ¼ d - Δi, c
Dði, c, s, d0 Þ þ
∑or d ∑e d r_orc, s, d , or - r_orc, s, d - 1, or 3 ð1 - lastc, s, d - 1, or Þ 0
0
0
0
" c, " sjs ∈ Sbp ∩ Sc , " d 4964
ð31Þ
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
∑ i∈I
s
∑ d e d - ðΔ 0
,
i c
þ 1Þ
Dði, c, s, d0 Þ þ
∑
∑
i ∈ Is d0 ¼ d - ðΔi, c - 1Þ
∑ i∈I
∑ d ¼d - Δ 0
s
,
ARTICLE
Dði, c, s, d0 Þ þ
i c
Dði, c, s, d0 Þ - Inv_bulkðc, s, dÞ e
∑or d ∑e d r_orc, s, d , or - r_orc, s, d - 1, or 3 ð1 - lastc, s, d - 1, or Þ 0
0
( Pr
∑
∑
∑
∑
i ∈ Is d0 ¼ d - ðΔi, c - 1Þ
∑
Dði, c, s, d0 Þ þ
i ∈ Is d0 e d - ðΔi, c þ 1Þ
∑
i ∈ Is d 0 ¼ d - Δ i, c
Dði, c, s, d0 Þ þ
Dði, c, s, d0 Þ - Inv_bulkðc, s, dÞ >
∑or d ∑e d ~r_orc, s, d , or - ~r_orc, s, d - 1, or 3
0
0
0
0
0
"c, " sjs ∈ Sbp ∩ Sc , " d
∑
∑
i ∈ I s d e df - ðΔ þ 1Þ i, c t
∑ i∈I
s
Dði, c, s, dÞ þ
f
∑
∑
∑
i ∈ Is d ¼ d f - Δ i, c t
ð32Þ Dði, c, s, dÞ þ
Y orjprob_orc, s, d, or > 0
∑
f
" c, " sjs ∈ Spp ∩ Sc , " t
∑
Dði, c, s, dÞ þ
∑
∑
∑
∑
ð33Þ Dði, c, s, dÞ þ
i ∈ Is d ¼ d f - Δ i, c m
i ∈ Is d ¼ df - ðΔ - 1Þ i, c m
∑
Dði, c, s, dÞ e
f
d e dm
r c, s , d
"c, " sjs ∈ Spp ∩ Sc , " m ð34Þ The modified demand constraints can then be transformed to probabilistic demand constraints, and eqs 35 and 36 are the bulk product underproduction and overproduction probabilistic demand constraints, respectively. Note that, for the bulk product demand constraints, the probability of the constraint being violated is equal to the product of the probability of the given order summation being realized and the probability of the given summation of shipment variables being realized. The positive variable del_sum(c, s, d, l) denotes the viability of a particular delivery summation combination, and section 4.3.5 describes how the aforementioned positive variable is constrained. ( Pr
∑
∑
i ∈ Is d0 e d - ðΔi, c þ 1Þ
∑
∑
i ∈ Is d0 ¼ d - ðΔi, c - 1Þ
Dði, c, s, d0 Þ þ
∑
Dði, c, s, d0 Þ þ
0
0
0
8 < Pr :i ∈ Is
∑
∑
1 - lastc, s, d0 - 1, or
)
0
Y
e
orjprob_orc, s, d, or > 0
∑l del_sumðc, s, d, lÞ 3 prob_lc, s, d, l
∑
"c, " sjs ∈ Sbp ∩ Sc , " d
∑
"
δ max 1, r
prob_orc, s, d, or 3
"c, " sjs ∈ Sbp ∩ Sc , " d
ð35Þ
∑
Dði, c, s, dÞ þ
i ∈ Is d ¼ df - ðΔ - 1Þ i, c t
∑
f
Dði, c, s, dÞ þ
Dði, c, s, dÞ þ Under_packðc, s, tÞ
39
rc, s, d 5
= ;
~r c, s, d -
d ¼ dt
" c, " sjs ∈ Spp ∩ Sc , " t
1 - lastc, s, d0 - 1, or -
0
f
∑
0
∑l del_sumðc, s, d, lÞ 3
d e dt - ðΔi, c þ 1Þ
i ∈ Is d ¼ df - ðΔ - 1Þ i, c m
∑or d ∑e d r_orc, s, d , or - r_orc, s, d - 1, or 3
e
ð36Þ The packed product probabilistic underproduction and overproduction demand constraints are shown in eqs 37 and 38, respectively. Since the packed products do not possess demand due date uncertainty, the probability of the packed product underproduction demand constraint (eq 37) being violated is equal to the product of the probability of the given left-hand side being realized with only the shipment arrivals associated with the binary variable del_sum(c, s, dft, l) having an arrival probability of 1 del_sum_aðc, s, d - 1, apÞ ¼
" c, " sjs ∈ Sp ∩ Sc , " d
ð51Þ 4.3.6. Objective Function. The robust multisite PPDM objective function shown in eq 52 penalizes bulk product underproduction and overproduction, packed product underproduction, transportation cost in the form of fixed and shipment proportional, and unit underutilization at each facility within the supply chain, as well as the maximization the gross profit generated from producing both the bulk and packed products. Each term within the objective function is assigned a weighting in accordance with its relative importance, as determined by the supply chain manager. 2
0
3
ð50Þ
Equation 51 mandates that only one delivery combination is possible for a given center, state, and day combination and, in turn, helps to tighten the LP relaxation of the problem.
∑ ∑ del_sumðc, s, d, lÞ ap l ¼ ðjAPj ðap - 1Þ þ ðap - 1Þ þ 1Þ 0
ð49Þ
del_sum_aðc, s, d, apÞ þ del_sum_aðc, s, d - 1, ap0 Þ - 1 e del_sumðc, s, d, lÞ
After the initial day within the time horizon, eqs 46 and 47 are enforced, to constrain the activation of del_sum(c, s, d, l). Since del_sum_ac, s, d, ap is activated when a given one-day-ahead addition possibility can be realized, eq 46 takes into account one-dayahead orders, while eq 47 covers nominal-day arrivals. For example, if there are two one-day-ahead shipment arrivals and three nominal-day shipment arrivals for a given center, state, and day combination, then delivery combination L12 =(AP3, AP4) should be activated, and eqs 46 and 47 accomplish this objective. Note that, since del_sum(c, s, d, l) is defined in terms of the binary variable del_sum_a(c, s, d, ap) through eqs 45-47, it can be declared as a positive variable with an upper bound equal to 1. del_sum_aðc, s, d, apÞ ¼
ð48Þ
0
3
0
ð47Þ An alternative to including eqs 46 and 47 would be to reformulate the product of del_sum_a(c, s, d, ap) and del_sum_a(c, s, d - 1, ap) through its equivalent convex and concave envelopes, as shown in eqs 48-50. However, the " ap, " c, " sjs ∈ Sp ∩ Sc , " djd > 1
4967
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
as value-at-risk, a loss function, f(x, u), first must be defined, where x represents a decision vector and u is a vector representing uncertainty with a probability distribution, p(u). Having generically defined the loss function, the value-at-risk problem can be formulated as seen in eq 53, where the objective is to minimize the value-at-risk (ξω), given that the probability of the value-atrisk being exceeded is less than or equal to (1 - ω), where ω is a user-specified confidence level that takes on a value between 0 and 1. ξω ¼ min ξ ξ ð53Þ R s:t: f ðx, uÞ e ξ pðuÞ du g ω After the ω value-at-risk (ξω) has been characterized by means of eq 53, the conditional expectation of the loss function being greater than ξω is expressed in eq 54: Z 1 Φω ðxÞ ¼ f ðx, uÞ 3 pðuÞ du ð54Þ ð1 - ωÞ f ðx, uÞ > ξω Following the work of Rockafellar and Uryasev,47,48 a special function (Fω(x, ξ), where ξ represents a given loss level) can be defined (see eq 55). Fω(x, ξ) has the properties illustrated within eq 56, and, because of these properties, it suffices to use Fω(x, ξ) instead of Φω(x) when applying the CVAR theory. Z 1 Fω ðx, ξÞ ¼ ξ þ maxð0, f ðx, uÞ - ξÞ 3 pðuÞ du ð1 - ωÞ u ∈ Rn ð55Þ
Figure 2. Sample average approximation (SAA) overview.
below some threshold that is defined by the parameter risk_aversion. Overall, eqs 59 and 60 are used to represent and bound the scenario-based CVAR 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 that represents 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. ξþ
Φω ðxÞ ¼ min Fω ðx, ξÞ ξ
ξω ¼ left end point of arg min Fω ðx, ξÞ ξ
Φω ðxÞ ¼ Fω ðx, ξω ðxÞÞ
ð56Þ
min Φω ðxÞ ¼ min Fω ðx, ξÞ ξ x, ξ Fω(x, ξ) can be approximated by having the integral replaced by a summation over many scenarios where sc represents a given scenario and πsc is the probability that a given scenario will occur (see eq 57). ∼
F ω ðx, ξÞ ¼ ξ þ
1 ð1 - ωÞ
∑sc πsc 3 maxð0, f ðx, scÞ - ξÞ
ð57Þ
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. Assuming that all scenarios have equal probabilities (i.e., 1/|SC|), eq 57 can be rewritten to the form shown in eq 58: ∼
F ω ðx, ξÞ ¼ ξ þ
1 ð1 - ωÞjSCj
∑sc maxð0, f ðx, sÞ - ξÞ
ð58Þ
where |Sc| is the number of scenarios. To include eq 58, which represents the scenario-based CVAR in a generic MILP problem, the inner maximization term must be reformulated as shown in eq 59, where z(sc) is a positive variable representing the adjusted loss function for a given scenario sc: zðscÞ g f ðx, scÞ - ξ
" sc
ð59Þ
With the aforementioned reformulation (eq 59), the scenariobased CVAR (F~ω(x, ξ)) can be restricted by means of eq 60 to be
1 ð1 - ωÞjSCj
∑sc zðscÞ e risk_aversion
ð60Þ
5.1. Sample Average Approximation. When applying the scenario-based CVAR framework, the issue of adequately sampling the given probability space and ensuring that the obtained solution is feasible when considering the complete uncertainty space arises. To maintain computational tractability, only a finite number of scenarios can be considered within the proposed optimization model, and therefore, an iterative approach must be adopted to guarantee that the final solution satisfies the specified CVAR constraints, when considering the entire uncertainty space. The sample average approximation (SAA), which was developed by Wang and Ahmed,51 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 given approach will be utilized when applying CVAR to the multisite planning under demand and transportation time uncertainty problem. The SAA entails solving the optimization model with a relatively small number of loss function scenarios and then applying a larger number of scenarios to calculate the average and standard deviation of the left-hand side of the encapsulated CVAR constraint (see eq 60). A z-transformation is then performed using the original right-hand side of the encapsulated CVAR constraint, along with the aforementioned average and standard deviation, and if the calculated z-value is large (i.e., z g 2), then it is unlikely, in a probabilistic sense, that the constraint would be violated when considering the entire uncertainty space. If the z-value is not sufficiently large, then the optimization model is resolved with a reduction in the allowable level of risk exposure to generate a more conservative solution that will hopefully satisfy the z-termination criterion. Figure 2 provides an overview of the SAA algorithm; for a detailed explanation of how 4968
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
The following additional indices, sets, parameters, and variables are required for the CVAR multisite PSPM. 6.2. Nomenclature
Index sc = scenario Set SC = scenario
Figure 3. Schematic of the CVAR Multisite production shifting planning model (PSPM).
to implement the SAA algorithm for planning under uncertainty problems, consult the works of Verderame and Floudas44 and AlQahtani and Elkamel.23
6. CVAR MULTISITE PRODUCTION SHIFTING PLANNING MODEL (CVAR MULTISITE PSPM) 6.1. Overview. In an effort to maintain model linearity, a scenario-based CVAR approach has been adopted. However, the computational complexity of replacing the nominal demand constraints with the set of CVAR demand constraints encapsulating the aforementioned demand and transportation time uncertainty within the multisite PPDM is prohibitive. To overcome this implementation issue, a novel production shifting planning model (PSPM) has been developed. The given linear programming model takes an existing production profile supplied from another operational planning model and shifts daily production and/or shipment levels to satisfy the CVAR demand constraints while concurrently attempting to minimize transportation cost and supplied production profile violation. Some model aggregation procedures also are undertaken to further mitigate the computational complexity brought on by implementing a scenario-based CVAR approach. Namely, the CVAR demand constraints are written on a product class basis as opposed to individual states, and the time indexing is aggregated as well. Since only a finite number of scenarios can be considered within the CVAR multisite PSPM, a SAA algorithm is applied to ensure that the obtained solution is a valid upper bound on the original minimization problem. Figure 3 provides a pictorial representation of how the CVAR multisite PSPM is applied. The quality of the final solution and the number of iterations of the SAA algorithm required to achieve convergence are affected by the initial supplied production profile. A valid candidate profile is the production profile that results from the application of the robust multisite PPDM. If the supplied production profile is obtained from the robust multisite PPDM, then the final production profile generated by the CVAR multisite PSPM represents a unified solution that has been constructed in a sequential fashion, using both the robust optimization framework (stage 1) and CVAR theory (stage 2). As a result, the proposed CVAR multisite PSPM can be characterized as an approach toward modeling demand and transportation time uncertainty which considers multiple uncertainty methodologies sequentially to generate a quality production profile with a high degree of uncertainty immunity.
Parameters dit = first day of week t dim = first day of month m r\_scc, s, d, sc = demand for state s at center c on day d related to scenario sc (intermediate demand due date scenario parameter) Δ_sci, c, sc = transportation time from facility i to center c related to scenario sc tot_pi, s, d = supplied planning production level at facility i for state s on day d ω = CVAR confidence level γ = CVAR exposure indicator Continuous Variables sla_del1(i, s, d), sla_del2(i, s, d), sla_del3(i, s, d) = delivery slack variables for facility i, state s, and day d ξbulk_under(c, t) = acceptable bulk product underproduction loss level for center c and week t ξbulk_over(c, m) = acceptable bulk product overproduction loss level for center c and month m ξpack_under(c, m) = acceptable packed product underproduction loss level for center c and month m ξpack_over(c, m) = acceptable packed product overproduction loss level for center c and month m zbulk_under(c, t, sc) = bulk underproduction loss function variable for center c, week t, and scenario sc zbulk_over(c, m, sc) = bulk overproduction loss function variable for center c, month m, and scenario sc zpack_under(c, m, sc) = packed underproduction loss function variable for center c, month m, and scenario sc zpack_over(c, m, sc) = packed overproduction loss function variable for center c, month m, and scenario sc. Using the above notation, the CVAR multisite PSPM is composed of the following constraints and objective function. 6.3. Model. 6.3.1. Planning Production Level Constraints. The primary objective of the CVAR multisite PSPM is to adapt an existing production profile originally generated from the application of another operational planning model by adjusting shipment and/or production totals, to satisfy the CVAR demand constraints that encapsulate the demand and transportation time uncertainty under consideration. Equation 61 provides the balance used to allocate the supplied daily production profile represented by the parameter tot_pi, s, d among the set of customer distributions centers. The presence of the slack variables sla_del1(i, s, d) and sla_del2(i, s, d) allows for alteration of the supplied production profile by permitting additional production or the placement of unnecessary product into inventory. tot_pi, s, d - sla_del1ði, s, dÞ þ sla_del2ði, s, dÞ ¼
∑
c ∈ Cs
Dði, c, s, dÞ
" sjs ∈ Sp , " iji ∈ Is , " d ð61Þ 4969
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Table 3. Robust Multisite PPDM Bulk Product Analysis
Equation 62 mandates that the level of overproduction in excess of the supplied production profile up to a given day d be less than or equal to the level of underconsumption denoted by the left-hand side summation. Therefore, eq 62 constrains the values that the slack variables sla_del1(i, s, d) and sla_del2(i, s, d) can undertake by attempting to enforce the logic condition that additional required production indicated by the activation of sla_del2(i, s, d) cannot be undertaken if it has not previously been accrued earlier in the time horizon, as determined by the activation of the slack variable sla_del1(i, s, d). The presence of an additional slack variable, sla_del3(i, s, d), which is penalized within the model’s objective function, allows for violation of the given logic condition, to maintain model feasibility when considering the CVAR demand constraints.
∑
d0 e d
sla_del1ði, s, d0 Þ g
∑
d0 e d
sla_del2ði, s, d0 Þ - sla_del3ði, s, dÞ
" sjs ∈ Sp , " iji ∈ Is , " d ð62Þ 6.3.2. CVAR Demand Constraints. The first step in the development of the CVAR demand constraints is the generation of the relevant underproduction and overproduction loss functions for both the bulk and packed product classes. Note that, when generating the given loss functions, the transportation time parameter for a given facility i, center c, and scenario sc (Δ_sci, c, sc) is determined by sampling a discrete, uniform, random variable with bounds of plus or minus one day from the nominal shipment time. The bulk product demand due date parameter for a given scenario (r\_scc, s, d) is calculated in part by sampling a discrete, uniform random variable with bounds of plus
Table 4. Robust Multisite PPDM Packed Product Analysis
or minus one day from the nominal due date, to take into account demand due date uncertainty. The demand amount uncertainty is incorporated by sampling a random variable with a mean value of 1 and a standard deviation that is defined by eq 29, and the sampled value is subsequently multiplied by the nominal amount to generate the realized demand amount. For the packed products, the demand due date is considered to be deterministic in nature; therefore, only the demand amount uncertainty must be considered, by means of sampling a normal random variable in the same fashion as that in the bulk product case. To maintain computational tractability, the loss functions are written on a per product class basis, as opposed to individual states. Similarly, the bulk product underproduction loss function is defined on a per week basis, while the overproduction bulk loss function, underproduction packed loss function, and the overproduction packed loss function are calculated on a monthly basis. The bulk product underproduction loss function is determined on a weekly basis, to stress the importance of satisfying individual bulk product orders. Equations 63 and 64 define the bulk product underproduction and overproduction loss functions, respectively: fbulk_under ðc, t, scÞ ¼ -
∑i
∑
∑
∑
s ∈ Sbp ∩ Sc di e d e df t t
prices 3 r_scc, s, d
∑
s ∈ Sbp ∩ Si ∩ Sc di - Δ_sci, c, sc e d e df - Δ_sci, c, sc t t
" c, " t, " sc 4970
prices 3 Dði, c, s, dÞ
ð63Þ
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Table 5. Robust Multisite PPDM and Scheduling Comparison Production Center I1 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
Production Levels 1
971
923
4.933
535
525
1.831
2
611
573
5.436
681
647
3.614
3
798
716
7.050
1655
1622
2.690
4
698
684
5.894
925
927
1.977
5
466
393
7.167
1977
1417
11.010
6
762
639
8.764
371
702
4.960
7
588
493
9.654
798
406
10.045
1 2
57194 35008
54595 33303
4.543 4.668
51040 69089
49624 65433
2.774 4.221
3
42206
37640
6.600
155120
151193
3.269
4
41657
41435
5.164
88524
88442
2.496
5
24267
20254
6.542
188676
129129
12.422
6
44768
36261
8.818
34740
66423
6.292
7
33197
25729
10.450
80701
38343
11.874
Profit Levels
Production Center I2 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
Production Levels 1 2
782 513
763 464
2.406 5.248
640 1542
630 1085
1.437 21.371
3
747
653
7.940
888
1264
2.936
4
486
501
5.819
721
730
2.141
5
443
356
7.897
1935
1490
9.191
6
588
608
6.026
328
734
1.979
7
381
346
6.329
507
363
4.011
1 2
43577 29514
42634 27052
2.163 4.659
59468 149970
58140 106549
2.233 21.366
3
43285
37290
8.078
83510
117227
3.766
4
27819
29849
5.111
68416
69316
2.804
5
22772
17955
7.298
187124
143005
9.891
6
35682
37123
5.302
31703
71058
2.567
7
20667
18089
5.966
49073
35038
4.598
Profit Levels
Production Center I3 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
Production Levels 1 2
919 638
855 590
7.001 7.196
613 723
600 693
3
748
649
9.162
1752
1677
3.799
4
705
719
6.527
897
919
2.393
4971
2.126 3.174
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Table 5. Continued Production Center I3 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
5 6
439 776
380 725
7.403 7.253
1904 319
1573 595
7.238 2.426
7
585
464
8.888
828
462
7.345
1
50960
47739
6.322
58335
56451
3.229
2
37301
34913
6.355
74826
71665
3.788
3
40004
33582
9.379
160989
152375
4.644
4
41341
43506
5.817
87923
90566
2.883
5
24650
21512
6.694
184268
148376
8.283
6 7
44480 33092
41574 23681
6.664 9.315
30840 82026
60414 45255
2.902 7.966
Profit Levels
Figure 4. Robust Gantt chart for 13 bottleneck reactors for three-month time horizon I1. fbulk_over ðc, m, scÞ ¼
∑i -
∑
∑
s ∈ Sbp ∩ Si ∩ Sc di - Δ_sci, c, sc e d e dmf - Δ_sci, c, sc m
∑
∑
s ∈ Sbp ∩ Sc di e d e dmf m
Dði, c, s, dÞ
r_scc, s, d
" c, " m, " sc ð64Þ The underproduction loss function (fbulk_under(c, t, sc)) takes on a nonzero value when the weighted total shipment for a given week t
at center c is less than the total weighted demand for the given combination of week t, center c, and scenario sc. Each shipment variable and demand due date scenario parameter is weighted by its respective price within eq 63, to ensure the production of moreprofitable products. In an analogous fashion, the bulk demand overproduction loss function takes on a positive value when the shipment level of bulk product for a given month m and center c is greater than the monthly demand level for bulk product for a given 4972
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Figure 5. Robust Gantt chart for 13 bottleneck reactors for three-month time horizon I2.
month m, center c, and scenario sc. The terms within eq 64 are not weighted by market prices, because the inventory cost of storing product is considered uniform. The packed product loss functions can be formulated in a similar fashion, and eqs 65 and 66 define the packed product underproduction and overproduction loss functions, respectively. fpack_under ðc, m, scÞ ¼ -
∑i
p
∑
∑
∑
p
s ∈ Sp ∩ Sc dmi e d e dmf
∑
-
p
∑
∑
s ∈ Sp ∩ Si ∩ Sc dmi - Δ_sci, c, sc e d e dmf - Δ_sci, c, sc
∑ p
∑
s ∈ Sp ∩ Sc dmi e d e dmf
-
∑i
∑
∑
∑
∑
b p
∑
i
∑
s ∈ Sbp ∩ Sc di e d e dmf m
-
Dði, c, s, dÞ
∑i
prices 3 Dði, c, s, dÞ
∑
c
Dði, c, s, dÞ f
dmi - Δ_sci, c, sc e d e dm - Δ_sci, c, sc
r_scc, s, d - ξbulk_over ðc, mÞ
ð68Þ
p
∑
∑
p
s ∈ Sp ∩ Sc dmi e d e dmf
∑
prices 3 r_scc, s, d
∑
s ∈ Sp ∩ Si ∩ Sc dmi - Δ_sci, c, sc e d e dmf - Δ_sci, c, sc
- ξpack_under ðc, mÞ zpack_over ðc, m, scÞ g
prices 3 r_scc, s, d
s ∈ Sbp ∩ Si ∩ Sc di - Δ_sci, c, sc e d e df - Δ_sci, c, sc t t
zpack_under ðc, m, scÞ g
ð65Þ
r_scc, s, d
s ∈ Sbp ∩ Sc di e d e df t t
∑i s ∈ S ∑∩ S ∩ S
ð67Þ
"c, " t, " sc
" c, " m, " sc
" c, " m, " sc ð66Þ Having defined the pertinent loss functions, the lower bounds on the bulk and packed loss function underproduction and overproduction variables can be defined as shown in eqs 6770, where the right-hand side of each of the aforementioned constraints is equivalent to the pertinent demand loss function. zbulk_under ðc, t, scÞ g
-
prices 3 Dði, c, s, dÞ
" c, " m, " sc
∑i
zbulk_over ðc, m, scÞ g
prices 3 r_scc, s, d
s ∈ Sp ∩ Si ∩ Sc dmi - Δ_sci, c, sc e d e dmf - Δ_sci, c, sc
fpack_over ðc, m, scÞ ¼
- ξbulk_under ðc, tÞ
-
∑ p
∑i
p
∑
prices 3 Dði, c, s, dÞ
" c, " m, " sc
∑
s ∈ Sp ∩ Si ∩ Sc dmi - Δ_sci, c, sc e d e dmf - Δ_sci, c, sc
∑
s ∈ Sp ∩ Sc dmi e d e dmf
ð69Þ Dði, c, s, dÞ
r_scc, s, d - ξpack_over ðc, mÞ
" c, " m, " sc ð70Þ The loss function variables are then further constrained through eqs 71-74, where the left-hand side of each of the constraints represents the CVAR and the right-hand side is the allowed level of risk exposure, which, for the underproduction case, is a userspecified fraction of the price-weighted nominal demand profile summation for either the bulk or packed product classes, whereas for the overproduction case, it is a fraction of the nominal demand 4973
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Figure 6. Robust Gantt chart for 13 bottleneck reactors for three-month time horizon I3.
profile summation for either the bulk or packed product classes. The decision to weight the right-hand side of the encapsulated CVAR constraints (eqs 71 and 73) by market prices is made in accordance with the units of the given loss function under consideration. 1 ∑t ξbulk_underðc, tÞ þ ð1 - ωÞjSCj ∑t ∑sc zbulk_underðc, t, scÞ eγ ∑ ∑d prices 3 rc, s, d " c ð71Þ s ∈ Sbp ∩ Sc
1 ∑m ξbulk_overðc, mÞ þ ð1 - ωÞjSCj ∑m ∑sc zbulk_overðc, m, scÞ eγ ∑ ∑d rc, s, d " c ð72Þ
products under consideration, relative to their historical demand profiles, eqs 75-79 are included within the CVAR multisite PSPM. Equations 75 and 76 constrain the level of production for each final product to be between the maximum and minimum sampled production totals on a per customer distribution center basis for the time horizon of interest. Equation 77 mandates that at least 80% of the nominal total demand at each center for a specific product be satisfied. Equation 78 ensures that the total level of bulk product shipped to a given center c will be g90% of the aggregate demand for the bulk product class at the given center, and eq 79 provides an analogous restriction on packed product shipment profiles.
∑i d e jDj -∑ðΔ
s ∈ Sbp ∩ Sc
1 ∑m ξpack_underðc, mÞ þ ð1 - ωÞjSCj ∑m ∑sc zpack_under ðc, m, scÞ eγ ∑ ∑d prices 3 rc, s, d " c ð73Þ
,
i c
∑i d e jDj -∑ðΔ
p
,
i c
∑m
eγ
∑m ∑sc zpack_over ðc, m, scÞ
∑ ∑d rc, s, d p
s ∈ Sp ∩ Sc
"c
Dði, c, s, dÞ e
∑d max r_scc, s, d, sc sc
ð75Þ
∑d min r_scc, s, d, sc sc
ð76Þ
" c, " sjs ∈ Sp ∩ Sc
s ∈ Sp ∩ Sc
1 ξpack_over ðc, mÞ þ ð1 - ωÞjSCj
- 1Þ
- 1Þ
Dði, c, s, dÞ g
" c, " sjs ∈ Sp ∩ Sc
ð74Þ
∑i d e jDj -∑ðΔ
6.3.3. Additional Production Bounds. To ensure an even distribution of production among all the bulk and packed
,
i c
- 1Þ
Dði, c, s, dÞ g
4 5
∑d rc, s, d
ð77Þ
" c, " sjs ∈ Sp ∩ Sc 4974
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Table 6. CVAR Multisite PSPM Bulk Product Analysis
∑i
∑
∑
s ∈ Sbp ∩ Si ∩ Sc d e jDj - ðΔi, c - 1Þ
9 g 10 s ∈ Sbp ∩ Sc
∑ ∑d rc, s, d
∑i
∑
∑
p s ∈ Sp ∩ Si ∩ Sc d e jDj - ðΔi, c - 1Þ
9 g 10 s ∈ Spp ∩ Sc
∑ ∑d rc, s, d
Table 7. CVAR Multisite PSPM Packed Product Analysis
Dði, c, s, dÞ "c
ð78Þ
facilities, a unit underutilization penalty term is not needed either. " min
∑i ∑c s ∈ S ∑∩ S ∩ S ∑d Dði, c, s, dÞ 3 σ_transi, c, s p
Dði, c, s, dÞ
þ
i
c
∑i s ∈ S∑∩ S ∑d sla_del3ði, s, dÞ p
"c
#
ð80Þ
i
ð79Þ
6.3.4. Objective Function. The objective function of the CVAR multisite PSPM is shown in eq 80, which minimizes the transportation cost of shipping product, to satisfy demand at the multiple customer distribution centers, and penalizes the violation of the supplied production profile. The penalization of the supplied production profile violation implicitly takes into account the maximization of gross profit, along with minimizing underproduction, overproduction, and unit underutilization, because the supplied production profile was generated with these objectives in mind. Also, the penalization of the supplied production profile violation ensures that the final production profile generated from the CVAR multisite PSPM explicitly considers pertinent parameter uncertainty by means of both the robust optimization framework and CVAR. Note that the explicit penalization of product underproduction and overproduction is not necessary, because of the presence of the CVAR demand constraints; also, because the CVAR multisite does not model any of the units present within the multiple production
7. COMPUTATIONAL STUDIES 7.1. Robust Multisite PPDM. The robust multisite PPDM was applied to the given supply chain planning under uncertainty problem with a three-month time horizon. All of the 1593 supplied demand due parameters (rc, s, d) were considered to be uncertain. The amount of each demand due date parameter was assumed to follow a normal distribution with a mean equal to the nominal value and a standard deviation that varied linearly with the time horizon. The time horizon was divided into seven planning periods, each having a time duration of approximately two weeks, to generate a tiering scheme for the standard deviation (σd), which was used to characterize each normal random variable, encapsulating the uncertainty of a demand due date parameter’s amount (eq 29). As previously stated, the demand due date for bulk products was assumed to follow a uniform, discrete distribution with bounds of plus or minus one day from the nominal, while the demand due date for packed 4975
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Figure 7. CVAR Gantt chart for 13 bottleneck reactors for three-month time horizon I1.
products was assumed to be deterministic in nature. Each transportation time parameter was assumed to follow a uniform discrete distribution with bounds of plus or minus one day from the nominal value. Note that all presented demand totals represent the nominal or mean demand totals for each customer distribution center. These mean values were used to provide a base level of comparison between the planning production level totals and the historical demand totals present at each customer distribution center. With regard to the objective function weighting, R was assigned a value of 1000, with the remaining coefficients all being given a value of 1, in accordance with each objective function term’s relative importance. The uncertainty level parameter ɛ was set equal to 1, to express the high degree of demand due date parameter uncertainty; the degree of constraint relaxation parameter (δr) was set equal to 0.01, to maintain a stringent level of constraint feasibility; and κ1 and κ2, which encapsulate the allowable degree of constraint violation for the robust packed product demand constraints, were set equal to 0.01 and 0.10, respectively. The parameter κ1 was given a lower value when compared to κ2 because it is considered more important to satisfy the weekly packed product underproduction demand constraints, as opposed to the monthly packed product overproduction demand constraints. Note that the given planning model was implemented using GAMS52 and solved with CPLEX.53 A relative optimality tolerance of 1% was used as the termination criterion, along with a time limit of 5 h and an integer solution limit of 200 for the given planning model. The executed model is large-scale and has 607 169 equations, 725 780 continuous
variables, and 41 490 binary variables. Tables 3 and 4 present the planning level production breakdown for the packed and bulk goods, respectively. Note that all production levels presented within this work are in terms of mass units, and all profit totals are based in currency units. The planning level results indicate that the robust multisite PPDM not only explicitly takes into account the various forms of supply chain uncertainty but also generates a production profile that reflects the relative transportation costs associated with shipping product from a given production facility to a particular customer distribution center. As seen from Table 4, a small degree of packed product overproduction has been permitted, even though the nominal monthly packed product demand overproduction constraints do not allow for the generation of inventory. The reason for this occurrence is that the nominal packed product overproduction demand constraints do not consider the possibility of day-ahead shipment arrivals, which are considered within the maximization of gross profit and the robust packed product overproduction demand constraints, and therefore, a relatively small amount of inventory is generated as a result. Note that all presented demand totals represent the nominal or mean demand totals for each customer distribution center. All demands were treated as uncertain, and these mean values were used to provide a base level of comparison between the planning production level totals and the historical demand totals present at each customer distribution center. To ascertain whether the generated production profile is a tight upper bound on the true production capacity of the supply chain, the medium-term scheduling framework that was 4976
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Table 8. CVAR Multisite PSPM and Scheduling Comparison Production Center I1 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
Production Levels 1
932
891
4.381
353
321
9.035
2
690
631
6.147
1496
1101
23.074
3
597
553
6.492
782
1169
1.499
4
852
816
5.888
1072
989
3.303
5
559
491
6.857
1346
1258
4.152
6
612
586
6.481
1056
1131
2.209
7
655
591
6.932
352
247
3.709
1 2
54170 41532
52126 38009
3.774 5.817
27433 162831
24241 121193
11.638 23.562
3
30257
28622
5.718
62400
102377
1.921
4
49685
47830
5.157
99955
94509
2.921
5
31215
27807
6.027
139504
125866
4.864
6
34662
33360
5.701
90927
97411
2.993
7
42917
36388
7.136
42287
30471
4.680
Profit Levels
Production Center I2 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
Production Levels 1 2
611 728
592 674
3.080 5.457
345 1723
327 1258
5.145 23.374
3
480
440
6.216
4
709
674
5.866
447
829
4.023
1275
1259
5
481
433
6.531
802
3.109
781
3.036
6
461
370
8.293
1116
895
6.305
7
455
491
6.398
597
448
8.076
1 2
34002 42787
33059 39835
2.772 5.072
27638 173742
26061 124636
5.703 25.167
3
26807
25038
5.468
39261
79154
4.484
4
40598
38812
5.167
120723
118880
3.496
5
26617
24224
5.763
70297
62376
4.762
6
28464
21714
8.327
109976
87449
7.954
7
23725
27021
5.963
67829
44484
10.899
Profit Levels
Production Center I3 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
Production Levels 1 2
808 769
744 714
7.983 7.546
286 1666
286 1244
0.000 21.596
3
639
574
4
789
782
8.312
671
1040
1.996
6.360
1015
870
5.434
4977
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Table 8. Continued Production Center I3 Bulk period
planning
scheduling
Packed cum diff (%)
planning
scheduling
cum diff (%)
5 6
526 626
437 576
7.933 7.955
1523 969
1104 927
11.951 10.756
7
626
539
8.727
451
385
11.024
1
45305
42068
7.145
22939
22939
0.000
2
44832
41757
7.002
178571
130470
23.870
3
32668
29348
7.843
55731
97111
2.612
4
46530
46588
5.653
96401
80562
6.379
5
30033
25126
7.263
142045
91395
14.769
6 7
35590 40356
32459 32362
7.495 9.300
93980 53228
90902 40940
12.937 13.778
Profit Levels
Figure 8. CVAR Gantt chart for 13 bottleneck reactors for three-month time horizon I2.
developed by Janak et al.,54,55 which is based on the continuoustime, unit-specific, event-based, short-term scheduling model that was developed by Floudas and co-workers,56-67 was used to determine the allocation of plant resources for each production facility within the supply chain. The given medium-term scheduling model rigorously models all of the unit operations occurring within each production facility, and, as a result, it
provides an accurate representation of the supply chain’s true production capacity. The application of the scheduling model to each scheduling subhorizon entailed using an optimality tolerance of 0.001% as the termination criterion, along with a time limit of 3 h and an integer solution limit of 40, and the given medium-term scheduling model was implemented using GAMS52 and solved with CPLEX.53 For a given scheduling 4978
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
Figure 9. CVAR Gantt chart for 13 bottleneck reactors for three-month time horizon I3.
Table 9. Average Percent Planning Period Cumulative Overproduction Comparison Bulk production center
avg (%)
Packed std (%)
avg (%)
std (%)
Robust Multisite PPDM I1
6.985
1.738
5.161
3.826
I2
5.952
1.869
6.152
7.206
I3 average
7.633 6.857
0.994 1.534
4.071 5.128
2.270 4.434
I1
6.168
0.869
6.712
7.613
I2
5.977
1.559
7.581
7.194
I3
7.831
0.745
8.965
7.272
average
6.659
1.058
7.753
7.360
CVAR Multisite PSPM
subhorizon of 3 days, the medium-term scheduling model contained 6198 binary variables, 35 769 continuous variables, and 344 320 constraints. Table 5 presents the planning period production and profit totals for the planning and scheduling levels, as well as the planning-period-based percent cumulative difference between the planning and scheduling level totals for each production facility. The percent cumulative difference can be defined as the difference between the planning production/ profit total up to and including the planning period of interest and the scheduling level production/profit total up to and
including the given planning period, divided by the aforementioned planning level production/profit total; this fractional quantity is multiplied by 100 to convert it to a percentage. It is preferable that the percent cumulative difference between the planning and scheduling levels for each production facility be positive and relatively small for each planning period, to demonstrate that the supplied planning level production profile is a tight upper bound on the true production capacity of the supply chain, and Table 5 exemplifies this characteristic. Figures 4-6 also indicate that the bottleneck units present within each production facility have been sufficiently utilized during the process of attempting to satisfy the supplied production profile. 7.2. CVAR Multisite PSPM. The CVAR multisite PSPM was applied to the same three-month supply chain planning under uncertainty problem, to validate the given approach. The proposed linear programming model used the production profile generated from the robust multisite PPDM as an initialization point. The CVAR confidence level ω was set equal to 0.90, along with the CVAR risk exposure indicator (δ) being given a value of 0.15, to generate a conservative planning level production profile. The SAA algorithm was utilized with a termination criterion of g2 for each z-value, related to the encapsulated CVAR demand constraints (see eqs 71-74). Each application of the CVAR multisite PSPM considered an independent set of 100 demand profile scenarios and calculated the aforementioned z-values, using 50 000 demand profile scenarios. The CVAR multisite PSPM was implemented using GAMS,52 solved by means of CPLEX,53 and contained 107 515 continuous variables and 4979
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research 40 006 constraints. Tables 6 and 7 provide the planning level production breakdown among the three production facilities. The outstanding packed product demand for the second customer distribution center (see Table 7) indicates the influence of the aggregated CVAR demand constraints, as well as the fact that the obtained solution is only guaranteed to be a feasible upper bound on the optimal solution as per the SAA algorithm. Despite this shortcoming, the generated production profile does, for the most part, efficiently allocate production responsibilities by implicitly taking into account the production capacity of each facility by mitigating the violation of the supplied production profile and explicitly considering the relevant transportation costs. Once again, the listed demand totals represent the nominal demand totals for each customer distribution center under consideration. After generating the planning level production profile through the application of the CVAR multisite PSPM within the SAA algorithm, the medium-term scheduling framework that was developed by Janak et al.54,55 was applied using the aforementioned production profile. Note that the planning-period- and production-center-specific percent cumulative difference is used to assess the ability of the planning level production profile to provide a tight upper bound on the production capacity of the supply chain in question, and Table 8 demonstrates that the CVAR multisite PSPM is capable of accomplishing this objective. Figures 7-9 indicate that the bottleneck units are sufficiently utilized, which is another indication that the supplied production profile allows for efficient utilization of the supply chain’s production capacity. 7.3. Comparison. Table 9 provides a basis of comparison between the robust multisite PPDM and CVAR multisite PSPM, where the average and standard deviation of the percent cumulative difference on a planning period basis are presented for each of the three production facilities, along with the averages of these quantities over all three production facilities. The results indicate that both the robust multisite PPDM and the CVAR multisite PSPM are capable of providing tight upper bounds on the supply chain’s production capacity, as well as being very similar to each other on an aggregate basis. However, note that the robust multisite PPDM does not require any modeling aggregation, the application of an iterative SAA algorithm, or an initial supplied production profile. The proposed CVAR multisite PSPM does have the advantage of modeling uncertainty through both the robust optimization framework and CVAR theory to generate a unified production profile. As a result, both the robust multisite PPDM and CVAR multisite PSPM are viable alternatives toward addressing the supply chain planning under an uncertainty problem, because both approaches explicitly consider the multiple forms of system uncertainty while still providing a tight upper bound on the supply chain’s true production capacity.
8. CONCLUSION The operational planning of a multisite production and distribution network entails determining the daily production and shipment profiles for the entire supply chain. To increase the reliability of the generated profiles, the multiple forms of system uncertainty must be explicitly taken into account. Demand due date and demand amount uncertainty, along with transportation time uncertainty, are considered within the proposed operational planning model. Both the robust optimization framework and conditional value-at-risk (CVAR) theory have been utilized to address the given supply chain planning under uncertainty problem. Based on the computational results on a large-scale
ARTICLE
industrial case study, we conclude that the robust multisite production disaggregation model (PPDM) and the CVAR multisite production shifting planning model (PSPM) are capable of providing tight upper bounds on the production capacity of the supply chain while concurrently taking the aforementioned forms of system uncertainty into account. Overall, this work has demonstrated that both the robust optimization framework and CVAR theory can be successfully applied in the development of optimization models aimed at addressing the operational planning of a multisite production and distribution network under demand and transportation time uncertainty.
’ AUTHOR INFORMATION Corresponding Author
*Tel.: 1-609-258-4595. E-mail address: fl
[email protected].
’ ACKNOWLEDGMENT The authors gratefully acknowledge financial support from the National Science Foundation (CMMI-0856021). ’ REFERENCES (1) Simchi-Levi, D.; Kaminsky, P.; Simchi-Levi, E. Managing the Supply Chain: The Definitive Guide for the Business Professional; McGraw-Hill: New York, 2004. (2) Shah, N. Pharmaceutical supply chains: Key issues and strategies for optimization. Comput. Chem. Eng. 2004, 28, 929–941. (3) Shah, N. Process industry supply chains: Advances and challenges. Comput. Chem. Eng. 2005, 29, 1225–1235. (4) Varma, V. A.; Reklaitis, G. V.; Blau, G. E.; Pekny, J. F. Enterprisewide modeling and optimization—An overview of emerging research challenges and opportunities. Comput. Chem. Eng. 2007, 31, 692–711. (5) Papageorgiou, L. G. Supply chain optimization for the process industries: Advances and opportunities. Comput. Chem. Eng. 2009, 33, 1931–1938. (6) Wilkinson, S. J.; Shah, N.; Pantelides, C. C. Aggregate modeling of multipurpose plant operation. Comput. Chem. Eng. 1995, 19, S583–S588. (7) Bassett, M. H.; Pekny, J. F.; Reklaitis, G. V. Decomposition techniques for the solution of large-scale scheduling problems. Process Syst. Eng. 1996, 42, 3373–3387. (8) Kallrath, J. Diskrete Optimierung in der Chemischen Industrie; Springer: Berlin-Heidelberg, 1995. (9) Kallrath, J.; Wilson, J. M. Business Optimization Using Mathematical Programming; Macmillan: Houdmills, Basingstoke, U.K., 1997. (10) Kallrath, J.; Timpe, C. H. Optimal operational planning in large multi-site production networks. Eur. J. Oper. Res. 2000, 126, 422–435. (11) Kallrath, J. Combined strategic and operational planning—An MILP success story in chemical industry. OR Spectrum 2002, 24, 315–341. (12) Kallrath, J. Solving planning and design problems in the process industry using mixed integer and global optimization. Ann. Oper. Res. 2005, 140, 339–373. (13) Kallrath, J.; Maindl, T. I. Real Optimization with SAP-APO; Springer: Heidelberg, Germany, 2006. (14) Sung, C.; Maravelias, C. T. An attainable region approach for production planning of multiproduct processes. AIChE J. 2007, 53, 1298–1315. (15) Verderame, P. M.; Floudas, C. A. Operational planning of a multisite production and distribution network. Comput. Chem. Eng. 2009, 33, 1036–1050. (16) Li, Z.; Ierapetritou, M. Process scheduling under uncertainty: Review and challenges. Comput. Chem. Eng. 2008, 32, 715–727. 4980
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research (17) Verderame, P. M.; Elia, J. A.; Li, J.; Floudas, C. A. Planning and Scheduling under Uncertainty: A Review Across Multiple Sectors. Ind. Eng. Chem. Res. 2010, 49, 381–389. (18) Gupta, A.; Maranas, C. D. A two-stage modeling and solution framework for multisite midterm planning under demand uncertainty. Ind. Eng. Chem. Res. 2000, 39, 3799–3813. (19) Gupta, A.; Maranas, C. D.; McDonald, C. M. Mid-term supply chain planning under demand uncertainty: Customer demand satisfaction and inventory management. Comput. Chem. Eng. 2000, 24, 2613–2621. (20) Levis, A. A.; Papageorgiou, L. G. A hierarchical solution approach for multi-site capacity planning under uncertainty in the pharamaceutical industry. Comput. Chem. Eng. 2004, 28, 707–725. (21) Wu, D.; Ierapetritou, M. Hierarchical approach for production planning and scheduling under uncertainty. Chem. Eng. Process. 2007, 46, 1129–1140. (22) Colvin, M.; Maravelias, C. T. A stochastic programming approach for clinical trial planning in new drug development. Comput. Chem. Eng. 2008, 32, 2626–2642. (23) Al-Qahtani, K.; Elkamel, A. Robust planning of multisite refinery networks: Optimization under uncertainty. Comput. Chem. Eng. 2010, 34, 985–995. (24) Clay, R. L.; Grossmann, I. E. A disaggregation algorithm for the optimization of stochastic planning models. Comput. Chem. Eng. 1997, 21, 751–774. (25) Pistikopoulos, E. N.; Dua, V. Planning under uncertainty: A parametric optimization approach. In Proceedings of Third International Conference on Foundations of Computer-Aided Process Operations; Pekny, J. F., Blau, G. E., Eds.; American Institute of Chemical Engineers (AIChE): New York, 1998; pp 164-169. (26) Ryu, J.; Dua, V.; Pistikopoulos, E. N. A bilevel programming framework for enterprise-wide process networks under uncertainty. Comput. Chem. Eng. 2004, 28, 1121–1129. (27) Charnes, A.; Cooper, W. W. Chance-constrained programming. Manage. Sci. 1959, 6, 73–79. (28) Charnes, A.; Cooper, W. W. Chance constraints and normal deviates. J. Am. Stat. Assoc. 1962, 57, 134–148. (29) Charnes, A.; Cooper, W. W. Deterministic equivalents for optimizing and satisficing under chance constraints. Oper. Res. 1963, 11, 18–39. (30) Petkov, S. B.; Maranas, C. D. Multiperiod planning and scheduling of multiproduct batch plants under demand uncertainty. Ind. Eng. Chem. Res. 1997, 36, 4864–4881. (31) Li, P.; Garcia, H. A.; Wozny, G. Chance constrained programming approach to process optimization under uncertainty. Comput. Chem. Eng. 2008, 32, 25–45. (32) Mitra, K.; Ravindra, G. D.; Patwardhan, S. C.; Sardar, G. Midterm supply chain planning under uncertainty: A multiobjective chance constrained programming framework. Ind. Eng. Chem. Res. 2008, 47, 5501–5511. (33) Ben-Tal, A.; Goryashko, A.; Guslitzer, E.; Nemirovski, A. Adjustable robust rolustions for uncertain linear programs. Math. Program. A 2004, 99, 351–376. (34) Ben-Tal, A.; Nemirovski, A. Robust solutions of uncertain linear programs. Oper. Res. Lett. 1999, 25, 1–13. (35) Ben-Tal, A.; Nemirovski, A. Robust solutions of linear programming problems contaminated with uncertain data. Math. Program. A 2000, 88, 411–425. (36) Ghaoui, L. H. E.; Oustry, F.; Lebret, H. Robust solutions to uncertain semidefinite programs. SIAM J. Opt. 1998, 9, 33–52. (37) Ghaoui, L. H. E. Robust solutions to least-square problems with uncertain data. SIAM J. Matrix Anal. Appl. 1997, 18, 1035–1064. (38) Bertsimas, D.; Pachamanova, D.; Sim, M. Robust linear optimization under general norms. Oper. Res. Lett. 2004, 32, 510–516. (39) Bertsimas, D.; Sim, M. Robust discrete optimization and network flows. Math. Program. B 2003, 98, 49–71. (40) Bertsimas, D.; Sim, M. The price of robustness. Oper. Res. 2004, 52, 35–53.
ARTICLE
(41) 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. (42) 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–195. (43) Verderame, P. M.; Floudas, C. A. Operational planning of largescale industrial batch plants under demand due date and amount uncertainty: I. Robust optimization. Ind. Eng. Chem. Res. 2009b, 48, 7214–7231. (44) 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. 2010, 49, 260–275. (45) Artzner, P.; Delbaen, F.; Eber, J. M.; Heath, D. Coherent measure of risk. Math. Finance 1999, 9, 203–228. (46) Ogryczak, W. Stochastic Dominance Relation and Linear Risk Measures. In Financial Modelling—Proceedings of the 23rd Meeting of the Euro Working Group; Progress and Business Publishers: Krakow, Poland, 2000. (47) Rockafellar, R.; Uryasev, S. Optimization of conditional valueat-risk. J. Risk 2000, 2, 21–43. (48) Rockafellar, R.; Uryasev, S. Conditional value-at-risk for general loss distributions. J. Banking Finance 2002, 26, 1443–1471. (49) Carneiro, M. C.; Ribas, G. P.; Hamacher, S. Risk management in the oil supply chain: A CVaR approach. Ind. Eng. Chem. Res. 2010, 49, 3286–3294. (50) 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. (51) Wang, W.; Ahmed, S. Sample average approximation of expected value constrained stochastic programs. Oper. Res. Lett. 2008, 36, 515–519. (52) Brooke, A.; Kendrick, D.; Meeraus, A.; Raman, R. GAMS: A User’s Guide; GAMS: San Francisco, CA, 2003. (53) ILOG CPLEX 9.0 User’s Manual; ILOG, Inc.: Mountain View, CA, 2005. (54) 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. (55) 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–8269. (56) 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. (57) 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–4374. (58) 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–3461. (59) Lin, X.; Floudas, C. Design, synthesis and scheduling of multipurpose batch plants via an effective continuous-time formulation. Comput. Chem. Eng. 2001, 25, 665–674. (60) 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. (61) 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. (62) 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. 4981
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982
Industrial & Engineering Chemistry Research
ARTICLE
(63) 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. (64) 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. (65) 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–686. (66) Shaik, M. A.; Floudas, C. A. Novel unified modeling approach for short-term scheduling. Ind. Eng. Chem. Res. 2009, 48, 2947–2964. (67) Li, J.; Floudas, C. A. Optimal Event Point Determination for Short-Term Scheduling of Multipurpose Batch Plants. Ind. Eng. Chem. Res. 2010, 49, 7446–7469.
4982
dx.doi.org/10.1021/ie101401k |Ind. Eng. Chem. Res. 2011, 50, 4959–4982