Stage-t scenario dominance for risk-averse multi-stage stochastic mixed-integer programs

This paper presents a new and general approach, named “Stage-t Scenario Dominance,” to solve the risk-averse multi-stage stochastic mixed-integer programs (M-SMIPs). Given a monotonic objective function, our method derives a partial ordering of scenarios by pairwise comparing the realization of uncertain parameters at each time stage under each scenario. Specifically, we derive bounds and implications from the “Stage-t Scenario Dominance” by using the partial ordering of scenarios and solving a subset of individual scenario sub-problems up to stage t. Using these inferences, we generate new cutting planes to tackle the computational difficulty of risk-averse M-SMIPs. We also derive results on the minimum number of scenario-dominance relations generated. We demonstrate the use of this methodology on a stochastic version of the mean-conditional value-at-risk (CVaR) dynamic knapsack problem. Our computational experiments address those instances that have uncertainty, which correspond to the objective, left-hand side, and right-hand side parameters. Computational results show that our “scenario dominance"-based method can reduce the solution time for mean-risk, stochastic, multi-stage, and multi-dimensional knapsack problems with both integer and continuous variables, whose structure is similar to the mean-risk M-SMIPs, with varying risk characteristics by one-to-two orders of magnitude for varying numbers of random variables in each stage. Computational results also demonstrate that strong dominance cuts perform well for those instances with ten random variables in each stage, and ninety random variables in total. The proposed scenario dominance framework is general and can be applied to a wide range of risk-averse and risk-neutral M-SMIP problems.


Introduction
Risk is a fundamental issue arising in finance, insurance, project management, and many other areas. The analysis of risk in stochastic programming has recently been quite popular and raised many research questions, from model formulation to algorithmic approaches to tackle the problem difficulty, see, for instance, Schultz (2005), Schultz and Tiedemann (2006), Eichhorn and Römisch (2005), Ruszczyński (2013), and Yin and Büyüktahtakin (2021), among others. This paper addresses the computational difficulty of solving multi-stage stochastic mixed-integer programming problems (M-SMIPs) involving risk-averse objectives.
Stochastic mixed-integer programs can be formulated as an extensive form of a mixedinteger program (MIP), where a discrete stochastic process is represented by a finite number of realizations, i.e., scenarios. The size of the MIP grows exponentially in the number of decision stages and scenarios, leading to large-scale optimization formulations. Incorporating the riskaversion criterion in addition to the expectation in the objective function further complicates the problem.
In this paper, we propose a new and general method of "stage-t scenario dominance" to effectively solve the risk-averse M-SMIPs. Our main idea is to derive cutting planes and bounds using the dominance characteristics of scenarios defined by a partial ordering and the solution of a new scenario sub-problem, which includes all the non-anticipativity and original constraints. Specifically, we solve the problem for a number of specific scenario sub-problems and use this information and the partial ordering of scenarios to generate cutting planes so that the general problem can be solved to optimality or close-to-optimality by a cut-and-branch algorithm. We demonstrate the computational benefit from using the proposed approaches on the risk-averse stochastic dynamic knapsack problem.

Literature review and paper contributions
Solution algorithms proposed for risk-averse multi-stage stochastic programs with integer variables are quite limited due to the non-convexity arisen by integrality. Most approaches proposed for risk-averse multi-stage stochastic optimization problems are extensions of the solution techniques proposed for their risk-neutral counterparts. Schultz and Tiedemann (2006) develop a Lagrangian decomposition algorithm that is based on the relaxation of nonanticipativity constraints to solve the scenario-based formulation of two-stage mixed-integer stochastic programming involving the conditional value-at-risk (CVaR). Heinze and Schultz (2008) later extended the scenario-decomposition algorithm of Schultz and Tiedemann (2006) combined with a branch and bound scheme to the case of multi-stage stochastic integer programs with CVaR risk objective. Yu et al. (2018) employ the Stochastic Dual Dynamic Integer Programming (SDDiP) approach and Lagrangian cuts of Zou et al. (2019) for solving the risk-averse multi-stage stochastic facility location problem using expected conditional risk measures (ECRMs). Dentcheva and Martinez (2012) introduce Second-order Stochastic Dominance (SSD) constraints and explore the mathematical properties in two-stage stochastic problems. Escudero et al. (2016) study the first-and second-order stochastic dominance constraints induced by mixed-integer recourse in risk-averse multi-stage stochastic programs. Our scenario dominance approach is fundamentally different than the stochastic dominance constraints introduced in the literature (Müller and Stoyan 2002) in that we pairwise compare the scenarios based on the problem parameters and derive implications using the solution of scenario sub-problems, while in stochastic dominance random variables are compared based on their probability distribution function.
There has also been a stream of research studying bounds and approximation for riskaverse multi-stage stochastic mixed-integer programming problems (Maggioni and Pflug 2019;Mahmutogulları et al. 2018;Sandıkçı and Özaltın 2017;Bakir et al. 2020;Maggioni et al. 2014;Maggioni and Pflug 2016;Sandıkçı et al. 2013). Our focus in this paper is different from those by providing cutting plane algorithms based on the scenario dominance concept to tackle the computational difficulty of solving risk-averse multi-stage stochastic programs. Particularly, our scenario dominance-based cutting planes provide bounds on a single-scenario objective rather than the full objective function as considered in previous studies on bounds and approximations [e.g., see Maggioni and Pflug (2019) and Mahmutogulları et al. (2018)]. Furthermore, our cuts are driven based on a new scenario sub-problem, which is different from the classical scenario sub-problem that have been widely used in prior studies [see, e.g., Madansky (1960), ]. Ruszczyński (2002) reformulates the two-stage stochastic programming problems with probabilistic constraints as large-scale mixed integer programming problems with knapsack constraints. The author uses the partial ordering of the scenarios to define precedence relationships between the variables of the knapsack problem and then defines the induced knapsack covers (Park and Park 1997;Boyd 1993) for the precedence-constrained knapsack polyhedra. Our definition of scenario dominance is similar to the ordering of scenarios that was discussed in Ruszczyński (2002); but the way we use the scenario ordering concept is fundamentally different than the approach in Ruszczyński (2002) because we compare the single scenario objectives of two scenarios based on a newly-defined scenario sub-problem to derive cuts for the risk-averse M-SMIPs, while Ruszczyński (2002) uses the ordering of scenarios to define induced knapsack covers. Furthermore, different than former work, our study, for the first time, defines the notion of "stage-t scenario dominance" and utilizes it for solving the risk-averse M-SMIPs more efficiently with tractable computational times.
Due to the limited number of available solution algorithms for the risk-averse multi-stage stochastic programs with integer variables, the main contribution of this paper is to provide a new approach named "scenario dominance" to solving risk-averse M-SMIPs. Specifically, we provide the following contributions: (i) presenting a new approach, named as "stage-t scenario dominance" that compares scenarios based on the realization of uncertain parameters under each scenario up to a certain time period t, derives a partial ordering of scenarios, and studies the implications of this ordering; (ii) defining a new scenario sub-problem, which involves all constraints of the original problem, including the non-anticipativity constraints, as opposed to the classical scenario sub-problem, which ignores the non-anticipativity relations among scenarios, as commonly used in the literature (Birge and Louveaux 2011;Lulli and Sen 2004); (iii) deriving and proving the scenario-based formulation for the mean-CVaR multi-stage stochastic mixed-integer program; (iv) demonstrating the use of "scenario dominance" theory to derive bounds and cutting planes for the risk-averse mean-CVaR M-SMIPs based on the solution of a new scenario sub-problem; (v) deriving results on the minimum number of stage-t scenario-dominance relations in multi-stage scenario trees, and thus cutting planes generated; (vi) presenting the computational benefit using the proposed "scenario dominance" cuts and bounds to solve the risk-averse dynamic stochastic knapsack problem by conducting computational experiments; and (vii) showing that "scenario dominance" cuts can significantly reduce the solution time for large-scale mean-risk stochastic multi-stage multi-dimensional knapsack problems with both integer and continuous variables, whose structure is similar to the mean-risk M-SMIPs. The proposed framework is general and can be applied to a wide range of risk-averse M-SMIP problems.

Risk-neutral multi-stage stochastic MIPs
We consider a finite number of stages t ∈ T := {1, . . . , T }. Let ξ 1 be known, or deterministic, and (ξ 2 , . . . , ξ T ) denotes a sequence of uncertain parameters where ξ [t] := {ξ 1 , ξ 2 , . . . , ξ t } is the information observed by stage t. Let n t , q t , m t , and k t represent the number of decision variables, the number of integer variables, the number of constraints, and the number of random variables, respectively, at time t ∈ T . Let R denote the set of real numbers. We assume that ξ 2 , . . . , ξ T is a discrete-time stochastic process with a finite probability space ( t , F t , P t ), where the support, t ∈ R k t , represents a finite sample space for t = 2, . . . , T , and = T , F t is the set of all events defined on t , and P t is the corresponding probability distribution. Let R + and Z + denote the set of positive real numbers and positive integers, respectively. Let The sequence of decisions and the realizations of uncertain parameters are represented as We define the multi-stage stochastic mixed-integer program (M-SMIP) to be of the nested form below: where c 1 ∈ R n 1 , A 1 ∈ R m 1 ×n 1 , and b 1 ∈ R m 1 are known; and the uncertain parameters that realize as the stochastic process ξ evolves are given by c t (ξ [t] ) ∈ R n t , A t (ξ [t] ) ∈ R m t ×n t , H t (ξ [t] ) ∈ R m t ×n t , and b t (ξ [t] ) ∈ R m t ; and E ξ t+1 |ξ [t] [•] is the expectation with respect to ξ t+1 conditioned on the realization ξ [t] , for t = 2, . . . , T . Note that for T = 2, the M-SMIP (1) is a two-stage SMIP.
Numerical optimization of multi-stage stochastic programs usually requires discrete probability distributions since otherwise, the multi-variate integrals involved become intractable (Heinze and Schultz 2008). So we consider a finite number of realizations of a discrete stochastic process ξ = (ξ 1 , . . . , ξ T ), which has a finite support = ξ 1 , . . . , ξ¯ such that | | =¯ for some positive integer¯ . Let ξ ω 1 be equal to the known parameter ξ 1 for each ω ∈ := 1, 2, . . . ,¯ . A realization of a sequence of random parameters in stages 1, . . . , T , ξ ω = ξ ω 1 , ξ ω 2 , ξ ω 3 , . . . , ξ ω T ∈ , is then referred as a scenario indexed by ω ∈ . Each scenario ξ ω ∈ is associated with probability p ω , which is computed as the product of the conditional probabilities of all nodes that belong to the scenario path ξ ω , such that ω∈ p ω = 1. The realizations of a scenario ξ ω up to stage t is denoted by ξ ω [t] := ξ ω 1 , ξ ω 2 , ξ ω 3 , . . . , ξ ω t . When there is a finite number of scenarios, the M-SMIP (1) can be equivalently formulated as a large-scale deterministic program by replicating the decision variable in stage t, x t , for each scenario ξ ω ∈ . Denoting x ω t as the decision variable for stage t = 1, . . . , T , the M-SMIP given in (1) can be equivalently formulated as the extensive form of the multi-stage stochastic MIP (E-MSMIP) (2a): where the random parameters become known with The non-anticipativity constraints (2d) ensure that x ω t and x ω t must attain the same values as long as the paths corresponding to scenarios ξ ω and ξ ω coincide up to stage t (i.e., ξ ω [t] = ξ ω [t] ). Those non-anticipativity constraints ensure the real-time implementability of solutions.

Time consistency in multi-stage stochastic MIPs
The classical stochastic programming approach given in E-MSMIP (2a) may perform poorly if certain realizations of the stochastic parameters lead to an excessively high cost to accept. Thus, a mean-risk optimization formulation, which considers a risk-averse measure in addition to the expectation in the objective function, performs better than the expectation objective when the decision-making is at high-stake. In this paper, we focus on an asymmetric risk measure called Conditional Value-at-Risk (CVaR), which has been widely used as a risk measure in decision-making problems under uncertainty because it is a coherent measure of risk (see, for instance, see Rockafellar et al. (2000), Pflug (2000), and Ogryczak and Ruszczynski (2002) for more detailed discussions).
The incorporation of risk measures, such as CVaR, into the two-stage risk-averse stochastic program is straightforward (Ahmed 2006;Schultz and Tiedemann 2006;Schultz 2005;Fábián 2008). However, the definition of risk aversion is not apparent in the multi-stage setting. Risk measures can be applied at every stage additively or to the complete scenario path, only defining risk at the end of the planning horizon or can be represented in a nested form. One of the most important considerations in modeling the risk-averse multi-stage stochastic programs is time consistency. Risk-neutral or risk-averse two-stage stochastic programs are time consistent, while time consistency for multi-stage risk-averse stochastic programs is not guaranteed. Based on the definition of Homem-de Mello and Pagnoncelli (2016), time consistency means that if you solve a multi-stage stochastic program today and find solutions for each node of a tree, you should find the same solutions if you re-solve the problem tomorrow given what was observed and decided today. Various researchers have studied how to address the time consistency issue in multi-stage risk-averse stochastic programs. For example, Shapiro et al. (2013) propose a nested conditional risk measure for multi-stage optimization problems and prove its time consistency. However, their nested conditional risk measure is formulated as Bellman equations (Bellman 1957;Büyüktahtakın 2011), which cannot be solved in explicit form. Homem-de Mello and Pagnoncelli (2016) address this drawback by proposing a class of expected conditional risk measures (ECRMs) that is proven to be time consistent.

Mean-CVaR multi-stage stochastic MIP
Next, we analyze the incorporation of the multi-period CVaR in the objective function of the E-MSMIP (2a) and present the mean-CVaR multi-stage SMIP formulation, similar to the ECRM approach of Homem-de Mello and Pagnoncelli (2016), by casting the resulting riskaverse problem with ECRM as an extensive form that includes some additional variables and risk linearization constraints. Our definition of the multi-period CVaR below is also consistent with other multi-period nested CVaR models in the literature (see, e.g., Heinze and Schultz 2008;Dupačová and Kozmík 2015;Alonso-Ayuso et al. 2018). (Eichhorn and Römisch 2005;Heinze and Schultz 2008)] The Conditional Value-at-Risk of a random variable c t (ξ [t] ) x t at confidence level α t ∈ [0, 1) at each time stage t ∈ T \{1} can be expressed as the optimal value of the following optimization problem:

Definition 1.1 [Multi-period CVaR
Proposition 1 Consider the following mean-CVaR problem: where X t (x [t−1] , ξ [t] ) denotes the feasibility set at stage t with decisions up to stage t − 1 (x [t−1] ) and the realization of the uncertain parameters up to stage t (ξ [t] ), and λ ∈ R + is a trade-off coefficient that represents the relative weight of the risk objective (C V a R α t [•]) with respect to the expectation function (E [•]). Given the above notations and definitions, we can equivalently formulate the mean-CVaR problem (4), as the following extensive mean-CVaR multi-stage stochastic MIP (CVaR M-SMIP) problem (5a): The proof of Proposition 1 is given in Appendix A.

Scenario-dominance
In this section, we introduce a new concept of "scenario dominance" and "dominance sets," considering the partial order relations of scenarios for the CVaR M-SMIP problem (5a). The partial ordering of scenarios is defined by comparing the scenario realizations of uncertain parameters. Using the partial ordering of scenarios and its implications based on the solution of a scenario sub-problem, we derive new bounds and cutting planes to improve the computational solvability of mean-CVaR M-SMIP formulation (5a).

Defining the scenario-dominance
Similar to the partial ordering of elements of a set, we define the partial ordering of scenarios in the scenario set . Because a scenario is defined as the realizations of a random variable over multiple stages t ∈ T , the partial ordering of the scenarios is obtained by the pointwise comparison of scenario realizations at each stage t. Thus, in order to obtain a partial ordering of two scenarios ξ ω ∈ and ξ ω ∈ , the realization of scenario ξ ω ∈ at each time period t ∈ T , ξ ω 1 , ξ ω 2 , ξ ω 3 , . . . , ξ ω T will be pairwise compared with the realization of scenario ξ ω ∈ for each t, ξ ω 1 , ξ ω 2 , ξ ω 3 , . . . , ξ ω T . Note that ξ ω 1 = ξ ω 1 because ξ 1 is deterministic. This comparison leads to the following definition for a partial order relation of the scenarios.
Using the definition of the "scenario dominance" above, we define the scenario-dominance sets as follows: . . , 12. The + 3,ξ 5 represents the index set of scenarios that are stage-3 dominated by ξ 5 , − 3,ξ 5 represents the index set of scenarios that stage-3 dominate ξ 5 , and N 3,ξ 5 represents scenarios neither stage-3 dominate nor are stage-3 dominated by ξ 5 Definition 2.3 Stage-t Scenario-dominance sets + t,ξ ω is the index set of scenarios, which are stage-t dominated by scenario ξ ω ∈ , − t,ξ ω is the index set of scenarios which stage-t dominate scenario ξ ω ∈ , and N t,ξ ω is the index set of scenarios, which neither stage-t dominate nor are stage-t dominated by ξ ω ∈ for t ∈ T as defined below: Example 1 Figure 1 presents an example realization of the objective cost parameter, c t , in two stages resulting in 12 different scenario realizations. To simplify the presentation, we omit c 1 in a scenario description and scenario-dominance comparisons since it takes the same value under all scenario realizations. Thus, here a scenario ξ i represents a realization of c t in stages 2 and 3, [c i 2 , c i 3 ]. In the example, the scenario ξ 5 realization is [49,49]. The index set of scenarios, which are stage-3 dominated by ξ 5 is defined as + 3,ξ 5 = {1, 2, 3}, while the index set of scenarios, which stage-3 dominate scenario ξ 5 is given by − 3,ξ 5 = {5, 9, 10, 11, 12}. All other scenarios that neither dominate nor are dominated by ξ 5 is given by N ξ 5 = {4, 6, 7, 8}.

The number of scenario-dominance relations
The number of scenario-dominance relations depends on a number of factors, including the structure of the scenario tree, how the random variables are generated, and the number of scenarios. Different methods have been proposed for generating a multi-period scenario tree (Heitsch and Römisch 2009;Pflug and Pichler 2015;Dupačová et al. 2000). Here we focus on the two most common types of multi-stage scenario trees based on the type of random variable generation, where random variables are completely independent and random variables are jointly distributed but stage-wise independent.
The number of scenario dominance cutting planes described in Sect. 2.4 relies on the number of scenario-dominance relations. Thus, next, we discuss the relations of two scenarios that share the same path up to a certain time period and exploit this result to compute the minimum number of scenario-dominance relations based on how the random variables are generated in a multi-stage scenario tree.

Proposition 2 Let ξ ω ∈
and ξ ω ∈ be two scenarios such that ξ ω Proof It follows from non-anticipativity constraints (5e) that are ensured in a multi-stage stochastic program and the definition of the stage-t scenario dominance.
By Proposition 2, the scenarios that share the same path up to some stage t dominate each other in a scenario tree where there is at least one random variable at each stage. Thus, the minimum number of scenario-dominance relations could be defined by considering those dominance relations between scenarios that share a common history, motivating the following propositions.

Proposition 3 Let (ξ 2 , . . . , ξ T ) be a sequence of independently identically distributed (i.i.d) random vectors of k t dimensions each. Each univariate component of a vector ξ t is jointly
distributed. In other words, components of ξ t are positively or negatively correlated. In a scenario tree with l branches at each node in each stage, the minimum total number of scenario-dominance relations, W , is equal to: Proof In this scenario tree, we have l T −1 scenarios. At stage t = 1, we have 1 node, and thus l T −1 /l 0 scenarios share the same node at stage 1. At stage t = 2, we have l nodes, and thus l T −1 /l 1 scenarios share the same node at stage 2. At stage t = 3, we have l 2 nodes, and thus l T −1 /l 2 scenarios share the same node at stage 3. Thus, at stage t = t, we have l t−1 nodes, and thus l T −1 /l (t−1) scenarios share the same node at stage t. Based on Proposition 2, for those scenarios that share the same path up to stage t dominate each other, there are l T −1 /l t−1 2 domination relations. Since l t−1 nodes are shared at stage t, there are at least a total of l t−1 l T −1 /l t−1 2 stage-t scenario-dominance relations. This results in a total number of T t=1 l (2T −t−1) scenario-dominance relations over all stages at the minimum. Proposition 4 Let (ξ 2 , . . . , ξ T ) be a sequence of independently identically distributed (i.i.d) random vectors of k t = k dimensions each. Each univariate component of a vector ξ t is independently distributed. In a scenario tree with l branches at each node in each stage, the minimum total number of scenario-dominance relations, W , is equal to: Proof In this scenario tree, we have l (T −1)k scenarios. At stage t = t, we have l (t−1)k nodes, and thus l (T −1)k t /l (t−1)k scenarios share the same node at stage t, resulting in a total of l (t−1)k l (T −1)k /l (t−1)k 2 domination relations at stage t. This results in at least a total number of T t=1 l (2T −t−1)k scenario-dominance relations over all stages.
Example 2 Consider a CVaR-SMKP instance in which T = 4, l = 2, and k t = 2 for t = 1, . . . , T . In the case where each univariate component of the random vector ξ t are jointly distributed, we have¯ = 8 scenarios and a minimum number of 64, 32, 16, and 8 scenario-dominance relations in stages one, two, three, and four, respectively, resulting in a total of at least W = 120 scenario-dominance relations. In the case where each univariate component of the random vector ξ t is independently distributed, we have¯ = 64 scenarios and a total of at least W = 5440 scenario-dominance relations at the minimum.

Scenario sub-problem and bounds
In this section, we define a new scenario sub-problem, which optimizes variables related to a subset of scenarios and assigns feasible solutions for all the other variables. This scenario sub-problem helps us derive new lower and upper bounds on the mean-CVaR multi-stage stochastic MIP problem. As a main contribution of the paper, using the solution of the scenario sub-problem and the scenario-dominance relations described in Sect. 2.1, we derive new cutting planes to chop off non-optimal points and find solutions faster to the highly-complex mean-CVaR multi-stage stochastic MIPs. To introduce the definitions and propositions below, let us define the feasible set for the original problem. The set of feasible solutions of the CVaR M-SMIP (5a) is denoted by X as given below: Definition 2.4 (Scenario group sub-problem) Let S ⊆ . The scenario group sub-problem involving each ω ∈ S, (P S ), is formulated as follows: Note that the problem (P S ) (8) is equivalent to the original problem P (5a), if S = .
Single scenario sub-problem The scenario group sub-problem for the set consisting of the single scenario, ξ ω , is denoted as the scenario-ξ ω problem (P {ω} ) and formulated as follows: Remark 1 The scenario-ξ ω problem (P {ω} ) is an MIP, whose feasible region comprises all the variables and all the constraints of the original problem (5a) while its objective function includes only one specific scenario with index ω. Consequently, the P {ω} leads to a complete solution (x, η, v), which is feasible for (5b)-(5g) but optimizes only for one scenario ω.
In the following proposition, we assume that the objective function for each scenario is non- implying that c ω t ≥ 0 for each ω ∈ and t ∈ T . Note that this assumption is not necessary when S = , since Z = Z .

Proposition 5 Let Z be the optimal objective value of the original problem P (5a). Then, P S is a relaxation of P, that is Z ≥ Z S ∀S ⊆ .
Proof Define a new variable H ∈ R. The following problem is equivalent to the original problem P (5a): Similarly, the scenario group sub-problem P S (8) can be equivalently written as in the following problem: Clearly,X ⊆X , which proves the proposition.

Definition 2.5 Stage-t single scenario sub-problem
The single scenario problem, including the variables and constraints up to stage t, for t ∈ T , is denoted as the stage-t scenario-ξ ω problem (P {ω} t ) and formulated as follows: Definition 2.6 (Relaxed scenario group sub-problem) The relaxed scenario group subproblem involving each ω ∈ S, (P S R ), is formulated as follows: where which has been widely studied in the literature [see, e.g., Madansky (1960), Ahmed (2013), and Sandıkçı and Özaltın (2017)] because in the P {ω} R the objective function and each constraint of the problem (5a) are defined only for a specific scenario with scenario index ω and non-anticipativity constraints (5e) are relaxed, while all variables and constraints are kept in the P {ω} . Clearly, the P {ω} R is a relaxation of P {ω} as stated in Proposition 6. Using the solution of the single-scenario sub-problems, we provide lower and upper bounds, which could then be used to strengthen the linear programming relaxation of the original problem, as demonstrated under the computational results under Sect. 4.5.

Proposition 8 (Lower bound) Suppose that the objective function for each scenario is non-
negative. Then Z S and ω∈S Z {ω} provides a lower bound on the optimal objective function value of the original problem P (5a), Z , as in the following inequality: Proof The first inequality follows by Proposition 5, while the second inequality is the generalization of the superadditivity result in Proposition 7. (13) is tightest when S = . In this case, the non-negativity assumption for each scenario's objective function value is not required.

Upper bound
Proposition 9 (Upper bound) Letẋ ω be optimal solutions to the scenario-ξ ω problem P {ω} ; and letZ ω be the objective function value of P at solutionẋ ω . Then,

Scenario dominance cuts
In this section, we derive new cutting planes based on the "scenario dominance" concept defined in Sect. 2.1. We first define Lemmas 1 and 2 to motivate the main result on the scenario dominance cutting planes (Theorem 1). Then we present a corollary (Corollary 1.1) based on Theorem 1 to generalize the scenario dominance cuts to multiple time stages and further improve their computational benefits.
In this section, we present results regarding the general uncertainty, that is, the uncertainty in the right-hand side b ω t , left-hand side (A ω t , H ω t ), and the objective coefficients c ω t . For notational brevity, we denote P {ω} as P ω and Z {ω} as Z ω for the rest of the paper. We give the following two definitions that are commonly used in the results presented in the rest of the paper.
Definition 2.8 Let (ẋ ω ,η ω ,v ω ) be a feasible solution for the scenario-ξ ω problem (9) and Z ω be the corresponding objective function value, i.e., Lemma 1 below states that the portion of the optimal objective function value for the original problem (5a) that corresponds to the scenario-ξ ω is bounded below by the optimal objective value of the scenario-ξ ω problem.
Lemma 1 Let (x * , η * , v * ) be the optimal solution to the original problem (5a) and (ẋ ω ,η ω ,v ω ) be the optimal solution of the scenario-ξ ω problem (9). Then, the portion of the optimal solution for the original problem (5a) that corresponds to scenario-ξ ω , (x ω * , η ω * , v ω * ), satisfies the following inequality: Proof This is a proof by contradiction. For some ω ∈ and any (x ω * , η ω * , v ω * ) assume that Then define a new solutionx to P as follows: By (16), the new solution x ω t ,ή ω t ,v ω t is feasible and reduces the value of Z ω by (16), contradicting to the optimality of (ẋ ω ,η ω ,v ω ) for the scenario-ξ ω problem.
Lemma 2 below defines the relations of the optimal objective function value of the subproblems that correspond to two scenarios that have a dominance relation.
Lemma 2 Let ξ ω be a scenario such that ω ∈ − T ,ξ ω ⊆ for a given scenario ξ ω . Let (ẍ ω ,η ω ,v ω ) and (ẋ ω ,η ω ,v ω ) be the optimal solution of the scenario-ξ ω and the scenarioξ ω problems, Z ω and Z ω be the corresponding objective values, respectively. Then the following inequality holds: Proof Add the following equalities to both scenario-ξ ω and scenario-ξ ω problems, P ω and P ω : Define the new problems asP ω andP ω , respectively. Define the following sets of scenarios: Note that the set 1 only includes scenariosω = {ω, ω } that share a common history with scenario ω up to stage t > 1 and = 1 ∪ 2 ∪ {ω} ∪ ω .
Using those two lemmas above, we show the main theorem below on the scenario dominance cut that provides a lower bound on the portion of the optimal objective value for the original problem (5a) that corresponds to the scenario-ξ ω problem. This lower bound is obtained by solving the single-scenario problems that are dominated by the scenario ξ ω . (5a) is observed in the objective function c ω t , in the right-hand side coefficient b ω t , and the left-hand side coefficient (A ω t , H ω t ). Let ξ ω and ξ ω be two scenarios such that ξ ω ξ ω , ω, ω ∈ and

Theorem 1 Scenario dominance cuts Assume that uncertainty in problem
and c ω t ≥ c ω t for all t ∈ T . Then, the portion of the optimal solution for the original problem (5a) that corresponds to scenario-ξ ω , (x ω * , η ω * , v ω * ), satisfies the following inequality: Proof The first inequality below follows from Lemma 1, while the second follows from Lemma 2.
In the corollary following Theorem 1, we generalize the cuts to multiple stages by solving each stage-t problem for t ∈ T .
be the optimal solution for the stage-t scenario-ξ ω problem (10) and Z ω t be the corresponding optimal objective value such that Corollary 1.1 Stage-t Scenario dominance cuts Let ξ ω and ξ ω be two scenarios l, ω ∈ such that ξ ω t ξ ω t for t ∈T = 1, . . . ,t ⊆ T , wheret ∈ T . Then, the portion of the optimal solution for the original problem (5a) that corresponds to scenario-ξ ω up to stage t, (x ω t * , η ω t * , v ω t * ), satisfies the following inequality: Proof The proof follows from the application of Theorem 1 to each stage t ∈ T .

Remark 4
The stage-t scenario dominance cuts (27) are generalizations of the scenario dominance cuts (26). Both (26) and (27) are satisfied by all optimal solutions to the original problem (5a); however they are not valid inequalities for the feasible region of the problem (5a), i.e., they both may cut-off a non-optimal feasible solution, while preserving the optimal solution. This type of inequalities are not uncommon (see, e.g., Bender's optimality cuts Benders (2005), supervalid inequalities of Israeli and Wood (2002), and the pre-processing cuts presented in Kibis et al. (2020)), and they are useful because they help to shrink the feasible region while preserving the optimal solution.
Remark 5 For problem instances, P, where it is easy to find a feasible solution but hard to find the optimal one, we expect that the single scenario sub-problem, P ω , could be solved much faster than P because the P ω is a relaxation of P, as shown in Proposition 5. Furthermore, in the P ω , we optimize the value of the variables related to a single scenario ω, while we only need to assign a feasible solution to all other variables in order to satisfy the constraint set, X . This observation is also justified by the computational results for the dynamic knapsack problems in Sect. 4. On the other hand, if finding a feasible solution is hard, the single scenario sub-problems could take a long time to solve. In such cases, we suggest solving the relaxed single scenario sub-problem P ω R rather than the P ω . Because the P ω R is merely a relaxation of the P ω , including only the constraints and variables of a single scenario, the lower bound (13), scenario dominance cuts (26), and stage-t scenario dominance cuts (27) would still hold when P ω R is used to generate those inequalities.

Risk-averse multi-stage stochastic knapsack
We consider a class of stochastic mean-CVaR multi-stage multi-dimensional mixed 0-1 knapsack problems (SMKPs) as a benchmark set for our computations. The classical deterministic, single-dimensional, and risk-neutral version of the problem is known as a maximization binary knapsack problem (KP): given a set of items with profits and sizes, and the capacity value of a knapsack, the objective is to maximize the total profit of selected items in the knapsack while satisfying the capacity constraint (Kellerer et al. 2004). Another widely used variant is the minimization knapsack problem (GüNtzer and Jungnickel 2000; Han and Makino 2009;Csirik 1991), that is, given a set of items associated with costs and sizes, our goal is to minimize the total cost of selected items while ensuring that the knapsack is covered, i.e., the sum of the items' sizes is larger than or equal to a constant. GüNtzer and Jungnickel (2000) consider a minimization equivalent to the maximization KP, where binary selection variables are negated in the formulation. A generalization of the KP is the multi-dimensional knapsack (MKP), which involves multiple constraints, each corresponding to a knapsack, whose capacity should be satisfied. The structure of the MKP is quite similar to the general integer programs with multiple constraints. When the costs of items are not fixed and are instead given by a probability distribution, the problem is broadly referred to as the stochastic knapsack problem (Morton and Wood 1998;Dean et al. 2008;Angulo et al. 2016). Angulo et al. (2016) and Guan et al. (2009) present cutting plane approaches to improve the solvability of the stochastic dynamic knapsack problem. In the computational experiments, we tackle solving the general stochastic multistage multi-dimensional knapsack problem (SMKP) with integer and binary variables and a mean-risk objective.

Mean-CVaR multi-stage stochastic knapsack formulation
We adapt the scenario-dominance cuts presented in Sect. 2.4 to solve a class of mean-CVaR stochastic multi-dimensional mixed 0-1 knapsack problems (CVaR-SMKPs) and demonstrate their computational benefit on this class of multi-stage stochastic MIPs. CVaR-SMKP represents the general mean-CVaR multi-stage stochastic mixed-integer programs with binary and continuous variables in the first and subsequent stages.
To formulate the CVaR-SMKP, we define the following notation: Indices: t: index for time period, where t ∈ T = {1, . . . , T }, and T is the final time period. i: index for item, where i ∈ I = {1, . . . , I }, and I is the maximum number of items. ω: index for scenario, where ω ∈ = 1, . . . ,¯ , and¯ is the maximum number of scenarios.

Variables:
• x ω it ∈ B: 0 − 1 binary variable to select an item i in period t under scenario ω ∈ , where B is the set of binary integers. • z ω t ∈ R + : continuous variable in period t under scenario ω ∈ . • y ω t ∈ B: 0 − 1 binary variable that relates to the variable x ω it in period t under scenario ω ∈ .
• η ω t ∈ R: VaR variable in period t under scenario ω ∈ . • v ω t ∈ R + : continuous variable that represents the excess cost over the VaR variable in period t under scenario ω ∈ .
All parameters defined below are real numbers, e.g., belong to the field of R: Stochastic parameters: q ω t : cost coefficient associated with variable y ω t in period t under scenario ω ∈ . r ω t : coefficient corresponding to the continuous variable z ω t in period t under scenario ω ∈ . b ω t : minimum total value of all selected items in period t under scenario ω ∈ . Deterministic parameters: c it : cost coefficient associated with variable x ω it for item i in period t. d t : cost coefficient associated with variable z ω t in period t. a it : value of item i in period t. κ it : size of item i in period t. w t : coefficient corresponding to the continuous variable y ω t in period t. h t : minimum total size of all selected items in period t.
All objective coefficients are assumed to be non-negative. Using the above notation, the CVaR-SMKP (28a) can be formulated as follows. x x z ω The mean-risk function (28a) minimizes the expected sum of knapsack costs and the weighted VaR and the weighted positive expected excess cost over the VaR over all items i ∈ I, periods t ∈ T , and scenarios ω ∈ . Constraints (28b) and (28c) define knapsack-related constraints, where constraint (28b) [(28c)] ensures that the sum of the value [size] of selected items should be larger than a given constant. Constraints (28d) define VaR and the excess cost under each scenario over the VaR variable η ω t . Constraints (28e)-(28g) represent the nonanticipativity constraints for the x, z, y, η, and v variables, respectively. Finally, (28h) represent binary integer restrictions on x and y variables, while constraints (28i) represent lower bounds on the z, η, and v variables, respectively. Note that CVaR-SMKP formulated above includes uncertainty in the objective function parameter q ω t , the left-hand side parameter r ω t , and the right-hand side parameter b ω t .

Bounds
Proposition 10 Let S ⊆ be a subset of scenarios. For a given ω ∈ S, let (ẋ ω ,ż ω ,ẏ ω ,η ω ,v ω ) and Z ω be the optimal solution and the corresponding objective function value of the scenario-ξ ω problem, respectively. Then an optimal solution to the CVaR-SMKP (28a), (x * , z * , y * , η * , v * ), satisfies the following inequality: Proof In the CVaR-SMKP (28a), we assume that all objective coefficients are non-negative, implying that the specific scenario objective function values are also non-negative. Then, the result follows immediately by applying inequality (13) in Proposition 8 for the CVaR-SMKP.
Example 3 Consider a CVaR-SMKP instance in which T = 4 and I = 2. This instance has = 8 scenarios, and the probability of each scenario is set as 1/8. The risk parameter λ is set as 1 and α t = 0.95 for each t ∈ T . We assume that all parameter values are set at 0 at period t = 1. The data pertaining the instance is: c 1t = (59, 35, 67  Now, let S = {1, 2, 3, 4}. Then Z 1 = 69.9, Z 2 = 80.6, Z 3 = 72.39, and Z 4 = 83.14. For S the lower-bound inequality (29) is: The minimum upper-bound is attained when solving the scenario problem with ω = 4, and thus the upper-bound inequality (30) is: The optimal objective value for this instance is 614.08. For S = {1, 2, 3, 4, 5, 6, 7, 8}, the lower bound of the inequality (31) can be improved to 614.08, while the upper bound in the inequality (32) stays the same.

Stage-t scenario-dominance cuts
Theorem 2 Let ξ ω be a scenario with ω ∈ S ⊆ . For each t ∈ T , let (ẋ t ω ,ż t ω ,ẏ t ω ,η t ω , v t ω ) represent a vector of optimal solutions to the stage-t scenario-ξ ω problem (10) and define the corresponding optimal objective value as: Then, the portion of the optimal solution up to stage t to the CVaR-SMKP (28a) corresponding to scenario ξ ω such that ξ t , is satisfied by the following inequality: Proof The result follows immediately by applying the inequality (27) in Corollary 1.1 for the mean-CVaR stochastic multi-stage knapsack problem (28a).

Corollary 2.1 The inequality (33) can be strengthened as in the following inequality:
Proof Because the solution (x ω * , z ω * , y ω * , η ω * , v ω * ) is a part of a feasible solution to the full problem defined in (28a), ∀ ω, ω ∈ such that ξ t ω ξ t ω we havē Because we have p ω ≥ p ω , we have C ≥C, and thus the following inequalities hold: Because the inequality (33) implies that A ≥ B, by the results above in (37) the inequality (34), which implies C ≥ B, strengthens the inequality (33).

Strong scenario-dominance cuts
Proposition 12 Let ξ ω and ξ ω be two scenarios that share a common history up to time represent a vector of optimal solutions to the stage-t scenario-ξ ω problem and define the objective value corresponding to the stage-t scenario-ξ ω problem at the solution

Then, there exists a feasible solution to the CVaR-SMKP (28a) corresponding to scenario-ξ
, that is satisfied by the following inequality: Remark 6 It is easy to see that there exists a feasible solution to the CVaR-SMKP (28a) that is satisfied by the inequalities (38) because those inequalities provide a lower bound on the x-, z-, y-, η-, and v-variables. Using the inequalities (38), we often get the optimal solutions, but sometimes they may cut-off the optimal solution. This is becauseẐ ω t is not optimal but a close-to-optimal objective value to the stage-t scenario-ξ ω problem. Therefore, one needs to carefully implement the strong scenario dominance cuts. In Sect. 4, we investigate the tradeoff between the quality of the best solution obtained by applying those inequalities versus the computational speed-up obtained with them. We, specifically, present the objective value gap defined by the percentage deviation between the best objective value found by CPLEX for the original problem and the best objective value obtained after applying the cuts (38) to analyze the impact of adding these aggressive inequalities that may cut-off the optimal solution. On the other hand, the stage-t scenario dominance inequalities (34) do not cut-off the optimal solution and reduce the optimality gap considerably, as shown in Section 4. However, the inequalities (38) may be preferable to the inequalities (34) if the user prefers a computational speed-up.

Implementation specifications
In our computational experiments, we demonstrate the effectiveness of the proposed bounds and cuts driven by "scenario dominance" in solving risk-averse multi-stage stochastic knapsack instances with different problem characteristics. In particular, we solve a variety of CVaR-SMKP instances with the following approaches: • cpx: direct solution of the model (28a) by CPLEX using its default settings.
The implementation of the cuts defined above involves appending them to the CVaR-SMKP formulations as CPLEX user cuts, and then solving the resulting model using CPLEX 12.71 with its default settings. The best feasible solution found within 50 CPU seconds (s) is fed into the optimizer before we run the model with the cuts in order to benefit from useful CPLEX cuts. The 50 s spent a prior to adding sdcut is also included in the cut generation time. We used a formula of 2 T −2 T −1 to determine the number of scenarios to be used for cut generation for solving the CVaR-SMKP. To drive this formula, we performed a preliminary experimental analysis on the fraction of the scenarios that are used for cut generation. After using various fractions, we observe that selecting the half of the scenarios divided by T − 1 strikes a good balance in terms of the cut generation time (which includes the solution of all selected scenario sub-problems) and the time benefit obtained from implementing the cuts. However, the number of scenarios to be used for cut generation could be set any number from 1 to¯ . We generated both scut and sdcut for periods t = {T /2, . . . , T } ∈T . The effectiveness of the proposed cuts is evaluated in terms of reducing the solution time, the integrality gap, the number of nodes solved in the branch and bound tree as well as the number of instances solved without memory problems relative to cpx, which is the direct solution of the model by CPLEX.
All computational testing was conducted on a desktop computer running with Intel i7 CPU and 64.0 GB of memory and using CPLEX 12.7.1. A time limit of 7200 CPU seconds was imposed for all test instances. We describe our test instance generation scheme for CVaR-SMKPs in Sect. 4.2, and present the detailed results in Sect. 4.3.

Instance generation and test problems
For each problem class, we randomly generate multiple instances using different seeds from a uniform distribution. All instances have the parameter T that defines the number of stages in the stochastic scenario tree with l = 2 branches emanating from each node. Thus the resulting scenario tree has a size of 2 T −1 scenarios. All CVaR-SMKP instances generated can be found at the Systems Optimization and Data Analytics Laboratory Mean-CVaR Stochastic Multi-Dimensional Mixed 0-1 Knapsack Test Problem Library (SODAL CVaRSMKP-LIB) (Büyüktahtakın 2021).

Knapsack instances
We generate test instances for CVaR-SMKP in a similar fashion to the two-stage binary knapsack instances of Angulo et al. (2016). Different than the study of Angulo et al. (2016), which solved a two-stage risk-neutral version of the multi-item knapsack problem with a maximum of 120 scenarios and 55 constraints, in our study, we consider a larger, risk-averse version of the problem with up to T = 10 stages, 512 scenarios, and 67,072 constraints. Also, rather than using only binary variables as in the formulation of Angulo et al. (2016), we consider a mix of binary and continuous variables, further complicating the multi-item knapsack problem considered.
We assume that all parameter values at period t = 1 are set at zero since uncertainty starts at stage two. The parameters c it , κ it , d t , a it , and w t are independent and identically distributed (i.i.d.) sampled from the uniform distribution over {1, . . . , 100}, e.g., U [1, R] where R = 100 and κ it ∈ U [1, R] for t = 2, . . . , T . We have considered uncertainty in the objective function parameter q ω t , the left-hand side parameter r ω t , and the right-hand side parameter b ω t . For the 2-branch problem, all uncertain parameters have two levels: low (L) and high (H), where q ω t , r ω t ∈ U [1, R/2] for the low level and q ω t , r ω t ∈ U [R/2 + 1, R] for the high level and b ω t = 3 4 t j=1 We set h t = 3 4 i∈I κ it + w t . By employing this scheme, we randomly generate instances for various number of stages (T ) and the number of items (I ). Ten random feasible instances for each stage and item combination (T , I ) ∈ {(5, 120), (6, 50), (7, 40), (8, 20), (9, 15), (10, 10)} are generated, resulting in a total number of 60 instances. As we focus on large scale problems of the form (5a), we assume p ω = 1/¯ , where¯ is the maximum number of scenarios, for the rest of the paper.
Those instances generated represent large-scale and hard CVaR-SMKP formulations with both the number of binary and continuous variables varying from 9656 to 71,140, and the number of constraints varying from 6096 to 67,072. The size and the difficulty of those instances increase with T and I .

Results for risk-averse multi-dimensional knapsack instances
In Tables 2, 3, 4, 5 and 6, we report results for a combination of columns defined as follows.
• I: number of items.
• time: CPU seconds required to solve the problem, excluding inequality generation time.
• Ctime: CPU seconds required to generate the cuts, including the solution of the single scenario sub-problems. • Ttime: total CPU seconds required to solve the problem, including inequality generation time.
• Tfac: time factor improvement by sdcut over cpx (Tfac= Time cpx /Time sdcut ), where Time cpx is the solution time (time) by cpx and Time sdcut is the total solution time (Ttime) by bc, sdc, or sdcut. • cut: number of cuts added to CPLEX as user inequalities.
• node (100K): number of nodes explored in the branch and bound tree in 100,000s.
• gap 1 : final optimality gap by cpx at the end of the time limit.
• gap 2 : percentage deviation between the best objective value found by cpx (obj cpx ) and the best objective value obtained by the considered cutting-plane method at the end of the time limit, e.g., gap 2 = 100 × (obj sdcut /obj cpx -1), where obj sdcut defines the best objective value obtained by sdcut. • IGap: percentage integrality gap of the formulation before inequalities are added (InitGap = 100 × (bestobj − relaxlb)/bestobj), where relaxlb and bestobj are objective function values of the initial LP relaxation and the best feasible solution by cpx, respectively. • RGap: percentage integrality gap of the formulation after inequalities are added (RootGap = 100 × (bestobj − rootlb)/bestobj), where rootlb is the objective function value of the initial LP relaxation after scut and bc are added. • GapImp: percentage improvement in the integrality gap at the root node (GapImp = 100 × (1-RGap/IGap). • unsol: the number of instances unsolved to optimality within the given time limit.
In all tables, the values under the time-related metrics (time, Ctime, Ttime, and Tfac) are rounded to the closest integer to improve the clarity of the presentation. Table 2 presents the initial gap (IGap), cut, CTime, the root gap (RGap), and the optimality gap improvement (GapImp) over cpx using the bc and scut together (bc+scut), where we solve both 2 T −2 and 2 T −1 number of scenario sub-problems for only one time period t = T . Here, we did not include the sdcut because the strong dominance cuts (38) are rather aggressive and may potentially cut-off the optimal solutions. The IGap, RGap, and GapImp columns are presented in Table 2. Each row of Table 2 presents results for an average of ten instances, while the Average row represents results for an average of 60 instances.

Results on bounds and Sdc
In Table 2, the percent initial gap shows an increasing trend as T increases, reaching 5.65% for instances with T = 10 and I = 10. From the table, the gap improvement due to bc+scut also shows an increasing trend as T increases. The average integrality gap of 2.37% is reduced to 1.25% and 0.10% by the implementation of bc+scut when 2 T −2 and 2 T −1 scenario sub-problems are solved, respectively. When cut generation is varied from solving 2 T −2 to solving 2 T −1 scenario sub-problems, the number of cuts generated as well as the cut generation time nearly double, as shown in Table 2. On the other hand, the average gap improvement due to scut and bc is 46.4%, when 2 T −2 scenarios are solved, while it increases to 94.3% when the number of scenario sub-problems solved increases to 2 T −1 . In other words, the gap improvement due to the bc + scut nearly doubles when the number of scenario sub-problems solved is doubled. Table 5 in Appendix B presents results that compare the performance of cpx, bc, scut, and sdcut instance-by-instance for the CVaR-SMKP test set with T = 6, I = 50, where λ = 1 and α t = 0.95. Each row of Table 5 represents an instance as specified under the Instance column, except the Average row block, which gives the average results over ten instances with T = 6 and I = 50. Here, bc and scut are implemented by solving 2 T −1 number of scenario sub-problems while sdcut is generated using 2 T −2 scenario sub-problems. In terms of the solution time performance, sdcut performs the best by reducing the cpx solution time from two hours to only 68 CPU seconds. The bc and scut do not help much reducing the solution time, except the instance number 4, where two hours of cpx solution time is reduced to 180 and 203 CPU seconds, respectively. The final MIP optimality gap of 0.08% by cpx is improved to 0.05 and 0.04% by bc and scut, respectively. The main benefit of bc and scut is observed in reducing the integrality gap by 93.2% on average, as shown under the GapImp column. Since the bc and scut would provide the highest benefit to close the gap from the linear programming relaxation of the problem, we recommend using both bc and scut to close the integrality gap in the state-of-the-art solvers, such as CPLEX, rather than adding those inequalities as cuts and solving the augmented problem for computational speed-up. Table 3 summarizes results using sdcut for instances with a variety of stages (T = 5 − 10) and items (I = 10 − 120) where λ = 1 and α t = 0.95. Each row of Table 3 presents results for an average of ten instances, while Average (Total) row reports the overall average (total) over sixty instances.

Results on Sdcut
Here, we do not include the bc because the initial experiments with the lower and upper bounding inequalities (29) and (30), shown in Section 4.5, do not provide any solution speedup over the sdcut, despite they close the initial integrality gap substantially. The difficulty of the instances presented in Table 3 is a factor of the time stage T as well as the number of items I . In all implementations, sdcut significantly reduces the solution time by cpx. For example, the solution time of 7287 is reduced to 68 CPU seconds for instances with T = 6 and I = 50, and for the T = 10 and I = 10 instances, the solution time of 5780 is reduced to 227 CPU seconds due to the sdcut. In particular, the average solution time is reduced from 6, 559 to 111 seconds, where 98 seconds are spent for cut generation, including the solution of the scenario sub-problems. When the total computation time is compared, a total of 39, 354 solution time of the cpx is reduced to 669 CPU seconds by the sdcut to solve all instances. While sdcut improves the solution time of the T = 5 and I = 120 instances with a factor of 110, it improves the solution speed of the T = 10 and I = 10 instances with a factor of 26. As shown in the overall averages in Table 3, sdcut improve results by a factor of 72, when all of the time factor improvement due to sdcut is averaged. The solution time reduction due to the sdcut is also consistent with the reduction in the number of nodes searched in the branch and bound tree. As the time period is increased from T = 5 to T = 10, the number of scenarios increases from 16 to 512. Because 2 T −2 T −1 * T −1 2 scenarios are solved as single scenario sub-problems, the time required for the cut generation (Ctime) shows an increasing trend as T increases. The average gap that defines the percentage deviation between obj obtained by cpx and obj by sdcut is 0.03%. This means that when all instances from T = 5 to T = 10 are averaged, the objective value provided by sdcut is 0.03% higher than the objective value provided by cpx. Out of 60 instances, cpx could not solve 53 of them to optimality within the 2-h of time limit, as shown under the unsol of Table 3, while the implementation with sdcut solves all of the 60 instances with only 0.03% deviation from cpx, within the 2-hour of time limit.
The Total row shows that a total time of 39,354 CPU seconds for solving 60 problems is reduced to only 669 CPU seconds (less than 11 minutes) by using sdcut, resulting in a time factor improvement of 432 for solving all instances presented in this table. Varying risk parameters does not change the performance of sdcut much (see detailed results for varying risk parameters in Table 6 and the related discussion under Appendix C).

Generalization to multiple random variables
We perform additional experiments to demonstrate the number of scenario-dominance relations and its impact on the performance of the sdcut as we vary the number of random variables. For those experiments, we have generated additional instances where we re-defined the parameter c it as an uncertain parameter c ω it , which represents the cost coefficient associated with variable x ω it for item i in period t under scenario ω, and assumed that all other parameters are deterministic. The random parameters c ω it are stage-wise independent but are jointly distributed over i for each t ∈ T . The c ω it could be positively or negatively correlated for each item i. A scenario ξ ω stage-t dominates a scenario ξ ω with respect to the problem (28a) if p ω ≥ p ω and c ω it ≥ c ω it hold for each i ∈ I and j = 2, . . . , t with t ∈ T \{1}. Since the case where each univariate component of the random vector is independently distributed could explode the size of the scenario tree, here we generated random variables as stage-wise independent but jointly distributed at each stage. However, results provide insights into the former case as well.
In Table 4, we compare the performance of sdcut with respect to cpx for a varying number of random variables in the scenario tree for instances with T = 10 and I = 10. In this table, I represents the number of randomly generated items in each stage t, RV shows the total number of random variables generated, Sd (t = 2) and Sd (t = T ) represent the number of scenario-dominance relations at stages two and T , and Sd (total) presents the total number of scenario-dominance relations over all time periods for each RV. For example, instances withÎ = 0 and RV = 0 are deterministic, whileÎ = 3 represents instances with three items randomly generated in each stage and thus that have a total ofÎ (T − 1) = 27 random variables since the values of the random parameters at stage 1 are assumed to be known.
Each row of Table 4 presents results for an average of ten instances. Since random variables are stage-wise independent but are dependent in each stage, the number of scenarios for each instance with a varying number of random variables is 512.
According to the average solution times in Table 4, sdcut consistently improves the solution times of cpx for instances with an increasing number of random variables by an overall average factor of 25, reducing the average solution time of 5,716 seconds to 231 seconds. As the number of random variables increases, the number of second-stage and last-stage scenario-dominance relations reduces. Note that the second-stage dominance relations are much larger than the last-stage dominance relations. When we reach 81 random variables, we also hit the minimum number scenario relations at stages two and T , and in total, as shown by Proposition 3. For the considered instances, having the minimum number scenario relations does not impact the performance of sdcut with respect to cpx. For example, when we haveÎ = 10 and RV = 90, we have the least total number of scenario-dominance relations (W = 261, 632), and 1705 cuts are generated under sdcut. While mipgap tends to increase as the number of random variables increase, the performance of sdcut is not impacted by an increased number of random variables. The results imply that when we increase the number of random variables in each stage, the algorithm still generates plenty of cutting planes due to an ample number of stage-t dominance relations.

Concluding remarks
In this paper, we introduced a new concept called "scenario dominance" and "stage-t scenario dominance" to effectively solve the general risk-averse multi-stage stochastic mixed-integer programs. This method derives the implications of one scenario on another by comparing the problem coefficients under each scenario and the optimal objective function values of the individual scenario sub-problems. We derive results on the minimum number of stage-t scenario-dominance relations in a multi-stage scenario tree. By studying those implications and the partial ordering of scenarios, we derive new bounds and cutting planes to tackle risk- averse SMIPs. We perform extensive computational experiments using the proposed cuts on the mean-CVaR stochastic dynamic knapsack problem. Our computations demonstrate that the proposed methodology reduces the solution time of the large-scale, risk-averse multi-stage stochastic mixed-integer optimization problems with varying risk characteristics by one-totwo orders of magnitude. Results also demonstrate that strong dominance cuts perform well for those instances with ten random variables in each stage and ninety random variables in total. They also have the potential to tackle a larger number of random variables because one can generate plenty of stage-t scenario-dominance relations even under the worst-case scenario. This study opens several exciting research directions worth investigating in the future. For example, "scenario dominance" approach could be investigated for other risk measures, such as the Expected Excess of a given target and the Excess Probability (Schultz and Tiedemann 2003). In this paper, we only used a small subset of the full scenario set to generate the scenario dominance cuts. A future experimental study could analyze the impact of using different fractions and subsets of scenarios and various scenario selection schemes on the quality of the bounds as well as the cuts proposed. Decomposition methods involving Lagrangian relaxation (CarøE and Schultz 1999;Nowak and Römisch 2000), column generation (Singh et al. 2009;Sen 2005), L-shaped decomposition (Angulo et al. 2016;Laporte and Louveaux 1993), progressive hedging (Watson and Woodruff 2011;Gade et al. 2016), scenario-decomposition (Birge 1985;Sandıkçı and Özaltın 2017), and stochastic dual dynamic integer programming (Zou et al. 2019) have been quite useful in solving various classes of multi-stage stochastic integer programs. Thus, another important future research direction is to integrate the proposed "scenario dominance" cut-and-bound generation scheme within a decomposition framework for solving mean-risk M-SMIPs.
In this paper, we focus on the type of stochastic multi-dimensional knapsack instances up to T = 10 stages, 10 random variables in each stage, and I = 120 items. The considered CVaR-SMKP deterministic counterparts, which involve the expected values of the uncertain parameters, are also computationally challenging. Future research could investigate a variant of CVaR-SMKP instances, which have easier deterministic formulations but a larger number of scenarios than considered in this work, and thus remain computationally hard to solve. For example, a future extension of this work could consider T = 5-stage I = 1-item CVaR-SMKPs withÎ = 10 independent random variables and 10 branches in each stage, with 10,000 scenarios. To tackle an instance like this, a parallel distributed platform would be more suitable for solving the individual scenario problems, reducing the dominance cut generation time, and thus improving the overall solution performance. In such a case, based on Proposition 4, we would have a minimum of over more than a hundred million scenario dominance relations, leading to an enormous number of scenario dominance cuts. Thus another important direction for research is to determine the most useful subset and the number of dominance cuts, which would not excessively increase the size of the problem formulation.
The concept of "scenario dominance" as defined in this paper is fundamentally different than the "stochastic dominance" relation defined in the literature (Müller and Stoyan 2002;Luedtke 2008) because two random variables are compared based on their probability distribution function in stochastic dominance, while in our scenario dominance approach we pairwise compare scenarios based on the realization of uncertain parameters at each time stage. Specifically, scenario dominance derives relations among scenarios by comparing possible outcomes of a single random variable or potential joint outcomes of more than one random variable. Yet, an interesting future direction would be to investigate the theoretical relations between the scenario dominance and the first (second)-order stochastic dominance.
The proposed approach in this paper does not require convexity or linearity, and thus can be extended to study the non-linear case and various types of risk-averse multi-stage stochastic MIPs. The generality of our method would require careful implementation of the proposed cuts for each specific problem, such as risk-averse multi-stage stochastic financial investment planning, capacity allocation, and network optimization problems.
Our approach also has some limitations. For example, problems with several uncorrelated random variables may exhibit less scenario-dominance relations and fewer dominance cuts than considered in this study. Therefore, an important future direction is to investigate further the boundaries of the scenario dominance cuts by exploring the breadth of problems that could benefit from our approach.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Table 5 presents results that compare the performance of cpx, bc, scut, and sdcut instanceby-instance for the CVaR-SMKP test set with T = 6, I = 50, where λ = 1 and α t = 0.95.   Table 6 summarizes results for (T , I ) = (5, 120) with varying values of risk parameters λ and α t . In particular, we consider 15 combinations of λ ranging from 0.001 to 1000 and α t ∈ {0.7, 0.95}. From each combination, three instances are generated with (T , I ) = (5, 120), resulting in 45 instances that differ in their λ and α t values, thus each row of Table 6 presents an average of the results for three instances. The Ttime (cpx) and Ttime (sdcut) columns of Table 6 gives the CPU times of the CPLEX performance for solving the CVaR-SMKP instances without cuts, and with sdcut, respectively, Ctime gives the sdcut generation time, gap 1 presents the cpx MIP gap, and gap 2 gives the percentage deviation between the best objective value found by cpx and sdcut. According to the total and average solution times, the solution times of cpx and sdcut are similar for varying λ and α t values. Thus we do not observe a significant change in instance difficulty and computational times with respect to the λ and α parameters. For all cases with varying λ and α t values, the sdcut drastically improves the performance with respect to cpx with a factor over 110 either when averaged or in total. Out of 45 cases, 45 of them cannot be solved to optimality by cpx within the 2-h limit, while all instances are solved within 2 minutes with the sdcut with an average optimality gap of 0.01%.