title: Multistage Robust Mixed-Integer Optimization Under Endogenous Uncertainty
authors: Wei Feng Yiping Feng Qi Zhang
text: 
  
  Buy Me That Look: An Approach for Recommending Similar Fashion Products
      
  Abhinav Ravi
  Data Sciences 
  
  Myntra
  
  Bengaluru, India 
  
  abhinav.ravi@myntra.com
  Sandeep Repakula
  Data Sciences 
  
  Myntra
  
  Bengaluru, India 
  
  sandeep.r@myntra.com
  Ujjal Kr Dutta
  Data Sciences 
  
  Myntra
  
  Bengaluru, India 
  
  ujjal.dutta@myntra.com
  Maulik Parmar
  Data Sciences 
  
  Myntra
  
  Bengaluru, India 
  
  parmar.m@myntra.com
  
  
      
  =====================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  Endogenous, i.e. decision-dependent, uncertainty has received increased interest in the stochastic programming community. In the robust optimization context, however, it has rarely been considered. This work addresses multistage robust mixed-integer optimization with decision-dependent uncertainty sets. The proposed framework allows us to consider both continuous and integer recourse, including recourse decisions that affect the uncertainty set. We derive a tractable reformulation of the problem by leveraging recent advances in the construction of nonlinear decision rules, and introduce discontinuous piecewise linear decision rules for continuous recourse. Computational experiments are performed to gain insights on the impact of endogenous uncertainty, the benefit of discrete recourse, and computational performance. Our results indicate that the level of conservatism in the solution can be significantly reduced if endogenous uncertainty and mixed-integer recourse are properly modeled.
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Robust optimization has proven to be an effective approach to decision making under uncertainty and has received considerable attention in recent years. While earlier works have only considered the case without recourse <cit.>, lending it a reputation of being overly conservative, effective means have been developed in the last decade to also address the case with recourse, in two- and multistage settings <cit.>. More recent efforts have focused on incorporating not only continuous but also discrete recourse <cit.>. For comprehensive reviews of the literature on robust optimization, we refer the reader to <cit.>, <cit.>, and <cit.>.
  
  The vast majority of existing works on robust optimization consider exogenous uncertainty, which is characterized by fixed uncertainty sets. Significantly fewer have addressed the case of endogenous, i.e. decision-dependent, uncertainty although, in practice, many uncertainties are inherently endogenous. Endogenous uncertainty has its origin in stochastic programming <cit.>. Mainly two types of endogenous uncertainty have been considered in the literature: (1) In the case of endogenous uncertainty of type 1, decisions affect the realization of the uncertain parameter by altering its underlying probability distribution. For example, a company may shift the probability distribution of the demand for its product toward higher values by lowering the selling price. (2) In the case of endogenous uncertainty of type 2, decisions affect the time at which an uncertain parameter materializes or its true value is revealed. The classical example for type-2 endogenous uncertainty is the size of an oilfield for which the true value cannot be determined until one starts drilling and extracting oil from it.
  
  The literature addressing type-1 endogenous uncertainty is relatively sparse. <cit.> considers single-stage stochastic network problems and uses Luce's choice axiom to develop an expression for the probability of routing along a path, which depends on the network design variables. <cit.> formulate a pre-disaster investment problem in which the failure probabilities of links in a transportation network can be altered by investment decisions related to strengthening those links. Discrete investment decisions are considered, which results in a two-stage stochastic program that allows the choice between a finite number of sets of failure probabilities. <cit.> apply a similar approach to a three-stage resource allocation planning problem for natural disaster relief under type-1 endogenous uncertainty. <cit.> propose two-stage models with probability distributions that are distorted through an affine transformation or a convex combination of multiple probability distributions. In addition, the authors consider parameterized distributions with the parameters being first-stage decision variables.
  
  A larger number of existing works focus on type-2 endogenous uncertainty. Applications include oil/gas field development <cit.>, capacity expansion in process networks <cit.>, open-pit mine production scheduling <cit.>, clinical trial planning <cit.>, R&D project portfolio management <cit.>, and vehicle routing <cit.>. Most commonly, the problem is formulated as a multistage stochastic program with discrete scenarios and explicit nonanticipativity constraints (NACs) that are active or inactive depending on the decisions related to the endogenous uncertainty. This results in a model that encodes a very large conditional scenario tree, which dramatically increases the computational complexity compared to the case with only exogenous uncertainty. Significant advances have been made in solving such stochastic programs, with approaches that focus on two general strategies: identifying redundant NACs that can be removed <cit.>, and applying tailored solution algorithms based on Lagrangean decomposition <cit.>, branch-and-cut <cit.>, Benders decomposition <cit.>, sequential scenario decomposition <cit.>, or heuristic knapsack decomposition <cit.>. <cit.> apply decision rules to obtain tractable conservative approximations for multistage stochastic programs with type-2 endogenous uncertain parameters that are continuously distributed. We refer to <cit.> for a comprehensive review of existing works in this area as well as a unifying framework that addresses problems with both exogenous and type-2 endogenous uncertainty.
  
  Only recently, endogenous uncertainty has also been considered in robust optimization. <cit.> and <cit.> address static robust optimization, i.e. without recourse, with decision-dependent polyhedral uncertainty sets. Similarly, <cit.> incorporate endogenous uncertainty into a multistage robust process scheduling framework; however, only first-stage decisions can affect the uncertainty set, and all recourse variables are continuous. <cit.> consider what the authors refer to as decision-dependent information discovery in two- and multistage robust optimization settings using a K-adaptability approach; here, the uncertainty set is fixed but the decision maker can decide whether or not to observe some uncertain parameters.
  
  In this work, we consider multistage robust optimization with mixed-integer recourse and decision-dependent uncertainty sets that can be altered in every stage. Specifically, we focus our discussion on polyhedral uncertainty sets that depend linearly on binary decision variables, which only affect the right-hand sides of the inequalities defining the uncertainty set. To derive tractable approximations of the resulting problem, we apply a decision rule approach based on the concept of lifted uncertainty <cit.>. This approach has recently been applied to model predictive control <cit.>, reservoir management <cit.>, and transmission expansion planning <cit.>. Here, we further expand the state of the art by introducing discontinuous piecewise linear decision rules for continuous recourse variables. The proposed framework significantly expands our capability to appropriately model endogenous uncertainty in robust optimization settings, with applicability in a variety of areas, such as network design, revenue management, and multiperiod planning. In our computational experiments, we demonstrate the benefits of considering endogenous uncertainty as well as both continuous and binary recourse, and provide results on the computational performance of the proposed models.
  
  The remainder of this paper is organized as follows. In Section <ref>, we present a two-stage robust mixed-integer optimization formulation with endogenous uncertainty, approximate it using decision rules in a lifted space, and derive a tractable mixed-integer linear programming (MILP) reformulation. The proposed approach is then extended to the multistage case in Section <ref>. In Section <ref>, we apply the proposed models to a two-stage design problem and a multistage production planning problem. Finally, in Section <ref>, we close with some concluding remarks.
  
  
  
    
  Notation We use lowercase and uppercase boldface letters to denote vectors and matrices, respectively, e.g. x∈ℝ^n and A∈ℝ^n × m. Scalar quantities are denoted by non-boldface letters. We use ∘ and 1(·) to denote the Hadamard multiplication operator and the indicator function, respectively. Furthermore, 0 and e denote the zero and all-ones vectors, respectively, while e_k is the standard basis vector whose kth element is 1.
  
  
  
  
  § THE TWO-STAGE CASE
  
  
  
  From a technical standpoint, the main novelty of this work lies in the incorporation of uncertainty sets that can be affected by binary recourse decisions across multiple stages. This requires a nontrivial integration of decision-dependent uncertainty sets and a means of modeling binary recourse, for which we apply a decision rule approach based on the concept of lifted uncertainty. However, primarily for the sake of clarity, we first discuss the two-stage case in this section, where the uncertainty set only depends on first-stage decisions. Here, we focus on the description of the decision-dependent uncertainty set and the lifting technique that also allows us to introduce discontinuous piecewise linear decisions rules for continuous recourse variables, which are significantly more flexible than traditional affine decision rules. Then, in Section <ref>, we turn to the general multistage case, which involves uncertainty sets that depend on binary recourse decisions.
  
  In the following, we present our approach to solving two-stage robust MILPs under endogenous uncertainty that can generally be formulated as follows:
  
  
  
      x_1 
           ξ^⊤( A_nx + D_ny) + ã_n^⊤x̃(ξ) + d̃_n^⊤ỹ(ξ) ≤ξ^⊤b_n  ∀  n ∈𝒩,  ξ∈Ξ(y) 
         x∈ℝ^P,  y∈{0,1}^Q
         x̃(ξ) ∈ℝ^P,  ỹ(ξ) ∈{0,1}^Q  ∀ ξ∈Ξ(y),
  
  
  where x and y are the first-stage variables, and x̃ and ỹ are the second-stage variables, which are functions of the uncertain parameters ξ∈ℝ^K. For ease of exposition, we assume that the objective function (<ref>) is certain, which is without loss of generality since x_1 can simply be the auxiliary variable introduced in an epigraph reformulation. The set of constraints is denoted by 𝒩, and as stated in (<ref>), they have to hold for all ξ in an uncertainty set Ξ. Variable domains are specified in constraints (<ref>) and (<ref>). Note that we assume fixed recourse and that all discrete variables are binary.
  
  Formulation (<ref>) implies that the uncertainty set Ξ depends on the binary decisions y. Indeed, similar to <cit.>, we consider decision-dependent uncertainty sets of the following form:
  
      Ξ(y) = {ξ∈ℝ^K: Wξ≤Uy,  ξ_1 = 1 },
  
  which is assumed to be a compact polyhedron for any feasible y, with W∈ℝ^M × K and U∈ℝ^M × Q. We assume that ξ_1 = 1, which eases the notation when introducing constant terms in the linear decision rules as shown in Section <ref>.
  
  
  Consider the following example of a decision-dependent uncertainty set:
  
      Ξ(y) = { ξ∈ℝ_+^3: 
        [                     ξ_2 ≤ 7 y_1 + 8 y_2;                            ξ_3 ≤ 13 y_2;                -8ξ_1 -ξ_2 +2ξ_3 ≤ 15y_2;       -13ξ_1 +ξ_2 + ξ_3 ≤ 7 y_1 + 2 y_2; 25 ξ_1 + 4 ξ_2 -7 ξ_3 ≤ 21 y_1 +  11y_2;              40 ξ_1 - 8 ξ_2 - 3 ξ_3 ≤ 0;                                 ξ_1 = 1 ]}.
  
  Here, ξ_3 can only be nonzero if y_2 = 1, and ξ_2 is forced to be 0 when y_1 and y_2 are both 0.
  
  Figure <ref> shows the projections of the uncertainty set Ξ onto the two-dimensional (ξ_2,ξ_3)-space for different values of y_1 and y_2. One can see that the choice of y affects the facets of Ξ, leading to different uncertainty sets. In this case, Ξ is largest when y_1=1 and y_2=1 (see Figure <ref>); in fact, Ξ(1,1) is a superset of Ξ(0,1) (Figure <ref>) and Ξ(1,0) (Figure <ref>), with the latter being merely the line for which 5 ≤ξ_2 ≤ 7. Note that it is not generally true that there always exists a Ξ(y^*) such that Ξ(y^*) is a superset of all Ξ(y) for y≠y^*. Furthermore, in this example, Ξ is an empty set for y_1=0 and y_2=0; hence, there has to be a constraint in problem (<ref>), e.g. y_1 + y_2 ≥ 1, that excludes this solution.
  
  
  
  Clearly, the decision-dependent uncertainty set given by (<ref>) captures type-1 endogenous uncertainty as the uncertain parameters' support changes with y. Depending on the constraints of the problem, it can also consider type-2 endogenous uncertainty. For example, if y_2 = 0 renders uncertain parameter ξ_3 irrelevant for the problem, e.g. if it is multiplied by y_2 everywhere it appears in the model, then the uncertainty set can encode the case in which ξ_3 only materializes if y_2 = 1. In general, we assume that an uncertain parameter whose materialization is decision-dependent can take the value zero, which is without loss of generality as it can always be achieved with a simple linear translation.
  
  
  In this work, we assume that the uncertainty set can only be affected by binary variables, which allows the final reformulation to be an MILP. The proposed approach can also be applied if y are continuous; however, in that case, we would arrive at a nonconvex mixed-integer nonlinear program (MINLP). Furthermore, the approach can be extended to consider random recourse and polyhedral uncertainty sets in which decisions also affect the left-hand sides; however, this would result in a substantial increase in computational complexity and is hence not discussed in this work.
  
  
  
  
  
  
  
  
  
   §.§ Lifted Uncertainty Set
  
  
  
  We apply the lifting technique proposed by <cit.> to derive tractable decision rules. The basic idea is to lift the original uncertain parameters onto a higher-dimensional space such that linear decision rules can be applied to the lifted parameters, allowing the construction of more flexible and binary decision rules.
  
  We first insert r_i-1 breakpoints into the marginal support of each uncertain parameter ξ_i such that
  
      ξ_i^min < p_i^1 < p_i^2 < ⋯ < p_i^r_i - 1 < ξ_i^max  ∀  i = 1, …, K,
  
  where ξ_i^min and ξ_i^max are the lower and upper bounds of ξ_i, respectively. As such, {ξ∈ℝ^K: ξ^min≤ξ≤ξ^max} is the smallest hyperrectangle that contains Ξ(y) for all feasible y. We now define a lifting operator L̅: ℝ^K↦ℝ^K that maps the original uncertain parameters ξ∈ℝ^K onto a K-dimensional space with K =∑_i=1^K r_i. The vector of lifted uncertain parameters is denoted by ξ̅ =  ( ξ̅_1, …, ξ̅_K  ) ∈ℝ^K with ξ̅_i =  ( ξ̅_i^1, …, ξ̅_i^r_i ) ∈ℝ^r_i, and the lifting operator L̅ = ( L̅_1, …, L̅_K ) with L̅_i = ( L̅_i^1, …, L̅_i^r_i) is defined as follows:
  
      ξ̅_i^j = L̅_i^j(ξ):= 
          {[                                 ξ_i                                  if                             r_i = 1;                    inf{ξ_i, p_i^j }                                  if                    r_i > 1,   j = 1; sup{inf{ξ_i, p_i^j } - p_i^j-1, 0 }                                  if        r_i > 1,   j = 2, …, r_i - 1;              sup{ξ_i - p_i^j-1, 0 }                                  if                   r_i > 1,   j=r_i. ] .
  
  As illustrated in Figure <ref>, L̅_i^j(ξ) is a piecewise linear function of ξ_i. By construction, the original and lifted uncertain parameters have the following relationship:
  
      ξ_i = e^⊤ξ̅_i   ∀  i = 1, …, K.
  
  
  
  
  
  
  
  
  Piecewise linear decision rules for continuous recourse variables can be derived using lifted uncertain parameters defined by L̅. To also allow binary recourse, we define a second lifting operator L: ℝ^K↦{0,1}^K as proposed by <cit.>. Here, we apply the same breakpoints as in (<ref>) and have K = ∑_i=1^K g_i with g_i = max{ 1,   r_i - 1 }. We introduce another vector of lifted uncertain parameters ξ̂ =  ( ξ̂_1, …, ξ̂_K  ) ∈{ 0, 1 }^K with ξ̂_i =  (ξ̂_i^1, …, ξ̂_i^g_i) ∈{ 0, 1 }^g_i, and define L =  (L_1, …, L_K  ) with L_i = ( L_i^1, …, L_i^g_i) as follows:
  
      ξ̂_i^j = L_i^j(ξ) := 
          {[                         1                        if                   r_i = 1;            1(ξ_i ≥ p_i^j)                        if r_i > 1,   j = 1, …, g_i, ] .
  
  which is illustrated in Figure <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We introduce a new vector of lifted uncertain parameters, ξ' =  (ξ̅,  ξ̂) ∈ℝ^K' with K' = ∑_i=1^K r_i + g_i. Note that ξ_1 = 1, hence r_1 = 1, L̅_1^1 (ξ) = 1, and L_1^1 (ξ) = 1. The lifted uncertainty set is then defined as follows:
  
      Ξ'(y) =  {ξ̃ = 
          ( ξ,  ξ') 
          = ( ξ,  ξ̅,  ξ̂) ∈ℝ^K + K':
          [     ξ∈Ξ(y); ξ̅ = L̅(ξ);  ξ̂= L( ξ) ]}.
  
  Due to the discontinuity of ξ̂ at each breakpoint (see Figure <ref>), Ξ' is an open set. Let Ξ̇ := cl(Ξ') be the closure of Ξ', Ξ̇_ξ' denote the projection of Ξ̇ onto the space of ξ', and the second-stage variables be functions of the lifted uncertain parameters ξ'. The two-stage problem then becomes
  
  
  
      x_1 
           ξ^⊤( A_nx + D_ny) + ã_n^⊤x̃(ξ') + d̃_n^⊤ỹ(ξ') ≤ξ^⊤b_n  ∀  n ∈𝒩,  ( ξ,  ξ') ∈Ξ̇( y) 
         x∈ℝ^P,  y∈{0,1}^Q
         x̃(ξ') ∈ℝ^P,  ỹ(ξ') ∈{0,1}^Q  ∀ ξ'∈Ξ̇_ξ'(y),
  
  
  which has the same optimal value as (<ref>), and there is a one-to-one mapping between feasible and optimal solutions of problems (<ref>) and (<ref>), as shown in <cit.>.
  
  However, because the uncertain parameters considered are correlated and decision-dependent, it is difficult to formulate the exact closed-form expression for Ξ̇(y). Hence, we introduce a tractable outer approximation Ξ(y), which stems from the convex hull of the marginal support of ξ̃_i for all i = 1,…,K. By replacing Ξ̇(y) with Ξ(y), the two-stage problem is subsequently transformed into (<ref>), whose solution is rather a conservative approximation of the one derived by (<ref>):
  
  
  
      x_1 
           ξ^⊤( A_nx + D_ny) + ã_n^⊤x̃(ξ') + d̃_n^⊤ỹ(ξ') ≤ξ^⊤b_n  ∀  n ∈𝒩,  ( ξ,  ξ') ∈Ξ( y) 
         x∈ℝ^P,  y∈{0,1}^Q
         x̃(ξ') ∈ℝ^P,  ỹ(ξ') ∈{0,1}^Q  ∀ ξ'∈Ξ_ξ'(y).
  
  
  
  To obtain Ξ(y), first notice that both L̅ and L can be considered piecewise mapping functions with each piece being a subinterval Ξ_i^j formed by two consecutive breakpoints:
  
      Ξ_i^j = {ξ_i ∈ℝ:   p_i^j-1≤ξ_i ≤ p_i^j }.
  
  For ease of exposition, we set p_i^0 = ξ_i^min and p_i^r_i = ξ_i^max. Let 𝒱^j_i = {p_i^j-1, p_i^j}, then for every i = 1,…,K, we can define the following sets of vertices in the lifted space:
  
      𝒱_i  =  ⋃_j=1^r_i𝒱_i^j 
         𝒱_i^j  = {ṽ_i =   (v_i,  v̅_i,  v̂_i ): 
           [                                     v_i ∈𝒱_i^j;                           v̅_i= L̅_i (e_i v_i); v̂_i = lim_ξ_i → v_i,  ξ_i ∈Ξ_i^jL_i (e_i ξ_i) ]}  ∀  j = 1, …, r_i,
  
  which allow us to formulate the following convex hull representation for the closure of the marginal support of ξ̃_i:
  
      Ξ_i = {ξ̃_i =  (ξ_i, ξ̅_i, ξ̂_i ):
              [           ∑_j=1^r_i∑_ṽ_i ∈𝒱_i^jλ_i^j(ṽ_i) = 1;     ξ_i = ∑_j=1^r_i∑_ṽ_i ∈𝒱_i^jλ_i^j(ṽ_i) v_i;   ξ̅_i = ∑_j=1^r_i∑_ṽ_i ∈𝒱_i^jλ_i^j(ṽ_i) v̅_i;   ξ̂_i = ∑_j=1^r_i∑_ṽ_i ∈𝒱_i^jλ_i^j(ṽ_i) v̂_i; λ_i^j(ṽ_i) ∈ℝ_+   ∀  j=1, …, r_i,  ṽ_i ∈𝒱_i^j ]},
  
  where λ_i^j(ṽ_i) denotes the coefficient associated with a particular vertex ṽ_i ∈𝒱_i^j. If the original uncertain parameters are independent and exogenous, Ξ̇ can be exactly represented as the Cartesian product of all Ξ_i. In the general and endogenous case, we have
  
      Ξ̇(y) ⊆Ξ(y) := {∏_i=1^K Ξ_i }⋂{ (ξ,  ξ' ) ∈ℝ^K + K':
        Wξ≤Uy,  ξ_1 = 1 }.
  
  Since Ξ(y) is generally a superset of Ξ̇(y), it may result in a more conservative solution. However, it is worth mentioning that the uncertainty set in the space of the original uncertain parameters, i.e. the projection of Ξ(y) onto the ξ-space, remains unchanged. This outer approximation only applies to the new lifted uncertain parameters ξ', which are used for the construction of the decision rules, as we will show in Section <ref>. Therefore, the use of this outer approximation should not be interpreted as expanding the uncertainty set to robustify against, but rather as further restricting the set of possible decision rules. As such, this increased conservatism could be compensated by adjusting the locations or increasing the number of breakpoints.
  
  
  In this small example, we illustrate the relationship between Ξ̇(y) and Ξ(y). Consider the following decision-dependent uncertainty set:
  
      Ξ(y) = {ξ∈ℝ^2:   0 ≤ξ_2 ≤ 5 ξ_1 - y,  ξ_1  = 1 }.
  
  We place one breakpoint, p_2^1, at the center of the marginal support of ξ_2, and analyze the resulting Ξ̇(y) and Ξ(y). For ease of visualization, we only show the projections of these two sets onto the two-dimensional (ξ_2,ξ̂_2^1)-space for different values of y. If y=0, Ξ̇(y) and Ξ(y) are in fact the same, which is depicted by the gray-shaded area in Figure <ref>. However, if y = 1, the projection of Ξ̇(y) is the red-shaded area in Figure <ref>, while that of Ξ(y) is the green-shaded area in Figure <ref>. Clearly, Ξ̇(y) is a subset of Ξ(y).
  
  
  
  
  
  
   §.§ Decision Rule Approximation and Reformulation
  
  
  
  Following the decision rule approach, we solve problem (<ref>) approximately by restricting the recourse decisions to adjust according to some decision rules, which are functions of the uncertain parameters. Here, we apply decision rules parameterized as follows:
  
  
  
      x̃=∑_i=1^K ( X_i ξ̅_i + X_i ξ̂_i )
          
         ỹ = ∑_i=1^K Y_i ξ̂_i
  
  
  where X_i ∈ℝ^P× r_i, X_i ∈ℝ^P× g_i, and Y_i ∈{-1, 0, 1}^Q× g_i for i=1, …, K. We emphasize that the structure of (<ref>) allows discontinuous piecewise linear decision rules; this is in contrast to most existing works in the literature, which use affine or continuous piecewise linear decisions rules for continuous recourse variables. Also note that, as shown in <cit.>, the given domain for Y_i is sufficient to allow binary ỹ using the decision rule in (<ref>). Before proceeding to the reformulation, we illustrate with the following simple example that, unlike in continuous optimization, it is often crucial for the continuous recourse variables to follow discontinuous decision rules in mixed-integer optimization.
  
  
  Consider the following MILP:
  
      2x̃_1 +  x̃_2 
           ỹ_1 ≤x̃_1 ≤ 3ỹ_1 
          3ỹ_2 ≤x̃_2 ≤ 5ỹ_2 
         x̃_1 + x̃_2 = ξ
         ỹ_1 + ỹ_2 ≤ 1 
         x̃_1,  x̃_2 ∈ℝ
         ỹ_1,  ỹ_2 ∈{0,1 },
  
  where the continuous variables x̃_1 and x̃_2 take nonzero values if and only if the respective binary variables ỹ_1 and ỹ_2 are equal to 1. For ξ∈Ξ = {ξ∈ℝ:   1 ≤ξ≤ 5}, the optimal solution to (<ref>) as a function of ξ is as follows:
  
      {[ ỹ_1 = 1,  x̃_1 = ξ ,  ỹ_2 = x̃_2 = 0                        if   1 ≤ξ < 3;  ỹ_1 = x̃_1 = 0,  ỹ_2 = 1,  x̃_2 = ξ                         if   3 ≤ξ≤ 5 ] .
  
  which is shown in Figure <ref>. One can see that x̃_1 and x̃_2 follow discontinuous piecewise linear functions. In fact, (<ref>) represents the only feasible solution when ξ≠ 3 and the only optimal solution when ξ = 3. As a result, restricting the continuous variables to be continuous piecewise linear functions of ξ would render the problem infeasible.
  
   
  
  
  By substituting the decisions rules (<ref>) into (<ref>), constraints (<ref>) become
  
      ξ^⊤f_n + ∑_i=1^K [
          ã_n^⊤X_i ξ̅_i + 
           (
          ã_n^⊤X_i + d̃_n^⊤Y_i
           ) ξ̂_i ]
          ≤ 0   ∀  n ∈𝒩,   (ξ,  ξ̅,  ξ̂) ∈Ξ(y)
  
  with f_n = A_n x + D_n y - b_n. Following standard robust optimization techniques, we first apply the worst-case reformulation:
  
      max_ (ξ,  ξ̅,  ξ̂) ∈Ξ(y) {ξ^⊤f_n + ∑_i=1^K [
          ã_n^⊤X_i ξ̅_i + 
           (
          ã_n^⊤X_i + d̃_n^⊤Y_i
           ) ξ̂_i ] }≤ 0   ∀  n ∈𝒩,
  
  which then leads to the following reformulation due to strong duality of the left-hand-side maximization problems:
  
      e^⊤ρ_n + μ_n^⊤Uy≤ 0   ∀  n ∈𝒩
         ρ_ni + μ_n^⊤w_i v_i ≥ f_ni v_i + ã_n^⊤X_i v̅_i + ( ã_n^⊤X_i + d̃_n^⊤Y_i ) v̂_i 
                             ∀  n ∈𝒩,   i=1, …, K,  (v_i,  v̅_i,  v̂_i ) ∈𝒱_i 
         ρ_n ∈ℝ^K,  μ_n ∈ℝ_+^M   ∀  n ∈𝒩,
  
  where ρ_n and μ_n are the dual variables of the maximization problem,  f_ni and ρ_ni denote the ith element of f_n and ρ_n respectively, and w_i is the ith column vector of matrix W.
  
  In addition, since the decision rule in (<ref>) is guaranteed to yield integer ỹ, the integratility constraints ỹ(ξ') ∈{ 0,1 }^Q in (<ref>) can be relaxed to
  
      0≤ỹ(ξ') ≤e  ∀ ξ'∈Ξ_ξ'(y),
  
  which, using similar arguments as above, can be reformulated into the following set of constraints:
  
      Ωe - ΨUy≥0
         Y_i v̂_i - ω_i + Ψw_i v_i ≥0  ∀  i=1, …, K,   (v_i,  v̅_i,  v̂_i ) ∈𝒱_i 
         Ωe + ΨUy≤e
         Y_i v̂_i - ω_i - Ψw_i v_i ≤0  ∀  i=1, …, K,  (v_i,  v̅_i,  v̂_i ) ∈𝒱_i 
         Ω,  Ω∈ℝ^Q× K,   Ψ,  Ψ∈ℝ_+^Q× M,
  
  where Ω, Ω, Ψ and Ψ are the matrices of dual variables, and ω_i and ω_i are the ith columns of Ω and Ω, respectively.
  
  Finally, we arrive at the following formulation:
  
      x_1 
           f_n = A_n x + D_n y- b_n   ∀  n ∈𝒩
         (<ref>),  (<ref>)
         x∈ℝ^P,  y∈{0,1}^Q
         X_i ∈ℝ^P× r_i,  X_i ∈ℝ^P× g_i,  Y_i ∈{-1, 0, 1}^Q× g_i  ∀  i=1, …, K.
  
  Note that (<ref>) and (<ref>) contain bilinear terms; however, each bilinear term consists of a continuous variable multiplied by a binary variable and can hence be exactly linearized <cit.>. As a result, problem (<ref>) can be solved as an MILP.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § THE MULTISTAGE CASE
  
  
  
  In this section, the proposed methodology is extended to multistage robust MILPs with endogenous uncertainty of the following form:
  
  
  
      e_1^⊤x_1 
           A^1 x_1 + D^1 y_1 ≤b^1 
         A^t(ξ^[t]) x_1 + D^t(ξ^[t]) y_1 + ∑_t'=2^t [ A_t'^t x_t'(ξ^[t']) + D_t'^t y_t'(ξ^[t']) ] ≤b^t (ξ^[t]) 
                                       ∀  t = 2, …, T,  ξ^[t]∈Ξ^[t](y^[t-1]) 
         x_1 ∈ℝ^P_1,  y_1 ∈{ 0, 1 }^Q_1
         x_t(ξ^[t]) ∈ℝ^P_t,  y_t(ξ^[t]) ∈{ 0, 1 }^Q_t  ∀  t = 2, …, T,  ξ^[t]∈Ξ^[t]( y^[t-1])
  
  
  where x_t and y_t are the continuous and binary variables, respectively, in stage t. The uncertain parameters whose true values are observed in stage t are denoted by ξ_t = (ξ_t1, …, ξ_tK_t) ∈ℝ^K_t. Again for notational convenience, K_1 = 1 and the only element of ξ_1, ξ_11, is assumed to be 1. We further define all uncertain parameters observed up to stage t as ξ^[t] = ( ξ_1, ξ_2, …, ξ_t ). Moreover, A^t(ξ^[t]) ∈ℝ^N_t × P_1, D^t(ξ^[t]) ∈ℝ^N_t × Q_1, and b^t (ξ^[t]) ∈ℝ^N_t are coefficient matrices or vectors that depend linearly on ξ^[t], i.e.
  
      A^t(ξ^[t]) = ∑_t'=1^t ∑_i=1^K_t'A_t'i^t ξ_t'i,   D^t(ξ^[t]) = ∑_t'=1^t ∑_i=1^K_t'D_t'i^t ξ_t'i,   b^t (ξ^[t]) = ∑_t'=1^t ∑_i=1^K_t'b_t'i^t ξ_t'i,
  
  where A_t'i^t ∈ℝ^N_t × P_1, D_t'i^t ∈ℝ^N_t × Q_1, and b_t'i^t ∈ℝ^N_t.
  
  Formulation (<ref>) indicates that the uncertainty set Ξ^[t] in stage t depends on binary decisions made in previous stages, y^[t-1] = (y_1, …, y_t-1), where y_t̂(ξ^[t̂]) are recourse variables for t̂≥ 2. This means that, in contrast to the two-stage case, the changing of the uncertainty set is now also a recourse decision. Here, we consider decision-dependent uncertainty sets of the following form:
  
      Ξ^[t](y^[t-1]) = {ξ^[t]∈ℝ^K^[t]:
          W^[t]ξ^[t]≤∑_t̂=1^t-1U_t̂^t y_t̂(ξ^[t̂]),  ξ_11 = 1
          },
  
  where K^[t] = ∑_t'=1^t K_t', W^[t]∈ℝ^M_t × K^[t] and U_t̂^t ∈ℝ^M_t × Q_t̂. We assume that Ξ^[t](y^[t-1]) is a compact polyhedron for any feasible y^[t-1]. Note that Ξ^[t](y^[t-1]) reduces to the uncertainty set in the two-stage problem shown in (<ref>) when t=2.
  
  
  
   §.§ Lifted Uncertainty Set
  
  
  With r_ti-1 breakpoints placed inside the marginal support of each uncertain parameter ξ_ti, the two lifting operators defined in Subsection <ref> can also be applied in the multistage case. We have L̅_t: ℝ^K_t↦ℝ^K_t with K_t = ∑_i=1^K_t r_ti and L_t: ℝ^K_t↦ℝ^K_t with K_t = ∑_i=1^K_t g_ti, where g_ti = max{ 1, r_ti-1 }, for all t=1, …, T. The lifted uncertain parameters ξ̅_t =  (ξ̅_t1, …, ξ̅_t K_t) and ξ̂_t =  (ξ̂_t1, …, ξ̂_t K_t) are then defined as follows:
  
  
  
      ξ̅_ti^j = L̅_ti^j ( e_tiξ_ti)   ∀  t=1, …, T,    i =1,…, K_t,   j=1, …, r_ti
         ξ̂_ti^j = L_ti^j  ( e_tiξ_ti)   ∀  t=1, …, T,   i =1,…, K_t,   j=1, …, g_ti.
  
  
  
  Let ξ̃_t = ( ξ_t, ξ̅_t, ξ̂_t ) ∈ℝ^K_t with K_t = K_t + K_t + K_t and ξ̃^[t] = ( ξ̃_1, …, ξ̃_t ). The lifted uncertainty set in stage t is then
  
      Ξ'^[t](y^[t-1]) = {ξ̃^[t]∈ℝ^K^[t]: 
          [ ξ^[t] =  ( ξ_1, ξ_2, …, ξ_t  ) ∈Ξ^[t](y^[t-1]);             ξ̅_t' = L̅_t'(ξ_t')   ∀  t'=1,…, t;              ξ̂_t' = L_t'(ξ_t')   ∀  t'=1,…, t ]},
  
  where K^[t] = ∑_t'=1^t K_t'.
  
  
  
  
  
  
  
  
  
  
  
  Since Ξ'^[t] is an open set due to the discontinuity in ξ̂_t', we also require the outer approximation of its closure, which we denote by Ξ^[t]. For every t=1,…,T and i=1, …, K_t, we define the following sets of vertices in the lifted space:
  
      [                                                                                                                                                                                                                 𝒱_ti = ⋃_j=1^r_ti𝒱_ti^j; 𝒱_ti^j = {ṽ_ti = ( v_ti, v̅_ti, v̂_ti):
               [                                         v_ti∈𝒱_ti^j;                           v̅_ti = L̅_ti( e_ti v_ti); v̂_ti = lim_ξ_ti→ v_ti,  ξ_ti∈Ξ_ti^jL_ti( e_tiξ_ti) ]}  ∀  j = 1, …, r_ti, ]
  
  where Ξ_ti^j = {ξ_ti∈ℝ: p_ti^j-1≤ξ_ti≤ p_ti^j } and 𝒱_ti^j = { p_ti^j-1,   p_ti^j }, assuming that p_ti^0 = ξ_ti^min and p_ti^r_ti = ξ_ti^max. This allows us to formulate the following convex hull representation of Ξ^[t]:
  
      [ Ξ^[t](y^[t-1]) = {ξ̃^[t]∈ℝ^K^[t]:
          [            ∑_j=1^r_t'i∑_ṽ_t'i∈𝒱_t'i^jλ_t'i^j(ṽ_t'i) = 1   ∀  t' = 1, …, t,   i = 1, …, K_t';   ξ_t'i = ∑_j=1^r_t'i∑_ṽ_t'i∈𝒱_t'i^jλ_t'i^j(ṽ_t'i) v_t'i  ∀  t' = 1, …, t,   i = 1, …, K_t'; ξ̅_t'i = ∑_j=1^r_t'i∑_ṽ_t'i∈𝒱_t'i^jλ_t'i^j(ṽ_t'i) v̅_t'i  ∀  t' = 1, …, t,   i = 1, …, K_t'; ξ̂_t'i = ∑_j=1^r_t'i∑_ṽ_t'i∈𝒱_t'i^jλ_t'i^j(ṽ_t'i) v̂_t'i  ∀  t' = 1, …, t,   i = 1, …, K_t';                                                    W^[t]ξ^[t]≤∑_t̂=1^t-1U_t̂^ty_t̂(ξ̃^[t̂]);                    λ_t'i^j(ṽ_t'i) ∈ℝ_+   ∀  t' = 1, …, t,   i = 1, …, K_t',   j=1, …, r_t'i ]}. ]
  
  
  
  
   §.§ Decision Rule Approximation and Reformulation
  
  
  
  Analogous to the decision rules (<ref>) introduced in the two-stage problem, we apply the following decision rules to the recourse variables x_t and y_t for t ≥ 2:
  
  
  
      x_t = ∑_t'=1^t ∑_i=1^K_t'X_t'i^t ξ̅_t'i + X_t'i^t ξ̂_t'i
         y_t = ∑_t'=1^t ∑_i=1^K_t'Y_t'i^t  ξ̂_t'i
  
  
  where X_t'i^t ∈ℝ^P_t × r_t'i, X_t'i^t ∈ℝ^P_t × g_t'i, and Y_t'i^t ∈{-1,  0,   1}^Q_t × g_t'i. Note that in the case of type-2 endogenous uncertainty, nonanticipativity has to be further enforced depending on which uncertain parameters materialize. This is encoded in the definition of the uncertainty set, where an unmaterialized uncertain parameter is forced to take the value zero. For linear decision rules like the ones in (<ref>), this is equivalent to forcing the decision rule coefficients corresponding to unmaterialized uncertain parameters to zero.
  
  By substituting (<ref>) and (<ref>) into formulation (<ref>) and replacing Ξ^[t](y^[t-1]) with Ξ^[t](y^[t-1]), constraints (<ref>) become
  
      ∑_t'=1^t ∑_i=1^K_t'f_t'i^t ξ_t'i +  ∑_t” = 2 ^t ∑_t'=1^t”∑_i=1^K_t'[ A_t”^t X_t'i^t”ξ̅_t'i + ( A_t”^t X_t'i^t” + D_t”^t Y_t'i^t”) ξ̂_t'i] ≤0
        ∀  t=2, …, T, ξ̃^[t]∈Ξ^[t](y^[t-1])
  
  with f_t'i^t = A_t'i^t x_1 + D_t'i^t y_1 - b_t'i^t. The worst-case reformulation of (<ref>) is then
  
      max_ξ̃^[t]∈Ξ^[t](y^[t-1]) {∑_t'=1^t ∑_i=1^K_t'( f_t'i^t ξ_t'i +  ∑_t” = max{2,   t' }^t [ A_t”^t X_t'i^t”ξ̅_t'i + ( A_t”^t X_t'i^t” + D_t”^t Y_t'i^t”) ξ̂_t'i] ) }≤0
      ∀  t=2, …, T,
  
  which in turn can be reformulated into the following set of constraints: 
  
      
      Φ_t U_1^t y_1 + ∑_t'=1^t∑_i=1^K_t'δ_t'i^t ≤0  ∀  t=2,…, T  
         (f_t'i^t  - Φ_tw_t'i^[t]) v_t'i - δ_t'i^t + ∑_t” = max{2,   t' }^t [ A_t”^t X_t'i^t”v̅_t'i + ( A_t”^t X_t'i^t” + D_t”^t Y_t'i^t” + Φ_t U_t”^t Y_t'i^t”) v̂_t'i] ≤0
                                 ∀  t=2,…,T,   t' = 1, …, t,   i = 1, …, K_t',  ( v_t'i,  v̅_t'i,  v̂_t'i) ∈𝒱_t'i
         Φ_t ∈ℝ_+^N_t × M_t,  δ_t'i^t ∈ℝ^N_t  ∀  t=2,…,T,   t'= 1, …, t,   i=1,…, K_t'
  
  
  where U_t^t = 0 is introduced for notational convenience. The dual variables associated with the maximization problems in (<ref>) are denoted by Φ_t and δ_t'i^t, and w_t'i^[t] is the column vector corresponding to ξ_t'i in matrix W^[t]. Constraints (<ref>) contain bilinear terms comprised of continuous and integer variables as Y_t'i^t∈{-1,   0,   1}^Q_t× g_t'i. To facilitate the linearization, we express Y_t'i^t as the difference between binary variables, i.e.
  
      Y_t'i^t = Ẏ_t'i^t -Ÿ_t'i^t,   Ẏ_t'i^t + Ÿ_t'i^t≤1
  
  where Ẏ_t'i^t,  Ÿ_t'i^t∈{0,1}^Q_t× g_t'i, and the latter inequalities can be added to eliminate symmetry. Constraints (<ref>) then become
  
      (f_t'i^t  - Φ_tw_t'i^[t]) v_t'i - δ_t'i^t 
          + ∑_t” = max{2, t' }^t ( A_t”^t X_t'i^t”v̅_t'i + 
          
          [ A_t”^t X_t'i^t” +  (D_t”^t  + Φ_t U_t”^t ) ( Ẏ_t'i^t” -Ÿ_t'i^t”)  ] v̂_t'i) ≤0
                         ∀  t=2,…,T,   t' = 1, …, t,   i = 1, …, K_t',  ( v_t'i,  v̅_t'i,  v̂_t'i) ∈𝒱_t'i
  
  
  The integrality constraints on y_t in (<ref>) can be relaxed, i.e. 0≤y_t(ξ^[t]) ≤e, and reformulated as follows:
  
      .
          [                                                                                ∑_t'=1^t ∑_i=1^K_t'π_t'i^t - Θ^t U_1^t y_1 ≥0; ( Ẏ_t'i^t -Ÿ_t'i^t) v̂_t'i - π_t'i^t + Θ^t w_t'i^[t] v_t'i - ∑_t” = max{2,   t' }^t Θ^t U_t”^t (Ẏ_t'i^t” -Ÿ_t'i^t”) v̂_t'i≥0;                                                         ∀  t' = 1, …, t,   i= 1, …, K_t',  ( v_t'i,  v̅_t'i,  v̂_t'i) ∈𝒱_t'i;                                                        Θ^t ∈ℝ_+^Q_t × M_t,  π_t'i^t ∈ℝ^Q_t  ∀  t' = 1, …, t,   i= 1, …, K_t' ]}∀  t=2, …, T  
          . [                                                                                ∑_t'=1^t ∑_i=1^K_t'π_t'i^t + Θ^t U_1^t y_1 ≤e; ( Ẏ_t'i^t -Ÿ_t'i^t) v̂_t'i - π_t'i^t - Θ^t w_t'i^[t] v_t'i + ∑_t” = max{2,   t' }^t Θ^t U_t”^t (Ẏ_t'i^t” -Ÿ_t'i^t”) v̂_t'i≤0;                                                         ∀  t' = 1, …, t,   i= 1, …, K_t',  ( v_t'i,  v̅_t'i,  v̂_t'i) ∈𝒱_t'i;                                                       Θ^t ∈ℝ_+^Q_t × M_t,  π_t'i^t  ∈ℝ^Q_t  ∀  t' = 1, …, t,   i= 1, …, K_t' ]}∀  t=2, …, T.
  
  
  
  Finally, we arrive at the following reformulation of the multistage problem:
  
  
  
      e_1^⊤x_1 
           A^1 x_1 + D^1 y_1 ≤b^1 
         f_t'i^t = A_t'i^t x_1 + D_t'i^t y_1 - b_t'i^t   ∀  t= 2, …, T,   t' = 1, …, t,   i=1, …, K_t'
         (<ref>),  (<ref>),  (<ref>),  (<ref>)
         x_1 ∈ℝ^P_1,  y_1 ∈{ 0, 1 }^Q_1
         X_t'i^t ∈ℝ^P_t × r_t'i,  X_t'i^t ∈ℝ^P_t × g_t'i,  Ẏ_t'i^t,  Ÿ_t'i^t∈{0,1}^Q_t× g_t'i
                             ∀  t=2, …, T,  t' = 1, …, t,   i=1, …, K_t' .
  
  
  Every bilinear term in (<ref>) is composed of a continuous and a binary variable; hence, an MILP formulation is readily obtained after exact linearization of the bilinear terms.
  
  Note that in the decision rules given by (<ref>), a recourse variable in stage t is a function of all uncertain parameters observed up to stage t, i.e. the decision rules make use of all available information. In practice, however, this may not be the best choice as the model size and hence the computational performance strongly depend on the number of parameters involved in the decision rules. It has been observed in several multistage applications <cit.> that the optimal decision rules usually only depend on a small subset of uncertain parameters; hence, a common strategy to reduce computation time is to restrict the decision rules to depend on a smaller set of uncertain parameters. An intuitive choice is to let a recourse variable in stage t only depend on the uncertain parameters observed in the previous Δ t and the current stages. The required change in the reformulation is shown in Appendix <ref>.
  
  
  
  
  § COMPUTATIONAL CASE STUDIES
  
  
  
  In this section, we apply the proposed approach to a two-stage design problem and to a multistage production planning problem. All model instances were implemented in Julia v1.2.0 using the modeling environment JuMP v0.20.0 <cit.> and solved to 1% optimality using Gurobi v8.1.1 on a Intel Core i7-8700 CPU at 3.20 GHz machine with 8 GB RAM.
  
  
  
   §.§ Design for Flexible Production
  
  
  We consider the design of a production system that manufactures a single product for which the required production amount can vary across a wide range. Such production flexibility is especially important in systems with little or no product inventory capacity, e.g. in electricity generation or if the product is highly volatile. The production system can consist of a set of production units that all produce the same product but differ in capacity and cost. The production cost of a unit is assumed to be an affine function of the production amount. Moreover, while the minimum production amount, c^min, is known, the maximum production amount (i.e. the production capacity) is uncertain and is only known after the unit is built. However, we do know that the capacity will be between c^max - ĉ^max and c^max.
  
  
  
  
  
  
  
  Given a set of production units ℐ, the objective is to decide which subset of units to build such that all product demand d within a range [d^min, d^max] can be met exactly and the worst-case total cost is minimized. We can formulate the problem as the following two-stage robust optimization problem:
  
  
  
      ∑_i ∈ℐα_i z_i + max_ĉ∈𝒞(z), d ∈𝒟∑_i ∈ℐβ_i ỹ_i(ĉ,d) + γ_i x̃_i(ĉ,d) 
            z_i ∈{0,1}  ∀  i ∈ℐ
         ∑_i ∈ℐx̃_i(ĉ,d) = d   ∀ ĉ∈𝒞(z),   d ∈𝒟
         .
      [               ỹ_i(ĉ,d) ≤ z_i; x̃_i(ĉ,d) ≥ c^min_i ỹ_i(ĉ,d); x̃_i(ĉ,d) ≤ c^max_i ỹ_i(ĉ,d);    x̃_i(ĉ,d) ≤ c^max_i - ĉ_i;               x̃_i(ĉ,d) ∈ℝ_+;              ỹ_i(ĉ,d) ∈{0,1} ]}  ∀  i ∈ℐ,  ĉ∈𝒞(z),   d ∈𝒟
  
  
  with the uncertainty sets 𝒞(z) and 𝒟 defined as follows:
  
      𝒞(z) = {ĉ∈ℝ_+^|𝒥|: ĉ≤ĉ^max∘z}
  
  
      𝒟 = { d ∈ℝ_+: d^min≤ d ≤ d^max}.
  
  
  In problem (<ref>), z are the first-stage design variables while x̃ and ỹ are the second-stage operational variables. Production unit i is built if z_i = 1; ỹ_i = 1 if unit i is used to manufacture the product, and x̃_i denotes the corresponding production amount. In the objective function, α_i, β_i, and γ_i denote the capital, fixed production, and variable production costs for unit i, respectively.
  
  
  
    §.§.§ The Benefit of Discrete Recourse
  
  
  Demand d is modeled as an exogenous uncertain parameter. Notice that because of the equality constraints (<ref>), the problem will be infeasible if we apply static robust optimization, i.e. if we assume that x̃ and ỹ are not adjustable. In the following, we further demonstrate the importance of discrete recourse with an illustrative example, for which all data are provided in Appendix <ref>. It involves three alternative production units with different capacity ranges and production cost functions as shown in Figure <ref>, where the dashed line segments indicate the regions of uncertainty. Consider two cases with different demand ranges: Case A with d^min = min_i ∈ℐ c_i^min and d^max = ∑_i ∈ℐ (c_i^max - ĉ_i^max), and Case B with d^min = c^min_1 and d^max = c_1^max - ĉ_1^max.
  
   
  
  
  
  
  
  
  
  
  
  
  
  
  In Case A, the maximal possible demand d^max can be fulfilled by installing and operating all three units, while only Unit 2, which has the smallest capacity, should be operated to meet the minimal demand d^min. The problem is infeasible if only continuous recourse is considered since the decision of which units to operate cannot be adjusted depending on the realization of the demand. However, if the binary variables ỹ are adjustable, any demand within the given range can be met exactly.
  
  In Case B, the demand range is such that it can be covered by Unit 1 but neither Unit 2 nor Unit 3 alone. However, with both Units 2 and 3 installed and switching between these two units depending on the realization of the demand, meeting demand over the entire range is feasible. Moreover, the production costs for Units 2 and 3 are lower than for Unit 1. Hence, assuming the capital costs are the same, selecting Units 2 and 3 is a better solution than selecting Unit 1. Yet, this solution is only feasible if we have both binary and continuous recourse, i.e. if x̃ and ỹ are both adjustable. If only x̃ are adjustable, the only robust feasible solution is to choose Unit 1.
  
  Table <ref> shows the optimal (worst-case) costs for Cases A and B, each solved once allowing only continuous recourse and then with both continuous and binary recourse. The decision rules are constructed using three equidistantly generated breakpoints for each uncertain parameter in both Case A and B. One can see that Case A is only feasible if we allow both continuous and binary recourse. In Case B, considering binary in addition to continuous recourse significantly improves the optimal value while ensuring the same level of robustness.
  
  
   
  
  
  
  
  
  
  
  
  
  
   
  
  
    §.§.§ The Impact of Endogenous Uncertainty
  
  
  In this problem, we model ĉ as endogenous uncertain parameters. However, given the uncertainty set 𝒞(z), it is actually not necessary to do so. The main reason is that since the uncertain parameters are independent, ĉ_i = ĉ^max_i is always the worst case if z_i = 1. In addition, the value of ĉ_i does not affect the problem if z_i = 0. Hence, we can simply replace ĉ in (<ref>) with ĉ^max. The situation is different if the uncertain parameters are correlated. For example, consider the following budget-based uncertainty set for ĉ:
  
      𝒞^1 = {ĉ∈ℝ_+^|ℐ|: ∑_i ∈ℐĉ_i ≤τ∑_i ∈ℐĉ_i^max,  ĉ≤ĉ^max},
  
  where the total deviation of ĉ from zero across all possible units is bounded from above by τ∑_i ∈ℐĉ_i^max. In this case, for τ < 1, the worst-case value of ĉ_i may not be ĉ^max_i and cannot be easily determined a priori. 
  
  The uncertain parameter ĉ_i is endogenous because it only materializes if unit i is built. Otherwise, the parameter is physically meaningless and should therefore be irrelevant for the problem; as a result, the budget uncertainty set should change accordingly. This endogenous nature of the uncertainty is not capture in 𝒞^1. A more appropriate decision-dependent uncertainty set is
  
      𝒞^2 (z) = {ĉ∈ℝ_+^|ℐ|: ∑_i ∈ℐĉ_i ≤τ∑_i ∈ℐĉ_i^max z_i,  ĉ≤ĉ^max∘z},
  
  where ĉ_i is fixed to zero if z_i=0 and the budget only considers materialized uncertain capacities. One can see that, assuming ĉ^max > 0, 𝒞^2 (z) ⊂𝒞^1 for any z≠e; hence, using 𝒞^1 as the uncertainty set is expected to lead to overly conservative solutions.
  
  Consider Case B for which we can obtain the following analytical optimal solution:
  
      z_1 = 0,   z_2 = z_3 = 1,  ỹ_1 = 0,  x̃_1 = 0 
         {[                             ỹ_2 = 1,  x̃_2 = d,  ỹ_3 = x̃_3 = 0,                                     if  d^min≤ d ≤ c_2^max - ĉ_2;             ỹ_2 = ỹ_3 = 1,  x̃_2 = d - c_3^min,  x̃_3 = c_3^min,                  if  c_2^max - ĉ_2 < d ≤ c_3^min + c_2^max - ĉ_2; ỹ_2 = ỹ_3 = 1,  x̃_2 = c_2^max - ĉ_2,  x̃_3 = d - c_2^max + ĉ_2,                         if  c_3^min + c_2^max - ĉ_2 < d ≤ d^max. ].
  
  Then, we can determine the worst case, which depends on the choice of uncertainty set:
  
      For ĉ∈𝒞^1:          d^worst-1 = d^max,  ĉ_2^worst-1 = min{ĉ_2^max,  τ (  ĉ_1^max + ĉ_2^max + ĉ_3^max ) }. 
         For ĉ∈𝒞^2(z):          d^worst-2 = d^max,  ĉ_2^worst-2 =min{ĉ_2^max,  τ ( ĉ_2^max + ĉ_3^max) }.
  
  The optimal values (in k$), i.e. the minimum worst-case costs, for these two cases can then be computed as follows:
  
      Z^worst-1 = -15 + 4d^max + 2 ĉ_2^worst-1. 
          Z^worst-2 = -15 + 4d^max + 2 ĉ_2^worst-2.
  
  Evidently, we have Z^worst-1≥ Z^worst-2, which shows that the decision-dependent uncertainty set 𝒞^2(z) is less conservative than the fixed 𝒞^1. We can further see that Z^worst-1 and Z^worst-2 both depend on τ, which is depicted in Figure <ref>. Here, we also compare the analytical solutions with solutions obtained from solving the two-stage robust optimization problem with the proposed decision rules. The breakpoints for demand d are chosen to be c_2^max - ĉ_2^worst and c_3^min + c_2^max - ĉ_2^worst with ĉ_2^worst = ĉ_2^worst-1 or ĉ_2^worst = ĉ_2^worst-2 depending on the choice of uncertainty set. By doing so, as shown in Figure <ref>, we can recover the analytical optimal solutions.
  
   
  
  
  
    §.§.§ On the Selection of Breakpoints
  
  
  The quality of the proposed decision rules strongly depends on the choice of breakpoints. In Case B, the optimal breakpoints could be determined a priori; however, this is not generally true in more complex instances. In practice, we have to apply some heuristic to generate the breakpoints. The most intuitive one is to simply choose the number of breakpoints for each uncertain parameter and place them equidistantly inside its marginal support. However, it is recommended to utilize problem-specific features to design improved breakpoint generation procedures. For example, in this problem, according to constraint (<ref>), the demand d has to be equal to the sum of all built units' production amounts. This insight motivates a tailored method that uses all c_i^min and c_i^max that are within the range [d^min, d^max] as breakpoints for d. More generally, let p_i^j be a breakpoint for ĉ_i, choose all c_i^min∈ [d^min, d^max] and c_i^max - p_i^j ∈ [d^min, d^max] for i ∈ℐ and j=0, …, r_i to be breakpoints for d.
  
  For a randomly generated case with eight alternative production units (data provided in Appendix <ref>), we compare the two heuristic breakpoint generation methods described above. In the case of equidistant construction of breakpoints, we apply the same number of breakpoints to each of the nine uncertain parameters and examine the performance as we increase the number of breakpoints. When using the tailored method, we only apply breakpoints to d. The computational results are obtained with d^min = 43.1, d^max = 406.5 and τ = 0.5, as shown in Table <ref>. One can see that the problem is infeasible if no breakpoints are used. Similarly, it is infeasible if one breakpoint is placed at the center of the marginal support of every uncertain parameter. As the number of breakpoints increases, equidistant generation of breakpoints leads to improved solutions, albeit at higher computational cost since the model size grows with the number of breakpoints. Note, however, that the optimal value does not improve from 27 to 36 breakpoints. The same optimal value is achieved with the 15 breakpoints generated using the tailored method, where the problem was solved in 12 seconds, which is in contrast to the 172 seconds required to solve the instance with a total of 27 equidistantly placed breakpoints. 
  
  
  
  
  
   §.§ Multiperiod Production Planning
  
  
  In the second case study, we consider a multiperiod production planning problem with exogenous uncertain demands and endogenous uncertain production capacities. Endogenous, especially type-2 endogenous, uncertainty is prevalent in planning and scheduling applications as many task-related uncertainties, such as production capacity, yield, and processing time, only materialize if one decides to perform the task <cit.>.
  
  The multistage sequential decision-making process is depicted in Figure <ref>, where we apply the convention that a time period t starts at time point t-1 and ends at time point t. Before the start of the planning horizon, which is given by the set of time periods 𝒯 := {1, …, T }, we have to decide whether each unit i ∈ℐ := {1, …, I} should be upgraded such that its capacity is increased or the uncertainty associated with the capacity is reduced. This first-stage binary decision is denoted by z_i and is associated with a fixed cost γ_i. We then observe the demand in time period 1, d_1, and decide on which units to run; hence, the binary variable y_1i, which is 1 if and only if unit i operates in time period 1, depends on the realization of d_1. The production capacity of unit i in time period t is c_i^max - ĉ_ti, where ĉ_ti is an uncertain parameter. The uncertainty in the capacity of a production unit only materializes if the unit is turned on; hence, ĉ_1 are only observed after y_1 are set. Once ĉ_1 are observed, the production amounts x_1, the purchasing amount p_1, and the resulting inventory level s_1 are determined; hence, these decisions depend on the realization of d_1 and ĉ_1. As indicated in Figure <ref>, this sequential decision-making process is carried out until the end of the planning horizon. As a result, we have a multistage problem with 2T+1 stages.
  
   
  
  The multistage robust production planning problem is formulated as follows:
  
      ∑_i ∈ℐγ_i z_i + max_d^[T]∈𝒟^[T],  ĉ^[T]∈𝒞^[T](z,  y^[T])∑_t ∈𝒯∑_i ∈ℐα_i x_ti + β_i y_ti + θ_t p_t + η_t s_t  
            z_i ∈{0,1}  ∀  i ∈ℐ
          y_ti∈{0,1}  ∀  i ∈ℐ,   t ∈𝒯,  d^[t]∈𝒟^[t],  ĉ^[t-1]∈𝒞^[t-1](z,  y^[t-1]) 
         .
      [ s_t = s_t-1 + ∑_i ∈ℐx_ti + p_t - d_t;        x_ti≤ c_i^max - ĉ_ti  ∀  i ∈ℐ;          x_ti≤ c_i^max y_ti  ∀  i ∈ℐ;          x_ti≥ c_i^min y_ti  ∀  i ∈ℐ;                          s_t ≤ s^max;                      s_t,   p_t ∈ℝ_+;                   x_ti∈ℝ_+   ∀  i ∈ℐ ]}  ∀  t ∈𝒯,  d^[t]∈𝒟^[t],  ĉ^[t]∈𝒞^[t](z,  y^[t])
  
  where α, β, θ, and η are cost parameters, and y_t = y_t(d^[t],  ĉ^[t-1]), x_t = x_t(d^[t],  ĉ^[t]), p_t = p_t(d^[t],  ĉ^[t]), and s_t = s_t(d^[t],  ĉ^[t]). We consider the following uncertainty sets:
  
  
  
      𝒟^[t] = {d^[t]∈ℝ_+^t:  d^min≤d≤d^max}
         𝒞^[t](z,  y^[t]) = {ĉ^[t]∈ℝ_+^tI: 
          [         ĉ_t'i≤ĉ_t'i^max-1 z_i y_t'i + ĉ_t'i^max-2 (1-z_i) y_t'i  ∀   t' = 1, …, t,   i ∈ℐ; ∑_i ∈ℐĉ_t'i≤τ_t'∑_i ∈ℐĉ_t'i^max-1 z_i y_t'i + ĉ_t'i^max-2 (1-z_i) y_t'i  ∀   t' = 1, …, t ]},
  
  
  where 𝒟^[t] is a simple box uncertainty set, while 𝒞^[t](z,  y^[t]) is a decision-dependent budget uncertainty set. The uncertain parameter ĉ_t'i is nonzero only when y_t'i = 1. Furthermore, in the case of y_t'i = 1, the upper bound of c_t'i is c_t'i^max-1 if z_i = 1 and c_t'i^max-2 otherwise. With c_t'i^max-1 < c_t'i^max-2, an equipment upgrade increases the minimum capacity and reduces the level of uncertainty in the capacity. Note that there are bilinear terms involving two binary variables in (<ref>), i.e. z_i y_t'i, which can be easily linearized.
  
  In the following, we consider multiple instances of a production planning problem with three units. The maximum inventory s^max is set to 5; all the other data are provided in Appendix <ref>. We construct decision rules using one breakpoint for each uncertain parameter and such that recourse variables only depend on the uncertain parameters from the current time period.
  
  
  
    §.§.§ The Benefit of Discrete Recourse
  
  
  We first consider the case with two time periods and investigate how the solution depends on the equipment upgrade costs γ. We set γ_1 = γ̅, γ_2 = 1.5 γ̅, and γ_3 = 3 γ̅. The results from solving multiple instances with different γ̅ are shown in Figure <ref>. The pie chart associated with each instance indicates which units are being upgraded; colored fill means that we decide to upgrade the corresponding unit. One can see that as γ̅ increases, it becomes less worthwhile to invest in equipment upgrades, up to a point where we leave all units unchanged. Figure <ref> also shows the benefit of discrete recourse as we solve each instance once considering only continuous recourse and another time with both continuous and binary recourse. Comparing the optimal values, we see that a cost reduction of more than 35 % can be achieved if binary recourse is considered in addition to continuous recourse. Note that in both cases, continuous recourse variables can follow discontinuous piecewise linear decision rules.
  
   
  
  We also solve the problem with different numbers of time periods. Table  <ref> shows the results for T equal to 2, 3, 4, and 5 with γ̅ = 100. Again, one can see that the cost is substantially reduced if in addition to continuous recourse, also binary recourse is considered. Moreover, the cost reduction increases with the number of time periods.
  
  
  
  
  
    §.§.§ Discussion on Computational Performance
  
  
  The computational results in Table <ref> also show that the benefit of binary recourse comes at the cost of significantly greater computational complexity. For example, for T=5, the model with only continuous recourse solves in 1,560 s, while the computation time for the model with both continuous and binary recourse is about two orders of magnitude longer. One reason for the higher computational complexity is obviously the increased model size. However, it turns out that the by far larger contributing factor is the “looseness” of the MILP formulation. Recall that the incorporation of binary recourse variables that affect the uncertainty set results in a formulation that involves bilinear terms. Each of these bilinear terms consists of a binary variable and a continuous variable representing a “dual” variable that stems from the reformulation. The linearization of these bilinear terms involves the lower bounds, which are zero, and the upper bounds of the dual variables. Generally, we cannot find tight upper bounds on these variables a priori such that very large values have to be chosen, which leads to a very weak LP relaxation of the MILP. This explanation is consistent with our observation, which is that the actual optimal solution is usually found fairly quickly but the lower bound only improves very slowly.
  
  We further confirm our suspicion with a small experiment. After solving each of the instances with continuous and binary recourse to optimality, we update the bounds on the dual variables based on the optimal solution and re-solve the problem. Let Λ be the value of a dual variable at the optimal solution, then the tightest update that we can apply is to set the upper bound M to Λ. In addition, we consider two other update rules: M = 2 Λ and M = 2 Λ + 0.01. The computation times for re-solving the four instances with the three different bound update rules are shown in Table <ref>. In all cases, the computation times are drastically shorter than the ones required to solve the problems with the original bounds. Moreover, one can see that the solution times in case of the third update rule are significantly longer although the resulting bounds are only minimally larger. The reason is that at the optimal solution, most dual variables are zero; hence, the first two update rules fix all these variables to zero, which makes the problem considerably easier to solve. Note that these updated bounds do not result in a rigorous reformulation of the problem although the same optimal value is achieved. The sole purpose of this experiment is to demonstrate the impact of these bounds on the computation time.
  
  
  
  
  
  
  § CONCLUSIONS
  
  
  
  In this work, we addressed multistage robust optimization with mixed-integer recourse and endogenous uncertainty, considering polyhedral uncertainty sets that are affected by binary variables. Applying a decision rule approach, which relies on the concept of lifted uncertainty sets, we derived tractable reformulations for the two- and multistage cases. The proposed framework has significant modeling flexibility as it can incorporate uncertainty sets affected by recourse decisions, binary recourse, and continuous recourse variables that follow discontinuous piecewise linear decision rules. The main advantage of appropriately modeling endogenous uncertainty and mixed-integer recourse is manifested in the significant reduction in solution conservatism, as demonstrated in our computational experiments.
  
  Although the proposed reformulations enjoy the favorable tractability properties of robust optimization, the computational case studies also show that they tend to result in large and rather loose MILP formulations. Hence, future work will focus on the development of solution strategies that improve the computational performance.
  
  
  
  
  § ACKNOWLEDGMENTS
  
  We gratefully acknowledge financial support from the National Key Research and Development Program of China (No. 2019YFB1705004), Science Fund for Creative Research Groups of NSFC (No. 61621002), and China Scholarship Council (CSC) (No. 201906320317).
  
  
  
  
  
  § RESTRICTED DECISION RULES IN THE MULTISTAGE CASE
  
  
  
  Consider decision rules for recourse variables in stage t that only depend on uncertain parameters observed in the previous Δ t and current stages:
  
  
  
      x_t = X_11^t ξ̅_11 + X_11^t ξ̂_11 + ∑_t'= max{2,   t - Δ t }^t ∑_i=1^K_t'X_t'i^t ξ̅_t'i + X_t'i^t ξ̂_t'i
         y_t = Y_11^t  ξ̂_11 + ∑_t'= max{2,   t - Δ t }^t ∑_i=1^K_t'Y_t'i^t  ξ̂_t'i,
  
  
  where ξ_11 is still included in the decision rule as it accounts for the constant term. Following the same procedure presented in Subsection <ref>, this results in the following reformulation of (<ref>) considering (<ref>):
  
  
  
      f_t'i^t = A_t'i^t x_1 + D_t'i^t y_1 - b_t'i^t   ∀  t= 2, …, T,   t' = 1, …, t,   i=1, …, K_t'
         Φ_t U_1^t y_1 + ∑_t'=1^t∑_i=1^K_t'δ_t'i^t ≤0  ∀  t=2,…, T 
         (f_11^t  - Φ_tw_11^[t]) v_11 - δ_11^t 
          + ∑_t” = 2^t ( A_t”^t X_11^t”v̅_11 + 
          
          [ A_t”^t X_11^t” +  (D_t”^t  + Φ_t U_t”^t ) ( Ẏ_11^t” -Ÿ_11^t”)  ] v̂_11) ≤0
                                                                 ∀  t=2,…,T 
         (f_t'i^t  - Φ_tw_t'i^[t]) v_t'i - δ_t'i^t 
          + ∑_t” = t'^min{t,   t'+Δ t }( A_t”^t X_t'i^t”v̅_t'i + 
          
          [ A_t”^t X_t'i^t” +  (D_t”^t  + Φ_t U_t”^t ) ( Ẏ_t'i^t” -Ÿ_t'i^t”)  ] v̂_t'i) ≤0
                     ∀  t=2,…,T,   t' = 2, …, t,   i = 1, …, K_t',  ( v_t'i,  v̅_t'i,  v̂_t'i) ∈𝒱_t'i.
  
  
  The integrality constraints can be reformulated in the same fashion.
  
  
  
  § DATA FOR CASE STUDIES
  
  
  
  
  
  
  
  
  
  
  
  
  newapa
  
  
  
  