3782
Ind. Eng. Chem. Res. 2004, 43, 3782-3791
Short-Term Scheduling under Uncertainty Using MILP Sensitivity Analysis Zhenya Jia and Marianthi G. Ierapetritou* Department of Chemical and Biochemical Engineering, Rutgers University, Piscataway, New Jersey 08854
The aim of this paper is to develop an integrated framework within which to address the issue of uncertainty in short-term scheduling. The idea of inference-based sensitivity analysis for MILP problems is utilized within a branch-and-bound solution framework to determine the importance of different parameters and constraints and to provide a set of alternative schedules for the range of uncertain parameters under consideration. Two case studies are presented to illustrate the application of the proposed approach and to provide a comparison with parametric programming and robust optimizaton studies. 1. Introduction In the literature, a significant amount of work has been done in the area of short-term scheduling, which involves the determination of the order in which tasks use units and various resources and the detailed timing of the execution of all tasks so as to achieve the desired performance. Most of the proposed approaches can be classified into two broad categories based on the time representation they adopt. Early attempts were mainly concentrated on the discrete-time or time-indexed formulation, where the time horizon is divided into a number of intervals of equal duration. Following this approach leads to large-scale mixed integer linear programming (MILP) problems involving large number of binary variables.1 In recent years, a number of publications have focused on developing efficient methods based on a continuous-time representation. Two alternative continuous-time approaches have been developed: the first representation requires that each task take place only between two time points, as proposed by Schilling and Pantelides2 and Mockus and Reklaitis,3 whereas the other representation requires only that each task start at a time point, as proposed by Castro et al.4 An eventpoint-based approach was proposed by Ierapetritou and Floudas5,6 and applied to the scheduling of multipurpose batch and continuous processes. Most recently, Maravelias and Grossmann7 presented a novel continuoustime MILP model for general multipurpose batch plants scheduling that accounts for batch mixing/splitting, resource constraints, variable batch sizes and processing times, and various storage policies. In all of these studies, the problem data are assumed to be deterministic. However, in real plants, parameters such as raw material availability, processing times, and market requirements vary with respect to time and are often subject to unexpected deviations. Therefore, the consideration of uncertainty in scheduling problems becomes of great importance in preserving plant feasibility and viability during operations. Although a large number of papers have addressed uncertainty in process design, much less attention has been devoted to the issue of uncertainty in process planning and scheduling, mainly because of the increased complexity of the deterministic problem. Among * To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: (732)445-2971. Fax: (732)445-2421.
the works on this subject appearing in the literature is that of Shah and Pantelides,8 which addresses the problem of design of multipurpose batch plants considering different schedules for different sets of production requirements using a scenario-based approach9 and an approximate solution strategy. Ierapetritou and Pistikopoulos10 presented a two-stage stochastic programming formulation for the problem of batch plant design and operations under uncertainty. Straub and Grossmann11 proposed the idea of the stochastic flexibility index to evaluate the effect of uncertainty quantitatively. Bansal et al.12 proposed a parametric programming framework for the flexibility analysis design of linear systems, and in their subsequent work, they generalized and unified this approach is to handle nonlinear systems. The multiperiod planning and scheduling of multiproduct plants under demand uncertainty was addressed by Petkov and Maranas.13 The stochastic elements in their proposed model are expressed with equivalent deterministic forms, resulting in a convex mixed integer nonlinear programming (MINLP) problem. Schmidt and Grossmann14 considered the optimal scheduling of new product testing tasks and reformulated the initial nonlinear, nonconvex disjunctive model as an MILP using different sets of simplifying assumptions that give rise to different models. The uncertainties in planning and scheduling problems are generally described through probabilistic models. During the past decade, fuzzy set theory has been applied to scheduling optimization using heuristic search techniques.15,16 Recently, Balasubramanian and Grossmann17 developed MILP models for flowshop scheduling and new product development processing scheduling, based on a fuzzy representation of uncertainty. Daniels and Carrillo18 addressed the problem of β-robust scheduling in single-stage production facilities with uncertain processing times. Vin and Ierapetritou19 proposed a multiperiod programming model to improve the schedule performance of batch plants under demand uncertainty. Acevedo and Pistikopoulos20 addressed linear process engineering problems under uncertainty using a branch-and-bound algorithm, based on the solution of multiparametric linear programs at each node of the tree and the evaluation of the uncertain parameters space for which a node must be considered. However, most of the existing approaches can handle only a certain type of uncertain parameters, mostly uncertainty in product demands, and more importantly,
10.1021/ie0306731 CCC: $27.50 © 2004 American Chemical Society Published on Web 02/26/2004
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3783
the additional complexity makes them infeasible for realistic applications. In this paper, a novel framework is proposed for uncertainty analysis of scheduling problems based on the ideas of sensitivity analysis of the corresponding MILP problem, which is addressed using a branch-and-bound solution methodology. The paper is organized as follows: The deterministic model for batch plant scheduling that is adopted in this work is briefly described in the next section. In section 3, the basic background of parametric programming, sensitivity analysis, and robust optimization is provided, followed by the proposed approach for scheduling under uncertainty. Section 4 is then utilized to present the results of two examples and compare them with robust and parametric optimization studies. Finally, the work is summarized in section 5, and some future directions are discussed.
is considered, U ) H in constraints 8-10. In general, the objective function is to minimize the makespan as shown in expression or to maximize the total profit. Allocation constraints 2 state that only one of the tasks can be performed in each unit at an event point (n). Constraints 3 represent the condition that the material balance for each material at each event point n is equal to that at event point n - 1, adjusted by any amounts produced and consumed between event points n - 1 and n and delivered to the market at event point n. The storage and capacity limitations of production units are expressed by constraints 4 and 5, respectively. Constraints 6 are written to satisfy the demands of final products. Constraints 7-14 represent time limitations due to task duration and sequence requirements in the same or different production units. Parameters R(i,j) and β(i,j) are defined as follows
2. Deterministic Scheduling Formulation In this paper, the mathematical model used for batch plant scheduling follows the main ideas of the continuous-time formulation proposed by Ierapetritou and Floudas.5 The model involves the following constraints
2 R(i,j) ) T h (i,j) 3
minimize H or maximize
∑s ∑n price(s) d(s,n)
(1)
subject to
wv(i,j,n) e 1 ∑ i∈I
(2)
j
st(s,n) ) st(s,n - 1) - d(s,n) + Fp(s,i) b(i,j,n - 1) +
∑ i∈I j
∑ j∈J
i
Fc(s,i) ∑ b(i,j,n) ∑ i∈I j∈J j
(3)
i
st(s,n) e stmax(s)
(4)
2 h (i,j)/[Vmax(i,j) - Vmin(i,j)] β(i,j) ) T 3 where T h (i,j) is the mean processing time of task i in unit j. These expressions are based on the assumption that there is 33% variability of the processing time around the mean value, although different processing times can be easily adapted. When wv(i,j,n) ) 0, the last two terms in constraints 7 are equal to zero because of capacity constraints. Otherwise, the last two terms are added to Ts(i,j,n). Therefore, the duration of task i in unit j at event point n depends on the amount of material being processed. Extensive testing of this model is reported in the literature7 that compares this model favorably to other existing discrete- and continuous-time formulations for batch plant scheduling. 3. Scheduling under Uncertainty
Ts(i,j,n + 1) g Tf(i,j,n) - U[1 - wv(i,j,n)]
(8)
Ts(i,j,n) g Tf(i′,j,n) - U[1 - wv(i′,j,n)]
(9)
Ts(i,j,n) g Tf(i′,j′,n) - U[1 - wv(i′,j′,n)]
(10)
Ts(i,j,n + 1) g Ts(i,j,n)
(11)
Although a large number of papers in the literature deal with uncertainty issues concerning process design and production planning, as reported in the Introduction, the issue of uncertainty is not well studied for scheduling problems mainly because of the high complexity of the deterministic case. 3.1. Parametric Programming. In the past few years, increasing interest in parametric programming has arisen, mainly focusing on integer programming problems including linear and nonlinear formulations.21,22 In this section, a brief description of the basic approach proposed by Pertsinidis et al.23 is provided that covers the right-hand side scalar variations for MILP problems. This approach was selected for its simplicity for comparison with the results obtained from the proposed uncertainty analysis framework. The MILP problem considered has the following form
Tf(i,j,n + 1) g Tf(i,j,n)
(12)
z(θ) ) min cTx + dTy
Tf(i,j,n) e H
(13)
Ts(i,j,n) e H
(14)
Vmin(i,j) wv(i,j,n) e b(i,j,n) e Vmax(i,j) wv(i,j,n) (5)
∑n d(s,n) g r(s)
(6)
Tf(i,j,n) ) Ts(i,j,n) + R(i,j) wv(i,j,n) + β(i,j) b(i,j,n) (7)
s.t.
where U denotes an upper bound of the makespan, for the cases where the objective is the minimization of the makespan. For the cases where maximization of profit
Ax + Dy e b(θ) xL e x e xU L
(15)
U
θ eθeθ
x ∈Rm, y∈(0, 1)t where b(θ) ) b0 + θr represents the vector of uncertain
3784 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
parameters, b0 and r are constants, and scalar θ varies in the interval [θL, θU]. First, the problem is solved at point b ) b0 + θLr, and the optimal solution is obtained to be (x*, y*) with objective value zL. The aim is then to identify a break point, θ′, beyond which the current integer solution is no longer optimal. By performing an LP sensitivity analysis on the above formulation with fixed binary variables y*, one can observe that the optimal value z varies with respect to θ as a linear function
z(θ) ) zL + λθ, ∀θL e θ e θU′ e θU
(16)
where [ θL, θU′] is the sensitivity range of the current basis. By including the above equality, which should hold at the break point, and an integer cut to exclude the current optimal solution, the following MILP problem is solved to obtain θ′ and a new solution (x′*, y′*).
z(θ) ) min cTx + dTy s.t.
Ax + Dy - θr e b0
∑ yi - ∑ yi e | F1| - 1
i∈F1
Bertsimas et al.26 Their method provides a natural generalization of the Farkas lemma for IP problems and leads to a method of performing sensitivity analysis. A method of sensitivity analysis for mixed integer linear programming based on the idea of inference duality was presented by Dawande and Hooker.27 It reveals that any perturbation that satisfies a certain system of linear inequalities will reduce the optimal value by no more than a prespecified amount. The inference-based sensitivity analysis consists of two parts: dual analysis, which determines how much the problem can be perturbed while keeping the objective function value in a certain range, and primal analysis, which gives an upper bound on how much the objective function value will increase if the problem is perturbed by a certain amount. The dual solution is obtained by using inference methods to generate constraints at every node that are violated by the branching cuts. The dual solution can be viewed as a proof of optimality and can be utilized to determine the parameter perturbations under which the dual solution still provides a valid proof. More specifically, the main results of the inference-based sensitivity analysis are summarized below. Consider the general mixed integer problem
i∈F0
cTx + dTy - zL - λθ ) 0
minimize (17)
xL e x e xU L
subject to Ax g a, 0 e x e h, xj integer, j ) 1, ..., k (18)
U′
θ eθeθ m
z ) cx
t
x∈R , y∈(0, 1)
where F1 and F0 denote the index sets of the binary variables with values equal to 1 and 0, respectively, at the current solution and |F1| is the cardinality of F1. The next break point and the successor optimal schedule can be obtained through the above procedure within the interval [θ′, θU]. 3.2. MILP Sensitivity Analysis. The formulation presented in section 2 corresponds to a mixed integer linear programming (MILP) problem where the binary variables wv(i,j,n) denote the assignment of tasks i to units j at event points n throughout the time horizon. Therefore, the effects of operating parameters on the plant performance can be investigated through the sensitivity analysis of the MILP model of the deterministic scheduling problem. Although sensitivity analysis theory is well developed in linear programming, efforts are still being made to handle the integer programming case, mainly because of the lack of optimality criteria for the integer optimization problems. Schrage and Wolsey24 examined the effect of a small perturbation on the right-hand side or objective function coefficients in an integer program by collecting dual information at each node of the branchand-bound tree while solving the original integer program and using a recursive scheme to obtain an upper bound on the objective function. Their results were extended to nonlinear integer programming problems by Skorin-Kapov and Granot.25 Pertsinidis et al.23 developed a sensitivity analysis algorithm for parametric MILP programs that also provides the results of the MILP master problem for MINLP case, while employing an algorithm that provides a sequence of improving parametric lower and upper bounds for parametric NLP subproblems. An algebraic geometry algorithm for solving integer programming problems was presented by
and assume a perturbation of all problem parameters to give
minimize
z ) (c + ∆c)x
subject to (A + ∆A)x g a + ∆a
(19)
0 e x e h, xj integer, j ) 1, ..., k The constraint z g z* - ∆z remains valid for the perturbations ∆A and ∆a in the parameters appearing on the left- and right-hand sides of the constraints if there exist sp1, ..., spn that satisfy the following set of inequalities n
λpi
n
j pj - u pj ) - λpi ∆ai e rp Aiu pj + ∑spj (u ∑ j)1 j)1 spj g λpi ∆Aij, spj g -qpj , j ) 1, ..., n
(20)
n
qpj u j pj + λpa - zp + ∆zp ∑ j)1
rp ) -
Likewise, this constraint remains valid for a perturbation ∆c of the coefficients of the objective function if there exist sp1, ..., spn that satisfy the following set of inequalities n
∆cj u pj ∑ j)1
- spj (u j pj - u pj ) g -rp (21)
spj g -∆cj, spj g -qpj , j ) 1, ..., n where qpj ) λpi Aij - λp0cj; p corresponds to the leaf node where the dual variable of the objective function (λp0) equals 1; u pj and u j pj denote the lower and upper bounds, respectively, of xj at node p; zp is the objective
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3785
value at node p; and ∆zp ) z* - zp. Leaf nodes are the nodes at which the branch-and-bound procedure terminates on the basis of standard fathoming criteria.28 Thus, by utilizing constraints 20 and 21 in the scheduling problem, the range of parameters where the objective remains within certain limits can be identified and used to evaluate alternative schedules at the branchand-bound tree. Morever, the importance of different constraints and parameters is obtained and can be utilized to improve future plant operability. More details on this subject are given in section 3.4. 3.3. Robust Optimization. To improve the schedule flexibility prior to its execution, it is important to measure the performance of a deterministic schedule under changing conditions due to uncertainty. The standard deviation (SD) is one of the most commonly used metrics for evaluating the robustness of a schedule. To evaluate the SD, the deterministic model with a fixed sequence of tasks [wv(i,j,n)] is solved for different realizations of uncertain parameters that define the set of scenarios k that results in different makespans Hk. The SD is then defined as
SDavg )
x
∑k
∑
Hk (Hk - Havg)2 k , Havg ) (22) ptot (ptot - 1)
where Havg is the average makespan over all scenarios and ptot denotes the total number of scenarios. A detailed discussion of different robustness metrics can be found in Samsatli et al.29 Vin and Ierapetritou19 proposed a robustness metric that takes the infeasible scenarios into consideration. In case of infeasibility, the problem is solved to meet the maximum demand possible by incorporating slack variables in the demand constraints. Then, the inventory of all raw materials and intermediates at the end of the schedule are used as initial conditions in a new problem with the same schedule to satisfy the unmet demand. The makespan under infeasibility (Hcorr) is determined as the sum of those two makespans. Vin and Ierapetritou’s proposed robustness metric is defined as19
SDcorr )
x∑ k
(Hact - Havg)2 (ptot - 1)
(23)
where Hact ) Hk if scenario k is feasible and Hact ) Hcorr if scenario k is infeasible. The idea of robust optimization (RO) was developed by Mulvey et al.30 to handle the trade-off associated with solution and model robustness. A solution to an optimization problem is considered to be solution robust if it remains close to the optimal solution for all scenarios; the solution is model robust if it remains feasible for most scenarios. To test the robust otimization approach, we utilize the robust metric proposed by Vin and Ierapetritou19 (SDcorr), together with the concept of the upper partial mean (UPM) introduced by Ahmed and Sahinidis.31 They define the variance measure δ h as follows
δ h)
∑k pkδk,
δk ) max(0, Hk -
∑k pkHk)
(24)
where δk corresponds to the positive deviation of the makespan under scenario k from the expected value. δ h
is used to penalize the variability of the objective function utilizing a penalty coefficient of λ. Following these ideas, the proposed robust optimization formulation has the form
minimize
∑k pkHk + λ∑k pkδk + ω∑k pk∑s slackk(s) subject to
wv(i,j,n) e 1 ∑ i∈I
(25) (26)
j
stk(s,n) ) stk(s,n - 1) - dk(s,n) + Fp(s,i) bk(i,j,n - 1) + Fc(s,i)
∑ i∈I j
∑ j∈J
∑ i∈I
i
j
bk(i,j,n) ∑ j∈J
(27)
i
stk(s,n) e stmax(s)
(28)
Vmin(i,j) wv(i,j,n) e bk(i,j,n) e Vmax(i,j) wv(i,j,n) (29)
∑n dk(s,n) + slackk(s) g r(s)
(30)
Tfk(i,j,n) ) Tsk(i,j,n) + R(i,j) wv(i,j,n) + β(i,j) bk(i,j,n) (31) Tsk(i,j,n + 1) g Tfk(i,j,n) - U[1 - wv(i,j,n)]
(32)
Tsk(i,j,n) g Tfk(i′,j,n) - U[1 - wv(i′,j,n)]
(33)
Tsk(i,j,n) g Tfk(i′,j′,n) - U[1 - wv(i′,j′,n)] (34) Tsk(i,j,n + 1) g Tsk(i,j,n)
(35)
Tfk(i,j,n + 1) g Tfk(i,j,n)
(36)
Tfk(i,j,n) e Hk
(37)
δk g Hk -
∑k pkHk,
δk g 0, Hk e U
(38)
The artificial variables in demand constraints 30 are introduced to measure demand infeasibility, which is penalized in objective function 25 using a penalty coefficient ω. In the case of uncertain demands or an objective that involves the maximization of profit, the infeasibility term is not included because, in such cases, the demand is always satisfied so as to maximize the profit and minimize the artificial variables. Note that the size of the problem, in terms of the number of continuous variables and constraints, increases in proportion to the number of scenarios considered. The implementation of this approach and a comparison of the results with those from the proposed approach are presented in section 4 for an example problem. 3.4. Proposed Uncertainty Analysis Approach. The basic idea of the proposed approach is to utilize the information obtained from the sensitivity analysis of the deterministic solution to determine (a) the importance of different parameters and constraints and (b) the range of parameters in which the optimal solution remains unchanged. The main steps of the proposed approach are shown in Figure 1. More specifically, the proposed analysis consists of two parts. In the first part, important information about the effects of different parameters is extracted following the sensitivity analy-
3786 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
Figure 1. Flowchart of proposed approach.
sis step, whereas in the second part, alternative schedules are determined and evaluated for different uncertainty ranges. First, the deterministic scheduling problem is solved at the nominal values using a branch-and-bound solution approach, and the dual multipliers λp are collected at each leaf node p. Then, the inference-based sensitivity analysis as described in the previous section is performed for all of the important scheduling parameters, including demands, prices, processing times, and capacities. Note that only the dual information of the nodes that correspond to nonzero dual variables is required. Using the results of this analysis, a number of very interesting questions can be answered regarding the plant robustness to parameter changes. In particular, these questions include: How do the capacities of the units affect the production objective? What is the range of product demand that can be covered, and how much the profit would be affected by such changes? What is the effect of a price change on the objective value? What is the significance of the constraints involved in the model? Are there any redundant sets of constraints? The first question can be answered by imposing the same perturbation to capacity constraints 5 for the different units involved in the production of specific products and determining the change in the objective function (∆z). The unit with the largest effect on the objective function value is also the most critical unit for the production of this product, and thus, a change in its capacity will result in the largest production change. Similarly, the rest of the questions can be answered by analyzing perturbations of the appropriate constraints together with the effects in the objective function. Results for two examples are given in the next section. In the second part of the analysis, the sensitivity information is used to define the range of uncertain parameters where the schedule is optimal and to identify alternative schedules in different uncertainty ranges. Sets of constraints 20 and 21 are used to determine the range of uncertain parameters for certain changes in the objective function. The branch-and-bound procedure is then continued on the nodes with the
Figure 2. State-task network representation for example 1.
objective value within the predicted limits to identify new optimal solutions. The alternative schedules are evaluated using the robustness metric (SDcorr) as defined in section 3.3 and the average and the nominal schedule performance in terms of the objective function. Because the entire analysis is based on a single tree among a large number of possible branch-and-bound trees that can be used to solve the MILP, it provides conservative sensitivity ranges. Theoretically, the exact sensitivity ranges can be obtained by investigating an exponential number of branch-and-bound trees. However, using the above analysis, useful information can be extracted regarding the approximate range of parameter change for a certain objective function change, the plant robustness to parameter changes, and the importance of different parameters, as illustrated in the next section. 4. Illustrative Examples 4.1. Proposed Approach. Example 1 involves a single-product production line consisting of a mixer, a reactor, and a purifier, as shown in Figure 2. The data for this example were obtained from Ierapetritou and Floudas.5 First, the problem is solved with the objective of maximizing the production within the time horizon of 8 h, and the dual informaton (λp) is collected at each leaf node. The pertubations of capacities, demands, and prices can be viewed as ∆A, ∆a, and ∆c, respectively, in constraint sets 20 and 21, while the production change corresponds to ∆z. For a given magnitude of production change ∆z, one of these two linear systems (depending on which parameter is examined) is solved at each leaf node. By checking whether states s1, s2,..., sn, that satisfy the linear system at each leaf node exist, we can determine the maximum allowable perturbation. The following important findings are revealed: First, by perturbing the capacity constraints, it was found that the most critical unit of the production line is the
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3787
Figure 3. Branch-and-bound tree for example 1 with uncertain demand.
purifier. If the purifier capacity is increased, production will increase accordingly. The same information can be used to determine the allowable purifier size reduction that will result in a limited production decrease. For example, it was found that the purifier capacity can decrease by up to 4.24 units without reducing the profit by more than 5%. Important sensitivity information can be obtained by a simple primal analysis of the MILP problem that results in an upper bound on the optimal value given a perturbation of the problem data. In this way, it can be determined that a 10% increase in purifier size from 50 to 55 units will result in a maximum of a 2.75% increase in production profit from 71.16 to 73.12 units (upper bound on the objective function). The results of the sensitivity analysis can also be used to analyze the effects of product price on the objective function to investigate the relative importance of different products, as well as to make decisions regarding the production levels. For this example, it was found that there is one-to-one correspondence between a price increase and an objective function increase. In particular, a 5% price increase would result in a corresponding 5% increase in the objective function. The question of how demand fluctuations affect the profit function and the production time was also analyzed. To consider this factor, the required dual information involves the dual variables associated with the demand constraint at the nodes at which this constraint is active. It was found that a demand change from 60 to 65 units (8.3% increase) can be accommodated within the production line, resulting in an equivalent increase of production profit but also in an increased (by 2.6%) makespan, to 11.59 from 11.29 h. For the second part of the analysis, first, the demand is considered to be the uncertain parameter varying within the range of [20, 100]. The problem is solved at a nominal demand (50) with the objective function of
Table 1. Comparison of Alternative Schedules for Example 1 Considering Uncertain Demand Hnom (h) Havg (h) SDcorr
schedule 1
schedule 2
schedule 3
9.83 14.20 5.52
10.77 11.56 1.61
10.91 11.79 2.17
the minimum makespan. A branch-and-bound tree is constructed as illustrated in Figure 3, where the value next to each node indicates the objective value of the relaxed LP and the binary variables that are branched at each level are shown at the right side. Applying the inference duality sensitivity analysis, the following expression is obtained regarding the range of demand change following a specific objective function change (∆H)
-0.0966∆d e ∆H which means that, if the demand is increased by ∆d, the new makespan becomes at most Hnom + 0.0966∆d. When the demand is increased to 80, the current schedule becomes infeasible. Using the previous inequality, we get H′ e 12.73. By solving the LP problem at the leaf nodes with the perturbed demand, the remaining nodes are nodes 1 and 2, with objective values 12.13 and 12.73, respectively. Therefore, we continue the branch-and-bound procedure on these two nodes and find optimal schedule 2 and feasible schedule 3. After the schedules are determined to cover the uncertainty range under consideration, the schedules are evaluated on the basis of the corrected SD robustness metric defined in eq 14 considering five scenarios in the demand interval.20,100 The average makespan, standard deviation, and nominal makespan are reported in Table 1. Note that, although schedule 1 provides the best solution at the nominal point, it has poor performance over the entire interval. Schedule 2 apparently
3788 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004
Figure 4. Branch-and-bound tree for example 1 with uncertain processing time.
has the highest robustness and provides the shortest average makespan. Schedule 3 provides a compromise solution with robustness between those of schedules 1 and 2 and average and nominal performances close to the best values obtained with schedule 2. Next, the processing time of mixing, which is denoted as R(i,j), is considered as uncertain parameter in the deterministic model. The branch-and-bound tree is built as depicted in Figure 4 at the nominal point R(i1,j1) ) 3.0 with the objective of maximizing the profit, and schedule 1 is obtained as the optimal solution. Note that negative values appear because, to apply the dual analysis, the maximization objective is transformed to a minimization problem. When R(i1,j1) is increased to 4.0, the current schedule becomes infeasible. The new objective value is estimated to be
-Profit′ e -Profitnom + 24.48∆R ) -47.04 By investigating the nodes that follow nodes 1-3, the optimal solution is found to be the equivalent schedules 2 and 4, with schedule 3 as a feasible solution. Equivalent schedules are schedules that have different values for binary variable sets but correspond to the same scheduling sequence. This is due to model redundancy, which enables the determination of a variety of scheduling alternatives. These schedules are then evaluated with respect to their robustness within the uncertainty range of R(i1, j1) [3.0, 4.0], and the results are listed in Table 2. The optimal schedule at the nominal point has the largest SD, whereas schedule 2 gives the best solution robustness. Schedule 3 provides a better average profit compared to schedule 2, so it provides an attractive alternative. Example 25 considers two different products produced through five processing stages: heating, reactions 1-3, and separation of product 2 from impure E, as il-
Table 2. Comparison of Alternative Schedules for Example 1 Considering Uncertain Processing Time -Profitnom -Profitavg SDcorr
schedule 1
schedule 2
schedule 3
-71.52 -66.98 26.9
-65.27 -64.61 9.33
-65.27 -65.17 10.49
lustrated in the state-task network (STN) representation in Figure 5. For the first part of the analysis, the problem is solved with the objective of maximizing the profit within the time horizon of 12 h. After the sensitivity analysis is performed, the following information is obtained: It is found that the most critical task of the production line is reaction 2. If the processing capacity of reaction 2 in reactor 1 or reactor 2 is decreased by 11 units, the profit will be reduced by 5%, In contrast, very little or no change is observed in the objective function upon a change in the reaction 1 or reaction 3 processing capacity in both reactors. The objective value is also not sensitive to the change of other parameters; for example, the processing capacity of separation in the separator can drop by up to 120 units without decreasing the profit. Another important modeling issue that can be addressed is the question of constraint redundancy. Here, the importance of storage constraints is investigated, and it is found that these constraints are redundant because they are not active in any of the solution branch-and-bound nodes. More interestingly, the duration constraints are also found to be redundant, which means that the maximum processing capacities are already reached with the current processing times, so that the profit cannot be improved even with zero processing times assuming a fixed number of event points. For the second part of the analysis, the demand of product 2 is considered to be the uncertain parameter, varying within the range of [20, 80], and the objective function is modified to minimize the makespan. A
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3789
Figure 5. State-task network representation for example 2. Table 3. Comparison of Alternative Schedules for Example 2 Hnom (h) Havg (h) SDcorr
schedule 1
schedule 2
schedule 3
7.00 8.15 2.63
7.14 7.24 0.29
7.40 7.40 0.27
branch-and-bound tree is constructed at the nominal point r(p2) ) 50, and the dual information is stored at each node. Applying the inference duality sensitivity analysis, the following expression is obtained regarding the range of demand change following specific objective change (∆H)
-0.0297∆d e ∆H which means that, if the demand is increased by ∆d, the new makespan becomes at most Hnom + 0.0297∆d. When r(p2) is increased from 50 to 80, schedule 1 becomes infeasible. Then, we solve the LP problem at each leaf node with the demand of 80 and check the leaf nodes with the objective function values below the 7.89 value that is obtained using this inequlaity in the branch-and-bound tree. The new optimal solution is found to be schedule 2, and schedule 3 is one feasible solution. The schedules are then evaluated with respect to the mean and nominal makespan and standard deviation within the demand range of [20, 80], and the values are shown in Table 3. Compared with schedule 2, schedule 3 has a larger mean makespan but a lower standard deviation, which means higher robustness; therefore, depending on the decision maker’s attitude toward risk and the expected demand growth, one can choose schedule 3 over schedule 2, whereas schedule 1 remains a valid alternative if the demand is expected to remain constant. Note that, although example 2 is larger and more complex than example 1, the proposed uncertainty analysis does not substantially increase the problem complexity. That is due to the fact that the required information is already obtained from the solution of the deterministic problem.
4.2. Comparisons with Robust Optimization and Parametric Programming. Example 1 with demand as the uncertain parameter is solved using the robust optimization model as introduced in section 3.3. Note that only schedule 1 in Figure 3 can be obtained because any feasible schedule for that formulation has to satisfy the worst-case scenario, which is demand ) 100. For the same example, considering uncertain processing time and using robust optimization, schedule 3 in Figure 4 is found with the penalty coefficient λ in the range of [0, 1.6], and schedule 2 is obtained for λ in the interval of [1.6, +∞], but the intermediate option that corresponds to schedule 1 is not found. The robust optimization model is also applied to example 2, and schedules 1 and 2 can be obtained with different values of the penalty coefficient λ. However, schedule 3 is not found. Example 2 is used as a test problem of applying parametric programming23 and compared with the proposed uncertainty analysis. For the case in which the demand of product 2 is considered uncertain, the problem is first solved at the nominal point r0(p2) ) 50, which gives rise to the following demand constraint
∑n d(p2,n) g r0(p2)
(39)
The optimal solution obtained is schedule 1, as reported in Table 4. The binary variables are then fixed, and the LP sensitivity analysis with respect to r(p2) results to the following constraints for time horizon H
H ) H0 + 0.0297r(p2)
for 50 e r(p2) e 61.36 (40)
H ) H′0 + 0.0463r(p2)
for 61.36 e r(p2) e 67.5 (41)
where H0 and H′0 correspond to the makespans at r(p2) ) 50 and 61.36, respectively. The LP problem becomes infeasible for r(p2) greater than 67.5. If there exists a point r*(p2) ∈[50, 61.36] at which a new schedule becomes optimal, the objective value at this point will be equal to H0 + 0.0297r*(p2). Thus, the
3790 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 Table 4. Values of Binary Variables of Optimal Schedules for Example 2 task, unit
n1
n2
n3
heating, heater reaction 1, reactor 1 reaction 1, reactor 2 reaction 2, reactor 1 reaction 2, reactor 2 reaction 3, reactor 1 reaction 3, reactor 2 separation, still
Schedule 1 1 0 1 0 0 0 0 0
n0
0 0 1 1 0 0 0 0
0 0 0 0 0 1 1 0
0 0 0 1 0 0 0 1
heating, heater reaction 1, reactor 1 reaction 1, reactor 2 reaction 2, reactor 1 reaction 2, reactor 2 reaction 3, reactor 1 reaction 3, reactor 2 separation, still
Schedule 2 1 1 1 0 0 0 0 0
1 0 0 1 1 0 0 0
0 0 0 0 0 1 1 0
0 0 0 1 0 0 0 1
heating, heater reaction 1, reactor 1 reaction 1, reactor 2 reaction 2, reactor 1 reaction 2, reactor 2 reaction 3, reactor 1 reaction 3, reactor 2 separation, still
Schedule 3 1 1 1 0 0 0 0 0
0 0 0 1 1 0 0 0
0 0 0 0 0 1 1 0
0 0 0 0 0 0 0 1
following problem is solved
minimize H subject to
constraints 2-11
H ) H0 + 0.0297r(p2)
(42)
50 e r(p2) e 61.36
∑
(i,j,n)∈F1
wv(i,j,n) -
∑
wv(i,j,n) e |F1|
(i,j,n)∈F0
The last constraint in problem 42 is added to exclude the current schedule from the solution. A new optimal schedule (schedule 2) is obtained at the point r(p2) ) 56.85, which indicates that the range [50, 56.85] is the optimality range for schedule 1. If the problem is infeasible, a new optimal solution will be sought within the range of r(p2) ∈[61.36, 67.5]. By repeating the above steps, schedule 3 is found to be the optimal solution at the point r(p2) ) 83.87. Using parametric programming, all of the optimal schedules within the uncertainty range under consideration can be determined with the corresponding break points. However, a number of disadvantages are apparent. First, in the case in which multiple equivalent schedules exist, the integer cut cannot differentiate among the schedules. Thus, one more integer cut has to be added each time to exclude the equivalent schedule, which will evidently increase the complexity of the problem. Second, the underlying complexity of the multiparametric MILP problems far exceeds the effort of the proposed approach, which becomes infeasible for large-scale problems where multiple parameters need to be investigated simultaneously. The examples appearing in the literature do not involve more than a few binary variables, although most practical scheduling problems will involve on the order of hundreds assignment decisions. GAMS/CPLEX 7.5 is used for the solution of the deterministic formulation of the nominal
value and relaxed LP problems at each leaf node, as well as the solution of the robust optimization and parametric programming model. 5. Future Work and Directions An integrated framework is developed in this work to handle uncertainty in short-term scheduling based on the idea of inference-based sensitivity analysis for MILP problems and the utilization of a branch-andbound solution methodology. The proposed method leads to the determination of the importance of different parameters and the constraints on the objective function and the generation and evaluation of a set of alternative schedules given the variability of the uncertain parameters. The main advantage of the proposed method is that no substantial complexity is added compared with the solution of the deterministic case because the only additional required information is the dual information at the leaf nodes of the branch-and-bound tree. Two illustrative examples are presented in this paper to highlight the information extracted by the proposed approach and the complexity involved compared with robust and parametric studies. Work is in progress to further develop the proposed approach by improving the retrieval and storage of the required dual information at the leaf nodes in the branch-and-bound tree. This will enable the utilization of the approach in large-scale problems and the consideration of multiple parameters. Additional work is directed toward the application of the proposed approach to the area of refinery scheduling. Acknowledgment The authors gratefully acknowledge financial support from the National Science Foundation under NSF CAREER Grant CTS-9983406. Note Added after ASAP Posting. This article was released ASAP on February 26, 2004. The correct version, with a modification to the second sentence of paragraph 3.1, was posted on March 4, 2004. Nomenclature Indices i ) tasks j ) units n ) event points s ) states k ) scenarios Parameters pk ) probability of scenario k price(s) ) price of state s Fp(s,i), Fc(s,i) ) proportion of state s produced, consumed, respectively, from task i r(s) ) market requirement for state s at the end of the time horizon Vmin(i,j) ) minimum capacity of unit j when processing task i Vmax(i,j) ) maximum capacity of unit j when processing task i stmax(s) ) maximum storage capacity of state s Variables wv(i,j,n) ) binary variables that assign the beginning of task i in unit j at event point n
Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3791 b(i,j,n) ) amount of material undertaking task i in unit j at event point n d(s,n) ) amount of state s being delivered to the market at event point n H ) time horizon st(s,n) ) amount of state s at event point n Ts(i,j,n) ) starting time of task i in unit j at event point n Tf(i,j,n) ) finishing time of task i in unit j at event point n
Literature Cited (1) Kondili, E.; Pantelides, C. C.; Sargent, R. W. H. A general algorithm for scheduling batch operations. Comput. Chem. Eng. 1993, 17, 211. (2) Schilling, G.; Pantelides, C. C. A simple continuous-time process scheduling formulation and a novel solution algorithm. Comput. Chem. Eng. 1996, 20, S1221. (3) Mockus, L.; Reklaitis, G. V. Continuous-time representation approach to batch and continuous process scheduling. 1. MINLP formulation. Ind. Eng. Chem. Res. 1999, 38, 1973. (4) Castro, P.; Barbosa-Povoa, A. P. F. D.; Matos, H. An improved RTN continuous-time formulation for the short-term scheduling of multipurpose batch plants. Ind. Eng. Chem. Res. 2001, 40, 2059. (5) 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. (6) 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. (7) Maravelias, C. T.; Grossmann, I. E. New general continuoustime state-task network formulation for short-term scheduling of multipurpose batch plants. Ind. Eng. Chem. Res. 2003, 42, 3056. (8) Shah, N.; Pantelides, C. C. Design of multipurpose batch plants with uncertain production requirements. Ind. Eng. Chem. Res. 1992, 31, 1325. (9) Dembo, R. Scenario Optimization. Ann. Oper. Res. 1991, 30, 63. (10) Pistikopoulos, E. N.; Ierapetritou, M. G. A novel approach for optimal process design under uncertainty. Comput. Chem. Eng. 1995, 19, 1089. (11) Straub, D. A.; Grossmann, I. E. Design optimization of stochastic flexibility. Comput. Chem. Eng. 1993, 17, 339. (12) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N. Flexibility analysis and design using a parametric programming framework. AIChE J. 2002, 48, 2851. (13) Petkov, S. B.; Maranas, C. D. Multiperiod planning and scheduling of multiproduct batch plants under demand uncertainty. Ind. Eng. Chem. Res. 1997, 36, 4864. (14) Schmidt, C.; Grossmann, I. E. Optimization models for the scheduling of testing tasks in new product development. Ind. Eng. Chem. Res. 1996, 35, 3498.
(15) Fortemps, P. Jobshop scheduling with imprecise durations: a fuzzy approach. IEEE Trans. Fuzzy Syst. 1997, 5, 557. (16) Ishibuchi, H.; Yamamoto, N.; Murata, T.; Tanaka, H. Genetic algorithms and neighborhood search algorithms for fuzzy flowshop scheduling problems. Fuzzy Sets Syst. 1994, 67, 81. (17) Balasubramanian, J.; Grossmann, I. E. Scheduling optimization under uncertaintysAn alternative approach. Comput. Chem. Eng. 2003, 27, 469. (18) Daniels, R. L.; Carrillo, J. E. β-robust scheduling for singlemachine systems with uncertain processing times. IIE Trans. 1997, 29, 977. (19) Vin, J. P.; Ierapetritou, M. G. Robust short-term scheduling of multiproduct batch plants under demand uncertainty. Ind. Eng. Chem. Res. 2001, 40, 4543. (20) Acevedo, J.; Pistikopoulos, E. N. A multiparametric programming approach for linear process engineering problems under uncertainty. Ind. Eng. Chem. Res. 1997, 36, 717. (21) Dua, V.; Pistikopoulos, E. N. Algorithms for the solution of multiparametric mixed integer nonlinear optimization problems. Ind. Eng. Chem. Res. 1999, 38, 3976. (22) Dua, V.; Pistikopoulos, E. N. Algorithms for the solution of multiparametric mixed integer linear optimization problems. Ann. Oper. Res. 2000, 99, 123. (23) Pertsinidis, A.; Grossmann, I. E.; McRae, G. J. Parametric optimization of MILP programs and a framework for the parametric optimization of MINLPs. Comput. Chem. Eng. 1998, 22, S205. (24) Schrage, L.; Wolsey, L. Sensitivity analysis for branch and bound integer programming. Oper. Res. 1985, 33, 1008. (25) Skorin-Kapov, J.; Granot, F. Nonlinear integer programming: Sensitivity analysis for branch and bound. Oper. Res. Lett. 1987, 6, 269. (26) Bertsimas, D.; Georgia, P.; Tayur, S. A new algebraic geometry algorithm for integer programming. Manage. Sci. 2000, 46, 999. (27) Dawande, M. W.; Hooker, J. N. Inference-based sensitivity analysis for mixed integer/linear programming. Oper. Res. 2000, 48, 623. (28) Floudas, C. A. Nonlinear and Mixed-Integer Optimization; Oxford University Press: New York, 1995. (29) Samsatli, N. J.; Papageorgiou, L. G.; Shah, N. Robustness metrics for dynamic optimization models under parameter uncertainty. AIChE J. 1998, 44, 1993. (30) Mulvey, J. M.; Vanderbei, R. J.; Zenios, S. A. Robust optimization of large-scale systems. Oper. Res. 1995, 43, 264. (31) Ahmed, S.; Sahinidis, N. Robust process planning under uncertainty. Ind. Eng. Chem. Res. 1998, 37, 1883.
Received for review August 15, 2003 Revised manuscript received December 1, 2003 Accepted December 3, 2003 IE0306731