Computational strategies for non-convex multistage MINLP models with decision-dependent uncertainty and gradual uncertainty resolution

In many planning problems under uncertainty the uncertainties are decision-dependent and resolve gradually depending on the decisions made. In this paper, we address a generic non-convex MINLP model for such planning problems where the uncertain parameters are assumed to follow discrete distributions and the decisions are made on a discrete time horizon. In order to account for the decision-dependent uncertainties and gradual uncertainty resolution, we propose a multistage stochastic programming model in which the non-anticipativity constraints in the model are not prespecified but change as a function of the decisions made. Furthermore, planning problems consist of several scenario subproblems where each subproblem is modeled as a nonconvex mixed-integer nonlinear program. We propose a solution strategy that combines global optimization and outer-approximation in order to optimize the planning decisions. We apply this generic problem structure and the proposed solution algorithm to several planning problems to illustrate the efficiency of the proposed method with respect to the method that uses only global optimization.


Introduction
Stochastic programming is a mathematical programming framework for modeling optimization problems that involve uncertainty in the data which are represented by probability distributions.There are two broad classes of stochastic programming problems: chance constrained (Charnes and Cooper 1963) and problems with recourse (Dantzig 1955;Beale 1955).In this work we focus on "multistage stochastic programming with recourse" and for the sake of simplicity, we will refer to it throughout the paper as "stochastic program".The fundamental idea behind stochastic programming is the concept of recourse, which is the ability of the decision-maker to take corrective action after a random event has taken place over a sequence of stages.
In this paper the uncertain parameters are assumed to follow discrete probability distributions and the planning horizon consists of a fixed number of decision points.Using these two assumptions, the stochastic process can be represented using scenario trees.Figure 1 is a standard representation of a scenario tree having two uncertain parameters (ξ 1 , ξ 2 ), each having two discrete possible values (π), which leads to four scenarios.Uncertain parameters ξ 1 and ξ 2 reveal at the end of first and second time periods respectively.
In a standard scenario tree, each node for time period t represents a possible state of the system at that time period.Each arc represents the possible transition from one state in time period t to another state in time period t + 1.A path from the root node to a leaf node represents a scenario.Thus, a scenario is a combination of possible uncertain parameters in each of the time periods, which in turn define a stage.
Figure 2 is an alternative representation of the scenario tree in Fig. 1 proposed by Ruszczyński (1997).In this representation each scenario is denoted by a set of unique nodes.The horizontal lines connecting nodes in time period t indicate that these nodes have the same information in that time period, and are therefore indistinguishable.These non-anticipativity constraints suggest that decisions cannot be based on knowledge that will be revealed in the future.Collapsing the horizontal lines reduce the tree in Fig. 2 to the one shown in Fig. 1.For modeling the problem, scenario trees presented in Fig. 2 will be considered in order to explicitly handle the non-anticipativity constraints, and incorporate decision dependent uncertainties and gradual uncertainty resolutions (Tarhan and Grossmann 2008).(SP) is a standard linear stochastic program with |T | time periods and |S| scenarios, modeled using the alternative scenario tree proposed by Ruszczyński (1997).The parameters p s represent the probability of scenario s, while the variables x s t represent variables for time period t in scenario s.Equation (1) corresponds to the objective of maximizing the expectation of an economic criterion (e.g.net present value).Equation (2) represents the multi-period constraints for each scenario s.Equation (3) represents the non-anticipativity constraints which state that decisions x s t and x s t must be identical if scenario pair (s, s ) is indistinguishable in time period t .These non-anticipativity constraints (denoted by set N) correspond to the horizontal lines connecting nodes in different scenarios in Fig. 2. Equation (4) represents the bounds and integrality constraints on the variables.A s τ x s τ ≤ a s t , ∀(t, s)

(SP) max
(2) x s t = x s t , ∀(t, s, s ) ∈ N (3) One important aspect of this paper is decision-dependent uncertainty.According to Jonsbraten (1998), uncertainty in stochastic programming problems can be divided into two classes: exogenous uncertainty and endogenous uncertainty.Most previous work in the literature deals with problems with exogenous uncertainty (e.g.demands) where the optimization decisions cannot influence the stochastic process.Reviews of previous work on problems with exogenous uncertainty can be found in Sahinidis (2004) and Schultz (2003).Problems where stochastic processes are affected by decisions are said to possess endogenous uncertainty.According to Goel and Grossmann (2006), decisions affect the stochastic process in at least two ways: decisions can alter the probability distributions (type 1), or decisions can act to discover more accurate information (type 2).In this paper, we focus on the latter type of endogenous uncertainty (type 2) where the decisions act to gain more information and eventually resolve uncertainty.
Previous works on type 1 decision-dependent uncertainty include Ahmed (2000), Viswanath et al. (2004), Held and Woodruff (2005).Ahmed (2000) presented examples on network design, server selection, and facility location problems with decision-dependent uncertainties, showing that these programs can be reformulated as MILP problems and solved by LP-based branch-and-bound algorithms.Viswanath et al. (2004) addressed a network problem having endogenous uncertainty in survival distributions.The problem was a two-stage stochastic program in which first-period investment decisions were made for changing the survival probability distribution of arcs after a disaster.The aim was to find the investments that would minimize the expected shortest path from source to destination after a disaster.Held and Woodruff (2005) worked on another instance in which the problem included endogenous uncertainty in the structure of a network.In each stage of this problem, an operator tries to find the shortest path from a source to a destination after some of the nodes in the network becomes blocked.The aim was to maximize the probability of stopping the flow of goods or information in the network.
Previous work on type 2 uncertainty can be classified into two subgroups: the first subgroup assumes that the discovery of accurate information (resolution of endogenous uncertainty) occurs instantaneously (Pflug 1990;Jonsbraten et al. 1998;Goel 2005;Goel and Grossmann 2004;Goel et al. 2006;Boland et al. 2008) and the second subgroup assumes uncertainty resolves gradually over time because of learning (Tarhan and Grossmann 2008;Solak 2007;Tarhan et al. 2009;Stensland and Tjøstheim 1991;Dias 2002;Jonsbraten 1998;Harrison 2007).Pflug (1990) addressed the problems in the context of discrete-event dynamic systems in which the underlying stochastic process depended on the optimization decisions.Jonsbraten et al. (1998) addressed endogenous uncertainty where project decisions reduced the uncertainty, and proposed an implicit enumeration algorithm in which decisions that affect the uncertain parameter values are made in the first stage.Goel and Grossmann (2004) considered the gas field development problem under uncertainty in the size and quality of reserves where decisions on the timing of field drilling were assumed to yield an immediate resolution of the uncertainty.Goel et al. (2006) later proposed a branch-and-bound algorithm for solving the corresponding disjunctive/mixed-integer programming model where the lower bounds are generated by Lagrangean duality.Boland et al. (2008) applied multistage stochastic programming to open pit mine production scheduling, which is modeled as a mixed-integer linear program.They considered endogenous uncertainty where the excavation decisions resolve uncertainty in geology immediately.They followed a similar approach as Goel and Grossmann (2006) for modeling the problem, with the exception of eliminating some of the binary variables used in the general formulation to represent conditional nonanticipativity constraints.Furthermore, they solved the model in full-space without using a decomposition algorithm.Tarhan and Grossmann (2008) considered the synthesis of process networks.They assumed that there were uncertainties in the yields of the processes and that these uncertainties resolved gradually over time depending on the investment and operating decisions.Compared to the previous work that involved with instantaneous uncertainty resolution, gradual resolution brought some challenges to the model, underlying scenario tree, and the nonanticipativity constraints.Solak (2007) considered the stochastic programming models in which times of uncertainty realizations are dependent on the decisions made and uncertainty resolves gradually.They used a similar modeling approach as the one presented in Tarhan and Grossmann (2008).Specifically, they considered the project portfolio optimization problem that deals with the selection of research and development projects and determination of optimal resource allocations.They used sample average approximation method (Kleywegt et al. 2002) for solving the problem, where the sample problems are solved through Lagrangian relaxation and heuristics.Tarhan et al. (2009) addressed the planning of offshore oil or gas field infrastructure under uncertainty.The main uncertainties considered are in the initial maximum oil or gas flowrate, recoverable oil or gas volume, and water breakthrough time of the reservoir, which are represented by discrete distributions.Furthermore, it is assumed that these uncertainties are not immediately realized, but are gradually revealed as a function of design and operating decisions.To account for these decision-dependent uncertainties, the authors proposed a multistage stochastic programming model that captures the complex economic objectives and nonlinear reservoir behavior, and simultaneously optimizes the investment and operating decisions over the entire planning horizon.Stensland and Tjøstheim (1991) have worked on a discrete time problem for finding optimal decisions with uncertainty reduction over time and applied their approach to oil production.These authors expressed the uncertainty in terms of a number of production scenarios.Their main contribution was combining production scenarios and uncertainty reduction effectively for making optimal decisions.Dias (2002) presented four propositions to characterize technical uncertainty and the concept of revelation towards the true value of the variable.These four propositions, based on the theory of conditional expectations, are employed to model technical uncertainty.Jonsbraten (1998) considered gradual uncertainty reduction where all Fig. 3 Scenario tree with different uncertainty resolution scheme uncertainty is assumed to resolve at the end of the project horizon.The author used a decision tree approach for modeling the problem where Bayesian statistics are applied to find the probabilities of branches in the decision tree that are decision dependent.The author also proposed an algorithm that relies on the prediction of upper and lower bounds.Harrison (2007) used a different approach for optimizing two-stage decision making problems under uncertainty.Some of the uncertainty was assumed to resolve after the observation of the outcome of the first stage decision.The author developed a new method, called Bayesian Programming, where the corresponding integrals were approximated using Markov Chain Monte Carlo simulations, and decisions were optimized using simulated annealing type of meta-heuristics.
It is more convenient to represent the endogenous type of uncertainty using the alternative scenario tree in Fig. 2 (compared to Fig. 1) where the non-anticipativity constraints are explicitly handled as functions of decisions.To illustrate the effect of decisions on the scenario tree in Fig. 2, we assume that all the decisions are delayed for one period and no decision that can resolve uncertainty is made in time period 1.This forces all scenario pairs stay as indistinguishable until time period 2 and some become distinguishable at time period 3.The resulting scenario tree will be as shown in Fig. 3.This tree is different from the one in Fig. 2 where scenarios 2 and 3 at time period 2, and all scenarios at time period 3 are distinguishable.Therefore, in case of endogenous uncertainty, the set of non-anticipativity constraints (N ) in (3) is not fixed but it is a function of the decisions (x s t ).In this work, we provide a generic stochastic MINLP programming model that has decision-dependent uncertainty of type 2. In type 2, the decisions act to discover more accurate information and do not alter the underlying probability distribution.The model proposed in Tarhan et al. (2009) was designed specifically for the problem of oil and gas fields, and did not provide any generic model structure or discussion of the broader class of problems that it belongs to.In this work, we generalize the part of the model that relates the decisions to set of non-anticipativity constraints (N ), which changes as a function of decisions.Furthermore, we improve the computational efficiency of the solution strategy proposed by Tarhan et al. (2009) by incorporating both global and local optimization into the proposed optimization algorithm.
The outline of this paper is as follows.In Sect. 2 we present the generic mathematical model for the class of problems under consideration.In Sect. 3 the proposed solution approach is explained.Section 4 compares the results found by the proposed solution approach, and discusses the results on synthesis of process networks and planning of offshore oil or gas field infrastructure problems.

Mathematical model
In this section, only generic variable names are used in the explanation of the model.The generic variable names are selected based on whether a variable is decision, state, or recourse variable.The decision, state, and recourse variables are denoted by vectors d, x and u with dimension n d , n x , and n u , respectively.The first n d ≤ n d , n x ≤ n x , n u ≤ n u elements of these vectors are restricted to be discrete variables.Also, w and z are logic/binary variables used for defining the conditional non-anticipativity constraints.
The planning horizon is divided into |T | planning periods.The sequence of decisions is as follows: the decision variables d t are implemented at the beginning of time period t , which is followed by the resolution of uncertainty.The state variables x t are automatically calculated when the decisions are fixed.At the end of time period t , some recourse decisions u t are implemented (see Fig. 4).
Based on the above definitions, the proposed multistage stochastic MINLP model is as follows: In ( 5), the objective is to maximize the expected net present value, which is a linear or non-linear function of decision (d), state (x), and recourse (u) variables.Equation ( 6) and ( 7) are generic multi-period linear and non-linear constraints, respectively, that relate the decision (d), state (x), and recourse variables (u) for every scenario.There are two types of non-anticipativity constraints: initial (N I ) and conditional (N C ). Equation ( 8) represents the initial non-anticipativity constraints that hold regardless of any decision.These initial nonanticipativity constraints include not only the first period non-anticipativity constraints, but also some other subsequent period non-anticipativity constraints that must hold because of gradual uncertainty resolution.For instance, if at least two periods are required to distinguish two scenarios, then we can include non-anticipativity constraints for such scenario pairs for the second period in the initial non-anticipativity constraints.Equations ( 5)-( 8) and ( 12) are similar to the conventional stochastic programs (( 1)-( 4)), which do not include decision dependent uncertainty.To illustrate these constraints, we provide specific examples from the synthesis of process networks problem (details in Sect.4.1 and Appendix).In the process networks example, (6) corresponds to the linear constraints such as the material balance constraints (A.7), (A.8), (A.9), (A.10) on every node and process.Equation ( 7) corresponds to the non-linear constraints such (A.2) (non-linearity due to term (d qe,s i,t ) α ) that are used for calculating net present value for every scenario.
Equations ( 9)-( 11) are the constraints required to account for decision-dependent uncertainty and gradual uncertainty resolution.Equation ( 9) relates the binary/logic variables (w) with the discrete and continuous decision (d), state (x), and recourse (u) variables.These binary/logic variables (w) are used in (10) to model problem specific uncertainty resolution rules for every level l (∀l ∈ L) of gradual resolution and determine if scenario pair (s, s ) is indistinguishable or not (for levels l ∈ L, see Fig. 5 ).Since this notation is not standard, below we provide a simple example in the context of synthesis of process networks (details presented in Sect.4.1 and Appendix).In the process networks example (9) corresponds to constraints such as (A.4) and (A.5) which state that logic variable w l,s i,t where level l = 1 will be true if and only if process i has not been expanded, or not run as a pilot plant until time period t − τ .Level l = 1 stands for the resolution of uncertainty in Fig. 5. Similarly, based on our assumptions about resolution of uncertainty, w l,s i,t where level l = 2 will be true if and only if process i has operated only one time period or run as a pilot plant until time period t − τ .Level l = 2 specifies the resolution of uncertainty at the second level, i.e. partial resolution in Fig. 5. Equation (10) corresponds to (A.25) where the notation Q stands for the set of uncertain parameters and M q,l stands for the set of scenario pairs that differ by one uncertain parameter q in level l of uncertainty resolution.The (A.25) states that a scenario pair that is an element of M q,l , ∀q ∈ Q, ∀l ∈ L will be indistinguishable if and only if for each process that distinguishes the scenario pair, w l,s i,t holds true.For specific subsets M q,l used for the example problem, please see Appendix, Tables 5, 6, and 7.
The conditional non-anticipativity constraints in the disjunction (( 11)) are included into the model if two scenarios are indistinguishable at the end of time period t (z s,s t is true).The conditional non-anticipativity constraints dictate that the decisions at the beginning of Fig. 6 Simple process converting raw material to final product time period t + 1 in scenarios s and s must be identical.Otherwise they are not imposed on the feasible space.As proved in Goel and Grossmann (2006), it is enough to consider subsets of scenario pairs (M q,l ) that differ in only one uncertain parameter.Therefore, ( 10) and ( 11) are defined for scenario pairs (s, s ) that are an element of exactly one of such sets M q,l , ∀q ∈ Q, ∀l ∈ L. Equations ( 12)-( 13) represent the variables properties and integrality requirements.
As explained in the previous paragraph, gradual uncertainty resolution brings more complication to the conventional stochastic programs by ( 9)-( 11).One major difference between gradual resolution and instantaneous resolution cases (e.g.Goel and Grossmann 2006), which were discussed in the Sect. 1, is in the conditional non-anticipativity constraints.As mentioned in Sect. 1, we focus on the type 2 decision-dependent uncertainty where decisions act to gain more information and resolve uncertainty.In gradual resolution case, since the observed data may not be the true value, the conditional non-anticipativity constraints are imposed only on the decision variables (d t ), but not by the recourse variables (u t ).The decision maker observes some output from the system, but he/she may not differentiate certain scenario pairs due to the partial resolution.
For instance, let us assume there is a process which takes a raw material and generates a final product (Fig. 6).The demand for final product has to be satisfied exactly and in case of shortage, the final product can be bought from the market.There is uncertainty in the yield of the process, which resolves gradually depending on the production decisions.The decision at the beginning of each time period is the amount of raw material to buy and the recourse at the end of each period is the amount of final product to buy from the market to satisfy the shortage.
Although it may be more reasonable to assume a continuous type of distribution for the yield, we restrict it to a discrete one with four values (θ 1 , θ 2 , θ 3 , θ 4 ) to more easily explain the concept.The process runs a full year and yield data is calculated regularly using output/input data of the process.We assume that it requires two years of production to gather enough data to find the true value of the yield, i.e. resolve the uncertainty.We translate this into discrete scenarios in the following way: after one year of production (l = 2) the decision maker observes that the yield is on the low side (θ 1 or θ 2 ) or on the high side (θ 3 or θ 4 ) (as shown in Fig. 5).At this level l = 2 the decision maker may observe any of the discrete yields but if he/she finds out the yield is θ 1 , he/she can conclude for sure that the actual yield is not θ 3 or θ 4 but cannot conclude if it is θ 1 or θ 2 .At this partial resolution stage, the decision maker has to make identical decisions for the amount of raw material to purchase for the second year for each indistinguishable scenario groups (i.e.scenarios (1-4) and (5-8)) as shown in the scenario tree in Fig. 5.Note that, the recourse decisions (amount of final product to buy from the market) at the end of first year have to be different among some of the indistinguishable scenarios (e.g.scenario 1 and 3).This is because the demand stays the same for two scenarios, but the observed production changes due to different yields (θ 1 and θ 2 ), leading to different recourse decisions.Therefore, in contrast to the immediate resolution case (Goel and Grossmann 2006) we impose only the decision variables (d t ), not the recourse variables (u t ), in the disjunction in (11) for conditional non-anticipativity constraints.

Solution strategy
The size of the multistage stochastic MINLP model (P ) increases quadratically with the number of scenarios and linearly with the number of time periods.Therefore, it is difficult to solve this nonconvex MINLP model in fullspace for real size problems given that these problems are known to be NP hard.The model (P ) is composed of nonconvex mixed-integer nonlinear subproblems.The subproblems are connected through initial (8) and conditional (11) non-anticipativity constraints.When the initial and conditional non-anticipativity constraints are relaxed, each subproblem can be solved independently.In order to take advantage of this special problem structure, we propose a duality-based branch and bound algorithm along the lines of Goel and Grossmann (2006).For the sake of completeness, we will present the upper and lower bounding procedures, branching scheme, the solution algorithm (SP-GO) developed in previous works of Tarhan et al. (2009), and the improved algorithm (SP-OA) proposed in this paper.

Upper bounding procedure
In the proposed algorithm, the upper bound at the root node of the branch and bound tree is found by optimizing model (P LR 0 ) in which the logic constraints ( 9)-( 10) and disjunction (11) have been removed and the initial non-anticipativity constraint (8) are dualized as follows: 6)-( 7), (12) (15) The parameters λ s,s d,t represent the Lagrange multipliers corresponding to constraint (8).In order to find the tightest upper bound generated by model (P LR 0 ) at the root node, we consider the Lagrangean dual problem (P LD 0 ) that minimizes the model (P LR 0 ) in the space of multipliers.
At any node n, other than the root node in the branch and bound tree, model (P n ) is optimized to calculate upper bounds.In model (P n ), not only the initial non-anticipativity constraints, but also some conditional non-anticipativity constraints are included.Initial nonanticipativity constraints are added to the model regardless of any decision, whereas conditional ones from the relaxed disjunction are included based on the branching cuts.The conditional non-anticipativity constraints that apply in node n are included in the dynamic set N n C .The selection of which conditional non-anticipativity constraint to include into set N n C , as well as some necessary cuts to be added to (P n ) will be discussed.The model (P n ), not including any such cuts, is given as follows: 6)-( 7), ( 12) The upper bound at node n is generated by optimizing model (P LR n ) in which the nonanticipativity constraint ( 18) is dualized as follows: Similarly as in the root node, in order to find the tightest upper bound, we again consider the Lagrangean dual problem (P LD n ) which is minimization of the model (P LR n ) in the space of multipliers, where λ s,s d,t represents the Lagrange multipliers corresponding to constraint (18).Both models (P LR 0 ) and (P LR n ) can be decomposed into independent nonconvex MINLP subproblems for fixed values of the multipliers.Both models are relaxations of the model (P ) for any fixed values of the Lagrange multipliers, and in order to obtain valid upper bounds, these nonconvex subproblems must be globally optimized (for proof see Tarhan and Grossmann 2008).
Minimization of the Lagrangean dual models (P LD 0 ) and (P LD n ) in the space of multipliers is performed by the subgradient method proposed by Fisher (1985).

Lower bounding procedure
At each node, a lower bound is found using a heuristic that converts the solution found by the upper bound procedure to a feasible solution.Usually the solution found by the upper bound generation does not satisfy the non-anticipativity constraints.Feasible solutions are generated using a rolling horizon approach (Dimitriadis et al. 1997).Decisions in the indistinguishable scenario pairs are found by calculating the probability weighted average of variables in such scenarios.After fixing these variables and considering the resolution of uncertainty depending on the fixed decisions, the next period decisions are found iteratively by calculating in a similar fashion the probability weighted average of the variables in indistinguishable scenarios.This procedure is terminated when all scenarios are distinguishable, or the end of the planning horizon is reached, whichever comes first.Then, having these decisions fixed, the MINLP model is solved in full space using an outer approximation algorithm (Duran and Grossmann 1986), which may yield a feasible solution.In both example problems, the feasible solution is guaranteed to be available because of complete recourse.In the process synthesis problem, there is always the flexibility to buy the final product from the market to satisfy the demand, and in the oilfield problem there is always a feasible operating decision given a feasible infrastructure decision.Note that a feasible solution is always available, but there is no guarantee that the outer-approximation algorithm is going to find it.During the search, due to the non-convex nature of the problem, the generated cuts can cutoff the feasible solution space and may result in infeasible solutions.In such cases, the bound is not updated and the search resumes.

Branching
In general, at any node in the branch and bound tree, the solution of the Lagrangean dual may not satisfy the dualized non-anticipativity constraint (8), or the conditional nonanticipativity constraints in the relaxed disjunction (11) inferred by the decisions.In this Fig.7 Branching on the disjunction in constraint (11) case, new branches are generated from the current node by considering the violations in the dualized non-anticipativity constraints or the relaxed disjunction.

Branching on the dualized non-anticipativity constraints
At any node of the branch-and-bound tree, the solution of model (P LD n ) may not satisfy the dualized initial or conditional non-anticipativity constraints (N I ∪ N n C ).In this case, branching is performed over the violated constraints.If the violated non-anticipativity constraint involves discrete variables, then the branching strategy divides the feasible region into two, where one of the regions is restricted to A special case occurs during the branching of the first period non-anticipativity constraints.As these constraints should hold for each scenario pair, we can rewrite these equations as one expression, d 1 .The branching can then be performed so that the entire expression is restricted to ≤ dt and to ≥ dt (in case of continuous variables the inequalities take the form ≤ dt − ε and ≥ dt + ε.

Branching on the disjunctions
When the dualized non-anticipativity constraints are satisfied, it is possible to continue searching by separating the feasible space into two using the disjunction (11).The branching is performed using the conditional non-anticipativity constraints in the relaxed disjunctions that are supposed to hold, but do not do so given the values of the variables at the previous time periods.For instance, assume the decisions satisfy the initial non-anticipativity constraints, then the indistinguishable scenario pairs are inferred by using ( 9)-( 10).Given the indistinguishable scenario pairs, the branching is performed on z s,s t as shown in Fig. 7 if the variable values do not satisfy the conditional non-anticipativity constraints in disjunction (11).
The two generated nodes in Fig. 7 (nodes 1 and 2) have different properties.At node 1, z s,s t is fixed to false/zero and at node 2, z s,s t is fixed to true/one.At node 1 due to fixing z s,s t to false, it is required to add some cuts, which will guarantee the distinguishability of the scenario pair (s, s ).These cuts will be generated using constraints (9)-(10).To illustrate these cuts, we will again use the synthesis of process networks problem (Appendix).As an example, if we fix z s,s t to false, where the scenario pair (s, s ) belongs to the set M 1 , then using Table 6 and (A.25) we find that logic constraints ( 21) and ( 22) must be included in scenario subproblems s and s , respectively.
Different from node 1, at node 2 some non-anticipativity constraints coming from the relaxed disjunction are added to the set N n C in model (P n ) due to the fixing of z s,s t to true.Similar to node 1, cuts are generated and added to model (P n ) where z s,s t is fixed to true.These cuts are necessary to make the solution algorithm convergent (for proof see Goel and Grossmann 2006).

Solution algorithm (SP-GO)
Based on the explanations about the upper, lower bounding, and branching procedures, the proposed algorithm for optimizing multistage stochastic MINLP model (P ) is presented below.P denotes the list of open nodes each having an upper bound φ UB n found by the Lagrangean dual problem at node n, while φ LB represents the objective value of the best feasible solution obtained so far.Given these variables, the major steps of the algorithm are given below.
Step 2 Termination: If P = ∅, stop.The current best solution is optimal.Otherwise, repeat steps 2 to 6.
Step 3 Node selection: Select and delete node n from P based on the best bound.
Step 4 Bound generation: For the root node, generate the upper bound (φ UB 0 ) by applying a global optimization algorithm and obtain a lower bound (φ LB ) solutions ( d, x, û) and ( d, x, ū), respectively.(Details of this step (steps a-g) will be explained below.) Step 5 Fathoming: Delete from P all problems P with φ(P ) ≤ φ LB .
Step 6 Branching: Branch on the dualized non-anticipativity constraints or disjunctions that are violated by the solution ( d, x, û) of the relaxed problem (P n ).Generate two children nodes, add them to P.
Note that step 1 is the initialization and the algorithm iterates between steps 2 and 6 until convergence is achieved.The details of the step 4 (steps a to g) as follows: Step a Set iteration i = 0.
Step b While i is less than or equal to max_iteration, i = i + 1, repeat steps b through h.
Step c Generate upper bound: For fixed multipliers, use global optimizer for each MINLP subproblem to solve (P LR n ) to obtain solution ( d, x, û) with objective function value φ.
Step d Update upper bound: Update the upper bound by φ UB n = min{φ UB n , φ}.
Step e Generate lower bound: Generate a feasible solution ( d, x, ū) with objective value φ for the model (P ) based on the solution ( d, x, û) generated at step c.
Step f Update lower bound: Update the lower bound by φ LB = max{φ LB , φ}.
Step g Update multipliers: Update multipliers using the subgradient method.

Solution algorithm (SP-OA)
The bottleneck of the algorithm (SP-GO) is the upper bounding procedure (step c) in step 4, where every scenario subproblem has to be globally optimized in order to have valid upper bounds during the subgradient iterations.This is computationally expensive.Since it is sufficient to generate a valid upper bound only at the final subgradient iteration, the proposed algorithm, SP-GO, can be modified such that the non-convex subproblems are solved Fig. 8 Graphical comparison of SP-GO and SP-OA iterations using an outer-approximation algorithm instead of global optimization at the intermediate subgradient iterations.The reason for such a modification is that, aside from the fact that the outer-approximation algorithm is much faster than a global solver, in practice during the subgradient optimization we do not need to find the optimal multipliers, but simply multipliers that update the upper bound.This modification intends to improve the Lagrange multipliers at each iteration, without solving each scenario problem globally.The drawback of the approach is that it does not determine valid bounds during these intermediate subgradient iterations.
Note that the algorithm (SP-OA) combines the global optimization and outer-approximation algorithm for expediting the solution process without violating the validity of the bounds.Furthermore, the outer-approximation algorithm can be replaced by any algorithm that assumes convexity (e.g.generalized Benders decomposition, extended cutting plane), that will find local solutions in short time.We modify the steps c and d in the proposed algorithm (SP-GO) to incorporate the outer-approximation algorithm into the algorithm SP-OA.The type of optimization (global or outer approximation) is based on the iteration i and maximum iteration limit (max_iteration).If iteration i is one or the maximum iteration limit, every subproblem in the Lagrangean dual of the problem (P LR n ) is solved using global optimization; otherwise an outer-approximation algorithm that relies on convexity assumption is used.The reason for using global optimization at the first iteration is for initializing the variable values close to the global optimal solution, which is taken as an input for the outer-approximation algorithm at the second iteration.This improves the chances of finding optimal or near optimal solutions during the outer-approximation algorithm.The reason for using global optimization at the last iteration is for generating valid upper bounds that will be used for pruning the parts of the branch and bound tree and calculating the duality gap.In step d, the upper bound is updated only after global optimization is used for calculating valid bounds.
Figure 8 compares the typical profiles of the Lagrangean dual as a function of the multipliers.The solid lines represent a profile generated when all the subproblems are solved by global optimization, whereas the dashed lines represent a profile generated using outer approximation at the intermediate iterations.In the first and last iteration, a solution is found on the profile generated by global optimization (solid lines), and during the intermediate iterations outer-approximation (dashed lines) may end up in one of the local optimal solutions (x OA i ).

Results
In this section, we present results for comparing the solutions obtained with the algorithms SP-GO and SP-OA for the modified version of synthesis of process networks by Tarhan and Grossmann (2008) and planning of offshore oil or gas field infrastructure problems by Tarhan et al. (2009).

Synthesis of process networks
This section explains briefly the synthesis of process networks example and Sect.4.1.1analyzes the solution time and quality of SP-GO and SP-OA.In the synthesis of process networks problem, we consider the selection and capacity expansion of processes over a planning horizon given that there is uncertainty in the yields of some processes in a network of processes.Uncertainty can be reduced through investment in pilot plants and production.
The trade-off is that pilot plants delay the introduction of processes but at a reduced uncertainty.The cost of capacity is assumed to be concave with respect to the equipment size since the cost of most process equipment obeys the economies of scale.These nonlinear cost functions are different than the linear cost functions used in Tarhan and Grossmann (2008).
Figure 9 shows the specific network used in the synthesis of process networks problem.The demand for the final product A over the planning horizon is known and the company must satisfy that demand.The current production takes place only in process 3, which consumes an intermediate product B from the market.In case of production shortage it is possible to buy final product A from the market at a higher cost to satisfy the demand.Inventory for both the intermediate and final product can be maintained.Two new technologies (process 1 and process 2) are available to produce the intermediate product B from two different raw materials C or D. These new technologies have uncertainty in their yields, which gradually resolve over a two year period either with investments in pilot plants or with plant operation.
In Fig. 9 process 3 is already operational with an existing capacity of 3000 tons/year and a known yield of 70%.The only difference between the two different technologies (Process 1 and 2) is the variance of the yield distributions.Although they possess the same mean value, 75%, process 2 has a higher variance than process 1 (see Table 5).Process 1 can realize after first year a yield of 69 or 81% and after the second year 69, 73, 77 and 81%.Similarly, yields of process 2 after first year are 60 or 90% and after second year 60, 70, 80 and 90%.It is assumed that the probability of each scenario is the same.The uncertainty resolution and decisions are assumed to be related in the following way.The uncertainty can be resolved if the decision maker builds the actual plant and operates for two time periods.The two steps of resolution correspond to the operation at each year.The decision maker can also build a pilot plant, operate one year and resolve the uncertainty partially.At that point, the actual plant has to be built for completely resolving the remaining uncertainty.The advantage of building a pilot plant is to hedge against greater uncertainty but the drawback is the delayed production start time.The detailed optimization model, scenarios, and data for this instance are presented in Appendix, Tables 8-12.

Fig. 9 Schematic representation of synthesis of process networks problem
The size of this instance in full space with 10 periods and 16 scenarios is given in Table 1.Although the problem seems to be trivial, combining 16 scenario subproblems over 10 years with initial and conditional non-anticipativity constraints, leads to a large scale problem.

Performance analysis of SP-GO and SP-OA
The results in this section have been obtained on a Pentium-IV, 3.20 GHz Windows machine.Also, we employed AIMMS 3.8.4 for implementing the solution algorithm using solvers CPLEX 11.0, CONOPT 3.14, SNOPT 6.1, BARON 7.5.3,AOA (AIMMS Outer Approximation Module).
The specific synthesis problem has been optimized using both algorithms SP-GO and SP-OA with 2% worst case gap.This gap is calculated by adding the specified gap of 1% for the global optimizer and specified gap of 1% for the branch and bound tree.The global optimization algorithms have been run with 1% optimality gap, and both algorithms were terminated when the gap in the branch and bound tree is less than 1%.The results are presented in Table 2.In this instance, although both algorithms find close solutions, SP-GO finds a slightly better (0.1%) feasible solution than SP-OA (6.636 vs. 6.629).However, while SP-GO requires 90.7 hours to complete the branch and bound search, SP-OA requires 34.8 hours, a reduction of nearly 60%.In SP-OA the optimum solution was first found in 32.1 hours while SP-GO found it in 71.1 hours.One can also compare the best upper bounds found by the two algorithms.SP-OA reduces the best upper bound more than SP-GO by searching more nodes in the branch and bound tree, but it cannot find a better feasible solution.

Planning of offshore oil or gas field infrastructure
In this section, we briefly describe the specific planning of offshore oil or gas field infrastructure problem.Section 4.2.1 compares the results and the solution times found by the two algorithms (SP-GO, SP-OA).
We consider a field consisting of a single reservoir (Fig. 10), where a number of wells can be drilled and exploited for oil in every reservoir during the planning horizon.The problem involves making investment and operating decisions over the planning horizon.Investment decisions are selection of the number, type and capacity of facilities, and installation schedule of these facilities, as well as selection of types of wells and drilling schedule of wells.
The operating decisions are the amount of oil production for each time period given the limitations of the reservoirs.The goal is to capture the complex economic tradeoffs that arise from the investment and operating decisions in order to maximize the expected net present value of the project.The details and data for the planning of oil or gas field infrastructure under uncertainty are reported in detail in Tarhan and Grossmann (2008).Briefly, the main uncertainties considered are in the initial maximum oil or gas flowrate, recoverable oil or gas volume, and water breakthrough time of the reservoir, which are represented by discrete distributions.Furthermore, it is assumed that these uncertainties are not immediately realized, but are gradually revealed as a function of well drilling and production decisions.The model optimizes the investment and operating decisions over the entire planning horizon.Two instances, 1 and 2, have been optimized to show the efficiency of the proposed algorithm.In both instances the reservoir behavior is nonlinear (because of the water flowrate equations) but in instance 1, the maximum oil flowrate is assumed to be a linear function of the cumulative oil production, whereas in instance 2 it is a nonlinear curve.In both cases the nonlinearities give rise to non-convexities in the model.The size of the two instances are the same and given in Table 3.

Performance analysis of SP-GO and SP-OA
The results of this section have been obtained on a Pentium-IV, 3.20 GHz Windows machine.Also, we employed AIMMS 3.8 for implementing the solution algorithm using solvers The stopping criteria for Table 4 has been set partly by the practical considerations and the experience from solving similar problem types.In both instances, the solution algorithm SP-GO has been terminated either after 120 hours (5 days) or 5% of worst case gap.Instance 1 has been terminated due to bounding, whereas instance 2 has been terminated due to time limitation.The termination criteria for SP-OA has been set by the business requirement of finding solutions overnight (8 hours).Therefore, the algorithm has been terminated either after 8 hours or 5% of worst case gap.Similar to SP-GO, instance 1 has been terminated due to bounding whereas instance 2 has been terminated due to time limitation.
Table 4 shows that the solution times of both instances using SP-GO are long for practical purposes (52 and 120 hours).As presented in Table 4, in instance 1 using algorithm SP-OA, which combines the global optimization and outer-approximation, reduces the solution time by 90% and obtains better feasible solution (6.54 vs. 6.38) compared to SP-GO.In instance 2, although both algorithms are terminated due to time restrictions, SP-OA found a better solution (4.84 vs. 4.65) compared to SP-GO in a shorter run time (8 vs. 120 hours).Similar to the results obtained in the synthesis of process networks problem, large reductions in solution time are achieved by combining local and global MINLP solvers.Moreover, in these problems, SP-OA also found better feasible solutions.

Conclusion
We have presented in this paper a generic non-convex multistage MINLP model with decision dependent uncertainties.We have proposed an improvement on the duality-based branch and bound algorithm for solving the large sized instances where global optimization and outer-approximation algorithms are combined.The performance of the new algorithm (SP-OA) has been compared with the previous approach (SP-GO) for two problems, the synthesis of process networks and planning of offshore oil or gas field infrastructure.In the synthesis of process networks problem, SP-OA obtained a solution that is within 0.1% of the one by SP-GO while reducing the solution time about 60%.In the oil field infrastructure problem, SP-OA improved not only the best feasible solution found by SP-GO by 2-5%, but also the solution time by nearly 90%.The improvement made by SP-OA can be accounted for by combining the global optimizer and the outer-approximation algorithm to reduce the solution time.SP-OA takes advantage of the strengths of both algorithms, leading to large reductions in solution time without necessarily sacrificing the quality of the solution.The results also show that a combination of local and global MINLP solvers can lead to better solutions faster than only using a global solver.The global optimizer BARON (Sahinidis 2004) was used to generate valid bounds in longer time, while the outer-approximation algorithm AOA was used to update Lagrange multipliers in a shorter time without guarantee of finding global optimal solutions.where α is a fractional exponent, 0 < α < 1, which gives rise to a concave cost function.Disjunction (A.3) represents the input-output relationships for the processes with uncertain yields, θ l,s i , at each period and scenario.Note that at each time period and scenario, the yield of a process must be in one of the possible resolution levels l ∈ L.  .4) and (A.5) relate the Boolean variable w l,s i,t to the previous time period decision variables so that the some set of constraints in (A.3) will be valid.Assuming that the uncertainty will resolve in two steps, there are three possibilities.If we build neither a pilot plant nor an actual plant, the set of constraints in (A.3) that has zero yield (since there is no production capacity) will be valid as expressed by (A.4).If we build/operate the actual plant one year or run pilot plant one year, the set of constraints in (A.3), that have still uncertain yields (at level 2), will be valid.This is reflected by (A.5).Finally, if the plant operates two or more years or at least one year after pilot plant operation, the set of constraints in (A.3) that have exact yield will be valid (level 3).This is captured by (A.6) since a process has to be in one of the levels of uncertainty at any point in time.Note that instead of using (A.6) we could define the rule for l = 3 but we simplify the constraints by using (A.6).Instead of specifying a rule for level l = 3 as in (A.4)-(A.5),we only specify that the process yield must in one of the levels l ∈ L.   The mass balance constraints at each node in the network are given by, k∈ON(j ) x rate,s k,t = k∈IN(j ) x rate,s k,t , ∀j ∈ J, ∀s ∈ S, ∀t ∈ T (A.8) The balance constraint that relates inventory, purchase, sales and production for final product at consecutive periods in the network is given by the following equation, Similarly, the balance constraint that relates inventory, purchase, production and consumption of intermediate product at consecutive periods in the network is given by the following equation,  .18)states that at most one pilot plant can be installed for each process having uncertain yield, for each scenario through the entire project life.The logic constraint (A.19) states that if there has been an expansion in the process i until period t in scenario s, then there is no need for a pilot plant for that process after period t .

Fig. 4
Fig. 4 Representation of planning horizon

Fig. 5
Fig. 5 Scenario tree representation of gradual uncertainty resolution in yield ∀i ∈ IU, ∀s ∈ S(A.18) t max dt , d s t ≤ dt and the other d s t ≥ dt , d s t ≥ dt .dt is calculated by the probability weighted average of the variables, dt = t ≤ dt − ε and d s t ≥ dt + ε, d s t ≥ dt + ε where dt is calculated as shown above.

Table 3
Model size for the oil

Mathematical model for synthesis of process networks Nomenclature
AcknowledgementsThe authors acknowledge the financial support from Exxon-Mobil Upstream Research Company and partial support of the National Science Foundation under grant CTS-0521769.Whether or not process i is operated in period t , scenario s d Whether or not the yield of process i is in level l in period t , scenario s Amount of final product to put into inventory at the end of period t , scenario s u inv−int,s Amount of purchases of final product in period t , scenario s x Flowrate of stream k ∈ SK in period t , scenario s Time delays In this instance, for simplicity, we assume that uncertainty resolves in two steps (i.e. in three levels, L = {1, 2, 3}).Based on the above definitions, the model for the synthesis of process networks is as follows: (A.1) represents the expected net present value which is to be maximized over the set of scenarios S. Appendix:

Table 8
Various cost parameters used in the example

Table 12
Demand for final product in the example TJ, ∀s ∈ S, ∀t ∈ T (A.10) Constraint (A.11) forces the sales to satisfy demand exactly for each period t and scenario s, δ t u sales,s t = D t , ∀s ∈ S, ∀t ∈ T (A.11)Constraint (A.12) restricts the input flow to each process by the capacity of the corresponding process.∀i∈I, ∀s ∈ S, ∀t ∈ T (A.12)The capacity of every process in period t and scenario s is computed by the available capacity in the previous time period and the capacity expansion made (A.12).Logical constraint (A.16) states that operating a process for a given scenario s at period t requires an expansion at any period τ before t − τ , i.e. τ = 1, . . ., t − τ .
Constraint (A.17) forces a process i to operate in period t if its capacity is expanded in that period., ∀i ∈ IU, ∀s ∈ S, ∀t ∈ T (A.17) Constraint (A Logic equation (A.25) states that scenarios s and s are indistinguishable if for each process unit that differentiates the scenarios (s, s ) the logic variables (w l,s i,t ) are true at the level l ∈ L(i, s, s ).Finally, the disjunctive constraint (A.26) includes certain non-anticipativity constraints into the model if scenarios s and s are indistinguishable at the end of period t .Note that the scenarios are generated assuming that during the partial resolution, the yield of process 1 is the lowest 0.69 or the highest 0.81.The actual yield at full resolution is either , ∀i ∈ IU, ∀(t, s, s )∈ N I (A.22) d qe,s i,t = d qe,s i,t , ∀i ∈ I, ∀(t, s, s ) ∈ N I (A.23) d rate,s k,t = d rate,s k,t , ∀k ∈ DK, ∀(t, s, s ) ∈ N I (A.24) t , ∀t ∈ T , ∀q ∈ Q, ∀l ∈ L, (s, s ) ∈ M q,l (A.25)