A Generalized Knapsack-Problem Based Decomposition Heuristic for

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

A Generalized Knapsack-Problem Based Decomposition Heuristic for Solving Multistage Stochastic Programs with Endogenous and/or Exogenous Uncertainties Zuo Zeng, Brianna Christian, and Selen Cremaschi* Department of Chemical Engineering, Auburn University, Auburn, Alabama 36849, United States

Downloaded via UNIV OF WINNIPEG on July 3, 2018 at 16:54:41 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Optimization problems with decision-dependent (endogenous) and/or decision-independent (exogenous) uncertainties are commonly observed in the process industry, and multistage stochastic programming (MSSP) is one approach for modeling such problems. However, MSSPs grow quickly and become computationally intractable for real-world size problems due to their space and time complexities. This paper presents a generalized knapsack-problem based decomposition algorithm (GKDA) to efficiently obtain feasible solutions for large-scale MSSPs under endogenous and/or exogenous uncertainties. The GKDA decomposes the original MSSP into a series of knapsack problems and solves these problems at appropriate decision points of the planning horizon. We applied GKDA to obtain feasible solutions for four planning problems, which include continuous and/or discrete decision variables, and endogenous and/or exogenous uncertain parameters. The comparison of solutions obtained by GKDA to the optimum solutions for these problems revealed that the GKDA, in most cases, yielded tight feasible solutions with solution times up to 6 orders of magnitude faster than that of the deterministic equivalents. clinical trial planning problem.3 Uncertainty in whether a drug successfully completes a clinical trial or not is only realized after the corresponding clinical trial is completed, and the decision to push a drug through a clinical trial does not impact its success or failure. On the other hand, in a facility protection problem, the likelihood of a facility failing to deliver goods or services after a disruptive event is an example of type I endogenous uncertain parameter in which the values of the decision variables impact its realized value, which depends on the level of resources allocated as protection to that facility.11 This paper restricts its scope to problems with type II endogenous uncertain parameters. Multistage stochastic programming models with endogenous uncertain parameters (in the presence or absence of exogenous uncertain parameters) consider the full horizon scenarios directly, for which the decision variables are defined independently for each scenario.8 Then, a set of constraints called nonanticipativity constraints (NACs) are introduced to the model to prevent the anticipation of unrealized future outcomes when making decisions in individual scenarios. Unfortunately, the sizes of variables and constraints grow exponentially in these models with the increases in the numbers of scenarios and time periods that define the

1. INTRODUCTION Multistage stochastic programming (MSSP) is one of the approaches used to model and solve optimization problems under endogenous (i.e., decision-dependent) and/or exogenous (i.e., decision-independent) uncertainties, and it has been applied to several planning problems including among others synthesis of process networks with uncertain process yields,1,2 R&D pipeline management,3−5 oil field infrastructure planning,6 and artificial lift infrastructure planning.7 Multistage stochastic programming is a scenario-based approach, and includes multiple stages where decisions are made followed by realizations of uncertainty and potential recourse actions. Scenarios in a MSSP model correspond to possible future states of the system characterized by the uncertain parameters, and they are generated using the allowable combinations of the uncertain parameter outcomes.8 Endogenous and exogenous uncertain parameters in MSSP are defined based on how and when they are realized. Exogenous uncertain parameters are realized at a known fixed time period in the planning horizon irrespective of the values of the decision variables.9 For example, demand is generally considered to be independent of any capacity expansion decisions in process industries, and hence, is regarded as an exogenous uncertain parameter.1,2 In contrast, the values of the decision variables impact the realization times (type II) and/or the realized values of endogenous uncertain parameters (type I).10 An example of a type II endogenous uncertain parameter is encountered in the © XXXX American Chemical Society

Received: Revised: Accepted: Published: A

February 20, 2018 June 8, 2018 June 13, 2018 June 13, 2018 DOI: 10.1021/acs.iecr.8b00822 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

included two technology portfolio examples for which 5 and 10 projects were considered, and the solutions were within 3.4% to 8.3% of the optimal solution. A solution strategy for MSSP models with endogenous uncertainty whose deterministic equivalent is a nonconvex mixed integer nonlinear programming (MINLP) model was presented in Tarhan et al.17 The strategy utilizes global optimization and outer-approximation, and it was used to solve two planning problems, a modified version of the process synthesis problem1,2 and an offshore oilfield planning problem.14 Gupta and Grossmann18 explored improved Lagrangean decomposition frameworks for both process synthesis and oilfield planning problems. One of these frameworks yielded solutions with less than 1% optimality gap in a few iterations. For MSSP models with both endogenous and exogenous uncertainties, Apap and Grossmann8 proposed a sequential scenario decomposition (SSD) approach. The approach selects one scenario from the exogenous uncertainty scenario group at the initial time period, solves a series of subproblems with endogenous uncertainties to determine binary investment decisions, fixes the binary investment decisions to satisfy the first-time period endogenous and exogenous NACs, and solves the resulting model to obtain a feasible solution. The authors used the SSD approach with a Lagrangean relaxation strategy (dual bound) to solve the full-space problem under endogenous and exogenous uncertainties. To test their algorithm, the authors solved the process synthesis problem1,2 and the offshore oilfield planning problem,14 and they concluded that the SSD approach found high-quality feasible solutions several orders of magnitude quicker than solving the full space model. In this paper, we present a generalized knapsack-problembased decomposition algorithm (GKDA) to obtain feasible solutions for MSSP models under endogenous and/or exogenous uncertainties. The algorithm builds upon our previous heuristic, knapsack-problem based decomposition algorithm (KDA), which obtains feasible solutions to MSSP models with endogenous uncertainty.19 Instead of enumerating all uncertain parameter space as scenarios, the KDA decomposes the original multiperiod MSSP into a series of knapsack problems and solves these problems at appropriate decision points of the planning horizon. Our previous results revealed that KDA provided tight primal bounds for the problems studied, and generated these bounds several orders of magnitude faster than solving the deterministic equivalent.19,20 We later developed a branch-and-bound algorithm combining KDA and two dual bounding approaches to solve MSSPs under endogenous uncertainty.21 The algorithm uses progressive hedging (PH), or the solution of the individual scenario problems for generating dual bounds. The primal bounds are generated using the KDA.19 The branch-and-bound algorithm successfully addressed the space complexity of MSSP problems, but required considerable time to converge. The KDA can only be used to obtain feasible solutions for MSSPs under endogenous uncertainty that contain only discrete state variables and recourse actions impacting only objective function values (not constraint satisfaction). This paper generalizes the KDA to locate feasible solutions and to generate primal bounds for MSSP models under endogenous and exogenous uncertainties with both continuous and discrete state variables and recourse actions that may also impact scenario specific constraint satisfaction. We applied the GKDA to four different planning problems under uncertainty from process systems engineering (PSE) literature: new technology investment planning,22 clinical trial planning,3 synthesis of process networks,1,2 and artificial lift infrastructure planning.7 The general MSSP formulation that is

planning horizon, which results in both space and time complexities for real-world size problems. In general, recent literature focuses on developing approximation and decomposition approaches to address the computational complexity of MSSP models under endogenous and/or exogenous uncertainties. Goel and Grossmann12 proposed an approximation approach, which searched a restricted feasible region of the MSSP model to find a “good” solution. Their results revealed these “good” solutions yielded improved objective function values when compared to the optimum solution of the deterministic problem. A Lagrangean duality based branch-andbound algorithm, which is guaranteed to give the optimal solution, was introduced to solve similar problems.13 The results suggested that the branch-and-bound algorithm obtained significantly better solutions and tighter optimality gaps than the heuristic presented in Goel and Grossmann.13 The process synthesis problem with decision-dependent uncertainty in process yields, which are revealed gradually as investment in processes occurred, was modeled as a MSSP model by Tarhan and Grossmann.2 A duality-based branch-and-bound approach was able to solve instances of this problem with 16 scenarios to within 3% optimality gap.2 Tarhan et al.14 incorporated nonlinear reservoir models and gradual realization of uncertainty to the offshore gas field infrastructure planning problem,12 and solved the resulting MSSP models using a duality-based branchand-bound algorithm. The expected net present value of the solutions located by the algorithm for different instances of the problem were up to 22% better than the optimum solutions of the deterministic problems that used the expected values of the uncertain parameters; however, the authors noted that the solution times were “rather long”. Mercier and Hentenryck15 proposed a multistep anticipatory algorithm for solving pharmaceutical R&D pipeline management problems. The algorithm uses a sample average approximation (SAA) approach to generate a Markov Decision Process (MDP), solves the resulting MDP, and returns the greedy solution. The algorithm was tested using 12 instances of the R&D pipeline management problem, and yielded solutions that were 10% better than the dynamic programming equivalents for all instances. The authors concluded that the algorithm was, nevertheless, computationally expensive. Colvin and Maravelias4 explored a rolling-horizon approximation approach to solve large instances of the pharmaceutical clinical trial planning problem. The approach constructs a relaxed MSSP model by only retaining NACs for the stages of the “detailed time block”, solves this relaxed MSSP, and implements its solution to these stages. Then, the uncertainty associated with the implemented decisions are realized. The approach continues to construct and solve relaxed MSSP models rolling the “detailed time block” at each iteration until the end of the planning horizon. The authors were able to successfully solve cases with more than 1000 scenarios. Another approach, a branch-and-cut algorithm,5 initially constructs and solves a relaxed MSSP model that includes a subset of the NACs, and then iteratively adds the NACs that are violated. The authors concluded that the algorithm reduces the number of NACs that should be included in the problem formulation significantly and that this method would be advantageous for any stochastic programming formulation in which the majority of constraints are NACs. Solak et al.16 proposed a SAA algorithm for solving R&D portfolio optimization problems. The algorithm solves SAA problems to generate candidate solutions, and evaluates the quality of first-stage decisions of these solutions using a larger sample set of scenarios. Computational studies B

DOI: 10.1021/acs.iecr.8b00822 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Zis,,ts′ ⇔ H(bis,1, bis,2 , . . . , bis, t , yis,1 , yis,2 , . . . , yis, t )

addressed by the proposed GKDA is given in section 2. Section 3 outlines the details of the GKDA. Section 4 presents the application of the GKDA to find feasible solutions and primal bounds for the planning problems considered. Section 5 presents the results, discusses the solution quality and solution times, and compares them to the optimum solutions and solution times of deterministic equivalents. The last section summarizes the conclusions and the future directions of this work.

∀ (i , s , s′) ∈ SE ,

bis,1

Zis,,ts′ ∈ {0, 1}

yis, t , γts

2. GENERAL MULTISTAGE STOCHASTIC PROGRAMMING FORMULATION UNDER ENDOGENOUS AND EXOGENOUS UNCERTAINTIES A general MSSP formulation under endogenous and exogenous uncertainties, MSSPEX, is given in eqs 1−13. The presented model is adapted from the original model of Goel and Grossman.1 The planning horizon is defined as t ∈ T = {1, 2, 3,.., |T|}. The uncertain parameter θi is associated with the endogenous uncertainty per source i ϵ I. The model has binary state variables denoted by bi,ts, continuous state variables represented as yi,ts, and recourse-action decision variables defined as γst. The realizable values of uncertain parameter θi are represented using a finite set {θ1i ,θ2i ,···,θRi i}. The exogenous uncertain parameter realized at time t is denoted by ξt and it has a finite set of possible realizations {ξ1t ,ξ2t ,···,ξRt }. The full scenario set, S, is constructed using Cartesian product and it is defined as S ≔ {×t ∈ Tξt}×{×i ∈ Iθi}. We define an exogenous scenario pair set (t,s,s′) ∈ SX for scenarios s and s′ that are indistinguishable at time period t because they have the same realizations of the exogenous and endogenous uncertain parameters, i.e., SX ≔ {(t,s,s′):t ∈ T,ξst = ξst′, θsi = θsi ′|i ∈ I }. Similarly, we define an endogenous scenario pair set (i,s,s′) ∈ SE for scenarios s and s′ that are indistinguishable in terms of the realization of endogenous uncertain parameter associated with resource i and that have the same realizations for the exogenous uncertain parameters, i.e., SE ≔ {(i,s,s′):i ∈ I, θsi = θsi ′,ξst = ξst′|t ∈ T }.

(1)

s

s.t. gi , t , s(bis, t , yis, t , γts , θis , ξts) ≤ 0 hi , t , s(bis, t , yis, t , γts , θis , ξts) = 0

∀ i ∈ I, t ∈ T, s ∈ S ∀ i ∈ I, t ∈ T, s ∈ S

(2)

∀ i ∈ I,

∀ (s , s′) ∈ S

(4)

yis,1

∀ i ∈ I,

∀ (s , s′) ∈ S

(5)

=

yis,1′

bis, t = bis,′t

∀ i ∈ I,

∀ t ∈ T,

∀ (t , s , s′) ∈ SX

(6)

yis, t = yis,′t

∀ i ∈ I,

∀ t ∈ T,

∀ (t , s , s′) ∈ SX

(7)

s

s

γt = γt ′

ÄÅ s , s ÉÑ ÅÅ Zi , t ′ ÑÑ ÅÅ ÑÑ ÅÅ ÑÑ ÅÅ s Ñ ÅÅbi , t = bis,′t ÑÑÑ ÅÅ ÑÑ ÅÅ s ÑÑ ∨ [ ¬ Zis,,ts′] ÅÅÅ yi , t = yis,′t ÑÑÑ ÅÅ ÑÑ ÅÅÅ s ÑÑ ÅÅ γ = γ s′ ÑÑÑ ÑÑ t Ö ÇÅÅ t

∀ t ∈ T,

∀ (t , s , s′) ∈ SX

∀s∈S

∀ (s , s′) ∈ S ,

∀ s ∈ S,

(11)

∀ i ∈ I,

∀ i ∈ I,

∀t∈T

∀t∈T

(12) (13)

3. GENERALIZED KNAPSACK-PROBLEM BASED DECOMPOSITION ALGORITHM The generalized knapsack-problem based decomposition algorithm (GKDA) decomposes time periods and scenarios of a MSSP model into a series of continuous or 0−1 knapsack subproblems (KSPs). The algorithm generates and solves these problems based on the realizations of uncertainty along the planning horizon. The general form of each individual KSP is given in eqs 14−17.

(3)

bis,1 = bis,1′

∈

∀ i ∈ I,

(10)

The model MSSPEX consists of five parts: (1) The objective function, eq 1, minimizes the expected value of the total cost associated with state variables bsi,t, ysi,t, and recourse-action decision variables γts, where ps is the probability of scenario s and the function Gs(bsi,t,ysi,t,γst,θsi ,ξst) calculates the total cost for scenario s. (2) Scenario specific constraints are represented via eq 2 for inequality and via eq 3 for equality constraints. These constraints may be used to model resource availabilities, sequencing constraints, material flow constraints, and capacity constraints for each scenario and time period. (3) Initial nonanticipativity constraints given in eqs4−5 ensure that the values of the decision variables associated with both endogenous and exogenous uncertain parameters are identical for the first time period. (4) The NACs for exogenous uncertain parameters are presented in eqs6−8. They are written for each time period for which scenario pairs (s,s′) remain indistinguishable due to exogenous uncertainty; that is, (t,s,s′) ∈ SX. (5) Conditional NACs, disjunction in eq 9, are added to the MSSP formulation for preventing anticipative decisions before two scenarios become distinguishable due to a realization of an endogenous uncertain parameter. The value of the indicator variable, Zi,ts,s′, in the disjunction is determined by eq 10, which is a linear function, represented by H(.), of binary and continuous state variables, bsi,t and ysi,t. If Zi,ts,s′ is true, the scenario pairs (i,s,s′) ∈ SE remain indistinguishable, and hence, the NACs for these scenarios are enforced. In MSSPEX, if functions Gs(.), gi,t,s(.), and hi,t,s(.) are linear, its deterministic equivalent can be formulated as a mixed-integer linear programing (MILP) model; otherwise the deterministic equivalent is a MINLP model.

MSSPEX: A general MSSP formulation under endogenous and exogenous uncertainties min ∑ ps Gs(bis, t , yis, t , γts , θis , ξts)

∈ {0, 1}

∀ t ∈ T, t > 1

KSP: General form of an individual knapsack subproblem

(8)

max: −Gest(bi , yi , γ , E[θi ] , E[ξ ])

(14)

s.t. ∀ (i , s , s′) ∈ SE , ∀ t ∈ T , t > 1 (9) C

gi(bi , yi , γ , E[θi ] , E[ξ ]) ≤ 0

∀i∈E

(15)

hi(bi , yi , γ , E[θi ] , E[ξ ]) = 0

∀i∈E

(16)


Article

Industrial & Engineering Chemistry Research bi ∈ {0, 1}, yi , γ ∈ ,

∀i∈E

and endogenous uncertain parameter values are represented by their respective expected values (E[θi] and E[ξ]) as there have been no uncertainty realizations. The GKDA constructs a single (n = 1) KSP at the root node (t = 1), KSP1,1. The eligible item set (E) for KSP1,1 is constructed by enforcing the initial time period NACs, sequencing constraints, and bounds on decision variables of MSSPEX at t = 1. The objective function and weight constraints of KSP1,1 are constructed. Then, the GKDA sets Nt, which is equal to the number of KSPs that should be constructed at time t, equal to 1. The optimal solution of the root-node KSP1,1 is used to determine the values of the decision variables at the first-time period of the planning horizon. The state variables in MSSPEX are set equal to their optimum values in KSP1,1. The endogenous uncertainties associated with state variables and exogenous uncertainties that are resolved at the first-time period are realized. If the scenario specific constraints of MSSPEX (gi,t,s(.) and hi,t,s(.)) do not contain any uncertain parameters, the values of the recourseaction decision variables determined based on the optimum solution of the KSP are guaranteed to be feasible. For those problems, the recourse-action decision variables are set equal to their values in the optimum solution of KSP1,1. If the constraints gi,t,s(.) and hi,t,s(.) include uncertain parameters, the values of the recourse-action decision variables are calculated using the state decision variable values determined by the optimum solution of the KSP for each unique realization of uncertainty. The new child KSPs for the next time period are generated based on the realized uncertainties. The number of new KSPs generated (Nt+1) depends on the number of realizations. For example, suppose that the optimal solution of KSP1,1 sets a binary state decision variable associated with an endogenous uncertain parameter equal to 1, and further suppose that that endogenous uncertain parameter has two outcomes (H and L), which are realized when the binary state decision variable becomes 1. For the sake of simplicity, assume that there are no realizations of the exogenous uncertain parameters at t = 1. The GKDA then generates two new child KSPs, KSP1,2, KSP2,2, estimates and updates eligible items set E, the objective function and weight constraints for each KSP based on the realizations. For endogenous uncertainties, the GKDA generates child KSPs when the outcomes of selected items are realized. For exogenous uncertain parameters, the GKDA generates child KSPs at their realization times, which are known a priori. After the construction of child KSPs, the GKDA increments the time, and repeats the process of generating child KSPs, solving the KSPs, realizing the uncertainties, and determining recourse actions. The process continues until either the eligible item sets are empty or the end of the planning horizon is reached. The GKDA is guaranteed to yield a feasible solution as it satisfies all constraints of MSSPEX. The KSP generation method combined with recourse-action determination scheme ensures that scenario specific constraints are satisfied. The NACs of MSSPEX are satisfied by how the child KSPs and the eligible item sets are generated. Therefore, at termination, the GKDA provides a feasible solution for MSSPEX. The evaluation of eq 1 at this feasible solution yields a primal bound for the MSSPEX. The quality of the GKDA solution depends on how accurate item values estimate the impact of the decision variables associated with these items on MSSPEX objective function value. The next section discusses how item values are calculated. 3.1. Calculation of Item Values in Knapsack Subproblems. An item value in a KSP approximates the contribution of the decision variable associated with that item to the objective

(17)

The objective function and scenario specific constraints eqs1−3 of MSSPEX are used to construct the objective function and constraints of the KSP (eqs14−16). The objective function of the KSP maximizes the total value of the packed items. The decision variables (bsi,t,ysi,t, and γst) of MSSPEX are represented by indivisible items (bi), which have binary decision variables, and divisible items (yi and γ), which are continuous decision variables in KSP. The GKDA uses the concept of eligible item set E, where items e ∈ E are created by enumerating all allowable decisions at time period t in MSSPEX. The allowable decisions at time period t are determined using the state of the system at time period t, the NACs, the logic of the sequencing constraints (if present in MSSPEX), and the bounds (if present in MSSPEX). The state of the system at time period t is defined using the decision variable values and the uncertainty realizations up to time period (t − 1). The objective function ∑spsGs(bi,ts,yi,ts,γts,θis,ξts) (eq 1) of MSSPEX is estimated in each KSP by −Gest(bi,yi,γ,E[θi],E[ξ]) via taking the expected values of uncertain parameters θsi and ξst as E[θi] and E[ξ]. Generally, we directly transfer the coefficients of objective function associated with item e ∈ E as item values, Ve, in each KSP. There are two conditions for which this approach does not accurately approximate the value of the eligible items in KSPs. We introduce two approaches to improve item values, which are discussed in section 3.1. The inequality constraints gi,t,s and equality constraints hi,t,s in eqs2 and 3 of MSSPEX are directly transformed to gi and hi, i.e., weight constraints, in KSPs by removing both t and s indices. The coefficients of the decision variables in the equality and inequality constraints (eqs2 and 3) of MSSPEX define the item weights We, e ∈ E. It should be noted that the equations corresponding to sequencing constraints and bounds that were used for determining the allowable items are not included in weight constraints of the KSPs. The GKDA is outlined in Figure 1, where the nth KSP at time t is denoted by KSPn,t. In the first KSP generated, both exogenous

Figure 1. Generalized knapsack-problem based decomposition algorithm (KDA). D


Article


s s s variable values of (bi,|T| ,yi,|T| ,γ|T| ). Suppose that there are sequencing constraints in the scenario-specific constraint set that restrict the values of (bsi,|T|,ysi,|T|,γs|T|) on the values of previous decision variables {bsi,t,yi,ts,γst|t ≤ |T| − 1}. Then, this problem satisfies Condition 1. The objective function of a KSP generated at time period t for this problem using the MPGA is given in eq 21, and the item values are calculated based on the maximum expected contribution of each item to the objective function value.

function value of MSSPEX at the stage and time period the KSP is generated. The straightforward approach for calculating an item value is to use the MSSPEX objective function coefficients of the decision variable associated with the item. Suppose the objective function of MSSPEX is given in eq 18. min ∑ ps s

∑ ∑ [(V1,i ,t + V1,s i ,t )bis,t + (V2,i ,t + V2,s i ,t )yis,t i

+ (V3, t +

t

V3,s t )γts]

(18)

max ∑ E[V1, i , |T |]bi + E[V2, i , |T |]yi + E[V3, |T |]γ

(bi,ts,ysi,t,γst),

In eq 18, the decision variables are and there are both deterministic (V1,i,t, V2,i,t, V3,t) and uncertain (Vs1,i,t, Vs2,i,t, Vs3,t) coefficients. Using the straightforward approach, the objective function of a KSP generated at time period t is given by eq 19.

i

In Condition 2, the problem involves capital investment decisions with revenue streams associated with these investments in the future. For example, in planning offshore gas field developments,12 the decision to invest on a dedicated well platform at a field for producing gas can be made at any time period and the revenue associated with the produced gas will contribute to the overall profit for the rest of the planning horizon. Similar to Condition 1, due to the decomposition of time periods in GKDA, the KSPs cannot consider the impact of future revenue streams associated with the capital investment decisions. For these problems, the item costs (or negative values) associated with capital investment decision variables are calculated by converting the capital investment to annual capital charge for the remaining planning horizon.

min ∑ (V1, i , t + E[V1,s i , t |D])bi + (V2, i , t + E[V2,s i , t |D])yi i

+ (V3, t + E[V3,s t |D])γ

(19)

The deterministic coefficients at time period t (V1,i,t, V2,i,t, V3,t) contribute directly to the item values. To estimate the contris butions of the uncertain coefficients (Vs1,i,t, V2,i,t , Vs3,t), we calculate their conditional expectations at time period t where D = s {bsi,u,yi,u ,γsu,u ≤ t−1} (eq 19). In other words, the expected values are calculated based on the decision variable values up to the previous time period and associated endogenous and exogenous uncertain parameter realizations. Under certain conditions, the straightforward approach for generating item values may not accurately approximate the expected contribution of some decision variables to the MSSPEX objective function value in KSPs. In such cases, the GKDA may yield a poor feasible solution. We define two such conditions and propose different approaches to approximate their item values. • Condition 1: IF gains/revenues are only realized at certain stages or are only associated with certain decision variables, THEN generate the item value using the maximum potential gain approach (MPGA). • Condition 2: IF an item represents a capital investment decision variable in MSSPEX, THEN generate the item value using annual capital charge. Condition 1 considers the class of planning problems for which the decision maker cannot claim revenue or gains until a certain stage or until a certain outcome is realized in the planning horizon. For example, in financial planning and control problem (asset-liability management),9 the decision maker sequentially makes investment decisions at every stage; however, the gains are only associated with last state variables. This problem satisfied Condition 1, where all investment decisions implicitly contribute to the overall gain. The clinical trial planning problem also satisfied Condition 1, where the revenue associated with a drug is only realized once the drug successfully completes the last (i.e., third) clinical trial. With the maximum potential gain approach (MPGA), we calculate an item’s value based on the maximum expected contribution of the associated decision to the MSSP objective function value considering the remaining planning horizon. Suppose the objective function of MSSPEX is given as eq 20. min ∑ ps s

∑ {−V1,i ,|T|bis,|T| − V2,i ,|T|yis,|T| − V3,|T|γ|Ts|} i

(21)

4. APPLICATION OF GKDA TO FOUR PLANNING PROBLEMS We applied GKDA to four different planning problems: new technology investment planning,22 clinical trial planning,3 process network synthesis,1,2 and artificial lift infrastructure planning.7 For each problem, the following subsections present in detail how to construct KSPs, to determine eligible item lists, objective functions and weight constraints, and to implement the algorithm. 4.1. New Technology Investment Planning (NTIP) Problem. Givens for the new technology investment planning (NTIP) problem are a fixed length planning horizon divided into equal time periods [t ∈ T], a set of chemical processing technologies [i ∈ I], a set of chemicals [n ∈ N], and a set of technology maturity stages [sg ∈ {laboratory(L), pilot plant (PP), commercialization (C)}]. The objective is to determine the technology investment plan that yields the expected minimum total cost. Demand at each time period should be met either via production or via purchasing at market price. The yield of a technology is not known until the technology reaches the commercialization stage. To reach a stage, the installed capacity of a technology must reach the minimum capacity threshold for that stage. A technology may be abandoned at the completion of any stage. Following Colvin and Maravelias,3 we define the realizable values of uncertain parameter associated with project abandonment with set Ψi = {L(A),PP(A),PP(P)}, where PP(A) corresponds to the outcome that technology i is abandoned at pilot plant stage, and PP(P) corresponds to the outcome that technology i reaches the commercialization stage. There are four uncertain parameters associated with each new technology: (1) project abandonment (ψi), (2) process yield (χi), (3) learning-by-searching elasticity (αi), and (4) learningby-doing elasticity (βi). The learning elasticities define how the capacity expansion cost may change with investments in a twofactor learning curve.22 Demand for products is uncertain for time periods beyond the first, and there is an uncertain parameter

(20)

In eq 20, the decision variables are (bsi,t,ysi,t,γst), and only at the end of the planning horizon, the revenue is realized based on the decision E


Article


variable Xi. As such, the equation is incorporated as a constraint to the KSPs (eq 23). In eq 23, P(ai,bi) is the joint probability distribution of uncertain parameters ai ∈ αi (learning-by-searching elasticity) and bi ∈ βi (learning-by-doing elasticity), and CXi,t−1 and RDi,t−1 represent the cumulative installed capacity and research expenditure at the previous time period. The GKDA calculates these values as the sums of Xi and RDXi, respectively, for all parent KSPs. The scenario specific constraints of NTIP problem are included in KSPs as eqs24−30. Similar to CXi,t−1, the value of Yi,sg,t−1 is calculated based on the solutions of the parent KSPs. The resulting KSPs are MINLP models. After solving the root node KSP, the GKDA determines Nt+1 based on the realized uncertain parameters associated with the decisions made in the root node KSP, and constructs the corresponding child KSPs according to the realizations. Because the constraints gi,t,s(.) and hi,t,s(.) of the NTIP problem MSSP include uncertain parameters, the GKDA determines the values of the recourse-action decision variables using the state decision variable values determined by the optimum solution of the KSP for each unique realization of uncertainty for the realized parameters. The GKDA repeats the process until the end of the planning horizon. 4.2. Clinical Trial Planning Problem. The goal of the clinical trial planning problem is to identify the schedule of clinical trials for a set of candidate drugs, i ∈ I = {1,2,···,|I|}, that maximizes the expected net present value (ENPV) of the clinical trial pipeline over a planning horizon (t ∈ T). A drug has to complete three clinical trials in order (j ∈ J = {PI,PII,PIII}), and it can either pass (P) or fail (F) them. The binary decision variable xi,j,t,s is equal to 1 if drug i clinical trial j is started at time period t under scenario s, and zero otherwise. The outcomes of the clinical trials are uncertain and decision-dependent because the outcome of a clinical trial is realized once the trial is completed. The clinical trial outcomes are defined with set Θi = {PI(F),PII(F),PIII(F),PIII(P)}, where PIII(P) corresponds to the outcome that drug i passes trials PI, PII, and PIII. Clinical trial j of drug i has a known cost Cij, required resources ρijr (where r ∈ R = {1,2,···,|R|}), and fixed duration τij. After successfully completing all clinical trials, the revenue from drug i is realized. Potential revenue is defined using the maximum revenue, revmax i , associated with drug i. The patent life of a drug is assumed to start shrinking once the drug enters the pipeline, and associated losses are represented by two penalty terms: γDi (loss of market share) and γLi (loss of patent life). A slightly revised version of the MSSP formulation of Colvin and Maravelias5 and the corresponding nomenclature is given in Supporting Information. The GKDA starts by generating the eligible item set E for the root node KSP. The items k ∈ E are created by enumerating all possible decisions under the logic of the sequencing constraints for the clinical trials defined by the MSSP model. In KSPs, the decision variable xi,j,t,s is identified by xk, where k corresponds to a specific drug-clinical trial pair. The objective function of the original MSSP is given in eq 31, which reveals that the revenue associated with drug i can only be realized once the drug completes its last clinical trial (PIII) successfully, that is, the revenue is only associated with certain decision variables, i.e., Xi,PIII,t,s, which satisfies Condition 1. Thus, the item values are generated using MPGA, viz eq 32. They are generated considering the potential revenue associated with drug i, which is calculated by taking the maximum revenue associated with drug i and deducting the losses due to loss of active patent life and clinical trial costs.

(Dn,t) representing the demand of product n at time t. Hence, the NTIP problem contains endogenous (technology development, yield, and learning elasticities) and exogenous (demand) uncertainties. The NTIP problem has one binary state decision variable, the technology maturity stage (Yi,sg,t,s), and two continuous state variables, the level of capacity expansion (Xi,t,s), and the level of research expenditure (RDXi,t,s) for technology i at time period t in scenario s. The recourse actions include production rate (Mi,n,t,s) of material n using technology i at time period t in scenario s. The MSSP model for the NTIP problem, the deterministic equivalent of which is a MINLP model, was developed by Christian and Cremaschi.22 It is provided in the Supporting Information along with its nomenclature for completeness. The first step of the GKDA is the identification of eligible items for the root-node KSP. The items for the NTIP problem include all current-stage decision variables and recourse actions, the technology maturity stage (Yi,sg), the level of capacity expansion (Xi) and of research expenditure (RDXi) for each technology, and the production rate of each material using each technology (Mi,n). In KSPs, the variable Yi,sg is a binary decision variable, while the variables Xi, RDXi, and Mi,n are continuous decision variables. The item costs (which can be seen as the negative of item values) are the cost of capacity expansion for technology i (CCi), and purchasing price of material n (MCstn). The KSP objective function is constructed using the straightforward approach, and the general form of the KSPs for the NTIP problem is given in eqs22−30. ij max:−∑ MCstn·jjjjE[Dn] − j n k −

yz z {

∑ γi ,n,PR ·E[χi ]·Mi ,nzzzz i

∑ E[CCi]·Xi − ∑ RDXi i

(22)

i

s.t. E[CCi ] =

i CXi , t − 1 + Xi yz i zz · zz CX i ,0 k {

∑ ∑ P(ai , bi)·CCOi ·jjjjj

ij RDi , t − 1 + RDXi yz jj zz jj zz RD i ,0 k { ai ∈ αi bi ∈ βi

ai

Mi , n ≤ Yi ,3(CXi , t − 1 + Xi)E[ψi ] CXi , t − 1 + Xi ≥

CXimin , sg · Yi , sg

CXi , t − 1 + Xi ≤

min ∑ (CXimin , sg − CXi , sg + 1) · Yi , sg sg < 3

b

(23) (24) (25)

(26)

Yi , sg ≤ Yi , sg − 1

(27)

Yi , sg ≤ Yi , sg , t − 1

(28)

Xi ≤

Ximax

RDXi ≤

RDXimax

(29) (30)

In eq 22, the parameter γi,n is the ratio of stoichiometric coefficient of the chemical n and the primary reactant/feedstock PR for technology i multiplied by the ratio of the molecular weights. The values of expected demand E[Dn] and expected yield E[χi] are straightforward to calculate. For calculating E[CCi], the values of the outcomes depend on the value of the F


Article


Figure 2. Schematic of a simple process network1

These four KSPs represent the four possible outcomes of trial (A, PI) and (B, PI): (1) both drugs pass the clinical trials, (2) trail (A, PI) fail while trial (B, PI) passes, (3) trail (A, PI) passes while trial (B, PI) fails, and (4) both drugs fail clinical trial PI. The algorithm then updates the eligible item sets, objective functions, and weight constraints for all four KSPs. Because scenario specific constraints of clinical trial planning MSSP model (gi,t,s(.) and hi,t,s(.)) do not contain any uncertain parameters, the values of the recourse-action decision variables determined based on the optimum solutions of the KSPs are guaranteed to be feasible. The algorithm continues until either the eligible item sets are empty or the end of the planning horizon is reached. A detailed discussion of KDA and its application to this problem can be found in Christian and Cremaschi.19,20

l o

ENPV = f (Xi , j , t , s) =

∑ ps omoo∑ ∑ [revimaxXi ,PIII,t ,s s

n

t ′≤ t jij jj−X jj i , j ,1, s + ∑ Xi , j − 1, t − τi ,j−1, s ′ j j = PII,PIII j t ′> τi , j − 1 k yz t ′≤ t z ∑ Xi ,j ,t′ ,szzzzz − γi L(t + τi ,PIII)Xi ,PIII, t , s] z t′ { ij open j ∑ ∑ revi ,j fi ,j jjjjj−Xi ,PIII,t ,s j i j ≠ PI k yz t ′≤|T | t ′≤|T | z ∑ Xi ,j − 1,t′ − τi,j−1,s − ∑ Xi ,j ,t′ ,szzzzz z t ′> τi , j − 1 t′ { t ′≤|T | ij yz jj1 − ∑ X zz ∑ reviopen , j fi ,PI j i ,PI, t , s z ′ zz jj ′ i t k {

−γi D

−

+

+

+

+∑

i

t

∑

∑

∑

| o o o − ∑ cdt Ci , jXi , j , t , s} o o o i ,j,t o ~ i

Vk(i , j), t

max:

∑ Wk ,rxk ≤ Wrmax

∀r∈R

k

xk ∈ {0, 1}

∀k∈E

(34) (35)

4.3. Process Network Synthesis under Uncertain Process Yields. The goal of the process network synthesis problem is to determine optimal capacity expansion investments and process operation decisions that will maximize the total expected profit for the process network shown in Figure 2 from the sales of product A. The decisions determine which (if any) of the processes (i ∈ {Process 1, Process 2, Process 3}) should be installed or expanded (yexp,s ∈ {0,1}), which processes i,t should be operated (yoper,s ∈ {0,1}), the capacity levels of the i,t ≥0 ,s ∈  ) expansions (wiQE , the flow rates throughout the ,t

revirun , j , t fi , j + 1 Xi , j , t , s

ÉÑ yz ÑÑÑ ∑ τi ,j′zzzzz − cdtCi ,jÑÑÑÑÑ z ÑÑ j ′> j ÑÖ {

(33)

s.t.

j ∈{PI,PII} t >|T |− τi , j

ÄÅ ÅÅ ij ÅÅ j = EÅÅÅrevi max − γi Ljjjt + jj ÅÅ ÅÅÇ k

∑ Vkxk k∈E

(31)

(32)

s ∈ ≥ 0 ), purchases (xtpurch, s ∈ ≥ 0 ), sales network (wkrate, ,t

Item weight Wk,r is equal to the resource requirement of the associated drug-clinical trial pair for resource type r. The |R|-dimensional capacity vector {Wmax r } contains the available resources, ρmax r . The general form of KSPs for the clinical trial planning problem are given in eqs 33−35, and they are |R|-dimensional 0−1 knapsack problems. For this problem, the GKDA is identical to the KDA.19,20 Equation 34 is the scenario specific constraints of the original MSSP other than sequencing constraints, which are used for generating eligible items along with NACs. Once the root node KSP is solved, the number of child KSPs generated and their generation times are based on the outcomes of the root node KSP solution. For example, suppose the solution of the root node KSP selects two items associated with first clinical trials (PI) of drugs A and B, i.e., (A, PI) and (B, PI), a total of 22 = 4 child KSPs are generated at t = max(τA,PI,τB,PI) .

(xtsales, s ∈ ≥ 0 ), and levels of inventory (wtinv, s ∈ ≥ 0 ) of product A for satisfying the demand (dt) at each time period (t ∈ T) under each scenario (s ∈ S). The endogenous uncertain parameters of the problem are the yields of Process 1 and Process 2, which are represented by (θsi ) and are only resolved once the process is installed. The demand for the planning horizon is assumed to be known, and hence there are no exogenous uncertain parameters. The original MSSP formulation and the corresponding nomenclature developed by Goel and Grossmann1 for this problem is given in the Supporting Information. The GKDA starts by identifying the eligible items for constructing the root-node KSP. The eligible item set E includes oper expansion decisions (yexp i ), operation decisions (yi ), capacity QE levels of expansions (wi ), flow rates throughout the network purch (wrate ), sales (xsales), and level of inventory k ), purchases (x inv oper (w ) for product A. In KSPs, the decision variables yexp i and yi G


Article

Industrial & Engineering Chemistry Research rate purch sales are binaries, while the variables wQE , x , and winv are i , wk , x continuous. Note that KSPs constructed only consider the current time period, and the index for time period is not needed in the model. The NACs in GKDA are enforced by the KSP generation scheme, and hence, the scenario index is also not necessary for the decision variables of KSPs. The objective function of the original MSSP is given in eq 36. s s ,s rate, s purch, s max: ENPV = f (yiexp, , yioper, , wiQE , xtsales , s , wtinv, s) , t , wk , t , xt ,t ,t ÄÅ ÅÅ Å s oper, s ,s = ∑ ps ÅÅÅÅ− ∑ ∑ (FEiyiexp, + FOy + VEiwiQE i i ,t ,t ) ,t ÅÅ t i s ÅÇ ÑÉÑ ÑÑ s − ∑ ∑ VOk wkrate, − ∑ (αt xtpurch, s − βt xtsales, s − γtwtinv, s)ÑÑÑÑ ,t ÑÑ t k t ÑÖ

(yexp i,t ,

Equations 38−45 are the equality constraints of the KSPs (eq 16), and they ensure mass balance throughout the network, calculate inventory level, and enforce demand satisfaction. Equations 46−54 are the inequality constraints of KSPs (eq 15), and they constrain flow rates throughout the network either to the capacity levels of the processes or to the upper and lower limits defined by the problem, and ensure that only installed capacities for processes are operated. The KSPs use the expected values of the uncertain parameters, E[θ1] and E[θ2], for the yields of Process 1 and Process 2 before they are realized. When the solution of a KSP recommends installing and operating Process 1 or Process 2, the yield of the installed process is realized. Because the constraints gi,t,s(.) and hi,t,s(.) of the process network synthesis MSSP model include uncertain parameters, the values of the recourse-action decision variables are calculated using the state decision variable values determined by the optimum solution of the KSP for each unique realized yield value outside of the KSPs (see Figure 1). 4.4. Artificial Lift Infrastructure Planning Problem. During the typical lifetime of a horizontal gas well, when the gas flow rate drops below a critical value, that is, below the loading flow rate, liquid accumulates at the bottom of the well, and gas production stops.23 For the remainder of its lifetime, the well should be operated using one or more of the artificial lift methods (ALMs) to maintain production. Multiple ALMs may be installed in horizontal wells for achieving desirable production performance. An artificial lift infrastructure planning problem (ALIP) determines which ALM(s) to install and their operating windows throughout the planning horizon. The artificial lift infrastructure plan depends on the production and reservoir conditions, the completion details, and field information, some of which change during production. The ALIP is further complicated due to the highly stochastic nature of well production: before the installation of an ALM, the production rate using that ALM is uncertain and can only be realized after its installation. Hence, the problem contains endogenous uncertainty. The ALIP MSSP model contains one binary decision variable yi,t,p,s. The variable yi,t,p,s is equal to 1 if ALM i ∈ I is installed and operating at time period p ∈ T and continues to operate until time period t ∈ T under scenario s ∈ S, and zero otherwise. The model includes nine candidate ALMs (|I| = 9), which are plunger lift (PL), foam lift (FL), well head compression (WHC), velocity strings (VS), gas lift (GL), electrical submersible pump (ESP), sucker rod pump (SRP), progressing cavity pump (PCP), and jet pump (JP). It assumes that the well is loaded at the first time period, and the first ALM is installed at the second time period. The original MSSP model for the ALIP was introduced in Zeng and Cremaschi,7 and is given in the Supporting Information along with its nomenclature. The first step of GKDA is to construct the root node KSP. The items i ∈ E for the root node KSP are created by enumerating all allowable decisions at the second time period of the planning horizon under the logic of the sequencing constraints of the ALIP MSSP model. Therefore, the eligible item set of the root node KSP contains all candidate ALMs. The decision variable of KSP, yi, represents ALM i that is installed and currently operating. The item value is calculated using the expected net present value of the potential revenue and costs of installing and operating the associated ALM i for the rest of the planning horizon. As this item includes a capital investment decision, it satisfies Condition 2. Therefore, in item value calculations, the capital cost CCi,t at time t for ALM i is converted to annual capital charge using annual capital charge ratio ε. Figure 3

(36)

wQE i,t )

The decision variables and the associated items (yexp i,t , QE wi ) deal with capital investment decisions, and the model falls under Condition 2. Hence, the item values for these decision variables are generated using the corresponding annual capital charge. For the rest of the items in the KSPs, the straightforward approach is used to generate the item values. The objective function (eq 37) and the general form of the KSPs constructed by the GKDA for the process network synthesis problem are given in eqs37−54. In eq 37, the parameter ε corresponds to the annual capital charge ratio. oper max:− ∑ (εFEiyiexp + FOy + εVEiwiQE) − i i i

∑ VOkwkrate k

− (αx

purch

− βx

sales

inv

+ γw )

(37)

s.t. w3rate = E[θ1]w1rate

(38)

w4rate

(39)

=

E[θ2]w2rate

w8rate, = θ3w7rate

(40)

w5rate = w3rate + w4rate

(41)

w7rate = w5rate, s + w6rate

(42)

inv w inv = wpervious + (w8rate + x purch − x sales)

(43)

x sales = d

(44)

wicap = wicap,previos + wiQE

(45)

w3rate ≤ w1cap

(46)

w4rate

w2cap

(47)

w8rate ≤ w3cap

(48)

LiQEyiexp ≤ wiQE ≤ UiQEyiexp

(49)

L1inflow y1oper ≤ w3rate ≤ U1inflowy1oper

(50)

L 2inflow y2oper ≤ w4rate ≤ U2inflowy2oper

(51)

L3inflow y3oper ≤ w8rate ≤ U3inflowy3oper

(52)

yiexp + y exp_previous ≥ yioper

(53)

yioper

(54)

≤

≥

yiexp

H


Article


Figure 3. Objective function of MSSP and item values in KSPs for artificial lift infrastructure planning problem (see Supporting Information for nomenclature).

summarizes how the item values are calculated for KSPs generated at time period t based on the original MSSP objective function. The weight of each item is equal to 1, and at most one ALM can be operational at any time period, which defines the capacity constraints of the KSPs (eq 56). The general form of the KSPs generated at time period t used for solving the ALIP using GKDA is given in eqs55−57. The KSPs for this problem are 0−1 knapsack problems.

max:

∑ Vi ,tyi i∈E

∑ yi ≤ 1

(56)

i∈E

yi ∈ {0, 1}

∀i∈E

(57)

The GKDA solves the root node KSP, which results in the realization of production rate uncertainty for the selected ALM. Next, the GKDA determines Nt+1, which is equal to the number of realizations of the production rate uncertainty for the selected ALM, and generates Nt+1 child KSPs. Because the scenario specific constraints of ALIP MSSP model do not contain any uncertain parameters, the values of the recourse-action decision variables determined based on the optimum solutions of the KSPs are guaranteed to be feasible.

(55)

s.t. I


Article


Figure 4. Network topologies of the NTIP problem instances.22

5. RESULTS AND DISCUSSION We apply GKDA to generate feasible solutions for different instances of the four planning problems. To assess the GKDA solution quality, deterministic equivalents of the MSSPs given in the Supporting Information S.1−S.4 are solved for the instances when possible. All models and algorithms were implemented in Pyomo and solved using CPLEX 12.6.3 on a standard node of Auburn University Hopper Cluster. Data for all instances are available upon request from the corresponding author. 5.1. New Technology Investment Planning Problem. We solve four instances of the NTIP problem with different network topologies (Figure 4).22 The first instance (Figure 4a) compares a new technology at a Laboratory stage, T1, and a mature technology, T2, at the commercialization stage. The second one (Figure 4b) compares two laboratory stage technologies (T1 and T2). Instance 3 (Figure 4c) considers a retrofitting problem, where two of the technologies (T2 and T4) are at the laboratory stage, and two are mature (T1 and T3). The last NTIP problem instance considers a segment of the chemical process industry, where the desired product is ethylene (Figure 4d).24 Ethylene can be produced from naphtha using mature cracking technology and from biomass using two routes, fermentation−catalytic dehydration or gasification−catalytic conversion−catalytic dehydration. We assume that fermentation and

catalytic dehydration are mature, and that gasification and catalytic conversion are at laboratory stage. The planning horizon for all instances is three years discretized into one-year time periods. The deterministic equivalent of the NTIP problem MSSP yields a large-scale nonconvex MINLP model. For example, the MINLP model of the first instance has 1728 binary and 9409 continuous variables, and 46 561 constraints, and the MINLP model of the last instance has 1.983 million variables, of which 89 856 are binaries, and 3.905 million are constraints. These problems cannot be solved to global optimality using commercial global solvers. A relaxed MILP model was constructed using piecewise linear relaxation25 and exact linearization,26 and solved to obtain dual bounds.22 The GKDA provides the primal bounds. The results are summarized in Table 1. For problem instance 1, the GKDA solved 217 MINLP KSPs to optimality generating a feasible solution with an objective value of $151.8 billion. Each KSP required an average of 0.01 CPUs to solve and resulted in a total algorithm time of 40 CPUs. The dual bound is $146.4 billion, yielding a relative gap of 3.4%. The GKDA solution recommends installing 6 megatons (Mt) per technology at the first time period despite the possibility of T2 failing to be developed. The GKDA determines investment decisions based on the expected outcomes. For this problem instance, the expected success is high enough to justify a full J


Article

Industrial & Engineering Chemistry Research Table 1. Primal (GKDA) and Dual (Relaxed MILP) Bounds, Solution Times, And Relative Gaps for NTIP Problem Instances primal bound, GKDA

instance 1 instance 2 instance 3 instance 4

Table 2. Main Characteristics, ENPVs, and MSSP Solution Times of the Clinical Trial Planning Problem Instances

dual bound, MILP relaxation

ENPV ($)

solution time (CPUs)

ENPV ($)


gap (%)

152 101 6.88 168

40.0 1.1 0.6 1.7

147 91.3 6.88 161

4 1074 138 2313

3.4 9 0 4.1

development of T2. At time period two, the solution obtained by the GKDA recommends maximum research and capacity expansion investments in T2 if the development of T2 is successful. If the development of T2 is unsuccessful, the algorithm suggests installing 6 Mt of T1 capacity. At time period three, for scenarios where the development of T2 is successful, the GKDA solution selects to install capacity for T2. The level of capacity expansion depends on the realized yield value but is approximately 2 Mt. For scenarios where T2 is not successful, the GKDA solution installs an additional capacity of 6 Mt of T1. The GKDA yielded a feasible solution in 1.1 CPUs for the second instance (Table 1). The algorithm solved four MINLP KSPs, and the relative gap is 9%. The GKDA solution does not recommend any investments in either technology. In the MILP model, the capacity expansion cost is underestimated thus the value of expanding the undeveloped technologies is higher. The probability of success is not high enough to offset the potential loss in the KSPs. The third instance had a primal objective function value of $6.884 billion, which is equal to the relaxed MILP solution, and hence the GKDA solution is optimal. The optimal investment strategy is to use the existing capacities to produce the chemicals. The GKDA found the optimal solution in 0.6 CPUs by solving eight MINLP KSPs. The difference in the solution times among the three instances demonstrate that the solution time of the GKDA depends on the number of KSPs solved, which is a function of the uncertainty realizations that occur during the planning horizon. When applied to the fourth instance, the GKDA solution has an objective value of $168 billion. The dual bound is $161 billion, which yields a percent gap of 4.1%. The GKDA solved 19 MINLP KSPs with a total runtime of 1.7 CPUs. The solutions are very similar. The MILP solution recommends producing ethylene using the cracking technology, and purchasing ethylene when necessary to satisfy demand. No capacity expansions are recommended. The GKDA solution also does not recommend installing any new capacity. However, it recommends using existing fermentation and dehydration technologies to produce ethylene from biomass. 5.2. Clinical Trial Planning Problem. In this section, we apply GKDA to find feasible solutions to instances of the clinical trial planning problem with three-, four-, five-, six-, seven-, and ten products, and three clinical trials with different planning horizon lengths. We also solved the deterministic equivalents of the MSSPs for these instances to 0.1% optimality gap when possible. Number of trials, time periods, and scenarios for the instances are summarized in Table 2 along with the optimum ENPVs and the MSSP solution times for the instances that were solved to optimality. Table 2 illustrates the exponential growth in the number of scenarios with the increases in the number of drugs. With this growth, the MSSP models quickly become computationally intractable due to their space and time com-

case name

number of trials

number of time periods

number of scenarios

ENPV ($M)

MSSP solution time (CPUs)

3-prod 4-prod 5-prod 6-prod 7-prod 10-prod 15-prod

3 3 3 3 3 3 3

12 6 6 6 8 10 10

64 256 1 024 4 096 16 384 1 048 576 1 073 741 824

1189 1696 2082 2450

4 7 185 29 357

plexities. As expected, the MSSP solution times grow exponentially with the number of scenarios. For example, CPLEX 12.6.3 required more than eight CPU hours to solve the six-product instance. The MSSP models of seven-, and ten-product problems cannot be generated in RAM, and hence cannot be solved using CPLEX 12.6.3 due to their space complexities. The ENPVs of the GKDA solutions, the percent gaps from the optimum solutions (when available), and the GKDA solution times are compiled in Table 3. The percent gaps ranged from 0.9% to 1.8%. The GKDA obtained these feasible solutions very Table 3. Objective Function Values of the GKDA Solutions, the GKDA Solution Times, and the Corresponding Relative Gaps from the Optimum Solutions (Where Available) for the Clinical Trial Planning Problem Instances case name

ENPV of the GKDA solution

relative gap


3-prod 4-prod 5-prod 6-prod 7-prod 10-prod

1178 1677 2052 2407 2874 4078

0.9% 1.1% 1.5% 1.8% NA NA

1.3 1.5 4.7 3.4 20.0 1579.3

quickly, in 1.3 to 3.4 CPUs. For seven-, and ten-product instances, whose deterministic equivalents cannot be generated and solved, the GKDA yielded feasible solutions in 20.0 and 1579.3 CPUs. 5.3. Process Network Synthesis Problem with Uncertain Process Yields. The problem instance solved is derived from Goel and Grossmann,1 and its schematic is given in Figure 2. The planning horizon is 10 years discretized into one-year time periods. The realizations for the yield of Process 1 are 0.69 or 0.81 with equal probabilities. The uncertain parameter for yield of Process 2 has two equally likely realizations as well, 0.65 and 0.85. Process 3 has an initial capacity of 3 × 106 tons/year and a fixed yield of 0.7. The deterministic equivalent of the MSSP contains 400 binary variables, 685 continuous variables, and 5753 constraints. This is a relatively small MILP and CPLEX 12.6.3 solves it in 1 CPUs. The optimum ENPV is $11.459 × 106. The solution recommends installing Process 2 with a capacity of 10 × 106 tons/year at the first time period and expanding the capacity of Process 3 to 8 × 106 tons/year at the second time period. The ENPV of the GKDA solution is $9.159 × 106, which yields a percent gap of 20.0%. The capital charge rate used is t calculated by ir(1 + tir) , where ir = 0.1. The GKDA solution (1 + ir ) − 1

recommends installing 2.85 × 106 tons/year of Process 2 capacity in the second time period. The capacity of Process 2 is further expanded to 5.71 × 106 tons/year in the third time period, and to 10 × 106 tons/year in fourth time period. In the K


Article


5.4. Artificial Lift Infrastructure Planning Problem. We constructed a hypothetical instance of an ALIP based on a real-world field application: Woodford shale horizontal well. The production before the loading indicated a gas rate of 800 thousand standard cubic feet per day (Mscf/D), a liquid rate of 520 barrels per day (BPD). The well is currently operated using gas lift with a well-site compressor.27 The hypothetical instance assumes that the planning was carried out at the first loading prior to the installation of gas lift. The uncertain production rate parameter for each ALM is assumed to have two realizations, a High outcome, which is assumed to be 20% above its nominal value, and a Low outcome that is 20% below with equal probabilities. For the hypothetical instance, the operating envelope of SRP is modified by setting the minimum allowable flow rate to 300 BPD. The problem is solved for planning horizon lengths of 12, 16, 24, and 36 months (time periods). The maximum CPU time is limited to 20 h. The MSSP models are solved to a 0.1% optimality gap except for 24- and 36-month planning horizons, which cannot be solved within the allowed time limitations. Table 4 summarizes the optimum ENPVs and the GKDA solution ENPVs along with the solution times. The decision trees are presented in Figure 5. The decision trees of the

Table 4. Optimum ENPVs, MSSP Solution Times, ENPVs of the GKDA Solutions and GKDA Solution Times for Hypothetical ALIP Instances with Different Planning Horizons MSSP

GKDA

planning horizon (months)

ENPV ($)


ENPV ($)


12 16 24 36

2,783,894 3,416,506 NA NA

21 4659 72000 72000

2,783,894 3,416,506 4,260,763 5,252,280

0.2 0.4 0.7 1.3

solution, the capacity of Process 3 is expanded twice, to 4 × 106 tons/year in the third time period and to 8 × 106 tons/year in the fourth time period. The GKDA solution never recommends holding inventory because their item values in KSPs are negative, and hence the current item generation rules of the GKDA fails to correctly assess the potential benefit of holding inventory for future time periods. The differences in inventory levels along with the capacity expansion decisions explain the relatively large percent gap between the GKDA and optimum solutions for this problem.

Figure 5. Decision trees of the optimum (a,b) and the GKDA solutions (a−d) for the ALIP instances. L


Article


avenues for improving the GKDA primal bounds, such as introducing additional weight constraints for directing the GKDA solutions away from the greedy one and improving item value estimates in KSPs. We also plan to investigate the potential for using the GKDA as the primal bound generation approach with different dual bounding methods and iteration tactics to solve large-scale MSSPs under endogenous and exogenous uncertainties.

optimum and the GKDA solutions are identical for 12-month (Figure 5a) and 16-month (Figure 5b) planning horizons. Both solutions recommend installing SRP at the second time period to unload the well initially. If the realized production rate using SRP is the High value, then the solution recommends operating SRP until the end of the planning horizon. On the other hand, if the realized production rate is the Low value for SRP, then it is replaced by ESP at t = 7. The solution recommends this switch for the Low value because the operational envelope of SRP is violated before the end of the planning horizon is reached for these scenarios. It should be noted that the optimality of the GKDA solutions for this problem cannot be guaranteed. The decision trees for the 24-month and 36-month planning horizons are given in Figure 5c,d. For longer planning horizons, even for scenarios with High production rate realization of SRP, its operational envelope is violated before the end of the planning horizon. For these problems, the solution recommends installing ESP for scenarios with High SRP production rate realization. For scenarios with Low SRP production rate realization, the ALM selected is Gas Lift. A comparison of the optimum ENPVs to the ENPVs of the GKDA solutions in Table 4 reveals that the GKDA yielded the optimum solutions for the 12- and 16-month planning horizons. The MSSP of the ALIP does not have any resource constraints and allows at most one ALM to be installed and operational at any given time period of the planning horizon. These constraints are easily enforced by the GKDA. Furthermore, the item values in KSPs are able to capture the impact of each ALM’s operation fairly accurately at any given time period. It is important to note that the GKDA obtains these solutions with solution times up to six orders of magnitude faster than solving the deterministic equivalent of the ALIP MSSP.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.8b00822. Multistage stochastic programming models and nomenclature of new technology investment planning,22 clinical trial planning,6 synthesis of process networks1 and artificial lift infrastructure planning7 problems (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Selen Cremaschi: 0000-0002-9333-344X Notes

The authors declare no competing financial interest. Biographies

6. CONCLUSIONS This paper introduces a new generalized knapsack-problem based decomposition algorithm, GKDA, to quickly obtain feasible solutions and primal bounds for large-scale multistage stochastic programming (MSSP) models under endogenous and/or exogenous uncertainties. The GKDA decomposes time periods and scenarios of an MSSP model into a series of continuous or 0−1 knapsack sub-problems (KSPs). The algorithm, then, generates and solves these KSPs based on the realizations of uncertainty. Therefore, the GKDA addresses both time and space complexities associated with solving large-scale MSSPs under endogenous and/or exogenous uncertainties. The GKDA extends our previously developed heuristic algorithm, knapsack-problem based decomposition algorithm (KDA), to obtain feasible solution for MSSPs under both endogenous and exogenous uncertainties with both continuous and discrete state variables and recourse actions that may also impact scenario specific constraint satisfaction. The GKDA is applied to solve four different planning problems under uncertainty: new technology investment planning,22 clinical trial planning,3 process network synthesis,1,2 and artificial lift infrastructure planning.7 For all problems considered, GKDA efficiently generated feasible solutions and primal bounds. The relative gap of the GKDA solutions from the optimum solutions ranged from 0.0% to 20.0% for all problem instances solved in this paper, and the GKDA generated these solutions with up to six orders of magnitude faster than solving the deterministic equivalents of the MSSPs. The GKDA is a greedy algorithm, and hence, fails to capture the potential impact of current time period decisions on future time periods. Future work will focus on investigating various

Zuo Zeng graduated from Southeast University with a BSs in Chemical Engineering and a minor in Business Administration, and earned his Masterʼs Degree at Carnegie Mellon University. He is currently a Doctoral Candidate at Auburn University, where he focuses on modeling of and developing solution strategies for large-scale multistage stochastic programs with endogenous or/and exogenous uncertainties.

Brianna Christian is a data scientist at Goodyear. She holds a BSc from the University of Kansas in Chemical Engineering and recently M


Article

Industrial & Engineering Chemistry Research completed her Ph.D. in Chemical Engineering at Auburn University. Her research focuses on optimal decision-making under decision dependent uncertainty. Of interest is addressing the computational challenges that arise in solving large-scale multi-stage stochastic programs.

VS = velocity strings GL = gas lift ESP = electrical submersible pump SRP = sucker rod pump PCP = progressing cavity pump JP = jet pump NTIP = new technology investment planning Mscf/D = thousand standard cubic feet per day BPD = barrels per day Mt = megaton

■

REFERENCES

(1) Goel, V.; Grossmann, I. E. A class of stochastic programs with decision dependent uncertainty. Math. Programming. 2006, 108, 355− 394. (2) Tarhan, B.; Grossmann, I. E. A multistage stochastic programming approach with strategies for uncertainty reduction in the synthesis of process networks with uncertain yields. Comput. Chem. Eng. 2008, 32, 766−788. (3) Colvin, M.; Maravelias, C. T. A stochastic programming approach for clinical trial planning in new drug development. Comput. Chem. Eng. 2008, 32, 2626−2642. (4) Colvin, M.; Maravelias, C. T. Scheduling of testing tasks and resource planning in new product development using stochastic programming. Comput. Chem. Eng. 2009, 33, 964−976. (5) Colvin, M.; Maravelias, C. T. Modeling methods and a branch and cut algorithm for pharmaceutical clinical trial planning using stochastic programming. Eur. J. Oper Res. 2010, 203, 205−215. (6) Gupta, V.; Grossmann, I. E. Solution strategies for multistage stochastic programming with endogenous uncertainties. Comput. Chem. Eng. 2011, 35, 2235−2247. (7) Zeng, Z.; Cremaschi, S. Artificial lift infrastructure planning for shale gas producing horizontal wells. Proc. FOCAPO/CPC 2017, 8−12. (8) Apap, R. M.; Grossmann, I. E. Models and computational strategies for multistage stochastic programming under endogenous and exogenous uncertainties. Comput. Chem. Eng. 2017, 103, 233−274. (9) Birge, J. R.; Louveaux, F. Introduction to Stochastic Programming; Springer Science & Business Media, 2011. (10) Jonsbraten, T. W.; Wets, R.J. B.; Woodruff, D. L. A class of stochastic programs with decision dependent random elements. Ann. Oper Res. 1998, 82, 83−106. (11) Medal, H. R.; Pohl, E. A.; Rossetti, M. D. Allocating Protection Resources to Facilities When the Effect of Protection is Uncertain. IIE Trans. 2016, 48, 220−234. (12) Goel, V.; Grossmann, I. E. A stochastic programming approach to planning of offshore gas field developments under uncertainty in reserves. Comput. Chem. Eng. 2004, 28, 1409−1429. (13) Goel, V.; Grossmann, I. E.; El-Bakry, A. S.; Mulkay, E. L. A novel branch and bound algorithm for optimal development of gas fields under uncertainty in reserves. Comput. Chem. Eng. 2006, 30, 1076− 1092. (14) Tarhan, B.; Grossmann, I. E.; Goel, V. Stochastic programming approach for the planning of offshore oil or gas field infrastructure under decision-dependent uncertainty. Ind. Eng. Chem. Res. 2009, 48, 3078−3097. (15) Mercier, L.; Van Hentenryck, P. Amsaa: A Multistep Anticipatory Algorithm for Online Stochastic Combinatorial Optimization. Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems. CPAIOR; Lecture Notes in Computer Science, Vol. 5015; Springer: Berlin, Heidelberg, 2008. (16) Solak, S.; Clarke, J. P. B.; Johnson, E. L.; Barnes, E. R. Optimization of R&D project portfolios under endogenous uncertainty. Eur. J. Oper Res. 2010, 207, 420−433. (17) Tarhan, B.; Grossmann, I. E.; Goel, V. Computational strategies for non-convex multistage MINLP models with decision-dependent uncertainty and gradual uncertainty resolution. Ann. Oper Res. 2013, 203, 141−166.

Selen Cremaschi is B. Redd Associate Professor in the Department of Chemical Engineering at Auburn University. She earned her Ph.D. degree from Purdue University, and both her BSc and MSc degrees from Bogazici University (Turkey), all in Chemical Engineering. Her research interests are risk management, optimization, process synthesis and planning under uncertainty. Her research group works at the intersection of operations research and chemical engineering, and develops systems analysis and decision support tools for complex systems, mainly focusing on pharmaceutical and energy industry. Prior to joining Auburn University, she was a faculty member of the Russell School of Chemical Engineering at The University of Tulsa. She was the recipient of a Tulsa Tau Beta Pi Teaching Excellence award (2010), an NSF CAREER award (2011), and a Zelimir Schmidt Award for Outstanding Researcher (2013) among others.

■

ACKNOWLEDGMENTS This invited contribution is part of the I&EC Research special issue for the 2018 Class of Influential Researchers. The authors gratefully acknowledge funding from the National Science Foundation (NSF) through the Auburn University Integrative Graduate Research and Education Traineeship (IGERT) program (Award No. 1069004).

■

ABBREVIATIONS MSSP = multistage stochastic programming GKDA = generalized knapsack-problem based decomposition algorithm NACs = nonanticipativity constraints SAA = sample average approximation SSD = sequential scenario decomposition KDA = knapsack-problem based decomposition algorithm PH = progressive hedging MSSPEX = a general MSSP formulation with endogenous and exogenous uncertainties MILP = mixed-integer linear programing MINLP = mixed-integer nonlinear programing KSP = knapsack subproblem ENPV = expected net present value ALMs = artificial lift methods ALIP = artificial lift infrastructure planning problem PL = plunger lift FL = foam lift WHC = well head compression N


Article

Industrial & Engineering Chemistry Research (18) Gupta, V.; Grossmann, I. E. A new decomposition algorithm for multistage stochastic programs with endogenous uncertainties. Comput. Chem. Eng. 2014, 62, 62−79. (19) Christian, B.; Cremaschi, S. Heuristic solution approaches to the pharmaceutical R&D pipeline management problem. Comput. Chem. Eng. 2015, 74, 34−47. (20) Christian, B.; Cremaschi, S. Variants to a knapsack decomposition heuristic for solving R&D pipeline management problems. Comput. Chem. Eng. 2017, 96, 18−32. (21) Christian, B.; Cremaschi, S. A branch and bound algorithm to solve large-scale multistage stochastic programs with endogenous uncertainty. AIChE J. 2018, 64, 1262−1271. (22) Christian, B.; Cremaschi, S. A Multistage Stochastic Programming Formulation to Evaluate Feedstock/Process Development for the Chemical Process Industry. Chem. Eng. Sci. 2018, 187, 223. (23) Lea, J. F.; Nickens, H. V. Gas Well Deliquification; Gulf Professional Publishing, 2011. (24) Fahmi, I.; Nuchitprasittichai, A.; Cremaschi, S. A new representation for modeling biomass to commodity chemicals development for chemical process industry. Comput. Chem. Eng. 2014, 61, 77−89. (25) Misener, R.; Thompson, J. P.; Floudas, C. A. APOGEE: Global optimization of standard, generalized, and extended pooling problems via linear and logarithmic partitioning schemes. Comput. Chem. Eng. 2011, 35, 876−892. (26) Oral, M.; Kettani, O. A linearization procedure for quadratic and cubic mixed-integer problems. Oper. Res. 1992, 40 (1-supplement-1), S109−S116. (27) Valbuena, J. D.; Defining the artificial lift system selection guidelines in horizontal wells. M.S. Thesis, The University of Tulsa, 2015.

O


A Generalized Knapsack-Problem Based Decomposition Heuristic for

Recommend Documents