Article pubs.acs.org/IECR
Capacity Planning under Clinical Trials Uncertainty in Continuous Pharmaceutical Manufacturing, 2: Solution Method Arul Sundaramoorthy,† Xiang Li,†,§ James M. B. Evans,‡ and Paul I. Barton*,† †
Process Systems Engineering Laboratory, Department of Chemical Engineering, and ‡Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, United States S Supporting Information *
ABSTRACT: In Part 1 of this paper, we presented a scenario-based multiperiod mixed-integer linear programming (MILP) formulation for a capacity planning problem in continuous pharmaceutical manufacturing under clinical trials uncertainty. The number of scenarios and, thus, the formulation size grows exponentially with the number of products. The model size easily becomes intractable for conventional algorithms for more than 8 products. However, industrial-scale problems often involve 10 or more products, and thus a scalable solution algorithm is essential to solve such large-scale problems in reasonable times. In this part of the paper, we develop a rigorous decomposition strategy that exploits the underlying problem structure. We demonstrate the effectiveness of the proposed algorithm using several examples containing up to 16 potential products and over 65 000 scenarios. With the proposed decomposition algorithm, the solution time scales linearly with the number of scenarios, whereby a 16-product example with over 65 million binary variables, nearly 240 million continuous variables, and over 250 million constraints was solved in less than 6 h of solver time.
1. INTRODUCTION The capacity planning problem in the pharmaceutical industry is characterized by large technical uncertainty in the clinical trials. Since time to market is the most important driver in the pharmaceutical industry, the launch-related decisions for a new product, such as when and how much capacity to build, when and how much to produce, store, etc., must be taken into consideration much before its anticipated launch date. This means that the product-launch decisions are made when the new product is still in clinical development. The product can either pass or fail the trial, but it is mandatory that enough capacity is ensured for the smooth launch of the new product into the market. In Part 1 of this paper, we addressed the above problem and developed a mathematical framework based on two-stage stochastic programming (TSSP). In the process systems engineering literature, there exist some TSSP-based formulations for capacity planning in the pharmaceutical industry under clinical trials uncertainty. Maravelias and Grossmann1 proposed an integrated framework for simultaneous drug development and design of batch facilities. They adopted a TSSP approach to account for the uncertainty in the outcome of clinical trials. To solve large-scale problem instances, they proposed an iterative heuristic solution method. Levis and Papageorgiou2 developed a large-scale MILP model based on TSSP for multisite capacity planning in the pharmaceutical industry under clinical trials uncertainty. To reduce the computational burden of large-scale models, they proposed a hierarchical solution algorithm. In the proposed formulation, we model the uncertainty in the final outcome of the development process, and, thus the number of scenarios grows exponentially with the number of products as 2I, where I is the total number of products (see Part 1 of this paper). In real-life capacity planning problems, the number of products can easily be 10 or more. For example, a problem © 2012 American Chemical Society
with 16 products contains over 65 000 scenarios. Thus, the problem size gets more challenging and becomes quickly intractable for conventional algorithms. To solve such massive problems in realistic time, an efficient solution method is necessary. The special structures of TSSP with recourse are wellexploited by duality-based decomposition methods.3 These methods are computationally attractive because the sizes of subproblems solved are independent of the number of scenarios, and decomposition enables solution of the subproblems in parallel. Two of the most widely studied duality-based decomposition methods are (i) Benders decomposition4 (BD), also known as the L-shaped method,5 and (ii) Lagrangian decomposition6 (LD). In BD or its extension to nonlinear problems, generalized Benders decomposition7 (GBD), the original problem is projected onto the space of first-stage variables and then reformulated into a dual problem. Because the dual problem contains an infinite number of constraints, it is relaxed into a lower bounding problem with a finite subset of these constraints. While the sequence of lower bounds on the original problem solution value is generated by a sequence of solutions of lower bounding problems, the sequence of upper bounds on the original problem solution value is obtained by a sequence of solutions of upper bounding problems, which are generated by fixing the first-stage decisions to the solutions of the lower bounding problems. When the sequences of upper and lower bounds converge, we obtain an optimal solution to the original problem. However, such convergence depends heavily on strong duality in the dual reformulation, which is often nonexistent for nonconvex problems. It should be noted that when Received: Revised: Accepted: Published: 13703
February 6, 2012 August 2, 2012 September 19, 2012 October 9, 2012 dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
Industrial & Engineering Chemistry Research
Article
the first-stage decisions are fixed, the upper bounding problem can be decomposed into much smaller subproblems for each of the scenarios, which can then be solved in parallel, thus making this method computationally very attractive. In LD, each scenario contains a duplicate of the first-stage variables, which are then linked across the scenarios through additional equality constraints. By dualizing these equality constraints into the objective function, a convex but nonsmooth dual problem is generated, which can be solved using a subgradient method8 that facilitates decomposition into smaller subproblems corresponding to each of the scenarios. In the LD method, convergence to a global optimum in the absence of strong duality properties is usually guaranteed by using a branch-and-bound framework,9,10 in which convergence of the upper and lower bounds can only be guaranteed by branching in the full space. When the number of scenarios is large, fullspace branching can easily become intractable posing challenges for these LD-based methods. Li et al.11−13 extended the GBD framework to a more general framework, called nonconvex generalized Benders decomposition (NGBD), to handle nonconvex recourse problems in TSSP rigorously. In NGBD, nonconvexity in the form of nonconvex functions participating in the objective and/or constraints is effectively handled using McCormick’s relaxation technique.14,15 Here, an additional convex lower bounding problem is added as a surrogate for the original problem, and the traditional GBD iteration is applied to the surrogate problem instead of the original problem. In this work, we extend the NGBD framework to address TSSP with integer recourse, where the nonconvexity in recourse problems arises from binary variables. In the following section, we discuss decomposition of the original problem into subproblems that are used in the NGBD algorithm and discuss their properties. In Section 3, we present a pseudo-code for the NGBD algorithm and discuss the finite convergence of the algorithm. In Section 4, we present computational experiments on several large-scale capacity planning problems in continuous pharmaceutical manufacturing, followed by conclusions in Section 5.
negative of the objective function in the original model for the ease of discussion of the algorithmic development. In problem P, binary variables y represent the first-stage decisions whereas the mixed (continuous and binary) variables xs represent the second-stage decisions. While set Y is defined by the bounds on the y variables and the constraints containing only y variables, set X0s is defined by the bounds on xs variables and the constraints containing only xs variables. Note that the first J components of xs are binary variables. The constraints containing both y and xs variables are represented by As xs + By ≤ 0
∀s
Furthermore, cTy and bTs xs represent the costs associated with the first-stage and second-stage decisions, respectively. In NGBD, the concepts of projection, dualization, restriction, and relaxation are employed. We extend the framework of NGBD to solve TSSP with binary variables in the first stage, and mixed variables in the second stage. NGBD requires an additional convex lower bounding problem, which is used as a surrogate for the original problem. The sequence of lower bounds on the original problem is obtained by applying the GBD iteration to the surrogate lower bounding problem. On the other hand, the sequence of upper bounds on the original problem is obtained by solving the primal problem, which is obtained by fixing the first-stage decisions in the original problem. For each possible integer realization y ∈ Y, if there exists a tractable convex hull representation for the feasible set of P (with fixed y), then the GBD iteration can be readily applied to the original problem without incorporating the surrogate lower bounding problem. However, it is often not a trivial task to obtain the convex hull, because it often requires the addition of innumerable constraints and/or variables. Alternatively, we employ continuous relaxation of the second-stage binary variables by treating them as [0,1] continuous variables. This relaxes the original problem into the following lower bounding problem: S
min c Ty +
2. DECOMPOSITION STRATEGY The multiscenario multiperiod MILP model presented in Part 1 of this paper can be compactly represented as follows:
x1 , ..., xS , y
x1 , ..., xS , y
∑ κsbsTxs s=1
As xs + By ≤ 0 xs ∈ Xs
(P)
(LBP)
s = 1 , ..., S
s = 1 , ..., S
y∈Y
subject to As xs + By ≤ 0 xs ∈ Xs0
s=1
subject to
S
min c Ty +
∑ κsbsTxs
where set Xs contains relaxed [0,1] continuous variables instead of binary variables. When LBP is solved for fixed y, the problem becomes a linear program (LP) defined on a compact set and, thus, has a finite optimal objective value or is infeasible. Proposition 2.1. Any feasible point of problem P is also feasible for LBP, and the optimal objective value of LBP represents a lower bound on the optimal objective value of problem P. Proof. By definition, X0s ⊂ Xs, ∀s = 1, ..., S, and, thus, the lower bounding principle follows. Now, the sequence of lower bounds on the surrogate problem LBP is obtained from a master problem, which is generated by projection of LBP from the space of both continuous and binary variables to the space of first-stage binary variables, followed by
s = 1, ..., S
s = 1, ..., S
y∈Y
where Y = {y ∈ {0,1}ny:Dy ≤ d} and X0s = {xs ∈ Xs:xsj ∈ {0,1}, j = 1, ...,J}. The set Xs is given by Xs = {xs ∈ nx :Esxs ≤ es,xs ≥ 0, xsj ≤ 1, j = 1, ..., J}. It is assumed that set Y is nonempty, and set X0s is nonempty and compact for s = 1, ..., S. Note that we use minimization of the objective function in problem P, as opposed to maximization in the original model of Part 1 of this paper. However, maximizing a function is equivalent to minimizing the negative of the same function. We minimize the 13704
dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
Industrial & Engineering Chemistry Research
Article
reformulation of any subproblem for a fixed-integer realization y into its dual. The master problem is given as shown below: min η
k
zsPBP = min κs(c Ty k + bsTxs) xs
subject to
(MP)
η,y
As xs + By k ≤ 0
subject to S
S
∑
η≥
s=1
(κsbsTxs
inf
+
xs ∈ Xs
λsTAs xs)
T
+ (c +
xs ∈ Xs
∑ λsTB)y s=1
k
where zPBP represents the optimal objective value of PBPsk. s
∀ λ1 , ..., λS ≥ 0
k
S s=1
s
s
∀ (μ1 , ..., μS ) ∈ M
s=1
y ∈ Y, η ∈
where λs, μs ∈ m, s = 1, ...,S, and M = {(μ1, ..., μS):μ1, ...,μS ≥ 0,∑Ss=1∑mi=1μsi > 0}. It should be noted that master problem MP is obtained by replacing S
M1 = {(μ1 , ..., μS ): μ1 , ..., μS ≥ 0,
k
= ∑Ss=1zPBP . Note that zPBP s s Proposition 2.3. For the kth integer realization, if λ*s are Lagrange multipliers for problem PBPsk, ∀s = 1, ..., S, then (λ1*, ..., λS*) are Lagrange multipliers for problem PBPk. Proof. The proof is trivial and is omitted. If, for scenario s, problem PBPsk is infeasible, then problem PBPk is infeasible, in which case the following feasibility problem is solved:
S
(μsT As xs) + ∑ μsT By ∑ xinf ∈X
0≥
(PBPks )
S
k
z FP =
m
∑ ∑ μsi = 1}
min
x1,..., xS , v1,..., vS
∑ κs || vs || s=1
(FPk)
s=1 i=1
subject to
with S
M = {(μ1 , ..., μS ): μ1 , ..., μS ≥ 0,
As xs + By k ≤ vs
m
∑ ∑ μsi > 0} s=1 i=1
Such replacement facilitates the generation of valid subproblems. However, both the master problems with sets M and M1, respectively, are equivalent.11 Proposition 2.2. Problems LBP and MP are equivalent in the sense that: (1) Problem LBP is feasible if and only if problem LBP is feasible; (2) The optimal objective values of problems LBP and MP are the same; (3) The optimal objective value of problem LBP is attained with an integer realization if and only if the optimal objective value of problem MP is attained with the same integer realization. Proof. For any fixed integer realization y, problem LBP becomes an LP, for which strong duality holds.16 Thus, the result follows immediately from Theorems 2.1, 2.2, and 2.3 in the work of Geoffrion.7 The sequence of upper bounds on problem LBP is obtained from the primal bounding problems, which are generated by restricting the y variables in LBP to an element yk in the set Y, where superscript k denotes the kth integer realization visited by the primal bounding problem:
x1 ,..., xS
∑ κsbsTxs s=1
xs ∈ Xs
s = 1 , ..., S
vs ∈ Vs
s = 1 , ..., S
k
k
zsFP = min κs || vs ||
(FPks )
xs , vs
subject to As xs + By k ≤ vs xs ∈ Xs vs ∈ Vs k
where zFP represents the optimal objective value of problem s k
k
FPsk, and zFP = ∑Ss=1zFP s . Proposition 2.4. For the kth integer realization, if μ*s are Lagrange multipliers for problem FPsk, ∀s = 1, ..., S, then (μ*1 , ..., μS*) are Lagrange multipliers for FPk. Proof. The proof is trivial and is omitted. Since the master problem contains infinitely many constraints, it is difficult to solve the problem directly. Thus, we solve the following relaxed master problem by enforcing only a finite number of constraints:
(PBPk)
subject to As xs + By k ≤ 0
xs ∈ Xs
where zFP represents the optimal objective value of problem FPk, ∥νs∥ represents an arbitrary norm of the slack variable vector vs for s = 1, ..., S, and set Vs ⊂ {vs ∈ m : vs ≥ 0}. The violations of the constraints are measured by elements of vs, the norm of which is then minimized for minimum violation of the constraints. Now, FPk can be decomposed into subproblems for each of the S scenarios:
S
k
z PBP = min c Ty k +
s = 1 , ..., S
s = 1, ..., S
s = 1, ..., S
min η
k
where zPBP represents the optimal objective value of problem PBP k. Now, problem PBP k can be decomposed into subproblems for each of the S scenarios:
η,y
(RMP1k)
subject to 13705
dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
Industrial & Engineering Chemistry Research S
η≥
∑ s=1
Article
S
inf (κsbsTxs + (λsj)T As xs) + (c T +
xs ∈ Xs
η ≥ z PBP + (c T +
s=1
0 ≥ z FP + (∑ (μsi )T B)(y − y i )
S
∑ xinf ∈X s
s=1
((μsi )T As xs)
+
s
∑
∑
(μsi )T By
∀i∈W
k
r ∈ {r : yrt = 1, r = 1,..., ny} ≤ |{r: yrt = 1}|
∑ r ∈ {r : yrt = 0, r = 1,..., ny} k
−1
∑
k
∑
y −
r r ∈ {r : yrt = 1, r = 1,..., ny} ≤ |{r: yrt = 1}|
yr
∀t∈T ∪W
∀ i ∈ Wk
s=1
s=1
yr −
∀ j ∈ Tk
S
i
S
∑ (λsj)T B)(y − y j ) s=1
∀ j ∈ Tk 0≥
S
j
∑ (λsj)T B)y
r ∈ {r : yrt = 0, r = 1,..., ny} k
−1
yr
∀ t ∈ T ∪ Wk
y ∈ Y, η ∈
y ∈ Y, η ∈
Similarly, problem FRMP1k is equivalent to the following single-level MILP:
where the index sets
ny
T k = {j ∈ {1 , ..., k}: problem PBP is feasible for y = y j }
min ∑ yl y
and
(FRMPk)
l=1
subject to k
i
W = {i ∈ {1 , ..., k}: problem PBP is infeasible for y = y }
S
i
0 ≥ z FP + (∑ (μsi )T B)(y − y i )
While Lagrange multipliers λjs for problem PBPjs form an optimality cut for iteration j ∈ Tk, Lagrange multipliers μis form a feasibility cut for iteration i ∈ Wk. To avoid re-examination of the same integer realization, a set of canonical integer cuts is added to problem RMP1k, as described by Balas and Jeroslow.17 Proposition 2.5. Problem RMP1k is a relaxation of MP when the latter is augmented with the relevant canonical integer cuts excluding the previously examined integer realizations. Proof. This follows from Proposition 3.8 in Li et al.11 If Tk is empty, then problem RMP1k is unbounded, in which case the following feasibility relaxed master problem is solved:
s=1
∑
y
S
S
inf ((μsi )T As xs) +
∑ (μsi )T By
s
s=1
∑ x ∈X s=1
∑
s
yr −
r ∈ {r : yrt = 1, r = 1 , ..., ny} ≤ |{r: yrt = 1}|
∑ r ∈ {r : yrt = 0, r = 1 , ..., ny} k
−1
∀ i ∈ Wk
z
∀t∈W
S T l
= min c y + x1,..., xS
∑ κsbsTxs s=1
(PPl)
subject to
∀t∈W
As xs + By l ≤ 0 k
xs ∈ Xs0
k
For the kth integer realization, after solving PBP or FP subproblems, RMP1k or FRMP1k must be solved to generate a new integer realization. However, RMP1k and FRMP1k are bilevel optimization problems, which can be reformulated into single-level optimization problems by exploiting the separability of the functions in the continuous and integer variables. Thus, problem RMP1k is equivalent to the following single-level MILP: η,y
PPl
yr
y∈Y
min η
−1
yr
Note that the sizes of problems RMPk and FRMPk are independent of the number of scenarios in the original problem. Proposition 2.6. Problems RMP1 k and RMP k are equivalent. Proof. This follows from the results in Section 2.5 in the work of Geoffrion.7 Finally, the sequence of upper bounds on the original problem is obtained by solving the primal problems, which are generated by fixing the y variables in problem P to an element yl in the set Y, where the superscript l denotes the lth integer realization visited by the primal problem:
subject to 0≥
r ∈ {r : yrt = 0, r = 1 , ..., ny} k
y∈Y
(FRMP1k)
l=1
∑
y −
r r ∈ {r : yrt = 1, r = 1,..., ny} ≤ |{r: yrt = 1}|
ny
min ∑ yl
∀ i ∈ Wk
s = 1 , ... , S
s = 1 , ... , S
l
where zPP represents the optimal objective value of problem PPk, which can be decomposed into subproblems for each of the S scenarios: l
zsPP = min κs(c Ty l + bsTxs) xs
(PPls)
subject to As xs + By l ≤ 0
(RMPk)
xs ∈ Xs0
subject to 13706
dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
Industrial & Engineering Chemistry Research
Article
Figure 1. Generation of subproblems in NGBD. Subproblems in the dashed boxes are not solved whereas those in the gray solid boxes are solved. l
tolerances. Note that standard algorithms for MILP are finite, so these tolerances are not strictly necessary. However, as for algorithms for the full model (P), such tolerances can accelerate the solution of hard instances. The pseudo-code for the NGBD algorithm is as follows: Initialize. 1. Counters k = 0, l = 1, and the index sets T0 = ⌀, W0 = ⌀, U0 = ⌀. 2. Upper bound on P, UBP = +∞; upper bound on LBP, UBLBP = +∞; lower bound on LBP, LBLBP = −∞. 3. Set tolerances ε′ and ε. Note that ∑Ss=1εs ≤ ε, where εs is the tolerance for the solution of PPs. 4. Specify initial integer realization, yk=1. repeat if k = 0 or (RMPk is feasible and LBLBP < UBLBP − ε′ and LBLBP < UBP − ε), then repeat
where zPP represents the optimal objective value of PPls, and s l l zPP = ∑Ss=1zPP s . Figure 1 illustrates the generation of subproblems through continuous relaxation, projection/dualization, restriction and relaxation. In Figure 1, the problems in the dashed boxes are not solved directly in NGBD. On the other hand, the subproblems in the gray boxes with solid lines, including RMP, FRMP, PBP subproblems, FP subproblems, and PP subproblems are solved in NGBD. Note that FRMP and FP subproblems are solved only when the RMP is unbounded and any of the PBP subproblems are infeasible, respectively. Furthermore, the sizes of these subproblems are independent of the number of scenarios in the original problem.
3. ALGORITHM Similar to the outer approximation algorithm of Kesavan et al.,18 the NGBD algorithm contains two loops: inner and outer. While the inner loop generates the sequence of lower bounds on the original problem by solving the surrogate lower bounding problems, the outer loop generates the sequence of upper bounds on the original problem by solving the primal problems. Assumption. The primal bounding problems are solved exactly in a finite number of steps because they are LP problems. In the inner loop, the traditional GBD iteration is applied to LBP. Thus, the sequence of upper bounds on LBP (UBLBP) is obtained by solving the primal bounding problems, which are LPs. On the other hand, the sequence of lower bounds on LBP (LBLBP) is obtained by solving the relaxed master problems, which are MILPs. We use a tolerance only for the relaxed master problems since the primal bounding problems are assumed to be solved exactly. Let ε′ ≥ 0 be the tolerance for the solution of RMP. In the outer loop, the sequence of upper bounds on P (UBP) is obtained by solving the primal problems, which are MILPs. Thus, we use a tolerance for the primal problems as well. Let ε ≥ 0 be the tolerance for the solution of PP. For the ease of presentation of the proposed algorithm, we use absolute tolerances for RMP and PP. However, one could also use relative tolerances or both absolute and relative
set k = k+1 1. Solve PBPks , for s = 1, ..., S sequentially. If PBPks is feasible for all the scenarios with Lagrange multipliers λks , then add an optimality cut to RMPk, according to λk1, ..., λkS, and set k k Tk = Tk−1∪k. If zPBP = ∑Ss=1zPBP < s k
UBLBP, then update UBLBP = zPBP , y* = yk, k* = k. 2. If problem PBPsk is infeasible for scenario s,̅ stop solving it for s = s ̅ + 1, ..., S, and set Wk = Wk−1∪k. Then set μks = 0, for s = 1, ..., s ̅ − 1, and solve FPsk, for s = s,̅ ..., S, and obtain the corresponding Lagrange multipliers μks . Add a feasibility cut to RMPk, according to μk1, ..., μkS. 3. If Tk = ⌀, solve FRMPk to ε′-optimality; otherwise, solve problem RMPk to ε′-optimality. If problem RMPk is 13707
dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
Industrial & Engineering Chemistry Research
Article
feasible, set LBLBP to the ε′-optimal solution value of problem RMPk. In either case, set yk+1 to the y value at the solution of the corresponding problem. until LBLBP ≥ UBLBP − ε′ or (RMPk or FRMPk is infeasible).
to consider building two greenfield continuous facilities (F1 and F2) and their possible expansions in the future. We are given the probability of the outcome of phase II and phase III clinical trials for each product. The decisions of interest to the company are when to build the greenfield facilities and when to expand them, and in each scenario, how much capacity, which product to produce in which facility, and in what amounts to produce, store, etc. In Part 1 of this paper, we studied five different examples to illustrate the features and performance of the proposed model. These examples varied in terms of the number of products in the portfolio and, thus, the number of scenarios addressed. Our preliminary study showed that when we increased the number of products to more than eight, the size of the full-space model becomes prohibitively large, so that even the state-of-the-art commercial MILP solver cannot solve it within a reasonable time limit and/or the solution requires very large memory resources. In this computational study, we show how such computational resource limitations can be effectively managed using the proposed decomposition algorithm for large-scale problems. To illustrate the computational performance of the NGBD algorithm, we consider eight different examples containing up to 16 products and over 65 000 scenarios. In all the examples, we consider a constant capacity size and fixed investment cost. Tables S1−S4 in the Supporting Information provide the data for all the examples. We used GAMS 23.7/CPLEX 12.3 on a 3.2 GHz Intel Xeon CPU with 12GB RAM on a Windows platform. GAMS was used to program the NGBD algorithm, and CPLEX was used to solve all LP and MILP subproblems. In the NGBD algorithm, both absolute and relative tolerances for PP were set to 10−3, while those for RMP and FRMP were set to 10−5. Furthermore, the integer realization (yk=1) was initialized to zero in all the examples. “Full model” refers to solution of problem P directly with CPLEX. The absolute and relative optimality criteria for the full model were set to 10−3. Table 1 provides a summary of model and solution statistics for all the examples. Note that the solution times reported in Table 1 are only the solver times of CPLEX in GAMS. The number of first-stage variables remains the same, whereas that of second-stage variables increases exponentially with the number of products. Note that Table 1 shows the number of second-stage variables, both binary and continuous, per scenario. To compute the actual number of second-stage variables, multiply the number of variables per scenario with the number of scenarios. For the biggest example (Ex8) with 16 products, the total number of second-stage variables is ∼305 million, which constitutes over 65 million binary variables and nearly 240 million continuous variables. However, note that, in the decomposition approach, it is never necessary to assemble all these variables and constraints in memory at one time. The full model was solved within tolerance and time limit for examples containing up to eight products (800%. Notice that the time taken for example Ex5 to get even a suboptimal solution with 800% gap is roughly 30 times more than the time taken to solve the previous example (Ex4) within tolerance. On the other hand, the NGBD algorithm finds an ε-optimal solution for example Ex5 in less than 5 min of CPU time. Thus, the proposed NGBD algorithm can serve as a supporting tool for decision makers by providing them a faster, better, and reliable solution.
■
Figure 5. Convergence of LBD and UBD in the NGBD algorithm for example Ex3.
ASSOCIATED CONTENT
S Supporting Information *
The data for all the examples are given as Supporting Information. Table S1 contains the probability of outcome of clinical trials for each product, while Table S2 provides the fixed investment costs, constant capacity sizes, minimum/maximum production rates, and maximum number of capacity expansion. Table S3 lists the price and cost details of products, and Table S4 contains the demand details of products from the launch to the peak periods. The descriptions for symbols used in Tables S1−S4 are given in Part 1 of this paper. This information is available free of charge via the Internet at http://pubs.acs.org/.
solving the LP relaxation of the full MILP model P, is very close to the optimal objective value of P in all of the examples. This potentially contributes to the faster convergence of both the branch-and-bound (full) and NGBD algorithms. For the NGBD algorithm, only a few integer realizations are required for convergence. For example, the number of integer realizations for the primal bounding problem is only two for examples Ex3−Ex8, whereas the number of integer realizations for the primal problem is just one (see Table 1). Interestingly, despite the increase in the problem size in examples Ex3−Ex8, the number of integer realizations required remains the same. This may not be the case for less-tight formulations as more iterations and integer realizations would be required for the branch-and-bound and NGBD algorithms, respectively.
■
AUTHOR INFORMATION
Corresponding Author
*Tel.: +1 617 253 6526. Fax: +1 617 253 5042. E-mail: pib@ mit.edu. Present Address
5. CONCLUSION In this two-part paper, we addressed an important strategic problem faced by major pharmaceutical companies. With the advent of continuous manufacturing in the pharmaceutical industry, there exist several opportunities to improve and challenges to address throughout the supply chain. In Part 1 of this paper, we developed a mathematical framework for capacity planning problem under clinical trials uncertainty based on TSSP. Since the proposed formulation resulted in a large-scale MILP model, whose size increases exponentially with the number of products, the existing state-of-the-art commercial solvers are not sufficient to solve such massive problems. This motivated us to develop an efficient decomposition algorithm in this part of the paper. The proposed NGBD algorithm efficiently solved problems with up to 16 products and over 65 000 scenarios in less than 6 h of solver time. The largest example (Ex8) contained over 250 million constraints, and about 305 million variables, of which over 65 million variables are binary and nearly 240 million variables are continuous. The solution time of the NGBD algorithm scales linearly with the number of scenarios in the problem, which makes it attractive to solve industrial-scale problems. Furthermore, one can achieve an increased speedup with the NGBD algorithm by employing massively parallel computing since the subproblems can be solved in parallel. It may be argued that, for strategic problems, such as the one considered in this work, the solver time can be much longer than a few hours. However, often, the decision maker would like to revise the formulation or rerun the model for some new parameters and analyze the results. If the solution time takes longer than a reasonable time (for example, days when solved
§
Department of Chemical Engineering, Queen’s University, 19 Division Street, Kingston, Ontario, K7L 3N6, Canada. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was funded by Novartis Pharmaceuticals as part of the Novartis-MIT Center for Continuous Manufacturing. The authors would like to thank the Executive Director of the Novartis-MIT Center for Continuous Manufacturing, Dr. Stephen R. Sofen, for insight into the problem.
■
REFERENCES
(1) Maravelias, C. T.; Grossmann, I. E. Simultaneous planning for new product development and batch manufacturing facilities. Ind. Eng. Chem. Res. 2001, 40, 6147−6164. (2) Levis, A. A.; Papageorgiou, L. G. A hierarchical solution approach for multi-site capacity planning under uncertainty in the pharmaceutical industry. Comput. Chem. Eng. 2004, 28, 707−725. (3) Birge, J. R.; Louveaux, F. Introduction to Stochastic Programming; Springer: New York, 1997. (4) Benders, J. F. Partitioning procedures for solving mixed-variables programming problems. Numer. Math. 1962, 4, 238−252. (5) Van Slyke, R. M.; Wets, R. L-shaped linear programs with applications to optimal control and stochastic programming. SIAM J. Appl. Math. 1969, 17, 638−663. (6) Fisher, M. L. The Lagrangian relaxation method for solving integer programming problems. Manage. Sci. 1981, 27, 1−18. (7) Geoffrion, A. M. Generalized Benders decomposition. J. Opt. Theory Appl. 1972, 10, 237−260. (8) Fisher, M. L. An applications oriented guide to Lagrangian relaxation. Interfaces 1985, 15, 10−21.
13710
dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
Industrial & Engineering Chemistry Research
Article
(9) Karuppiah, R.; Grossmann, I. E. A Lagrangean based branch-andcut algorithm for global optimization of nonconvex mixed-integer nonlinear programs with decomposable structures. J. Global Opt. 2008, 41, 163−186. (10) Khajavirad, A.; Michalek, J. J. A deterministic Lagrangian-based global optimization approach for quasiseparable nonconvex mixedinteger nonlinear programs. J. Mech. Des. 2009, 131, 051009-1− 051009-8. (11) Li, X.; Tomasgard, A.; Barton, P. I. Nonconvex generalized Benders decomposition for stochastic separable mixed-integer nonlinear programs. J. Opt. Theory Appl. 2011, 151, 425−454. (12) Li, X.; Armagan, E.; Tomasgard, A.; Barton, P. I. Stochastic pooling problem for natural gas production network design and operation under uncertainty. AIChE J. 2011, 57, 2120−2135. (13) Li, X.; Tomasgard, A.; Barton, P. I. Decomposition strategy for the stochastic pooling problem. J. Global Opt. 2011, DOI: 10.1007/ s10898-011-9792-0. (14) McCormick, G. P. Computability of global solutions to factorable nonconvex programs: Part IConvex underestimating problems. Math. Program. 1976, 10, 147−175. (15) Scott, J. K.; Stuber, M. D.; Barton, P. I. Generalized McCormick relaxations. J. Global Opt. 2011, 51, 569−606. (16) Bertsekas, D. P. Nonlinear Programming, Second Edition; Athena Scientific: Cambridge, MA, 1999. (17) Balas, E.; Jeroslow, R. Canonical cuts on the unit hypercube. SIAM J. Appl. Math. 1972, 23, 61−69. (18) Kesavan, P.; Allgor, R. J.; Gatzke, E. P.; Barton, P. I. Outer approximation algorithms for separable nonconvex mixed-integer nonlinear programs. Math. Prog., Ser. A 2004, 100, 517−535.
13711
dx.doi.org/10.1021/ie3003254 | Ind. Eng. Chem. Res. 2012, 51, 13703−13711
