Multisite Planning under Demand and ... - ACS Publications

Jump to Robust Optimization Framework - Before presenting the robust multisite planning model, an overview of the robust optimization framework is ...
0 downloads 0 Views 4MB Size
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



∑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